auto_differentiate.cpp

auto_differentiate.cpp

#include <math.h>
#include <stdio.h>

template <class Number>
Number f(Number x) {
    Number z = x;
    for (int i = 0; i < 4; ++i) {
        z = z*z - 0.1 + sin(z*3);
    }
    if (z < 1) {
        return z;
    } else {
        return z / (z + 1);
    }
}


struct RealDiff {
    double x;
    double dx;

    RealDiff(double x, double dx) : x(x), dx(dx) {}

    typedef RealDiff Me;

    Me operator-(double y) { return Me(x - y, dx); }
    Me operator+(double y) { return Me(x + y, dx); }
    Me operator*(double y) { return Me(x * y, dx * y); }
    Me operator/(double y) { return Me(x / y, dx / y); }
    Me operator+(Me y) { return Me(x + y.x, dx + y.dx); }
    Me operator*(Me y) { return Me(x * y.x, dx*y.x + y.dx*x); }
    Me operator/(Me y) { return Me(x / y.x, (dx*y.x - y.dx*x) / (y.x*y.x)); }
    bool operator<(double y) { return x < y; }

    friend Me operator/(double x, Me y) { return Me(x,0)/y; }

    friend Me sin(Me y) { return Me(sin(y.x), cos(y.x)*y.dx); }
};

int main() {
    double a = -3, b = +3;
    for (int i = 0; i < 10000; ++i) {
        double t = i  / (10000.0 - 1);
        double x = a + (b-a)*t;
        RealDiff r = f(RealDiff(x, 1.0));
        double numeric_diff = (f(x + 1e-5) - f(x)) / 1e-5;
        printf("%g %g %g %g\n", x, r.x, r.dx, numeric_diff);
    }
}