// Romberg Integration

#include <stdio.h>
#include <math.h>

#define true 1
#define false 0
#define PI  3.1415926535

double func(double x)
{
   return sin(x);
}

void main()
{
   double romberg[2][15];  // 각 table 값을 저장   
   double h, X, sum;
   double left_end = 0.0;  // 좌측 경계값
   double right_end = PI;  // 우측 경계값
   int rows = 6;           // 원하는 열의 갯수
   int i, j, k, L, M, OK = true;

   if (OK) { 
      h = right_end - left_end;   
      romberg[0][0] = (func(left_end) + func(right_end)) / 2.0 * h;
      
      printf("Initial Data:\n"); 
      printf("Limits of integration = [%12.8f, %12.8f]\n", left_end, right_end);
      printf("Number of rows = %3d\n", rows);
      printf("\nRomberg Integration Table:\n");
      printf("\n%12.8f\n\n", romberg[0][0]);
       
      for ( i=2; i<=rows; i++) { 
         /* approximation from Trapezoidal method */ 
	 sum = 0.0;
	 M =  exp((i - 2) * log(2.0)) + 0.5;
	 for (k=1; k<=M; k++) sum = sum + func(left_end + (k - 0.5) * h);
	 romberg[1][0] = (romberg[0][0] + h * sum) / 2.0;

	 /* extrapolation */    
	 for (j=2; j<=i; j++) {
	    L = exp(2* (j - 1) * log(2.0)) + 0.5;
	    romberg[1][j-1] = romberg[1][j-2]+(romberg[1][j-2]-romberg[0][j-2])/(L - 1.0);
	 }

         for (k=1; k<=i; k++) printf(" %11.8f", romberg[1][k-1]); 
	 printf("\n\n");
         h = h / 2.0;   
         for (j=1; j<=i; j++) romberg[0][j-1] = romberg[1][j-1];
      }
   }
}