complex.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_COMPLEX_H
#define THEORETICA_COMPLEX_H
 
#ifndef THEORETICA_NO_PRINT
#include <sstream>
#include <ostream>
#endif
 
#include <array>
 
#include "../core/error.h"
#include "../core/real_analysis.h"
 
 
namespace theoretica {
 
 
    template<typename Type = real>
    class complex {
        public:
 
            Type a;
 
            Type b;
 
 
            complex() : a(0), b(0) {}
 
 
            complex(Type real_part, Type imag_part)
                : a(real_part), b(imag_part) {}
 
 
            complex(Type real_part) : a(real_part), b(0) {}
 
 
            complex(const complex& z) : a(z.a), b(z.b) {}
 
 
            inline complex& operator=(const complex& z) {
                a = z.a;
                b = z.b;
                return *this;
            }
 
 
            template<typename T = Type>
            inline complex& operator=(const std::array<T, 2>& v) {
                a = v[0];
                b = v[1];
                return *this;
            }
 
 
            inline Type Re() const {
                return a;
            }
 
 
            inline Type& Re() {
                return a;
            }
 
 
            inline friend Type Re(const complex& z) {
                return z.a;
            }
 
 
            inline friend Type& Re(complex& z) {
                return z.a;
            }
 
 
            inline Type Im() const {
                return b;
            }
 
 
            inline Type& Im() {
                return b;
            }
 
 
            inline friend Type Im(const complex& z) {
                return z.b;
            }
 
 
            inline friend Type& Im(complex& z) {
                return z.b;
            }
 
 
            inline complex conjugate() const {
                return complex(a, -b);
            }
 
 
            inline Type sqr_norm() const {
                return a * a + b * b;
            }
 
 
            inline Type norm() const {
                return sqrt(sqr_norm());
            }
 
 
            inline complex inverse() const {
 
                const Type n = sqr_norm();
 
                if(n < MACH_EPSILON) {
                    TH_MATH_ERROR("complex::inverse", n, DIV_BY_ZERO);
                    return complex(nan(), nan());
                }
 
                return conjugate() / n;
            }
 
 
            inline complex& invert() {
 
                const Type n = sqr_norm();
 
                if(n < MACH_EPSILON) {
                    TH_MATH_ERROR("complex::invert", n, DIV_BY_ZERO);
                    a = (Type) nan();
                    b = (Type) nan();
                    return *this;
                }
 
                a = a / n;
                b = -b / n;
 
                return *this;
            }
 
 
            inline Type arg() const {
 
                // Real number case
                if(abs(b) < MACH_EPSILON) {
                    return (a >= 0) ? ((Type) 0) : ((Type) PI);
                }
 
                // Pure imaginary number case
                if(abs(a) < MACH_EPSILON) {
                    return (b >= 0) ? ((Type) PI / 2.0) : ((Type) -PI / 2.0);
                }
 
                // Use the 2-parameter arctangent in the general case
                return atan2(b, a);
            }
 
 
            inline complex operator+() const {
                return complex(a, b);
            }
 
 
            inline complex operator+(const complex& z) const {
                return complex(a + z.a, b + z.b);
            }
 
 
            inline complex operator-() const {
                return complex(-a, -b);
            }
 
 
            inline complex operator-(const complex& z) const {
                return complex(a - z.a, b - z.b);
            }
 
 
            inline complex operator*(const complex& z) const {
                return complex(
                    a * z.a - b * z.b,
                    a * z.b + b * z.a
                );
            }
 
 
            inline complex operator/(const complex& z) const {
                return operator*(z.inverse());
            }
 
 
            inline complex& operator+=(const complex& z) {
                return (*this = complex(a + z.a, b + z.b));
            }
 
 
            inline complex& operator-=(const complex& z) {
                return (*this = complex(a - z.a, b - z.b));
            }
 
 
            inline complex& operator*=(const complex& z) {
                return (*this = complex(
                    a * z.a - b * z.b,
                    a * z.b + b * z.a
                ));
            }
 
 
            inline complex& operator/=(const complex& z) {
                return (*this = operator*(z.inverse()));
            }
 
 
            inline complex operator+(Type k) const {
                return complex(a + k, b);
            }
 
 
            inline complex operator-(Type k) const {
                return complex(a - k, b);
            }
 
 
            inline complex operator*(Type k) const {
                return complex(a * k, b * k);
            }
 
 
            inline complex operator/(Type k) const {
 
                if(abs(k) < MACH_EPSILON) {
                    TH_MATH_ERROR("complex::operator/", k, DIV_BY_ZERO);
                    return complex(nan(), nan());
                }
 
                return complex(a / k, b / k);
            }
 
 
            inline complex& operator+=(Type k) {
                return (*this = complex(a + k, b));
            }
 
 
            inline complex& operator-=(Type k) {
                return (*this = complex(a - k, b));
            }
 
 
            inline complex& operator*=(Type k) {
                return (*this = complex(a * k, b * k));
            }
 
 
            inline complex& operator/=(Type k) {
 
                if(abs(k) < MACH_EPSILON) {
                    TH_MATH_ERROR("complex::operator/=", k, DIV_BY_ZERO);
                    return complex(nan(), nan());
                }
 
                return (*this = complex(a / k, b / k));
            }
 
 
            inline bool operator==(const complex& z) const {
                return (a == z.a) && (b == z.b);
            }
 
 
            inline bool operator!=(const complex& z) const {
                return !(*this == z);
            }
 
 
            inline static complex rotor(Type rad) {
                return complex(cos(rad), sin(rad));
            }
 
 
            inline static complex i() {
                return complex(0, 1);
            }
 
 
            // Friend operators to enable equations of the form
            // (real) op. (complex)
 
            inline friend complex operator+(Type r, const complex& z) {
                return z + r;
            }
 
            inline friend complex operator-(Type r, const complex& z) {
                return -z + r;
            }
 
            inline friend complex operator*(Type r, const complex& z) {
                return z * r;
            }
 
            inline friend complex operator/(Type r, const complex& z) {
                return complex(r, 0) / z;
            }
 
 
            inline operator real () {
 
                if (abs(b) >= MACH_EPSILON) {
                    return nan();
                }
 
                return a;
            }
 
 
#ifndef THEORETICA_NO_PRINT
 
            inline std::string to_string() const {
 
                std::stringstream res;
 
                // TO-DO Check whether Type has comparison
                // operators (bicomplex with Type=complex !)
 
                res << a;
                res << (b >= 0 ? " + " : " - ");
                
                if(abs(b) != 1)
                    res << abs(b);
 
                res << "i";
 
                return res.str();
            }
 
 
            inline operator std::string() {
                return to_string();
            }
 
 
            inline friend std::ostream& operator<<(std::ostream& out, const complex& obj) {
                return out << obj.to_string();
            }
 
#endif
 
    };
 
}
 
#endif