Instantly share code, notes, and snippets.

# safiire/dual.h

Created October 11, 2016 02:08
Show Gist options
• Save safiire/b424b96481cecd596a65d35428501f01 to your computer and use it in GitHub Desktop.
My Dual Number implementation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 #pragma once #include #include #include #include "saf_math.h" //// Some more information for adding more functionality here: //// http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/other/dualNumbers/functions/ template class Dual { private: T _real = 0; T _dual = 0; public: Dual(T real) : _real(real), _dual(T(0)) {} Dual(T real, T dual) : _real(real), _dual(dual) { } //// Get real part T real() const{ return _real; } //// Get real part T dual() const{ return _dual; } //// Addition Dual operator+=(const Dual &other){ _real += other._real; _dual += other._dual; return *this; } Dual operator+(const Dual &other) const { auto copy = *this; return copy += other; } //// Subtraction Dual operator-=(const Dual &other){ _real -= other._real; _dual -= other._dual; return *this; } Dual operator-(const Dual &other) const { auto copy = *this; return copy -= other; } //// Multiplication Dual operator*=(const Dual &other){ T real = _real; T dual = _dual; _real = real * other._real; _dual = real * other._dual + dual * other._real; return *this; } Dual operator*(const Dual &other) const { auto copy = *this; return copy *= other; } //// Division Dual operator/=(const Dual &other){ if(_real == 0 && other._real == 0 && other._dual != 0){ // When both the terms are purely dual, and it's not a division by zero // Basically the ε factors cancel out and can give you a real answer b/d // Solving for x,y in (0 + bε) / (0 + dε) = (x + yε) // (0 + bε) = (x + yε) * (0 + dε) // (0 + bε) = 0 + xdε // x = b/d, y = absolutely anything because it gets annihilated by εε term. // So why not make it 0? _real = _dual / other._dual; _dual = 0; }else{ // Otherwise we divide the same way as a regular complex number // d * d.conjugate() is always Real(d)^2 (*this) *= other.conjugate(); auto reciprocal = 1 / (other._real * other._real); _real *= reciprocal; _dual *= reciprocal; } return *this; } Dual operator/(const Dual &other) const { auto copy = *this; return copy /= other; } //// Ratio between real and dual T ratio() const { return _real / _dual; } //// Conjugate Dual conjugate() const{ return Dual {_real, -_dual}; } //// Raise to an integer exponent //// TODO this can probably be more easily written //// as _real^exponent + _dual^(exponent-1)*exponent //// Since the dual part is the derivitive, to remove the loop Dual pow(size_t exponent) const { auto result = Dual(1); for(size_t i = 0; i < exponent; ++i){ result *= (*this); } return result; } //// Dual exponential //// signum(a) * ||z|| * e^(b/a), a !=0 Dual exp() const { if(_real == 0) return std::numeric_limits::quiet_NaN(); T real_exponent = std::exp(_real); return Dual {real_exponent, real_exponent * _dual}; } //// Magnitude T magnitude() const { return std::abs(_real); } //// Argument / Angle T argument() const { return _dual / _real; } }; //// Printing Dual numbers template std::ostream& operator<<(std::ostream &os, const Dual &d) { T real = d.real(); T dual = d.dual(); const char *sign; switch(saf_math::signum(d.dual())) { case -1: sign = " - "; dual *= T(-1); break; case 0: case 1: sign = " + "; } os << "(" << real << sign << dual << "ε)"; return os; }

### safiire commented May 7, 2018

This choice:

``````//  x = b/d, y = absolutely anything because it gets annihilated by εε term.
//  So why not make it 0?
``````

Is very mathematically wrong, but it's fun to play with because it forces dual numbers to have inverses most of the time.