/*
  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;
}