Calculus

#define NDEPS (value)
typedef double(*math_func)(double x);
double nthDerivative(math_func f, double x, int N){
    //perfectly fine derivative function
    double derivative(math_func f, double x);
    if(N<0) return NAN; //bogus value of N
    if(N==0) return f(x);
    if(N==1) return derivative(f,x);
    double* vals=calloc(2*N+9,sizeof(double)); //buffer region around the real values
    if(vals==NULL){ //oops! no memory
        return NAN;
    }
    int i,j;
    //don't take too small a finite difference
    double h=max(sqrt(DBL_EPSILON)*x,NDEPS);
    for(i=-(N+4);i<=N+4;i++){
        vals[i+N+4]=derivative(f,x+h*i);
    }
    //for(i=0;i<2*N+9;i++){printf("%.1e ",vals[i]);}putchar('\n');
    for(j=1;j<N;j++){
        double *vals2=calloc(2*N+9,sizeof(double));
        for(i=2;i<2*N+7;i++){
            vals2[i]=(-vals[i+2]+8*vals[i+1]-8*vals[i-1]+vals[i-2])/(12*h);
        }
        free(vals);
        vals=vals2;
        //for(i=0;i<2*N+9;i++){printf("%.1e ",vals[i]);}putchar('\n');
    }
    double result=vals[N+4];
    free(vals);
    return result;
}