dual2.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_DUAL2_H
#define THEORETICA_DUAL2_H
 
#ifndef THEORETICA_NO_PRINT
#include <sstream>
#include <ostream>
#endif
 
#include "../core/error.h"
#include "../core/constants.h"
#include "../algebra/algebra_types.h"
 
 
namespace theoretica {
 
 
    class dual2 {
        public:
 
            real a; // Real part
            real b; // First order dual part
            real c; // Second order dual part
 
            dual2() : a(0), b(0), c(0) {}
 
            dual2(real real_part, real dual1_part, real dual2_part)
                : a(real_part), b(dual1_part), c(dual2_part) {}
 
            dual2(real real_part, real dual1_part)
                : a(real_part), b(dual1_part), c(0) {}
 
            dual2(real real_part)
                : a(real_part), b(0), c(0) {}
 
            ~dual2() = default;
        
            dual2(const vec3& v) {
                a = v[0];
                b = v[1];
                c = v[2];
            }
 
            inline dual2& operator=(const vec3& v) {
                a = v[0];
                b = v[1];
                c = v[2];
                return *this;
            }
 
            inline dual2& operator=(real x) {
                a = x;
                b = 0;
                c = 0;
                return *this;
            }
 
            inline dual2& operator=(const std::array<real, 3>& v) {
                a = v[0];
                b = v[1];
                c = v[2];
                return *this;
            }
 
            inline real Re() const {
                return a;
            }
 
            inline real Dual1() const {
                return b;
            }
 
            inline real Dual2() const {
                return c;
            }
 
            inline dual2 conjugate() const {
                return dual2(a, -b, -c);
            }
 
            inline dual2 inverse() const {
 
                if(a == 0) {
                    TH_MATH_ERROR("dual2::inverse", 0, MathError::DivByZero);
                    return dual2(nan(), nan(), nan());
                }
 
                return dual2(1.0 / a, -b / square(a), 2.0 * square(b) / cube(a) - c / square(a));
            }
 
            inline dual2 operator+() const {
                return dual2(a, b, c);
            }
 
            inline dual2 operator+(const dual2& other) const {
                return dual2(a + other.a, b + other.b, c + other.c);
            }
 
            inline dual2 operator+(real r) const {
                return dual2(a + r, b, c);
            }
 
            inline dual2 operator-() const {
                return dual2(-a, -b, -c);
            }
 
            inline dual2 operator-(const dual2& other) const {
                return dual2(a - other.a, b - other.b, c - other.c);
            }
 
            inline dual2 operator-(real r) const {
                return dual2(a - r, b, c);
            }
 
            inline dual2 operator*(const dual2& other) const {
                return dual2(a * other.a,
                    a * other.b + b * other.a,
                    a * other.c + 2 * b * other.b + c * other.a);
            }
 
            inline dual2 operator*(real r) const {
                return dual2(a * r, b * r, c * r);
            }
 
            inline dual2 operator/(const dual2& other) const {
                return operator*(other.inverse());
            }
 
            inline dual2 operator/(real r) const {
 
                if(r == 0) {
                    TH_MATH_ERROR("dual2::operator/", r, MathError::DivByZero);
                    return dual2(nan(), nan(), nan());
                }
 
                return dual2(a / r, b / r, c / r);
            }
 
 
            inline dual2& operator+=(const dual2& other) {
 
                a += other.a;
                b += other.b;
                c += other.c;
                return *this;
            }
 
 
            inline dual2& operator+=(real r) {
 
                a += r;
                return *this;
            }
 
            inline dual2& operator-=(const dual2& other) {
 
                a -= other.a;
                b -= other.b;
                c -= other.c;
                return *this;
            }
 
            inline dual2& operator-=(real r) {
 
                a -= r;
                return *this;
            }
 
            inline dual2& operator*=(const dual2& other) {
 
                const real a_old = a;
                const real b_old = b;
                const real c_old = c;
                
                a = (a_old * other.a);
                b = (a_old * other.b) + (b_old * other.a);
                c = (a_old * other.c) + (2 * b_old * other.b) + (c_old * other.a);
                return *this;
            }
 
            inline dual2& operator*=(real r) {
                a *= r;
                b *= r;
                c *= r;
                return *this;
            }
 
            inline dual2& operator/=(const dual2& other) {
                *this = operator*(other.inverse());
                return *this;
            }
 
            inline dual2& operator/=(real r) {
 
                if(r == 0) {
                    TH_MATH_ERROR("dual::operator/=", 0, MathError::DivByZero);
                    a = nan();
                    b = nan();
                    c = nan();
                    return *this;
                }
 
                a /= r;
                b /= r;
                c /= r;
 
                return *this;
            }
 
 
            inline bool operator==(const dual2& other) const {
                return (a == other.a) && (b == other.b) && (c == other.c);
            }
 
 
            inline vec3 to_vec() const {
                vec3 res;
                res[0] = a;
                res[1] = b;
                res[2] = c;
                return res;
            }
 
            inline void from_vec(const vec3& v) {
                a = v[0];
                b = v[1];
                c = v[2];
            }
 
 
            // Friend operators to enable equations of the form
            // (real) op. (dual2)
 
            inline friend dual2 operator+(real a, const dual2& d) {
                return d + a;
            }
 
            inline friend dual2 operator-(real a, const dual2& d) {
                return -d + a;
            }
 
            inline friend dual2 operator*(real a, const dual2& d) {
                return d * a;
            }
 
            inline friend dual2 operator/(real a, const dual2& d) {
                return dual2(a, 0, 0) / d;
            }
 
 
#ifndef THEORETICA_NO_PRINT
 
            inline std::string to_string(const std::string& epsilon1 = "ε",
                const std::string& epsilon2 = "ε²") const {
 
                std::stringstream res;
 
                res << a;
                res << (b >= 0 ? " + " : " - ");
                res << abs(b);
 
                res << epsilon1;
 
                res << (c >= 0 ? " + " : " - ");
                res << abs(c);
 
                res << epsilon2;
 
                return res.str();
            }
 
 
            inline operator std::string() {
                return to_string();
            }
 
 
            inline friend std::ostream& operator<<(std::ostream& out, const dual2& obj) {
                return out << obj.to_string();
            }
 
#endif
 
    };
 
}
 
 
#endif