Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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);
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.