Automatic differentiation

dx.cpp

#include <cmath>
#include <iostream>

template <typename T>
struct dFloat {
    T x, dx;
    dFloat(T x, T dx = 0) : x(x), dx(dx) {}
};

template <typename T>
inline dFloat<T> operator + (const dFloat<T>& l, const dFloat<T>& r) {
    return dFloat<T> { l.x + r.x, l.dx + r.dx };
}

template <typename T>
inline dFloat<T> operator * (const dFloat<T>& l, const dFloat<T>& r) {
    return dFloat<T> { l.x * r.x, l.dx * r.x + r.dx * l.x };
}

template <typename T>
inline dFloat<T> operator - (const dFloat<T>& l, const dFloat<T>& r) {
    return dFloat<T> { l.x - r.x, l.dx - r.dx };
}

template <typename T>
inline dFloat<T> operator / (const dFloat<T>& l, const dFloat<T>& r) {
    return dFloat<T> { l.x / r.x, (l.dx * r.x - l.x * r.dx) / (r.x * r.x) };
}

template <typename T>
inline dFloat<T> sin(const dFloat<T>& f) { return dFloat<T> { std::sin(f.x),  f.dx * std::cos(f.x) }; }

template <typename T>
inline dFloat<T> cos(const dFloat<T>& f) { return dFloat<T> { std::cos(f.x), -f.dx * std::sin(f.x) }; }

template <typename T>
inline dFloat<T> tan(const dFloat<T>& f) {
    auto val = std::tan(f.x);
    return dFloat<T> { std::tan(f.x), f.dx * (val * val + 1) };
}

template <typename T>
inline dFloat<T> pow(const dFloat<T>& f, const dFloat<T>& e) {
    auto val = std::pow(f.x, e.x);
    return dFloat<T> { val, (val / f.x) * (e.x * f.dx + f.x * std::log(f.x) * e.dx) };
}

template <typename T>
inline dFloat<T> log(const dFloat<T>& f) {
    return dFloat<T> { std::log(f.x), f.dx / f.x };
}

template <typename T>
inline dFloat<T> log10(const dFloat<T>& f) {
    return dFloat<T> { std::log10(f.x), f.dx / (std::log(10.0f) * f.x) };
}

template <typename T>
inline dFloat<T> f(const dFloat<T>& x) {
    return pow(dFloat<T> { 2.0f } * sin(x) / x, log10(tan(x)) + cos(x));
}

dFloat<float> f_float(const dFloat<float>& x) {
    return f(x);
}

float df(float x) {
    using namespace std;
    auto a = cos(x) + log(tan(x))/log(10);
    auto b = 1 / (sin(x) * cos(x) * log(10)) - sin(x);
    auto c = sin(x) / x;
    return
        pow(2 * c, a) *
        ((cos(x) / x - sin(x) / (x * x)) * a / c + log(c) * b) +
        log(2) * pow(2 * c, a) * b;
}

int main() {
    size_t N = 100;
    for (size_t i = 0; i < N; ++i) {
        auto x = dFloat<float> { i * 0.1f, 1 };
        auto y = f(x);
        std::cout
            << "f(" << x.x << ") = " << y.x
            << ", df(" << x.x << ") = " << y.dx
            << ", ref = " << df(x.x)
            << std::endl;
    }
    return 0;
}