multidual_functions.h

Go to the documentation of this file.

 
 
 
#ifndef THEORETICA_MULTIDUAL_FUNCTIONS_H
#define THEORETICA_MULTIDUAL_FUNCTIONS_H
 
#include "./multidual.h"
#include "../core/real_analysis.h"
 
 
namespace theoretica {
 
 
    template<unsigned int N>
    multidual<N> square(multidual<N> x) {
        return x * x;
    }
 
 
    template<unsigned int N>
    multidual<N> cube(multidual<N> x) {
        return x * x * x;
    }
 
 
    template<unsigned int N>
    multidual<N> conjugate(multidual<N> x) {
        return x.conjugate();
    }
 
 
    template<unsigned int N>
    multidual<N> pow(multidual<N> x, int n) {
 
        const real pow_n_1_x = pow(x.Re(), n - 1);
        return multidual<N>(pow_n_1_x * x.Re(), x.Dual() * pow_n_1_x * n);
    }
 
 
    template<unsigned int N>
    multidual<N> sqrt(multidual<N> x) {
 
        real sqrt_x = sqrt(x.Re());
 
        if(sqrt_x == 0) {
            TH_MATH_ERROR("sqrt(multidual)", sqrt_x, MathError::DivByZero);
            return multidual<N>(nan(), vec<real, N>(nan()));
        }
 
        return multidual<N>(sqrt_x, x.Dual() * 0.5 / sqrt_x);
    }
 
 
    template<unsigned int N>
    multidual<N> sin(multidual<N> x) {
        return multidual<N>(sin(x.Re()), x.Dual() * cos(x.Re()));
    }
 
 
    template<unsigned int N>
    multidual<N> cos(multidual<N> x) {
        return multidual<N>(cos(x.Re()), x.Dual() * -sin(x.Re()));
    }
 
 
    template<unsigned int N>
    multidual<N> tan(multidual<N> x) {
 
        real cos_x = cos(x.Re());
 
        if(cos_x == 0) {
            TH_MATH_ERROR("tan(multidual)", cos_x, MathError::DivByZero);
            return multidual<N>(nan(), vec<real, N>(nan()));
        }
 
        return multidual<N>(tan(x.Re()), x.Dual() / square(cos_x));
    }
 
 
    template<unsigned int N>
    multidual<N> cot(multidual<N> x) {
 
        real sin_x = sin(x.Re());
 
        if(sin_x == 0) {
            TH_MATH_ERROR("cot(multidual)", sin_x, MathError::DivByZero);
            return multidual<N>(nan(), vec<real, N>(nan()));
        }
 
        return multidual<N>(cot(x.Re()), x.Dual() * (-1 / square(sin_x)));
    }
 
 
    template<unsigned int N>
    multidual<N> exp(multidual<N> x) {
        real exp_x = exp(x.Re());
        return multidual<N>(exp_x, x.Dual() * exp_x);
    }
 
 
    template<unsigned int N>
    multidual<N> ln(multidual<N> x) {
 
        if(x.Re() <= 0) {
            TH_MATH_ERROR("ln(multidual)", x.Re(), MathError::OutOfDomain);
            return multidual<N>(nan(), vec<real, N>(nan()));
        }
 
        return multidual<N>(ln(x.Re()), x.Dual() / x.Re());
    }
 
 
    template<unsigned int N>
    multidual<N> log2(multidual<N> x) {
 
        if(x.Re() <= 0) {
            TH_MATH_ERROR("log2(multidual)", x.Re(), MathError::OutOfDomain);
            return multidual<N>(nan(), vec<real, N>(nan()));
        }
 
        return multidual<N>(log2(x.Re()), x.Dual() * LOG2E / x.Re());
    }
 
 
    template<unsigned int N>
    multidual<N> log10(multidual<N> x) {
 
        if(x.Re() <= 0) {
            TH_MATH_ERROR("log10(multidual)", x.Re(), MathError::OutOfDomain);
            return multidual<N>(nan(), vec<real, N>(nan()));
        }
 
        return multidual<N>(log10(x.Re()), x.Dual() * LOG10E / x.Re());
    }
 
 
    template<unsigned int N>
    multidual<N> abs(multidual<N> x) {
        return multidual<N>(abs(x.Re()), x.Dual() * sgn(x.Re()));
    }
 
 
    template<unsigned int N>
    multidual<N> asin(multidual<N> x) {
 
        if(x.Re() >= 1) {
            TH_MATH_ERROR("asin(multidual)", x.Re(), MathError::OutOfDomain);
            return multidual<N>(nan(), vec<real, N>(nan()));
        }
 
        return multidual<N>(asin(x.Re()), x.Dual() / sqrt(1 - square(x.Re())));
    }
 
 
    template<unsigned int N>
    multidual<N> acos(multidual<N> x) {
 
        if(x.Re() >= 1) {
            TH_MATH_ERROR("acos(multidual)", x.Re(), MathError::OutOfDomain);
            return multidual<N>(nan(), vec<real, N>(nan()));
        }
 
        return multidual<N>(acos(x.Re()), x.Dual() * (-1 / sqrt(1 - square(x.Re()))));
    }
 
    template<unsigned int N>
    multidual<N> atan(multidual<N> x) {
        return multidual<N>(atan(x.Re()), x.Dual() / (1 + square(x.Re())));
    }
 
 
    template<unsigned int N>
    multidual<N> sinh(multidual<N> x) {
 
        real exp_x = exp(x.Re());
        return multidual<N>((exp_x - 1.0 / exp_x) / 2.0, x.Dual() * (exp_x + 1.0 / exp_x) / 2.0);
    }
 
 
    template<unsigned int N>
    multidual<N> cosh(multidual<N> x) {
 
        real exp_x = exp(x.Re());
        return multidual<N>((exp_x + 1.0 / exp_x) / 2.0, x.Dual() * (exp_x - 1.0 / exp_x) / 2.0);
    }
 
 
    template<unsigned int N>
    multidual<N> tanh(multidual<N> x) {
        
        const real exp_2x = exp(-2.0 * abs(x.Re()));
 
        real t;
        if (x.Re() >= 0.0)
            t = (1.0 - exp_2x) / (1.0 + exp_2x);
        else
            t = (exp_2x - 1.0) / (1.0 + exp_2x);
 
        const real dt = (4.0 * exp_2x) / square(1.0 + exp_2x);
 
        return multidual<N>(t, x.Dual() * dt);
    }
 
}
 
 
#endif