/* * Copyright (c) 2021, Cesar Torres * * SPDX-License-Identifier: BSD-2-Clause */ #pragma once #include #include namespace AK { template class [[gnu::packed]] Complex { public: constexpr Complex() : m_real(0) , m_imag(0) { } constexpr Complex(T real) : m_real(real) , m_imag((T)0) { } constexpr Complex(T real, T imaginary) : m_real(real) , m_imag(imaginary) { } constexpr T real() const noexcept { return m_real; } constexpr T imag() const noexcept { return m_imag; } constexpr T magnitude_squared() const noexcept { return m_real * m_real + m_imag * m_imag; } constexpr T magnitude() const noexcept { return hypot(m_real, m_imag); } constexpr T phase() const noexcept { return atan2(m_imag, m_real); } template static constexpr Complex from_polar(U magnitude, V phase) { V s, c; sincos(phase, s, c); return Complex(magnitude * c, magnitude * s); } template constexpr Complex& operator=(Complex const& other) { m_real = other.real(); m_imag = other.imag(); return *this; } template constexpr Complex& operator=(U const& x) { m_real = x; m_imag = 0; return *this; } template constexpr Complex operator+=(Complex const& x) { m_real += x.real(); m_imag += x.imag(); return *this; } template constexpr Complex operator+=(U const& x) { m_real += x; return *this; } template constexpr Complex operator-=(Complex const& x) { m_real -= x.real(); m_imag -= x.imag(); return *this; } template constexpr Complex operator-=(U const& x) { m_real -= x; return *this; } template constexpr Complex operator*=(Complex const& x) { T const real = m_real; m_real = real * x.real() - m_imag * x.imag(); m_imag = real * x.imag() + m_imag * x.real(); return *this; } template constexpr Complex operator*=(U const& x) { m_real *= x; m_imag *= x; return *this; } template constexpr Complex operator/=(Complex const& x) { T const real = m_real; T const divisor = x.real() * x.real() + x.imag() * x.imag(); m_real = (real * x.real() + m_imag * x.imag()) / divisor; m_imag = (m_imag * x.real() - x.real() * x.imag()) / divisor; return *this; } template constexpr Complex operator/=(U const& x) { m_real /= x; m_imag /= x; return *this; } template constexpr Complex operator+(Complex const& a) { Complex x = *this; x += a; return x; } template constexpr Complex operator+(U const& a) { Complex x = *this; x += a; return x; } template constexpr Complex operator-(Complex const& a) { Complex x = *this; x -= a; return x; } template constexpr Complex operator-(U const& a) { Complex x = *this; x -= a; return x; } template constexpr Complex operator*(Complex const& a) { Complex x = *this; x *= a; return x; } template constexpr Complex operator*(U const& a) { Complex x = *this; x *= a; return x; } template constexpr Complex operator/(Complex const& a) { Complex x = *this; x /= a; return x; } template constexpr Complex operator/(U const& a) { Complex x = *this; x /= a; return x; } template constexpr bool operator==(Complex const& a) const { return (this->real() == a.real()) && (this->imag() == a.imag()); } constexpr Complex operator+() { return *this; } constexpr Complex operator-() { return Complex(-m_real, -m_imag); } private: T m_real; T m_imag; }; // reverse associativity operators for scalars template constexpr Complex operator+(U const& a, Complex const& b) { Complex x = a; x += b; return x; } template constexpr Complex operator-(U const& a, Complex const& b) { Complex x = a; x -= b; return x; } template constexpr Complex operator*(U const& a, Complex const& b) { Complex x = a; x *= b; return x; } template constexpr Complex operator/(U const& a, Complex const& b) { Complex x = a; x /= b; return x; } // some identities template static constinit Complex complex_real_unit = Complex((T)1, (T)0); template static constinit Complex complex_imag_unit = Complex((T)0, (T)1); template static constexpr bool approx_eq(Complex const& a, Complex const& b, double const margin = 0.000001) { auto const x = const_cast&>(a) - const_cast&>(b); return x.magnitude() <= margin; } // complex version of exp() template static constexpr Complex cexp(Complex const& a) { // FIXME: this can probably be faster and not use so many "expensive" trigonometric functions return exp(a.real()) * Complex(cos(a.imag()), sin(a.imag())); } } template struct AK::Formatter> : Formatter { ErrorOr format(FormatBuilder& builder, AK::Complex c) { return Formatter::format(builder, TRY(String::formatted("{}{:+}i", c.real(), c.imag()))); } }; #if USING_AK_GLOBALLY using AK::approx_eq; using AK::cexp; using AK::Complex; using AK::complex_imag_unit; using AK::complex_real_unit; #endif