/********************************************************
*  Calculate a limited development of a real function   *
*  f(x) at point x0 with step h (at order 5).           *
* ----------------------------------------------------- *
* SAMPLE RUN:                                           *
*                                                       *
* Limited development of f(x) at x0:                    *
*                                                       *
* Function to develop: 1/(1-x)                          *
*                                                       *
* At point x0= 0                                        *
* Value of step h= 0.1                                  *
*                                                       *
* a1 =  1.0000416250                                    *
* a2 =  1.0000416250                                    *
* a3 =  1.0011207928                                    *
* a4 =  1.0011207928                                    *
* a5 =  1.0136468470                                    *
*                                                       *
* (The correct answer is a1=a2=a3=a4=a5=1)              *         
*                                                       *
* Function to develop: log(x)                           *
*                                                       *
* At point x0= 1                                        *
* Value of step h= 0.1                                  *
*                                                       *
* a1 =  1.0000057552                                    *
* a2 = -0.5000050520                                    *
* a3 =  0.3334508142                                    *
* a4 = -0.2501062334                                    *
* a5 =  0.2011301593                                    *
*                                                       *
* (the correct answer is a1=1, a2=-0.5, a3=1/3,         *
*  a4=-1/4 and a5=1/5)                                  *   
* ----------------------------------------------------- *
* Ref.: "Mathematiques en Turbo Pascal By Alain         *
*        Reverchon and Marc Ducamp, Editions Eyrolles,  *
*        Paris, 1991" [BIBLI 03].                       *
*                                                       *
*                            C++ version by J-P Moreau. *
*                                (www.jpmoreau.fr)      *
*********************************************************
Note: as you can see this direct method is not suitable to
      calculate a limited development in an accurate way.
      The next programs (ltddev1 to ltddev3) will allow us
      to calculate limited developments at x=0 up to a high
      order with a good accuracy for functions f*g, f/g or
      gof, knowing the limited developments of f(x) and g(x)
      at x=0.
----------------------------------------------------------*/
#include <stdio.h>
#include <math.h>

double  d,h,r,s,x0;
int n;


//sample function f(x)=1/(1-x)
double f(double x) {
  if (fabs(x-1) >1e-12)
    return 1.0/(1.0-x);
  else
    return 1e12;
}

/*sample function f(x)=Ln(x)
double f(double x) {
  if (x >1e-12)
    return log(x);
  else
    return -1e12;
}                         */

//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()  {
  
  printf("\n Limited development of a real function f(x):\n\n");
  printf(" Function to develop: 1/(1-x)\n");
  printf("\n At point x0= "); scanf("%lf", &x0);
  printf("\n Value of step h= "); scanf("%lf", &h);
  printf("\n\n");

  s=1;

  for (n=1; n<6; n++) {
    s *= n;
    r = ar_nderiv(x0,n,h) / s;
    printf(" a%d = %13.10f\n", n, r);
  }
  printf("\n");

}

// end of file ltddev.cpp