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

#ifndef D
#define FLOAT float
#else
#define FLOAT double
#endif

FLOAT function(FLOAT x) {
#ifndef NEG
    return exp(x);
#else
    return exp(-x);
#endif
}

/* 3 points Lagrange */
FLOAT deriv1(FLOAT (*f) (FLOAT), FLOAT x, FLOAT h) {
    FLOAT xl = x - h;
    FLOAT xr = x + h;
    double fl = (*f)(xl);
    double fr = (*f)(xr);

    return (fr - fl)/(2 * h);
}

/* 5 points Lagrange */
FLOAT deriv2(FLOAT (*f) (FLOAT), FLOAT x, FLOAT h) {
    FLOAT x0 = x - h * 2;
    FLOAT x1 = x - h;
    FLOAT x3 = x + h;
    FLOAT x4 = x + h * 2;
    double f0 = (*f)(x0);
    double f1 = (*f)(x1);
    double f3 = (*f)(x3);
    double f4 = (*f)(x4);

    return (f0 - f1 * 8 + f3 * 8 - f4)/(12 * h);
}

/* 9 points Lagrange */
FLOAT deriv3(FLOAT (*f) (FLOAT), FLOAT x, FLOAT h) {
    FLOAT x0 = x - h * 4;
    FLOAT x1 = x - h * 3;
    FLOAT x2 = x - h * 2;
    FLOAT x3 = x - h;
    FLOAT x5 = x + h;
    FLOAT x6 = x + h * 2;
    FLOAT x7 = x + h * 3;
    FLOAT x8 = x + h * 4;
    double f0 = (*f)(x0);
    double f1 = (*f)(x1);
    double f2 = (*f)(x2);
    double f3 = (*f)(x3);
    double f5 = (*f)(x5);
    double f6 = (*f)(x6);
    double f7 = (*f)(x7);
    double f8 = (*f)(x8);

    return (f0 * 3 - f1 * 32 + f2 * 168 - f3 * 672 +
	    f5 * 672 - f6 * 168 + f7 * 32 - f8 * 3)/(840 * h);
}

/* n points HdB */
/* http://huizen.dto.tudelft.nl/deBruijn/programs/shannon.htm */
/* n should be odd and larger than 4 */
FLOAT derivn(FLOAT (*f) (FLOAT), FLOAT x, FLOAT h, int n) {
    FLOAT *xk;
    double *fk;
    double sample = h;
    double sigma = sample * 2;
    double norm = 1 / (sqrt(2 * M_PI) * 2);
    int i;
    double result = 0.0;
    double d;
    int hn;

    if(n < 5 || (n & 1) != 1) {
	fprintf(stderr, "HdB called with wrong n: %d\n", n);
	exit(1);
    }
    xk = malloc(n * sizeof(*xk));
    fk = malloc(n * sizeof(*fk));
    hn = (n - 1) / 2;

    for(i = - hn; i <=  hn; i++) {
	fk[i + hn] = (*f)(xk[i + hn] = x + i * sample);
    }
    for(i = 0; i < n; i++) {
	d = (x - xk[i])/sigma;
	result -= fk[i] * d * exp(- d * d / 2) / sigma;
    }
    free(xk);
    free(fk);
    return norm * result;
}

main() {
    FLOAT h = 1;
    double e;
    FLOAT result1, result2, result3, result4, result5, result6, result7;
    FLOAT result8;
    int i;

#ifndef NEG
    e = exp(1.0);
#else
    e = - exp(- 1.0);
#endif
    for(i = 0; i < 23; i++) {
	result1 = deriv1(&function, 1.0, h);
	result2 = deriv2(&function, 1.0, h);
	result3 = deriv3(&function, 1.0, h);
	result4 = derivn(&function, 1.0, h, 5);
	result5 = derivn(&function, 1.0, h, 9);
	result6 = derivn(&function, 1.0, h, 17);
	result7 = derivn(&function, 1.0, h, 33);
	result8 = derivn(&function, 1.0, h, 65);
	printf("%2d %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e\n", i,
	    fabs((result1 - e) / e), fabs((result2 - e) / e),
	    fabs((result3 - e) / e), fabs((result4 - e) / e),
	    fabs((result5 - e) / e), fabs((result6 - e) / e),
	    fabs((result7 - e) / e), fabs((result8 - e) / e));
	h = h / 2;
    }
}