// https://github.com/autodiff/autodiff/issues/179
#define FORWARD
#ifdef FORWARD
#include <autodiff/forward/dual.hpp>
#else
#include <autodiff/reverse/var.hpp>
#endif
#include <complex>

template <typename T> T f(const T & x)
{
  // "clever" branch to avoid 0/0 at x=0
  // Breaks autodiff::reverse::var
  if(x == 0.0)
    return 1.0;
  else
    //return (sin(x)/x + x*(x-1.0);
    // Causes complex-step to require h>1e-161
    return (sin(x) + x*x*(x-1.0))/x;
}



int main(int argc, char *argv[])
{
  for( double x : {1.5, 1.0, 0.0} )
  {
    double y = f(x);

#ifdef FORWARD
    // dual number version
    autodiff::dual x_ad = x;
    autodiff::dual y_ad = f(x);
    double dydx_ad = autodiff::derivative( f<autodiff::dual>, autodiff::wrt(x_ad), autodiff::at(x_ad));
#else
    autodiff::var x_ad = x;
    autodiff::var y_ad = f(x_ad);
    autodiff::var dydx_ad = autodiff::derivatives(y_ad, wrt(x_ad))[0];
#endif

    std::complex<double> x_cs;
    x_cs.real(x);
    // https://mdolab.engin.umich.edu/wiki/guide-complex-step-derivative-approximation
    // recommends 1e-200. That's fun but a bad idea.
    x_cs.imag(1e-161);
    double dydx_cs = f(x_cs).imag()/1e-161;

    double dydx_fd = (f(x+1e-8) - f(x-1e-8))/(2*1e-8);

    printf("f(%g) = %g\n  f'_ad = %0.17f\n  f'_cs = %0.17f\n  f'_fd = %0.17f\n",
        (double)x,y,(double)dydx_ad,dydx_cs,(double)dydx_fd);
  }
}