/*
  Gauss-Legendre数値積分公式
 */
#include <stdio.h>
#include <math.h>
#define NPNT 101
//
/*
　Legendre多項式
*/
double Pn_Legendre(int n, double x)
{
  double p0, p1, p2, ri, pp;
  int r;
  //
  if (n==0) pp = 1.0;
  if (n==1) pp = x;
  if (n>=2)
    {
      p0 = 1.0;
      p1 = x;
      for (r=2; r<=n; ++r)
	{
	  ri = 1.0 / (double)r;
	  p2 = (2.0-ri)*x*p1 - (1.0-ri)*p0;
	  p0 = p1;
	  p1 = p2;
	}
      pp = p2;
    }
  return pp;
}
/*
  Legendre多項式の導関数
*/
double dPn_Legendre(int n, double x)
{
  double dpn, pn, pn_1;
  //
  if (n==0) dpn = 0.0;
  if (n>=1)
    {
      pn = Pn_Legendre(n, x);
      pn_1 = Pn_Legendre(n-1, x);
      dpn = n * (pn_1 - x * pn) / (1.0 - x*x);
    }
  return dpn;
}
/* 

数値積分公式の標本点（Legendre多項式の零点），重みを計算．
nmax : 点数の最大値
xk[n][k] : n点公式のk番目の標本点
wk[n][k] : n点公式のk番目の重み

*/
void point_weight_GaussLe(int nmax, double xk[NPNT+1][NPNT+1], double wk[NPNT+1][NPNT+1])
{
  int n, k;
  int j, jmax=50;
  double xk0, xk1, pn, dpn;
  double x, ww, pn_1;
  double err_tol = 1.0e-15, err_func;
  int l_conv;
  //
  for (n=1; n<=nmax; ++n)
    {
      for (k=1; k<=n; ++k)
	{
	  l_conv = 0;
	  /* initial values for Newton's method */
	  xk0 = cos(M_PI*(k-0.25)/(n+0.5));
	  for (j=1; j<=jmax; ++j)
	    {
	      pn = Pn_Legendre(n, xk0);
	      dpn = dPn_Legendre(n, xk0);
	      xk1 = xk0 - pn / dpn;
	      pn = Pn_Legendre(n, xk1);
	      dpn = dPn_Legendre(n, xk1);
	      err_func = err_tol * fabs(dpn);
	      if (fabs(pn) < err_func) 
		{
		  l_conv = 1;
		  break;
		}
	      xk0 = xk1;
	    }
	  xk[n][k] = xk1;
	  if (l_conv == 0) printf("k=%d : not convergent\n", k);
	}
      //
      for (k=1; k<=n; ++k)
	{
	  x = xk[n][k];
	  pn_1 = Pn_Legendre(n-1, x);
	  ww = 2 * (1.0 - x*x) / (n*n*pn_1*pn_1);
	  wk[n][k] = ww;
	}
    }
}
/*
Gauss-Legendre公式による数値積分
func : 被積分関数
(xmin, xmax) : 積分区間
n : 数値積分の点数
xk, wk : Gauss-Legendre数値積分公式の標本点，重み（point_weight_GaussLeで計算）
*/
double integral(double (*func)(double), double xmin, double xmax, int n, double xk[NPNT+1][NPNT+1], double wk[NPNT+1][NPNT+1])
{
  double c1, c2, s;
  int k;
  //
  c1 = 0.5 * (xmax + xmin);
  c2 = 0.5 * (xmax - xmin);
  s = 0.0;
  for (k=1; k<=n; ++k)
    s += c2 * wk[n][k] * func(c1 + c2 * xk[n][k]);
  //
  return s;
}
//
//--------------------------------------------------------------------------
//
// 被積分関数
//
double func(double x)
{
  return 1.0 / (1.0 + x * x);
}
/*
メインプログラム
 */
int main(void)
{
  double xk[NPNT+1][NPNT+1], wk[NPNT+1][NPNT+1];
  double xmin, xmax, s, se, error;
  int n, nmax = 16; // nmax : Gauss-Legendre公式の点数の最大値
  //
  point_weight_GaussLe(nmax, xk, wk); // Gauss-Legendre公式の標本点・重みを計算
  xmin = 0.2; // (xmin, xmax) : 積分区間
  xmax = 2.0;
  printf("xmin = %22.15e, xmax = %22.15e\n", xmin, xmax); 
  se = atan(xmax) - atan(xmin); // 積分の厳密値
  printf("I(exact) = %22.15e\n", se);
  printf("  n I(Gauss-Legendre rule) relative error\n");
  printf("-----------------------------------------\n");
  //
  for (n=4; n<=nmax; n+=2)
    {
      s = integral(func, xmin, xmax, n, xk, wk); // n点Gauss-Legendre公式による数値積分値
      error = fabs(s - se) / fabs(se);           // 数値積分値の相対誤差
      printf("%3d %22.15e %10.3e\n", n, s, error);
    }
  //
  return 0;
}