Skip to content

Instantly share code, notes, and snippets.

@safiire

safiire/dual.h

Created Oct 11, 2016
Embed
What would you like to do?
My Dual Number implementation
#pragma once
#include <iostream>
#include <cmath>
#include <limits>
#include "saf_math.h"
//// Some more information for adding more functionality here:
//// http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/other/dualNumbers/functions/
template <typename T>
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<T> operator+=(const Dual<T> &other){
_real += other._real;
_dual += other._dual;
return *this;
}
Dual<T> operator+(const Dual<T> &other) const {
auto copy = *this;
return copy += other;
}
//// Subtraction
Dual<T> operator-=(const Dual<T> &other){
_real -= other._real;
_dual -= other._dual;
return *this;
}
Dual<T> operator-(const Dual<T> &other) const {
auto copy = *this;
return copy -= other;
}
//// Multiplication
Dual<T> operator*=(const Dual<T> &other){
T real = _real;
T dual = _dual;
_real = real * other._real;
_dual = real * other._dual + dual * other._real;
return *this;
}
Dual<T> operator*(const Dual<T> &other) const {
auto copy = *this;
return copy *= other;
}
//// Division
Dual<T> operator/=(const Dual<T> &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<T> operator/(const Dual<T> &other) const {
auto copy = *this;
return copy /= other;
}
//// Ratio between real and dual
T ratio() const {
return _real / _dual;
}
//// Conjugate
Dual<T> conjugate() const{
return Dual<T> {_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<T> pow(size_t exponent) const {
auto result = Dual<T>(1);
for(size_t i = 0; i < exponent; ++i){
result *= (*this);
}
return result;
}
//// Dual exponential
//// signum(a) * ||z|| * e^(b/a), a !=0
Dual<T> exp() const {
if(_real == 0)
return std::numeric_limits<T>::quiet_NaN();
T real_exponent = std::exp(_real);
return Dual<T> {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 <typename T>
std::ostream& operator<<(std::ostream &os, const Dual<T> &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

This comment has been minimized.

Copy link
Owner Author

@safiire 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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.