#ifndef DEF_AUTOMATIC_DIFFERENTIATION
#define DEF_AUTOMATIC_DIFFERENTIATION

#include <cmath>
#include <ostream>

/// Implementation of dual numbers for automatic differentiation
/// reference : http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.89.7749&rep=rep1&type=pdf
template<typename Scalar>
class Dual
{
	public:
		Dual(const Scalar & _a, const Scalar & _b = Scalar(0.0))
			: a(_a),
			  b(_b)
			  {}
		
		static Dual d() {
			return Dual(Scalar(0.), Scalar(1.));
		}
		
		Dual & operator+=(const Dual & x) {
			a += x.a;
			b += x.b;
			return *this;
		}
		
		Dual & operator-=(const Dual & x) {
			a -= x.a;
			b -= x.b;
			return *this;
		}
		
		Dual & operator*=(const Dual & x) {
			b = a*x.b + b*x.a;
			a *= x.a;
			return *this;
		}
		
		Dual & operator/=(const Dual & x) {
			b = (x.a*b - a*x.b)/(x.a*x.a);
			a /= x.a;
			return *this;
		}
		
		Dual & operator++() { // ++x
			return ((*this) += Scalar(1.));
		}
		
		Dual & operator--() { // --x
			return ((*this) -= Scalar(1.));
		}
		
		Dual operator++(int) { // x++
			Dual copy = *this;
			(*this) += Scalar(1.);
			return copy;
		}
		
		Dual operator--(int) { // x--
			Dual copy = *this;
			(*this) -= Scalar(1.);
			return copy;
		}
		
		Dual operator+(const Dual & x) const {
			Dual res(*this);
			return (res += x);
		}
		
		Dual operator+(void) const { // +x
			return (*this);
		}
		
		Dual operator-(const Dual & x) const {
			Dual res(*this);
			return (res -= x);
		}
		
		Dual operator-(void) const { // -x
			return Dual(-a, -b);
		}
		
		Dual operator*(const Dual & x) const {
			Dual res(*this);
			return (res *= x);
		}
		
		Dual operator/(const Dual & x) const {
			Dual res(*this);
			return (res /= x);
		}
		
		bool operator==(const Dual & x) const {
			return (a == x.a);
		}
		
		bool operator!=(const Dual & x) const {
			return (a != x.a);
		}
		
		bool operator<(const Dual & x) const {
			return (a < x.a);
		}
		
		bool operator<=(const Dual & x) const {
			return (a <= x.a);
		}
		
		bool operator>(const Dual & x) const {
			return (a > x.a);
		}
		
		bool operator>=(const Dual & x) const {
			return (a >= x.a);
		}
		
		Scalar a;	/// Real part
		Scalar b;	/// Infinitesimal part
};

template<typename A, typename B, int N>
Dual<B, N> operator+(A const& v, Dual<B, N> const& x) {
	return (Dual<B, N>(v) + x);
}
template<typename A, typename B, int N>
Dual<B, N> operator-(A const& v, Dual<B, N> const& x) {
	return (Dual<B, N>(v) - x);
}
template<typename A, typename B, int N>
Dual<B, N> operator*(A const& v, Dual<B, N> const& x) {
	return (Dual<B, N>(v) * x);
}
template<typename A, typename B, int N>
Dual<B, N> operator/(A const& v, Dual<B, N> const& x) {
	return (Dual<B, N>(v) / x);
}

// Basic mathematical functions for Scalar numbers

// Trigonometric functions
template<typename Scalar> Scalar sec(const Scalar & x) {
	return Scalar(1.)/cos(x);
}

template<typename Scalar> Scalar cot(const Scalar & x) {
	return cos(x)/sin(x);
}

template<typename Scalar> Scalar csc(const Scalar & x) {
	return Scalar(1.)/sin(x);
}

// Inverse trigonometric functions
template<typename Scalar> Scalar asec(const Scalar & x) {
	return acos(Scalar(1.)/x);
}

template<typename Scalar> Scalar acot(const Scalar & x) {
	return atan(Scalar(1.)/x);
}

template<typename Scalar> Scalar acsc(const Scalar & x) {
	return asin(Scalar(1.)/x);
}

// Hyperbolic trigonometric functions
template<typename Scalar> Scalar sech(const Scalar & x) {
	return Scalar(1.)/cosh(x);
}

template<typename Scalar> Scalar coth(const Scalar & x) {
	return cosh(x)/sinh(x);
}

template<typename Scalar> Scalar csch(const Scalar & x) {
	return Scalar(1.)/sinh(x);
}

// Inverse hyperbolic trigonometric functions
template<typename Scalar> Scalar asech(const Scalar & x) {
	return log((Scalar(1.) + sqrt(Scalar(1.) - x*x))/x);
}

template<typename Scalar> Scalar acoth(const Scalar & x) {
	return Scalar(0.5)*log((x + Scalar(1.))/(x - Scalar(1.)));
}

template<typename Scalar> Scalar acsch(const Scalar & x) {
	return (x >= Scalar(0.)) ? log((Scalar(1.) + sqrt(Scalar(1.) + x*x))/x) : log((Scalar(1.) - sqrt(Scalar(1.) + x*x))/x);
}

// Other functions
template<typename Scalar> Scalar exp10(const Scalar & x) {
	return exp(x*log(Scalar(10.)));
}

template<typename Scalar> Scalar sign(const Scalar & x) {
	return (x >= Scalar(0.)) ? ((x > Scalar(0.)) ? Scalar(1.) : Scalar(0.)) : Scalar(-1.);
}

template<typename Scalar> Scalar heaviside(const Scalar & x) {
	return Scalar(x >= Scalar(0.));
}

template<typename Scalar> Scalar abs(const Scalar & x) {
	return (x >= Scalar(0.)) ? x : -x;
}

// Basic mathematical functions for Dual numbers
// f(a + b*d) = f(a) + b*f'(a)*d

// Trigonometric functions
template<typename Scalar> Dual<Scalar> cos(const Dual<Scalar> & x) {
	return Dual<Scalar>(cos(x.a), -x.b*sin(x.a));
}

template<typename Scalar> Dual<Scalar> sin(const Dual<Scalar> & x) {
	return Dual<Scalar>(sin(x.a), x.b*cos(x.a));
}

template<typename Scalar> Dual<Scalar> tan(const Dual<Scalar> & x) {
	return Dual<Scalar>(tan(x.a), x.b*sec(x.a)*sec(x.a));
}

template<typename Scalar> Dual<Scalar> sec(const Dual<Scalar> & x) {
	return Dual<Scalar>(sec(x.a), x.b*sec(x.a)*tan(x.a));
}

template<typename Scalar> Dual<Scalar> cot(const Dual<Scalar> & x) {
	return Dual<Scalar>(cot(x.a), x.b*(-csc(x.a)*csc(x.a)));
}

template<typename Scalar> Dual<Scalar> csc(const Dual<Scalar> & x) {
	return Dual<Scalar>(csc(x.a), x.b*(-cot(x.a)*csc(x.a)));
}

// Inverse trigonometric functions
template<typename Scalar> Dual<Scalar> acos(const Dual<Scalar> & x) {
	return Dual<Scalar>(acos(x.a), x.b*(-Scalar(1.)/sqrt(Scalar(1.)-x.a*x.a)));
}

template<typename Scalar> Dual<Scalar> asin(const Dual<Scalar> & x) {
	return Dual<Scalar>(asin(x.a), x.b*(Scalar(1.)/sqrt(Scalar(1.)-x.a*x.a)));
}

template<typename Scalar> Dual<Scalar> atan(const Dual<Scalar> & x) {
	return Dual<Scalar>(atan(x.a), x.b*(Scalar(1.)/(x.a*x.a + Scalar(1.))));
}

template<typename Scalar> Dual<Scalar> asec(const Dual<Scalar> & x) {
	return Dual<Scalar>(asec(x.a), x.b*(Scalar(1.)/(sqrt(Scalar(1.)-Scalar(1.)/(x.a*x.a))*(x.a*x.a))));
}

template<typename Scalar> Dual<Scalar> acot(const Dual<Scalar> & x) {
	return Dual<Scalar>(acot(x.a), x.b*(-Scalar(1.)/((x.a*x.a)+Scalar(1.))));
}

template<typename Scalar> Dual<Scalar> acsc(const Dual<Scalar> & x) {
	return Dual<Scalar>(acsc(x.a), x.b*(-Scalar(1.)/(sqrt(Scalar(1.)-Scalar(1.)/(x.a*x.a))*(x.a*x.a))));
}

// Hyperbolic trigonometric functions
template<typename Scalar> Dual<Scalar> cosh(const Dual<Scalar> & x) {
	return Dual<Scalar>(cosh(x.a), x.b*sinh(x.a));
}

template<typename Scalar> Dual<Scalar> sinh(const Dual<Scalar> & x) {
	return Dual<Scalar>(sinh(x.a), x.b*cosh(x.a));
}

template<typename Scalar> Dual<Scalar> tanh(const Dual<Scalar> & x) {
	return Dual<Scalar>(tanh(x.a), x.b*sech(x.a)*sech(x.a));
}

template<typename Scalar> Dual<Scalar> sech(const Dual<Scalar> & x) {
	return Dual<Scalar>(sech(x.a), x.b*(-sech(x.a)*tanh(x.a)));
}

template<typename Scalar> Dual<Scalar> coth(const Dual<Scalar> & x) {
	return Dual<Scalar>(coth(x.a), x.b*(-csch(x.a)*csch(x.a)));
}

template<typename Scalar> Dual<Scalar> csch(const Dual<Scalar> & x) {
	return Dual<Scalar>(csch(x.a), x.b*(-coth(x.a)*csch(x.a)));
}

// Inverse hyperbolic trigonometric functions
template<typename Scalar> Dual<Scalar> acosh(const Dual<Scalar> & x) {
	return Dual<Scalar>(acosh(x.a), x.b*(Scalar(1.)/sqrt((x.a*x.a)-Scalar(1.))));
}

template<typename Scalar> Dual<Scalar> asinh(const Dual<Scalar> & x) {
	return Dual<Scalar>(asinh(x.a), x.b*(Scalar(1.)/sqrt((x.a*x.a)+Scalar(1.))));
}

template<typename Scalar> Dual<Scalar> atanh(const Dual<Scalar> & x) {
	return Dual<Scalar>(atanh(x.a), x.b*(Scalar(1.)/(Scalar(1.)-(x.a*x.a))));
}

template<typename Scalar> Dual<Scalar> asech(const Dual<Scalar> & x) {
	return Dual<Scalar>(asech(x.a), x.b*(Scalar(-1.)/(sqrt(Scalar(1.)/(x.a*x.a)-Scalar(1.))*(x.a*x.a))));
}

template<typename Scalar> Dual<Scalar> acoth(const Dual<Scalar> & x) {
	return Dual<Scalar>(acoth(x.a), x.b*(-Scalar(1.)/((x.a*x.a)-Scalar(1.))));
}

template<typename Scalar> Dual<Scalar> acsch(const Dual<Scalar> & x) {
	return Dual<Scalar>(acsch(x.a), x.b*(-Scalar(1.)/(sqrt(Scalar(1.)/(x.a*x.a)+Scalar(1.))*(x.a*x.a))));
}

// Exponential functions
template<typename Scalar> Dual<Scalar> exp(const Dual<Scalar> & x) {
	return Dual<Scalar>(exp(x.a), x.b*exp(x.a));
}

template<typename Scalar> Dual<Scalar> log(const Dual<Scalar> & x) {
	return Dual<Scalar>(log(x.a), x.b/x.a);
}

template<typename Scalar> Dual<Scalar> exp10(const Dual<Scalar> & x) {
	return Dual<Scalar>(exp10(x.a), x.b*(log(Scalar(10.))*exp10(x.a)));
}

template<typename Scalar> Dual<Scalar> log10(const Dual<Scalar> & x) {
	return Dual<Scalar>(log10(x.a), x.b/(log(Scalar(10.))*x.a));
}

template<typename Scalar> Dual<Scalar> exp2(const Dual<Scalar> & x) {
	return Dual<Scalar>(exp2(x.a), x.b*(log(Scalar(2.))*exp2(x.a)));
}

template<typename Scalar> Dual<Scalar> log2(const Dual<Scalar> & x) {
	return Dual<Scalar>(log2(x.a), x.b/(log(Scalar(2.))*x.a));
}

template<typename Scalar> Dual<Scalar> pow(const Dual<Scalar> & x, const Dual<Scalar> & n) {
	return exp(n*log(x));
}

// Other functions
template<typename Scalar> Dual<Scalar> sqrt(const Dual<Scalar> & x) {
	return Dual<Scalar>(sqrt(x.a), x.b/(Scalar(2.)*sqrt(x.a)));
}

template<typename Scalar> Dual<Scalar> sign(const Dual<Scalar> & x) {
	return Dual<Scalar>(sign(x.a), Scalar(0.));
}

template<typename Scalar> Dual<Scalar> abs(const Dual<Scalar> & x) {
	return Dual<Scalar>(abs(x.a), x.b*sign(x.a));
}

template<typename Scalar> Dual<Scalar> fabs(const Dual<Scalar> & x) {
	return Dual<Scalar>(fabs(x.a), x.b*sign(x.a));
}

template<typename Scalar> Dual<Scalar> heaviside(const Dual<Scalar> & x) {
	return Dual<Scalar>(heaviside(x.a), Scalar(0.));
}

template<typename Scalar> Dual<Scalar> floor(const Dual<Scalar> & x) {
	return Dual<Scalar>(floor(x.a), Scalar(0.));
}

template<typename Scalar> Dual<Scalar> ceil(const Dual<Scalar> & x) {
	return Dual<Scalar>(ceil(x.a), Scalar(0.));
}

template<typename Scalar> Dual<Scalar> round(const Dual<Scalar> & x) {
	return Dual<Scalar>(round(x.a), Scalar(0.));
}

template<typename Scalar> std::ostream & operator<<(std::ostream & s, const Dual<Scalar> & x)
{
	return (s << x.a);
}

#endif