/*

  Class to define basic properties of Dual numbers. The dual number w is defined as  w = x + ey
  where x,y are reals (represented as doubles), and e is an infinitesimal number with the properties e^2=0 and e!=0 
  Hence we can think of e as somewhat similar to the imaginary number i with i^2 = -1.
  
*/

#include <cmath>
#include <iostream>
using namespace std;

// This is the header defining DualNumber class.
#include "DualNumber.hpp"

/*
  Below we define a function where the type of the argument the return value is
  chosen automatically by the compiler at compile time. This allows the same function
  declaration to be used for doubles, complex numbers, and dual numbers.
*/

template <typename Number>
Number f(Number x)
/*
    f(x) = 3x^3    df/dx = 6x
*/
{
   return 3.0*x*x;
}

template <typename Number>
Number g(Number x)
{
   return sqrt(x)*exp(x*x)*x*x*x;
}


int main()
{
   // Now lets unleash the power of automatic differentiation:
   double x=0.5;
   double delta=1e-5;                          // arbitrary step length.
   double dfdx=(f(x+delta)-f(x))/delta;        // Different quotient approximation for df/dx
   double dgdx=(g(x+delta)-g(x))/delta;        // Different quotient approximation for dg/dx

   cout << "Evaluating the function with a double:\n";
   cout << "x=" << x << "\n";
   cout << "f(x)=" << f(x) << "\n";
   cout << "df/dx=" << dfdx << " (approximately)\n";
   cout << "g(x)=" << g(x) << "\n";
   cout << "dg/dx=" << dgdx << " (approximately)\n";

   cout << "\n\nDual number algebra:\n";
   DualNumber w(x,1.0);
   cout << "w=" << w << "\n";   
   cout << "f(w)=" << f(w) << "\n";   
   cout << "g(w)=" << g(w) << "\n";   

   return 0;
}