/* Aitken Δ2乗法を用いた振動積分 */
/* oscillate_aitken.c */
#include <stdio.h>
#include <math.h>
#define NUMAX 20
#define EPS 1e-15
double func(double x);
double gauss_legendre(double a,double b,int p);

int main(void)
{
  int n,pannel,i,nu,k,kend,nuend;
  double half_period = M_PI; /* 周期の半分 */
  double s,w,x,xa,xb;
  double delta,temp; /* 差分 */
  //  double integral = M_PI_2; /* 真値 p.193*/
  //  double integral = M_PI_2/M_E; /* 真値 p.194*/
  double t[NUMAX+1][NUMAX+1];
  
  n = 8;
  pannel = 1; /* 複合ガウス-ルジャンドル則の半周期の積分のパネル数 */
  nu = 0;
  delta = 1.0;
  //  xa = xinf;
  //  xb = xa + half_period;
  xa = 0.0;
  xb = xa + half_period/((double )pannel);
  w = gauss_legendre(xa,xb,n);
  for (i = 1; i< pannel;i++){
    xa = xb;
    xb = xa + half_period/((double )pannel);
    w += gauss_legendre(xa,xb,n);
  }
  s = w;
  t[nu][0] = s;
  for(nu = 1;(nu <= NUMAX)&&(delta>EPS);nu++){
    xa = xb;
    xb = xa + half_period / ((double)pannel);
    w = gauss_legendre(xa,xb,n);
    w = gauss_legendre(xa,xb,n);
    for (i = 1; i< pannel;i++){
      xa = xb;
      xb = xa + half_period/((double )pannel);
      w += gauss_legendre(xa,xb,n);
    }    
    s += w;
    t[nu][0] = s;
    kend = nu/2;
    for(k = 1;k <= kend;k++){
      t[nu-2*k][k] = t[nu-2*k][k-1]
	-(t[nu-2*k+1][k-1]-t[nu-2*k][k-1])*(t[nu-2*k+1][k-1]-t[nu-2*k][k-1])/
	(t[nu-2*k+2][k-1]-2*t[nu-2*k+1][k-1]+t[nu-2*k][k-1]);
    }
    delta = 1.0;
    for(k = 2;k <= kend;k++){
      temp = fabs(t[nu-2*k][k] - t[nu-2*k+2][k-1]);
      if (delta > temp) delta = temp;
    }
  }
  nuend = nu -1;

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

  //  printf("真値%20.15f\n",integral);
  printf("零点間の積分は%dパネルGauss-Legendre %d 点則\n",pannel,n);
  for(nu = 2;nu <= nuend;nu++){
    printf("%2d ",nu);
    kend = nu/2;
    for(k = 1;k <= kend;k++){
      printf("%10.2e", t[nu-2*k][k]-t[nu-2*k+2][k-1]);
    }
    //    for(k = 0;k <= kend;k++){
    //      printf("%10.2e", t[nu-2*k][k]-integral);
    //    }
    printf("\n");
  }

  return 0;
}

double func(double x)
{
  double y;

  //  y = x * sin(x) / (x * x + 1.0); /* p.194 */
  y = sin(x) / x; /* 例7.14 p.193 */
  
  return y;
}

double gauss_legendre(double a,double b,int n)
/* ガウス-ルジャンドルn点則 */
{
  int j;
  double s;
  double xz[n+1][n+1];
  double w[n+1][n+1];
  double c,d,xj,wj;
  
  switch(n){
  case 4:
    xz[4][1] =-sqrt(3.0/7.0 + (2.0/7.0)*sqrt(6.0/5.0));
    xz[4][2] =-sqrt(3.0/7.0 - (2.0/7.0)*sqrt(6.0/5.0));
    xz[4][3] =sqrt(3.0/7.0 - (2.0/7.0)*sqrt(6.0/5.0));
    xz[4][4] =sqrt(3.0/7.0 + (2.0/7.0)*sqrt(6.0/5.0));
    w[4][1] = (18 - sqrt(30))/36.0;
    w[4][2] = (18 + sqrt(30))/36.0;
    w[4][3] = (18 + sqrt(30))/36.0;
    w[4][4] = (18 - sqrt(30))/36.0;
    break;
  case 5:
    xz[5][1] = -(1.0/3.0) * sqrt(5 + 2 * sqrt(10.0/7.0));
    xz[5][2] = -(1.0/3.0) * sqrt(5 - 2 * sqrt(10.0/7.0));
    xz[5][3] = 0.0;
    xz[5][4] = (1.0/3.0) * sqrt(5 - 2 * sqrt(10.0/7.0));
    xz[5][5] = (1.0/3.0) * sqrt(5 + 2 * sqrt(10.0/7.0));
    w[5][1] = (322 - 13 * sqrt(70))/900;
    w[5][2] = (322 + 13 * sqrt(70))/900;
    w[5][3] = 128.0 /225.0;
    w[5][4] = (322 + 13 * sqrt(70))/900;
    w[5][5] = (322 - 13 * sqrt(70))/900;
    break;
  case 6:
    xz[6][1] = -0.9324695142031521;
    xz[6][2] = -0.6612093864662645;
    xz[6][3] = -0.2386191860831969;
    xz[6][4] = 0.2386191860831969;
    xz[6][5] = 0.6612093864662645;
    xz[6][6] = 0.9324695142031521;
    w[6][1] = 0.1713244923791704;
    w[6][2] = 0.3607615730481386;
    w[6][3] = 0.4679139345726910;
    w[6][4] = w[6][3];
    w[6][5] = w[6][2];
    w[6][6] = w[6][1];
    break;
  case 8:
    xz[8][1] = -0.9602898564975363;
    xz[8][2] = -0.7966664774136267;
    xz[8][3] = -0.5255324099163290;
    xz[8][4] = -0.1834346424956498;
    xz[8][5] = 0.1834346424956498;
    xz[8][6] = 0.5255324099163290;
    xz[8][7] = 0.7966664774136267;
    xz[8][8] = 0.9602898564975363;
    w[8][1] = 0.1012285362903763;
    w[8][2] = 0.2223810344533745;
    w[8][3] = 0.3137066458778873;
    w[8][4] = 0.3626837833783620;
    w[8][5] = w[8][4];
    w[8][6] = w[8][3];
    w[8][7] = w[8][2];
    w[8][8] = w[8][1];
    break;
  }

  c = (b - a)/2.0;
  d = (b + a)/2.0;
  s = 0.0;
   
  for (j = 1;j <= n; j++){
    xj = xz[n][j];
    wj = w[n][j];
    s += wj * c * func(c * xj + d);
  }

  return s;
}