/* 複合シンプソン則 */
/* simpson.c */
#include <stdio.h>
#include <math.h>
#define MMAX 1024
double func(double x);

int main(void)
{
  int m,k;
  double a,b,h,x,s,j,j0,I = exp(2.0) - 1.0;//真値Iがわかっているとき、誤差を出力するため
  a = 0.0; /* 積分区間の下限*/
  b = 2.0; /* 積分区間の上限*/  
  h = (b-a) /2.0;
  j = h * func(a + h) / 3.0;
  s = h * (func(a) + func(b)) /3.0 + 4 * j;
  m = 1;
  printf("%4s %6s %20s %10s\n","m","h","S_m ", "err");
  printf("%4d %9.5f %20.15f %9.3e\n",m,h,s,s-I);
  
  for (m = 2; m <= MMAX; m*=2){
    h *= 0.5;
    j0 = j;
    j = 0.0;
    for (k=1;k<=m;k++){
      x = a + (2 * k - 1) * h;
      j += h * func(x) / 3.0;
    }
    s = 0.5 * s + 4 * j -j0;
    printf("%4d %9.5f %20.15f %9.3e\n",m,h,s,s-I);
  }

 return 0;
}

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