complex_analysis.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_COMPLEX_FUNCTIONS
#define THEORETICA_COMPLEX_FUNCTIONS
 
#include "./complex.h"
#include "../core/real_analysis.h"
 
 
namespace theoretica {
 
 
    template<typename T>
    inline complex<T> identity(complex<T> z) {
        return z;
    }
 
 
    template<typename T>
    inline complex<T> conjugate(complex<T> z) {
        return complex<T>(z.a, -z.b);
    }
 
 
    template<typename T>
    inline complex<T> inverse(complex<T> z) {
        return conjugate(z) / z.sqr_norm();
    }
 
 
    template<typename T>
    inline complex<T> square(complex<T> z) {
        return complex<T>(
            z.Re() * z.Re() - z.Im() * z.Im(),
            2 * z.Re() * z.Im()
        );
    }
 
 
    template<typename T>
    inline complex<T> cube(complex<T> z) {
 
        // Use the algebraic identity minimizing multiplications
        // (a + bi)^3 = (a^3 - 3ab^2) + (3a^2b - b^3)i
        const T a2 = z.Re() * z.Re();
        const T b2 = z.Im() * z.Im();
 
        const T real_factor = a2 - T(3) * b2;
        const T imag_factor = T(3) * a2 - b2;
 
        return complex<T>(z.Re() * real_factor, z.Im() * imag_factor);
    }
 
 
    template<typename T>
    inline complex<T> powf(complex<T> z, real p) {
 
        if(abs(z.Re()) < MACH_EPSILON && abs(z.Im()) < MACH_EPSILON) {
 
            if(abs(p) < MACH_EPSILON) {
                TH_MATH_ERROR("powf(complex, real)", 0, MathError::ImpossibleOperation);
                return complex<T>(nan(), nan());
            }
 
            return complex<T>(0, 0);
        }
 
        const real rho = powf(z.norm(), p);
        const real theta = p * z.arg();
 
        return complex<T>(rho * th::cos(theta), rho * th::sin(theta));
    }
    
 
    template<typename T>
    inline complex<T> exp(complex<T> z) {
        return complex<T>(cos(z.Im()), sin(z.Im())) * exp(z.Re());
    }
 
 
    template<typename T>
    inline real abs(complex<T> z) {
        return z.norm();
    }
 
 
    template<typename T>
    inline complex<T> sin(complex<T> z) {
 
        const complex<T> t = z * complex<>::i();
        return (exp(t) - exp(-t)) * complex<T>(0, -0.5);
    }
 
 
    template<typename T>
    inline complex<T> cos(complex<T> z) {
 
        const complex<T> t = z * complex<>::i();
        return (exp(t) + exp(-t)) / 2.0;
    }
 
 
    template<typename T>
    inline complex<T> tan(complex<T> z) {
 
        const complex<T> t = exp(z * complex<T>(0, 2));
        return (t - complex<T>(1, 0)) / (t + complex<T>(1, 0)) * complex<T>(0, -1);
    }
 
 
    template<typename T>
    inline complex<T> sqrt(complex<T> z) {
 
        if(abs(z.a) < MACH_EPSILON && abs(z.b) < MACH_EPSILON)
            return complex<T>(0);
 
        return complex<T>(
            INVSQR2 * sqrt((z.norm() + z.Re())),
            INVSQR2 * sqrt((z.norm() - z.Re()))  * sgn(z.b));
    }
 
 
    template<typename T>
    inline complex<T> ln(complex<T> z) {
        return complex<T>(ln(z.norm()), z.arg());
    }
 
 
    template<typename T>
    inline complex<T> asin(complex<T> z) {
 
        // For real z in [-1, 1], use real asin
        if(abs(z.Im()) < MACH_EPSILON && abs(z.Re()) <= 1.0 + MACH_EPSILON) {
            T real_part = (z.Re() > 1.0) ? 1.0 : ((z.Re() < -1.0) ? -1.0 : z.Re());
            return complex<T>(asin(real_part), 0);
        }
 
        // General formula with better stability
        return ln(complex<>::i() * z + sqrt(complex<T>(1, 0) - square(z))) * complex<T>(0, -1);
    }
 
 
    template<typename T>
    inline complex<T> acos(complex<T> z) {
 
        // For real z in [-1, 1], use real acos
        if(abs(z.Im()) < MACH_EPSILON && abs(z.Re()) <= 1.0 + MACH_EPSILON) {
            const T real_part = (z.Re() > 1.0) ? 1.0 : ((z. Re() < -1.0) ? -1.0 : z.Re());
            return complex<T>(acos(real_part), 0);
        }
        
        // General formula
        return ln(z + sqrt(square(z) - complex<T>(1, 0))) * complex<T>(0, -1);
    }
 
 
    template<typename T>
    inline complex<T> atan(complex<T> z) {
 
        return ln((complex<>::i() - z) / (complex<>::i() + z)) * complex<T>(0, -0.5);
    }
 
}
 
#endif