duals/dual is a single-header Dual number C++ template library that implements numbers of the type \((a + b \cdot \epsilon)\) where \( \epsilon \ne 0\), and \(\epsilon^2 = 0\).
duals::dual<> can be used for "automatic" differentiation, it can also recursively nest with itself for higher orders of differentiation (ie for the second derivative use duals::dual<duals::dual<T>>). It can also be used to differentiate parts of complex functions as std::complex<duals::dual<T>>. This file can be used stand-alone for dual number support, but for Eigen vectorization support rather the file duals/dual_eigen should be included.
using namespace duals::literals;
template <class T> T f(T x) { return pow(x,pow(x,x)); }
template <
class T> T df(T x) {
return pow(x,-1. + x + pow(x,x)) * (1. + x*
log(x) + x*pow(
log(x),2.)); }
template <
class T> T ddf(T x) {
return (pow(x,pow(x,x)) * pow(pow(x,x - 1.) + pow(x,x)*
log(x)*(
log(x) + 1.), 2.) +
pow(x,pow(x,x)) * (pow(x,x - 1.) *
log(x) +
pow(x,x - 1.) * (
log(x) + 1.) +
pow(x,x - 1.) * ((x - 1.)/x +
log(x)) +
pow(x,x) *
log(x) * pow(
log(x) + 1., 2.) )); }
int main()
{
std::cout << " f(2.) = " << f(2.) << "\n";
std::cout << " df(2.) = " << df(2.) << "\n";
std::cout << "ddf(2.) = " << ddf(2.) << "\n";
std::cout << " f(2+1_e) = " << f(2+1_e) << "\n";
std::cout << " f(2+1_e).dpart() = " << f(2+1_e).dpart() << "\n";
std::cout << " f((2+1_e) + (1+0_e)_e).dpart().dpart() = " << f(x).dpart().dpart() << "\n";
}
Produces (notice the derivative in the dual-part):
f(2.) = 16
df(2.) = 107.11
ddf(2.) = 958.755
f(2+1_e) = (16+107.11_e)
f(2+1_e).dpart() = 107.11
f((2+1_e) + (1+0_e)_e).dpart().dpart() = 958.755
How this works can be seen by inspecting the infinite Taylor series expansion of a function \( f(x) \) at \( a \) with \( x = a + b\epsilon \). The series truncates itself due to the property \( \epsilon^2 = 0 \), leaving the function's derivative at \( a \) in the "dual part" of the result (when \( b = 1 \)):
\[ \begin{split} f(a + b \epsilon) &= f(a) + f'(a)(b \epsilon) + \frac{f''(a)}{2!}(b \epsilon)^2 + \frac{f'''(a)}{3!}(b \epsilon)^3 + \ldots \\ \ &= f(a) + f'(a)(b \epsilon) \end{split} \]
The class is contained in a single c++11 header #include <duals/dual>, for extended Eigen support, instead include the header #include <duals/dual_eigen>.
Type X in the templates below can be any value which can be assigned to value_type.
Type X also indicates a limitation to dual numbers of the same depth but (possibly) different value_type as duals::dual<T>. For example, you can assign (or add/sub/mul/div) duals::dual<float> and duals::dual<double>, but you can not assign duals::dual<duals::dual<float>> to duals::dual<float>.
Here is a synopsis of the class:
namespace duals {
template<class T> class dual {
dual(
const & re = T(),
const & du = T());
template<
class X>
dual(
const dual<X> &);
template<
class X> dual<T> &
operator= (
const dual<X> &);
template<
class X> dual<T> &
operator+=(
const dual<X> &);
template<
class X> dual<T> &
operator-=(
const dual<X> &);
template<
class X> dual<T> &
operator*=(
const dual<X> &);
template<
class X> dual<T> &
operator/=(
const dual<X> &);
bool operator ==(const X &b) const;
bool operator !=(const X &b) const;
bool operator <(const X &b) const;
bool operator >(const X &b) const;
bool operator <=(const X &b) const;
bool operator >=(const X &b) const;
bool operator <(const dual<X> &b) const;
bool operator >(const dual<X> &b) const;
bool operator <=(const dual<X> &b) const;
bool operator >=(const dual<X> &b) const;
};
dual<T> random(dual<T> low = {0,0}, dual<T> high = {1,0})
dual<T> log10(dual<T>)
dual<T> pow(dual<T>, U)
dual<T> pow(U, dual<T>)
dual<T> pow(dual<T>, dual<T>)
dual<T> sqrt(dual<T>)
dual<T> cbrt(dual<T>)
dual<T> sin(dual<T>)
dual<T> cos(dual<T>)
dual<T> tan(dual<T>)
dual<T> asin(dual<T>)
dual<T> acos(dual<T>)
dual<T> atan(dual<T>)
dual<T> sinh(dual<T>)
dual<T> cosh(dual<T>)
dual<T> tanh(dual<T>)
dual<T> asinh(dual<T>)
dual<T> acosh(dual<T>)
dual<T> atanh(dual<T>)
template<T>
operator>>(basic_istream<charT, traits> &, dual<T> &);
template<T>
operator<<(basic_ostream<charT, traits> &,
const dual<T> &);
}
Some useful typedefs:
typedef dual<float> dualf;
typedef dual<double> duald;
typedef dual<long double> dualld;
typedef dual<dualf> hyperdualf;
typedef dual<duald> hyperduald;
typedef dual<dualld> hyperdualld;
typedef std::complex<dualf> cdualf;
typedef std::complex<duald> cduald;
typedef std::complex<dualld> cdualld;
There are also literals for dual parts defined in the inline namespace duals::literals. They are available with using namespace duals; or using namespace duals::literals.
using namespace duals::literals;
dualf x = 3 + 4_ef;
duald y = 3 + 4_e;
dualld z = 3 + 4_el;
And in case you dislike iostreams, there are some formatters for the {fmt} formatting library. These are disabled by default, but can be enabled by #define CPPDUALS_LIBFMT for the dual<> formatter, and/or #define CPPDUALS_LIBFMT_COMPLEX for the std::complex<> formatter. There are three custom formatting flags that control how the dual numbers are printed, these must come before the normal {fmt} formatting spec: '$', '*', ','. For me these are about 3x faster than iostreams.
#define CPPDUALS_LIBFMT
#define CPPDUALS_LIBFMT_COMPLEX
#inlcude <duals/dual>
using namespace duals::literals;
...
string s = fmt::format("{:}", 1 + 2_e);
string s = fmt::format("{:g}", 1 + 2_e);
string s = fmt::format("{:*}", 1 + 2_e);
string s = fmt::format("{:,}", 1 + 2_e);
string s = fmt::format("{:*,g}", complexd(1,2_e));
The "archaic Greek epsilon" logo is from Wikimedia commons
Some casual reading material:
