Dual numbers let you compute derivatives to machine precision.
The basic idea is that if there were a number ε ≠ 0 with ε2 = 0 then f(x + ε) = f(x) + f'(x)ε.
For example, (x + ε)2 = x2 + 2xε so the derivative of x2 is 2x.
Likewise, if there were a number with εj ≠ 0 for j = 1, ..., n-1 and εn = 0 then f(x + ε) = f(x) + f'(x)ε + ... + f(n-1)(x) εn-1/(n-1)!.
There are such "numbers". Let ε be a square matrix having upper diagonal entries equal to 1
and 0 everywhere else. Use
operator() in this library to extract derivatives. This is much simpler than the
machinery involved with Automatic Differentiation.