/* gauss_legendre.c */
/* ガウス・ルジャンドル則 */
#include <stdio.h>
#include <math.h>
#define NMAX 8
double func(double xx);

int main(void)
{
  int n,j;
  double a,b,c,d, xj, wj,s = 0.0;
  double I = exp(2.0) -1.0;//真値がわかっている場合にIを用いる
  double xz[NMAX+1][NMAX+1];
  double w[NMAX+1][NMAX+1];

  a = 0.0; /* 積分区間の下限*/
  b = 2.0; //b = 2 * M_PI; /* 積分区間の上限*/
  c = (b - a)/2.0;
  d = (b + a)/2.0;
  
  xz[2][1] = -sqrt(1.0/3.0);
  xz[2][2] = sqrt(1.0/3.0);  
  xz[3][1] = -sqrt(3.0/5.0);
  xz[3][2] = 0;
  xz[3][3] = sqrt(3.0/5.0);
  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));
  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[2][1] = 1.0;
  w[2][2] = 1.0;  
  w[3][1] = 5.0/9.0;
  w[3][2] = 8.0/9.0;
  w[3][3] = 5.0/9.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;
  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;

  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];
  
  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];
  
  printf("%4s %20s %10s\n","n","s","err");
  for (n = 2; n <= NMAX; n++){
    if (n == 7) continue;
    s = 0.0;
   
    for (j = 1;j <= n; j++){
      xj = xz[n][j];
      wj = w[n][j];
      s += wj * c * func(c * xj + d);
     }
    printf("%4d %20.15f %9.3e\n",n,s,s-I); /* 真値がわかっている場合 */
    // printf("%4d %20.15f\n",n,s); /* 真値わからない場合 */
  }

 return 0;
}

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