/* Romberg積分法 */
/* romberg.c */
#include <stdio.h>
#include <math.h>
#define NUMAX 8
#define EPS 1e-12
double func(double x);

int main(void)
{
  int k,kend,nu,nuend,pow4k;
  double s,w,x,xa,xb,h,delta,integral;
  double t[NUMAX+1][NUMAX+1];

  xa = 0.0; /* 積分区間の下限 */
  xb = 2.0; /* 積分区間の上限 */
  //  integral = exp(2.0) - 1.0; /* 真値 */

  nu = 0;
  h = xb - xa;
  s = 0.5 * h * (func(xa) + func(xb));
  t[nu][0] = s;
  delta = 1.0;
  for(nu = 1;(nu <= NUMAX)&&(delta > EPS);nu++){
    h *= 0.5;
    x = xa + h;
    w = 0.0;
    while(x < xb){
      w += func(x);
      x += 2 * h;
    }
    s = 0.5 * s + h * w;
    t[nu][0] =s;
    pow4k = 1;
    for(k = 1;k <= nu;k++){
      pow4k *= 4;
      t[nu-k][k] = t[nu-k+1][k-1]+(t[nu-k+1][k-1]-t[nu-k][k-1])/(pow4k-1.0);
    }
    delta = fabs(t[0][nu-1]-t[0][nu]);
    //    printf("nu = %d delta = %6.2e\n",nu,delta);
  }
  nuend = nu - 1;

  
  for(nu = 0;nu <= nuend;nu++){
    printf("%2d ",nu);
    if (nu <= 6)
      kend = nu;
    else
      kend = 6;
    for(k = 0;k <= kend;k++){
      printf("%18.15f", t[nu-k][k]);
    }
    printf("\n");
  }

  printf("差分\n");
  for(nu = 1;nu <= nuend;nu++){
    printf("%2d ",nu);
    for(k = 1;k <= nu;k++){
      printf("%10.2e", t[nu-k][k]-t[nu-k+1][k-1]);
    }
    printf("\n");
  }
  
  return 0;
}

double func(double x)
{
  double y;

  y = exp(x);
   
  return y;
}