#include <iostream>
#include <iomanip>
#include <math.h>

typedef double (*Function) (double);

/*
  This is a formula for approximating the derivative of
  f at the point x. (h - approximation parameter)
 */
double derivative(Function f, double x, double h)
{
   return (f(x+h)-f(x-h))/(2*h);
}

/* This is the Richardson extrapolation formula from page 509
   of your textbooks.
*/   
double D(int n, int k, Function f, double x, double h = 1)
{

   //make sure the parameters are valid.
   if (n<0) n=0;
   if (k<0) k=0;
   if (k>n) k=n;

   if (k==0)
      return derivative(f,x,h/pow(2,n));

   double p4 = pow(4,k);
   return p4/(p4-1)*D(n,k-1,f,x,h) - 1.0/(p4-1)*D(n-1,k-1,f,x,h);
}

int main()
{
   /* Using the above funtions to approximate the value of the
      derivative of 'sin' at the point 'M_PI/2'.
   */
   for (int n=0; n<5; n++)
   {
      std::cout << n << " ";
      std::cout << std::setiosflags(std::ios::fixed) << std::setprecision(8);
      
      for (int k=0; k<=n; k++)
         std::cout << D(n,k, sin, M_PI/3) << " ";
      std::cout << std::endl;
   }
}