dual.h

Go to the documentation of this file.

 
 
#ifndef THEORETICA_DUAL_H
#define THEORETICA_DUAL_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 dual {
        public:
 
            real a; // Real part
            real b; // "Dual" part
 
            dual() : a(0), b(0) {}
 
            dual(real real_part, real dual_part)
                : a(real_part), b(dual_part) {}
 
            dual(real real_part)
                : a(real_part), b(0) {}
 
            ~dual() = default;
        
            dual(const vec2& v) {
                a = v[0];
                b = v[1];
            }
 
            inline dual& operator=(const vec2& v) {
                a = v[0];
                b = v[1];
                return *this;
            }
 
            inline dual& operator=(real x) {
                a = x;
                b = 0;
                return *this;
            }
 
            inline dual& operator=(const std::array<real, 2>& v) {
                a = v[0];
                b = v[1];
                return *this;
            }
 
            inline const real& Re() const {
                return a;
            }
 
            inline real& Re() {
                return a;
            }
 
            inline friend real Re(const dual& d) {
                return d.a;
            }
 
 
            inline friend real& Re(dual& d) {
                return d.a;
            }
 
            inline const real& Dual() const {
                return b;
            }
 
            inline real& Dual() {
                return b;
            }
 
            inline friend real Dual(const dual& d) {
                return d.b;
            }
 
 
            inline friend real& Dual(dual& d) {
                return d.b;
            }
 
            inline dual conjugate() const {
                return dual(a, -b);
            }
 
            inline dual inverse() const {
 
                if(a == 0) {
                    TH_MATH_ERROR("dual::inverse", 0, DIV_BY_ZERO);
                    return dual(nan(), nan());
                }
 
                return dual(1.0 / a, -b / square(a));
            }
 
            inline dual operator+() const {
                return dual(a, b);
            }
 
            inline dual operator+(const dual& other) const {
                return dual(a + other.a, b + other.b);
            }
 
            inline dual operator+(real r) const {
                return dual(a + r, b);
            }
 
            inline dual operator-() const {
                return dual(-a, -b);
            }
 
            inline dual operator-(const dual& other) const {
                return dual(a - other.a, b - other.b);
            }
 
            inline dual operator-(real r) const {
                return dual(a - r, b);
            }
 
            inline dual operator*(const dual& other) const {
                return dual(a * other.a, a * other.b + b * other.a);
            }
 
            inline dual operator*(real r) const {
                return dual(a * r, b * r);
            }
 
            inline dual operator/(const dual& other) const {
                return dual(a / other.a,
                    (b * other.a - a * other.b) / square(other.a));
            }
 
            inline dual operator/(real r) const {
                return dual(a / r, b / r);
            }
 
 
            inline dual& operator+=(const dual& other) {
 
                a += other.a;
                b += other.b;
                return *this;
            }
 
 
            inline dual& operator+=(real r) {
 
                a += r;
                return *this;
            }
 
            inline dual& operator-=(const dual& other) {
 
                a -= other.a;
                b -= other.b;
                return *this;
            }
 
            inline dual& operator-=(real r) {
 
                a -= r;
                return *this;
            }
 
            inline dual& operator*=(const dual& other) {
 
                a = (a * other.a);
                b = (a * other.b) + (b * other.a);
                return *this;
            }
 
            inline dual& operator*=(real r) {
 
                a *= r;
                b *= r;
                return *this;
            }
 
            inline dual& operator/=(const dual& other) {
 
                a = (a / other.a);
                b = (b * other.a - a * other.b) / square(other.a);
                return *this;
            }
 
            inline dual& operator/=(real r) {
 
                if(r == 0) {
                    TH_MATH_ERROR("dual::operator/=", 0, DIV_BY_ZERO);
                    a = nan();
                    b = nan();
                    return *this;
                }
 
                a /= r;
                b /= r;
 
                return *this;
            }
 
 
            inline bool operator==(const dual& other) {
                return (a == other.a) && (b == other.b);
            }
 
 
            inline vec2 to_vec() const {
                vec2 res;
                res[0] = a;
                res[1] = b;
                return res;
            }
 
            inline void from_vec(const vec2& v) {
                a = v[0];
                b = v[1];
            }
 
 
            inline mat2 to_mat() const {
 
                mat2 m;
                m.at(0, 0) = a;
                m.at(1, 0) = 0;
                m.at(0, 1) = b;
                m.at(1, 1) = a;
                return m;
            }
 
 
            // Friend operators to enable equations of the form
            // (real) op. (dual)
 
            inline friend dual operator+(real a, const dual& d) {
                return d + a;
            }
 
            inline friend dual operator-(real a, const dual& d) {
                return -d + a;
            }
 
            inline friend dual operator*(real a, const dual& d) {
                return d * a;
            }
 
            inline friend dual operator/(real a, const dual& d) {
                return dual(a, 0) / d;
            }
 
 
#ifndef THEORETICA_NO_PRINT
 
            inline std::string to_string(const std::string& epsilon = "e") const {
 
                std::stringstream res;
 
                res << a;
                res << (b >= 0 ? " + " : " - ");
                res << abs(b) << epsilon;
 
                return res.str();
            }
 
 
            inline operator std::string() {
                return to_string();
            }
 
 
            inline friend std::ostream& operator<<(std::ostream& out, const dual& obj) {
                return out << obj.to_string();
            }
 
#endif
 
    };
 
}
 
 
#endif