/********************************************************
*  Estimate the nth derivative of a real function f(x)  *
*  at point x with step h (n from 1 to 5).              *
* ----------------------------------------------------- *
* SAMPLE RUN:                                           *
*                                                       *
* Nth derivative of a real function f(x):               *
*                                                       *
* Function to derivate: x^3                             *
*                                                       *
* At point x0= 10                                       *
* Order of derivative (1 to 5): 1                       *
* Value of step h= 0.1                                  *
*                                                       *
* f'(x0)=300.000000                                     *  
*                                                       *
* At point x0= 10                                       *
* Order of derivative (1 to 5): 2                       *
* Value of step h= 0.1                                  *
*                                                       *
* f''(x0)=60.000000                                     *
*                                                       *
* ----------------------------------------------------- *
* Ref.: "Analyse numérique en C By Alain Reverchon and  *
*        Marc Ducamp, Armand Colin Editeur, Paris,      *
*        1990" [BIBLI 14].                              *
*                                                       *
*               C++ (simplified) version by J-P Moreau. *
*                        (www.jpmoreau.fr)              *
********************************************************/
#include <stdio.h>
#include <math.h>


//sample function f(x)=x^3
double f(double x) {
  return x*x*x;
}

//returns the value of nth derivative of function
//f(x) at point x0 with increment step h
//Note that for polynomial functions h must not be
//too small. Here h=1 gives best results than h=0.01
double ar_nderiv(double x,int n,double h) {
  static int table[5][7] = {
	{3, 45, -9, 1, 0, 0, 60},
    {3, 270, -27, 2, 0, 0, 180},
	{4, -488, 338, -72, 7, 0, 240},
	{4, -1952, 676, -96, 7, 0, 240},
	{5, 1938, -1872, 783, -152, 13, 288}
  };
  int i, *coeff;
  double fa,fb,fc,t;
  if(n>5) return 0;
  coeff = table[n-1];
  fa = f(x);
  for (i=1, t=0; i<=coeff[0]; i++) {
    fb = f(x-i*h);
	fc = f(x+i*h);
    t = t + coeff[i]*((n%2) ?
		(fc-fb) : ((fc-fa) + (fb-fa)));
  }
  t /= coeff[6];
  for (i=1; i<=n; i++)   t /= h;
  return t;
}


void main() {
  double d,h,x; int n;
  printf("\n Nth derivative of a real function f(x):\n\n");
  printf(" Function to derivate: x^3\n\n");
  printf(" At point x0= "); scanf("%lf",&x);
  printf(" Order of derivative (1 to 5): "); scanf("%d",&n);
  printf(" Value of step h= "); scanf("%lf",&h);
  printf("\n");

  d = ar_nderiv(x,n,h);

  switch(n) {
    case 1: printf(" f'(x0)=%f\n\n",d); break;
    case 2: printf(" f''(x0)=%f\n\n",d); break;
    case 3: printf(" f'''(x0)=%f\n\n",d); break;
    case 4: printf(" f''''(x0)=%f\n\n",d); break;
    case 5: printf(" f'''''(x0)=%f\n\n",d); break;
    default: printf(" n must be between 1 and 5 !\n\n");
  }
}

//end of file nderiv.cpp