From f4f5a1c0e748ba4681e001d44a59211abbea3261 Mon Sep 17 00:00:00 2001 From: Cesar Torres Date: Fri, 19 Mar 2021 22:23:48 +0100 Subject: [PATCH] AK: Add complex number library Useful for diverse algorithms. Also added some tests for it. --- AK/Complex.h | 342 +++++++++++++++++++++++++++++++++++++++ AK/Concepts.h | 4 + AK/Tests/CMakeLists.txt | 1 + AK/Tests/TestComplex.cpp | 65 ++++++++ 4 files changed, 412 insertions(+) create mode 100644 AK/Complex.h create mode 100644 AK/Tests/TestComplex.cpp diff --git a/AK/Complex.h b/AK/Complex.h new file mode 100644 index 00000000000..832ead8b962 --- /dev/null +++ b/AK/Complex.h @@ -0,0 +1,342 @@ +/* + * Copyright (c) 2021, Cesar Torres + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#pragma once + +#include +#if __has_include() +# define AKCOMPLEX_CAN_USE_MATH_H +# include +#endif + +#ifdef __cplusplus +# if __cplusplus >= 201103L +# define COMPLEX_NOEXCEPT noexcept +# endif +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 COMPLEX_NOEXCEPT { return m_real; } + + constexpr T imag() const COMPLEX_NOEXCEPT { return m_imag; } + + constexpr T magnitude_squared() const COMPLEX_NOEXCEPT { return m_real * m_real + m_imag * m_imag; } + +# ifdef AKCOMPLEX_CAN_USE_MATH_H + constexpr T magnitude() const COMPLEX_NOEXCEPT + { + // for numbers 32 or under bit long we don't need the extra precision of sqrt + // although it may return values with a considerable error if real and imag are too big? + if constexpr (sizeof(T) <= sizeof(float)) { + return sqrtf(m_real * m_real + m_imag * m_imag); + } else if constexpr (sizeof(T) <= sizeof(double)) { + return sqrt(m_real * m_real + m_imag * m_imag); + } else { + return sqrtl(m_real * m_real + m_imag * m_imag); + } + } + + constexpr T phase() const COMPLEX_NOEXCEPT + { + return atan2(m_imag, m_real); + } + + template + static constexpr Complex from_polar(U magnitude, V phase) + { + if constexpr (sizeof(T) <= sizeof(float)) { + return Complex(magnitude * cosf(phase), magnitude * sinf(phase)); + } else if constexpr (sizeof(T) <= sizeof(double)) { + return Complex(magnitude * cos(phase), magnitude * sin(phase)); + } else { + return Complex(magnitude * cosl(phase), magnitude * sinl(phase)); + } + } +# endif + + template + constexpr Complex& operator=(const Complex& other) + { + m_real = other.real(); + m_imag = other.imag(); + return *this; + } + + template + constexpr Complex& operator=(const U& x) + { + m_real = x; + m_imag = 0; + return *this; + } + + template + constexpr Complex operator+=(const Complex& x) + { + m_real += x.real(); + m_imag += x.imag(); + return *this; + } + + template + constexpr Complex operator+=(const U& x) + { + m_real += x.real(); + return *this; + } + + template + constexpr Complex operator-=(const Complex& x) + { + m_real -= x.real(); + m_imag -= x.imag(); + return *this; + } + + template + constexpr Complex operator-=(const U& x) + { + m_real -= x.real(); + return *this; + } + + template + constexpr Complex operator*=(const Complex& x) + { + const T 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*=(const U& x) + { + m_real *= x; + m_imag *= x; + return *this; + } + + template + constexpr Complex operator/=(const Complex& x) + { + const T real = m_real; + const T 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/=(const U& x) + { + m_real /= x; + m_imag /= x; + return *this; + } + + template + constexpr Complex operator+(const Complex& a) + { + Complex x = *this; + x += a; + return x; + } + + template + constexpr Complex operator+(const U& a) + { + Complex x = *this; + x += a; + return x; + } + + template + constexpr Complex operator-(const Complex& a) + { + Complex x = *this; + x -= a; + return x; + } + + template + constexpr Complex operator-(const U& a) + { + Complex x = *this; + x -= a; + return x; + } + + template + constexpr Complex operator*(const Complex& a) + { + Complex x = *this; + x *= a; + return x; + } + + template + constexpr Complex operator*(const U& a) + { + Complex x = *this; + x *= a; + return x; + } + + template + constexpr Complex operator/(const Complex& a) + { + Complex x = *this; + x /= a; + return x; + } + + template + constexpr Complex operator/(const U& a) + { + Complex x = *this; + x /= a; + return x; + } + + template + constexpr bool operator==(const Complex& a) + { + return (this->real() == a.real()) && (this->imag() == a.imag()); + } + + template + constexpr bool operator!=(const Complex& a) + { + return !(*this == a); + } + + 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+(const U& b, const Complex& a) +{ + Complex x = a; + x += b; + return x; +} + +template +constexpr Complex operator-(const U& b, const Complex& a) +{ + Complex x = a; + x -= b; + return x; +} + +template +constexpr Complex operator*(const U& b, const Complex& a) +{ + Complex x = a; + x *= b; + return x; +} + +template +constexpr Complex operator/(const U& b, const Complex& a) +{ + 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); + +# ifdef AKCOMPLEX_CAN_USE_MATH_H +template +static constexpr bool approx_eq(const Complex& a, const Complex& b, const double margin = 0.000001) +{ + const auto x = const_cast&>(a) - const_cast&>(b); + return x.magnitude() <= margin; +} + +//complex version of exp() +template +static constexpr Complex cexp(const Complex& a) +{ + //FIXME: this can probably be faster and not use so many expensive trigonometric functions + if constexpr (sizeof(T) <= sizeof(float)) { + return expf(a.real()) * Complex(cosf(a.imag()), sinf(a.imag())); + } else if constexpr (sizeof(T) <= sizeof(double)) { + return exp(a.real()) * Complex(cos(a.imag()), sin(a.imag())); + } else { + return expl(a.real()) * Complex(cosl(a.imag()), sinl(a.imag())); + } +} +} +# endif + +using AK::Complex; +using AK::complex_imag_unit; +using AK::complex_real_unit; +# ifdef AKCOMPLEX_CAN_USE_MATH_H +using AK::approx_eq; +using AK::cexp; +# endif +#endif diff --git a/AK/Concepts.h b/AK/Concepts.h index 603ed2d367d..a7ef8edb014 100644 --- a/AK/Concepts.h +++ b/AK/Concepts.h @@ -38,12 +38,16 @@ concept Integral = IsIntegral::value; template concept FloatingPoint = IsFloatingPoint::value; +template +concept Arithmetic = IsArithmetic::value; + #endif } #if defined(__cpp_concepts) && !defined(__COVERITY__) +using AK::IsArithmetic; using AK::IsFloatingPoint; using AK::IsIntegral; diff --git a/AK/Tests/CMakeLists.txt b/AK/Tests/CMakeLists.txt index 391bfbd0e06..5fc624697c1 100644 --- a/AK/Tests/CMakeLists.txt +++ b/AK/Tests/CMakeLists.txt @@ -16,6 +16,7 @@ set(AK_TEST_SOURCES TestCircularDeque.cpp TestCircularDuplexStream.cpp TestCircularQueue.cpp + TestComplex.cpp TestDistinctNumeric.cpp TestDoublyLinkedList.cpp TestEndian.cpp diff --git a/AK/Tests/TestComplex.cpp b/AK/Tests/TestComplex.cpp new file mode 100644 index 00000000000..9cccd525f1e --- /dev/null +++ b/AK/Tests/TestComplex.cpp @@ -0,0 +1,65 @@ +/* + * Copyright (c) 2021, Cesar Torres + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#include +#include + +TEST_CASE(Complex) +{ + auto a = Complex { 1.f, 1.f }; + auto b = complex_real_unit + Complex { 0, 1 } * 1; + EXPECT_APPROXIMATE(a.real(), b.real()); + EXPECT_APPROXIMATE(a.imag(), b.imag()); + +#ifdef AKCOMPLEX_CAN_USE_MATH_H + EXPECT_APPROXIMATE((complex_imag_unit - complex_imag_unit).magnitude(), 0); + EXPECT_APPROXIMATE((complex_imag_unit + complex_real_unit).magnitude(), sqrt(2)); + + auto c = Complex { 0., 1. }; + auto d = Complex::from_polar(1., M_PI / 2.); + EXPECT_APPROXIMATE(c.real(), d.real()); + EXPECT_APPROXIMATE(c.imag(), d.imag()); + + c = Complex { -1., 1. }; + d = Complex::from_polar(sqrt(2.), 3. * M_PI / 4.); + EXPECT_APPROXIMATE(c.real(), d.real()); + EXPECT_APPROXIMATE(c.imag(), d.imag()); + EXPECT_APPROXIMATE(d.phase(), 3. * M_PI / 4.); + EXPECT_APPROXIMATE(c.magnitude(), d.magnitude()); + EXPECT_APPROXIMATE(c.magnitude(), sqrt(2.)); +#endif + EXPECT_EQ((complex_imag_unit * complex_imag_unit).real(), -1.); + EXPECT_EQ((complex_imag_unit / complex_imag_unit).real(), 1.); + + EXPECT_EQ(Complex(1., 10.) == (Complex(1., 0.) + Complex(0., 10.)), true); + EXPECT_EQ(Complex(1., 10.) != (Complex(1., 1.) + Complex(0., 10.)), true); +#ifdef AKCOMPLEX_CAN_USE_MATH_H + EXPECT_EQ(approx_eq(Complex(1), Complex(1.0000004f)), true); + EXPECT_APPROXIMATE(cexp(Complex(0., 1.) * M_PI).real(), -1.); +#endif +} + +TEST_MAIN(Complex)