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) {
        return z * z * z;
    }
 
 
    template<typename T>
    inline complex<T> powf(complex<T> z, real p) {
        return exp(p * ln(z));
    }
    
 
    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);
 
 
        // Left-sided sign function for continuity along the negative real axis
        const T imag_sign = (z.Im() < 0) ? T(-1.0) : T(1.0);
 
        return complex<T>(
            INVSQR2 * sqrt((z.norm() + z.Re())),
            INVSQR2 * sqrt((z.norm() - z.Re())) * imag_sign);
    }
 
 
    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