My Dual Number implementation

dual.h

#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;
}

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.