complex.hpp

Go to the documentation of this file.

//=======================================================================
// Copyright (c) 2014-2023 Baptiste Wicht
// Distributed under the terms of the MIT License.
// (See accompanying file LICENSE or copy at
//  http://opensource.org/licenses/MIT)
//=======================================================================

#pragma once

namespace etl {

template <typename T>
struct complex {
    using value_type = T; 

    value_type real; 
    value_type imag; 

    constexpr complex(const T& re = T(), const T& im = T()) : real(re), imag(im) {
        // Nothing else to init
    }

    constexpr complex(const complex& rhs) = default;

    template <typename X>
    constexpr complex(const complex<X>& rhs) : real(rhs.real), imag(rhs.imag) {
        // Nothing else to init
    }

    complex& operator=(const T& rhs) noexcept {
        real = rhs;
        imag = 0.0;

        return *this;
    }

    complex& operator=(const complex& rhs) noexcept = default;

    complex& operator=(const std::complex<T>& rhs) noexcept {
        real = rhs.real();
        imag = rhs.imag();

        return *this;
    }

    template <typename X>
    complex& operator+=(const complex<X>& rhs) {
        real += rhs.real;
        imag += rhs.imag;

        return *this;
    }

    template <typename X>
    complex& operator-=(const complex<X>& rhs) {
        real -= rhs.real;
        imag -= rhs.imag;

        return *this;
    }

    template <typename X>
    complex& operator*=(const complex<X>& rhs) {
        T ac = real * rhs.real;
        T bd = imag * rhs.imag;

        T bc = imag * rhs.real;
        T ad = real * rhs.imag;

        real = ac - bd;
        imag = bc + ad;

        return *this;
    }

    template <typename X>
    complex& operator/=(const complex<X>& rhs) {
        T ac = real * rhs.real;
        T bd = imag * rhs.imag;

        T bc = imag * rhs.real;
        T ad = real * rhs.imag;

        T frac = rhs.real * rhs.real + rhs.imag * rhs.imag;

        real = (ac + bd) / frac;
        imag = (bc - ad) / frac;

        return *this;
    }
};

template <typename T>
inline bool operator==(const complex<T>& lhs, const complex<T>& rhs) {
    return lhs.real == rhs.real && lhs.imag == rhs.imag;
}

template <typename T>
inline bool operator!=(const complex<T>& lhs, const complex<T>& rhs) {
    return !(lhs == rhs);
}

template <typename T>
inline complex<T> operator-(complex<T> rhs) {
    return {-rhs.real, -rhs.imag};
}

template <typename T>
inline complex<T> operator+(const complex<T>& lhs, const complex<T>& rhs) {
    return {lhs.real + rhs.real, lhs.imag + rhs.imag};
}

template <typename T>
inline complex<T> operator-(const complex<T>& lhs, const complex<T>& rhs) {
    return {lhs.real - rhs.real, lhs.imag - rhs.imag};
}

template <typename T>
inline complex<T> operator*(const complex<T>& lhs, const complex<T>& rhs) {
    T ac = lhs.real * rhs.real;
    T bd = lhs.imag * rhs.imag;

    T bc = lhs.imag * rhs.real;
    T ad = lhs.real * rhs.imag;

    return {ac - bd, bc + ad};
}

template <typename T>
inline complex<T> operator*(const complex<T>& lhs, T rhs) {
    return {lhs.real * rhs, lhs.imag * rhs};
}

template <typename T>
inline complex<T> operator*(T lhs, const complex<T>& rhs) {
    return {lhs * rhs.real, lhs * rhs.imag};
}

template <typename T>
inline complex<T> operator/(const complex<T>& lhs, T rhs) {
    return {lhs.real / rhs, lhs.imag / rhs};
}

template <typename T>
inline complex<T> operator/(const complex<T>& lhs, const complex<T>& rhs) {
    T ac = lhs.real * rhs.real;
    T bd = lhs.imag * rhs.imag;

    T bc = lhs.imag * rhs.real;
    T ad = lhs.real * rhs.imag;

    T frac = rhs.real * rhs.real + rhs.imag * rhs.imag;

    return {(ac + bd) / frac, (bc - ad) / frac};
}

template <typename T>
T abs(complex<T> z) {
    auto x = z.real;
    auto y = z.imag;
    auto s = std::max(std::abs(x), std::abs(y));

    if (s == T()) {
        return s;
    }

    x = x / s;
    y = y / s;

    return s * std::sqrt(x * x + y * y);
}

template <typename T>
T arg(complex<T> z) {
    auto x = z.real;
    auto y = z.imag;

    return atan2(y, x);
}

template <typename T>
complex<T> sqrt(complex<T> z) {
    auto x = z.real;
    auto y = z.imag;

    if (x == T()) {
        auto t = std::sqrt(std::abs(y) / 2);
        return {t, y < T() ? -t : t};
    } else {
        auto t = std::sqrt(2 * (abs(z) + std::abs(x)));
        auto u = t / 2;

        if (x > T()) {
            return {u, y / t};
        } else {
            return {std::abs(y) / t, y < T() ? -u : u};
        }
    }
}

template <typename T>
complex<T> invsqrt(complex<T> z) {
    return complex<T>(T(1)) / sqrt(z);
}

template <typename T>
complex<T> cbrt(complex<T> z) {
    auto z_abs = etl::abs(z);
    auto z_arg = etl::arg(z);

    auto new_abs = std::cbrt(z_abs);
    auto new_arg = z_arg / 3.0f;

    return {new_abs * std::cos(new_arg), new_abs * std::sin(new_arg)};
}

template <typename T>
complex<T> invcbrt(complex<T> z) {
    return complex<T>(T(1)) / cbrt(z);
}

template <typename T>
complex<T> log(complex<T> z) {
    return {std::log(etl::abs(z)), etl::arg(z)};
}

template <typename T>
complex<T> log2(complex<T> z) {
    return etl::log(z) / etl::log(etl::complex<T>{T(2)});
}

template <typename T>
complex<T> log10(complex<T> z) {
    return etl::log(z) / etl::log(etl::complex<T>{T(10)});
}

template <typename T>
complex<T> sin(complex<T> z) {
    return {std::sin(z.real) * std::cosh(z.imag), std::cos(z.real) * std::sinh(z.imag)};
}

template <typename T>
complex<T> cos(complex<T> z) {
    return {std::cos(z.real) * std::cosh(z.imag), -std::sin(z.real) * std::sinh(z.imag)};
}

template <typename T>
complex<T> tan(complex<T> z) {
    return sin(z) / cos(z);
}

template <typename T>
complex<T> cosh(complex<T> z) {
    return {std::cosh(z.real) * std::cos(z.imag), std::sinh(z.real) * std::sin(z.imag)};
}

template <typename T>
complex<T> sinh(complex<T> z) {
    return {std::sinh(z.real) * std::cos(z.imag), std::cosh(z.real) * std::sin(z.imag)};
}

template <typename T>
complex<T> tanh(complex<T> z) {
    return sinh(z) / cosh(z);
}

template <typename T>
inline complex<T> inverse(complex<T> x) {
    return {x.imag, x.real};
}

template <typename T>
inline complex<T> inverse_conj(complex<T> x) {
    return {-x.imag, x.real};
}

template <typename T>
inline complex<T> conj_inverse(complex<T> x) {
    return {x.imag, -x.real};
}

template <typename T>
inline complex<T> conj(const complex<T>& c) {
    return {c.real, -c.imag};
}

template <typename T>
inline T get_imag(const std::complex<T>& c) {
    return c.imag();
}

template <typename T>
inline T get_imag(const etl::complex<T>& c) {
    return c.imag;
}

template <typename T>
inline T get_real(const std::complex<T>& c) {
    return c.real();
}

template <typename T>
inline T get_real(const etl::complex<T>& c) {
    return c.real;
}

template <typename T>
inline std::complex<T> get_conj(const std::complex<T>& c) {
    return std::conj(c);
}

template <typename T>
inline etl::complex<T> get_conj(const etl::complex<T>& c) {
    return {c.real, -c.imag};
}

template <typename T>
std::ostream& operator<<(std::ostream& os, const etl::complex<T>& c) {
    return os << "C(" << c.real << "," << c.imag << ")";
}

} //end of namespace etl