/*************************************************************
* Program to demonstrate integration of a real function F(x) *
*                    by Romberg's method                     *
* ---------------------------------------------------------- *
* REFERENCE: "Mathematiques en Turbo-Pascal (Part 1) by      *
*             Marc Ducamp and Alain Reverchon, Eyrolles,     *
*             Paris, 1987" [BIBLI 03].                       *
*                             C++ version by J-P Moreau      *
*                                 (www.jpmoreau.fr)          *
* ---------------------------------------------------------- *
* SAMPLE RUN:                                                *
* (Integrate sin(x) from x=0 to x=1)                         *
*                                                            *
* Integral of a function F(X) by Romberg's method            *
*                                                            *
* Input begin and end values for x variable:                 *
*                                                            *
*         X0 = 0                                             *
*         X1 = 1                                             *
*                                                            *
* Input desired precision: 1e-10                             *
*                                                            *
*                                                            *
* Value of integral   :  4.596977e-001                       *
*                                                            *
* Obtained precision  :  9.599083e-011                       *
*                                                            *
* Number of iterations:  4                                   *
*                                                            *
*************************************************************/
#include <stdio.h>
#include <math.h>

#define MAXITER  15

double  x0,       // begin x value
        x1,       // end   x value
        prec,     // desired precision
        integral, // result of integral
        obtprec;  // obtained precision

int     niter;    // number of actual iterations


  // Given function to integrate
  double FUNC(double x) {
    return sin(x);
  }


  /******************************************************
  * Integral of a function FUNC(X) by Romberg's method  *
  * --------------------------------------------------- *
  * INPUTS:                                             *
  *          a      begin value of x variable           *
  *          b      end value of x variable             *
  *       prec      desired precision                   *
  *    itermin      minimum number of ietrations        *
  *    itermax      maximum number of iterations        *
  *                                                     *
  * OUTPUTS:                                            *
  *    obtprec      obtained precision for integral     *
  *          n      number of iterations done           *
  *                                                     *
  * RETURNED VALUE  the integral of FUNC(X) from a to b *
  *                                                     *
  ******************************************************/ 
  double RombergIntegral(double a,double b,double prec, double *obtprec,
                                           int *n, int itermin, int itermax)  {
    int    i,j;
    double pas,r,s,ta;
    double t[MAXITER][MAXITER];
    if (itermax>MAXITER) itermax=MAXITER;
    r = FUNC(a);
    ta = (r + FUNC(b)) / 2;
    *n=0;
    pas=b-a;
    t[0][0]=ta*pas;
	do {
      *n = *n + 1;
      pas = pas/2;
      s=ta;
	  for (i=1; i<pow(2,*n); i++)
        s += FUNC(a+pas*i);
      t[0][*n]=s*pas;
      r=1;
	  for (i=1; i<*n+1; i++) {
        r=r*4;
        j=*n-i;
        t[i][j]=(r*t[i-1][j+1] - t[i-1][j])/(r-1);
	  }
      *obtprec = fabs(t[*n][0] - t[*n-1][0]);
      if (*n>=itermax) return t[*n][0];
	} while (*obtprec >prec || *n<itermin); 
    return t[*n][0];
  }

          
void main()  {
  
  printf("\n Integral of a function F(X) by Romberg's method\n\n");
  printf(" Input begin and end values for x variable:\n\n");

  printf("         X0 = "); scanf("%lf",&x0);
  printf("         X1 = "); scanf("%lf",&x1);
  
  printf("\n Input desired precision: "); scanf("%lf",&prec);

  integral = RombergIntegral(x0,x1,prec,&obtprec,&niter,1,50);

  printf("\n\n\n Value of integral   : %e\n\n", integral);
  printf(" Obtained precision  : %e\n\n", obtprec);
  
  printf(" Number of iterations: %d\n\n", niter);

}

// End of file tromberg.cpp