dual2_functions.h

Go to the documentation of this file.

 
 
 
#ifndef THEORETICA_DUAL2_FUNCTIONS_H
#define THEORETICA_DUAL2_FUNCTIONS_H
 
#include "./dual2.h"
#include "../core/real_analysis.h"
 
 
namespace theoretica {
 
 
    // Second derivative of the composite function:
    // (f(g))'' = f''(g) * g'' + f'(g) * g''
 
 
    inline dual2 square(dual2 x) {
        return x * x;
    }
 
 
    inline dual2 cube(dual2 x) {
        return x * x * x;
    }
 
 
    inline dual2 conjugate(dual2 x) {
        return x.conjugate();
    }
 
 
    inline dual2 pow(dual2 x, int n) {
 
        const real pow_n_2_x = pow(x.Re(), n - 2);
 
        real f = pow_n_2_x * square(x.Re());
        real df = pow_n_2_x * x.Re() * n * x.Dual1();
        real d2f = (pow_n_2_x * n * (n - 1)) * square(x.Dual1()) + df * x.Dual2();
 
        return dual2(f, df, d2f);
    }
 
 
    inline dual2 sqrt(dual2 x) {
 
        const real sqrt_x = sqrt(x.Re());
 
        if(sqrt_x == 0) {
            TH_MATH_ERROR("sqrt(dual2)", sqrt_x, MathError::DivByZero);
            return dual2(nan(), nan());
        }
 
        return dual2(
            sqrt_x,
            0.5 / sqrt_x * x.Dual1(),
            -0.25 / x.Re() / sqrt_x * square(x.Dual1()) + 0.5 / sqrt_x * x.Dual2()
        );
    }
 
 
    inline dual2 sin(dual2 x) {
 
        const real sin_x = sin(x.Re());
        const real cos_x = cos(x.Re());
 
        return dual2(
            sin_x,
            cos_x * x.Dual1(),
            -sin_x * square(x.Dual1()) + cos_x * x.Dual2()
        );
    }
 
 
    inline dual2 cos(dual2 x) {
 
        const real sin_x = sin(x.Re());
        const real cos_x = cos(x.Re());
 
        return dual2(
            cos_x,
            -sin_x * x.Dual1(),
            -cos_x * square(x.Dual1()) - sin_x * x.Dual2()
        );
    }
 
 
    inline dual2 tan(dual2 x) {
 
        const real sin_x = sin(x.Re());
        const real cos_x = cos(x.Re());
 
        if(cos_x == 0) {
            TH_MATH_ERROR("tan(dual2)", cos_x, MathError::DivByZero);
            return dual2(nan(), nan());
        }
 
        return dual2(
            tan(x.Re()),
            x.Dual1() / square(cos_x),
            2.0 * sin_x / cube(cos_x) * square(x.Dual1()) + x.Dual2() / square(cos_x)
        );
    }
 
 
    inline dual2 cot(dual2 x) {
 
        const real sin_x = sin(x.Re());
        const real cos_x = cos(x.Re());
 
        if(sin_x == 0) {
            TH_MATH_ERROR("cot(dual2)", sin_x, MathError::DivByZero);
            return dual2(nan(), nan());
        }
 
        return dual2(
            cot(x.Re()),
            -x.Dual1() / square(sin_x),
            -2.0 * cos_x / cube(sin_x) * square(x.Dual1()) - x.Dual2() / square(sin_x)
        );
    }
 
 
    inline dual2 exp(dual2 x) {
 
        const real exp_x = exp(x.Re());
 
        return dual2(
            exp_x,
            x.Dual1() * exp_x,
            square(x.Dual1()) * exp_x + x.Dual2() * exp_x
        );
    }
 
 
    inline dual2 sinh(dual2 x) {
 
        const real sinh_x = sinh(x.Re());
        const real cosh_x = cosh(x.Re());
 
        return dual2(
            sinh_x,
            cosh_x * x.Dual1(),
            cosh_x * square(x.Dual1()) + sinh_x * x.Dual2()
        );
    }
 
 
    inline dual2 cosh(dual2 x) {
 
        const real sinh_x = sinh(x.Re());
        const real cosh_x = cosh(x.Re());
 
        return dual2(
            cosh_x,
            sinh_x * x.Dual1(),
            sinh_x * square(x.Dual1()) + cosh_x * x.Dual2()
        );
    }
 
 
    inline dual2 tanh(dual2 x) {
        
        const real sinh_x = sinh(x.Re());
        const real cosh_x = cosh(x.Re());
 
        const real tanh_x = sinh_x / cosh_x;
 
        return dual2(
            tanh_x,
            x.Dual1() / square(cosh_x),
            2.0 * sinh_x / cube(cosh_x) * square(x.Dual1()) + x.Dual2() / square(cosh_x)
        );
    }
 
 
    inline dual2 ln(dual2 x) {
 
        if(x.Re() <= 0) {
            TH_MATH_ERROR("ln(dual2)", x.Re(), MathError::OutOfDomain);
            return dual2(nan(), nan());
        }
 
        return dual2(ln(x.Re()),
            x.Dual1() / x.Re(),
            -square(x.Dual1()) / square(x.Re()) + x.Dual2() / x.Re()
        );
    }
 
 
    inline dual2 log2(dual2 x) {
 
        if(x.Re() <= 0) {
            TH_MATH_ERROR("log2(dual2)", x.Re(), MathError::OutOfDomain);
            return dual2(nan(), nan());
        }
 
        return dual2(
            log2(x.Re()),
            x.Dual1() * LOG2E / x.Re(),
            -square(x.Dual1()) * LOG2E / square(x.Re()) + x.Dual2() * LOG2E / x.Re()
        );
    }
 
 
    inline dual2 log10(dual2 x) {
 
        if(x.Re() <= 0) {
            TH_MATH_ERROR("log10(dual2)", x.Re(), MathError::OutOfDomain);
            return dual2(nan(), nan());
        }
 
        return dual2(
            log10(x.Re()),
            x.Dual1() * LOG10E / x.Re(),
            -square(x.Dual1()) * LOG10E / square(x.Re()) + x.Dual2() * LOG10E / x.Re()
        );
    }
 
 
    inline dual2 abs(dual2 x) {
        return dual2(abs(x.Re()), x.Dual1() * sgn(x.Re()), 0);
    }
 
 
    inline dual2 asin(dual2 x) {
 
        if(x.Re() >= 1) {
            TH_MATH_ERROR("asin(dual2)", x.Re(), MathError::OutOfDomain);
            return dual2(nan(), nan());
        }
 
        const real sqrt_1mx2 = sqrt(1 - square(x.Re()));
 
        return dual2(
            asin(x.Re()),
            x.Dual1() / sqrt_1mx2,
            square(x.Dual1()) * x.Re() / ((1 - square(x.Re())) * sqrt_1mx2)
                + x.Dual2() / sqrt_1mx2
        );
    } 
 
 
    inline dual2 acos(dual2 x) {
 
        if(x.Re() >= 1) {
            TH_MATH_ERROR("acos(dual2)", x.Re(), MathError::OutOfDomain);
            return dual2(nan(), nan());
        }
 
        const real sqrt_1mx2 = sqrt(1 - square(x.Re()));
 
        return dual2(
            acos(x.Re()),
            -x.Dual1() / sqrt_1mx2,
            -square(x.Dual1()) * x.Re() / ((1 - square(x.Re())) * sqrt_1mx2)
                - x.Dual2() / sqrt_1mx2
        );
    }
 
 
    inline dual2 atan(dual2 x) {
 
        return dual2(
            atan(x.Re()),
            x.Dual1() / (1 + square(x.Re())),
            -square(x.Dual1()) * 2 * x.Re() / square(1 + square(x.Re()))
                + x.Dual2() / (1 + square(x.Re()))
        );
    }
 
}
 
 
#endif