/* グッラグ法(修正中点法にリチャードソン補外) */
/* gragg.c */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define KMAX 8
#define NUMAX 8
#define EPS 1e-14
double func(double t,double x);
double euler(double a,double b,double x_ini,int n);
double mod_midpoint(double a,double b,double x_ini,int n);
int count;

int main(void)
{
  int k,kend,nu,pow4k;//Richardson extrapolation
  int n;
  double h,t,x,x0,t_ini,t_end,x_ini,mx;
  //  double tv; /* 真値 */
  double delta; /* 差分 */
  double T[NUMAX+1][NUMAX+1];

 /* 初期条件 */
  t_ini = 0.0;
  t_end = 1.0;
  x_ini = 1.0;
  //  tv = exp(-t_end * t_end); /* 真値 */
 /* 初期条件 */

  for(nu = 0;nu <= NUMAX;nu++){
    if (nu == 0){
      n = 1;
      x = x_ini;
      x = x + (t_end - t_ini) * func(t_ini,x);
      T[nu][0] = x;
    }
    else{
      n *= 2;
      mx = mod_midpoint(t_ini,t_end,x_ini,n);
      T[nu][0] = mx;
    }
  }

  for(nu = 1;nu <= NUMAX;nu++){
    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);
    }
  }    

  for(nu = 0;nu <= NUMAX;nu++){
    printf("%2d ",nu);
    kend = nu;
    for(k = 0;k <= kend;k++){
      printf("%18.15f", T[nu-k][k]);
    }
    printf("\n");
  }

  printf("差分\n");
  for(nu = 0;nu <= NUMAX;nu++){
    printf("%2d ",nu);
    if (nu <= 4)
      kend = nu;
    else
      kend = 4;
    for(k = 1;k <= nu;k++){
      printf("%10.2e", T[nu-k][k]-T[nu-k+1][k-1]);
    }    
    //    for(k = 0;k <= nu;k++){
    //      printf("%10.2e", T[nu-k][k]-tv);
    //    }
    printf("\n");
  }

  printf("関数値計算回数=%d\n",count);
 
  return 0;
}

double func(double t, double x)
{
  count++;
  return -2 * t * x;
  //  return -2 * x;
}

double euler(double a,double b,double x_ini,int n)
{
  int i;
  double h,x,t;

  h = (b - a) /n;
  i = 0;
  t = a;   /* 初期条件 */
  x = x_ini; /* 初期条件 */
  for (i = 1;i <= n;i++){
    x += h * func(t,x);
    t += h;
  }

  return x;
}

double mod_midpoint(double a,double b,double x_ini,int n)
{
  int nu;
  double h,x,x0,x00,mx,t;

  h = (b - a) / n;
  t = a;
  x = x_ini;
  
  nu = 1;
  x0 = x_ini;
  x = x + h * func(t,x); //Euler
  t += h;
  nu = 2;
  x00 = x0;
  x0 = x;
  x = x00 + 2 * h * func(t,x);
  for (nu = 3;nu <= n;nu++){
      x00 = x0;
      x0 = x;
      t += h;
      x = x00 + 2 * h * func(t, x);
  }
  mx = 0.5 * (x0 + x + h * func(t+h,x));

  return mx;
}