#include #include /* Do a friends of one of FHK's problems. Sum(n=0..infinity) 1/(n+pi). */ /* First forwards (finding out where no futher change happens, and then */ /* backwards. Note the maths says the answer is infinite... */ double g(int i) { return 1.0/(i + M_PI); } int main() { double d, dd; int n = 0, i; float f=0, ff; do { ff = f; f += 1.0/(n + M_PI); n++; } while (f != ff); printf("Forward (float) sum until no change: %11.7f at n=%d\n", f, n); f = 0; for (i = n-1; i>=0; i--) f += 1.0/(i + M_PI); printf("Backward (float) sum: %11.7f\n", f); d = 0; for (i = 0; i=0; i--) d += g(i); printf("Backward (double) sum: %20.16f\n", d); /* now use Kahan float summation, using ff as guard digits... */ ff = 0; f = g(0); for (i = 1; i=0; i--) { float y = g(i) - ff; float t = f + y; ff = (t - f) - y; f = t; } printf("Backward (float) Kahan sum: %11.7f\n", f); /* now use Kahan summation, using dd as guard digits... */ dd = 0; d = g(0); for (i = 1; i=0; i--) { double y = g(i) - dd; double t = d + y; dd = (t - d) - y; d = t; } printf("Backward (double) Kahan sum: %20.16f\n", d); }