/*
  compute an integral over a finite interval 
  by the Legendre-Gauss quadrature formula
  
  copyright : Hidenori Ogata, 29 March 2004, ver.1
  
  (usage)
  Write 
  "point_weight_LeGauss(n, xk, wk);
  .....
  LeGauss(func, xmin, xmax, n, xk, wk, &itgl); " 
  in the main program. The parameters are as follows.
  
  func : the integrand function, which is double-precision valued 
  with one variable in double precision
  xmin double : the lower bound of integration
  xmax double : the upper bound of integration
  n    int    : the number of the abscissae of the Legendre-Gauss quadrature formula
  
  xk   double[] : the abscissae of the quadrature formula 
  wk   double[] : the weights of the quadrature formula
  itgl double * : the computed value of the integral
*/
// #include <stdio.h>
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
#define NPNT 101
//
// Gauss-Legendre数値積分公式のクラス
//
class GLquad {
  double xk[NPNT+1][NPNT+1], wk[NPNT+1][NPNT+1];
  int nmax;
public:
  GLquad(int nmax0);
  // 数値積分
  // func : 被積分関数
  // (xmin, xmax) : 積分区間
  // n : Gauss-Legendre公式の点数
  //
  double integral(double (*func)(double), double xmin, double xmax, int n){
    return LeGauss(func, xmin, xmax, n);
  }
  double LeGauss(double (*func)(double), double xmin, double xmax, int n);
private:
  void point_weight_LeGauss(double xk[NPNT+1][NPNT+1], double wk[NPNT+1][NPNT+1]);
  double Pn_Legendre(int n, double x);
  double dPn_Legendre(int n, double x);
};
//
// コンストラクタ
// nmax0 : Gauss-Legendre公式の点数の最大値
//
GLquad::GLquad(int nmax0)
{
  nmax = nmax0;
  point_weight_LeGauss(xk, wk);
}
//
//void point_weight_LeGauss(int, double double *xk, double *wk);
//void LeGauss(double (*func)(double), double, double, int, 
//             double [], double [], double *);
//
/* 

　数値積分公式の標本点（Legendre多項式の零点），重みを計算．

  compute the abscissae and the weights of the Legendre functions 
  by Newton's method
  (input)
  n  int      : the order n of the Legendre polynomial Pn(x)
                n must be <= NPNT-1.
  (output)
  xk double[] : the zeros of Pn(x)
  wk double[] : the weights of the Legendre-Gauss quadrature
*/
void GLquad::point_weight_LeGauss(double xk[NPNT+1][NPNT+1], double wk[NPNT+1][NPNT+1])
{
  int 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 (int 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) 
	    {
	      cout << "k=" << k << " : not convergent" << endl;
	      //	      printf("k=%3d: not convergent\n", k);
	      //           exit(0);
	    }
	}
      //
      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;
	}
    }
}
//
/*
  compute the integral over a finite interval [xmin, xmax]
  by the (n-point) Legendre-Gauss quadrature formula
  (input)
  func : the integrand function, which is double-precision valued 
         with one variable in double precision
  xmin double : the lower bound of integration
  xmax double : the upper bound of integration
  n    int    : the number of the abscissae of the Legendre-Gauss quadrature formula
  xk   double[] : the abscissae of the quadrature formula 
  wk   double[] : the weights of the quadrature formula
  (output)
  itgl double * : the computed value of the integral
*/
double GLquad::LeGauss(double (*func)(double), double xmin, double xmax, int n)
{
    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;
}
//
/*
　Legendre多項式

  Legendre polynomial of order n
  (input)
  n int : the order n of the Legendre polynomial Pn(x)
  x int : the variable x 
*/
double GLquad::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多項式の導関数

  derivative of the Legendre polynomial of order n: Pn'(x)
  (input)
  n int : the order n of the Legendre polynomial
  x int : the variable x 
*/
double GLquad::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;
}
//
//--------------------------------------------------------------------------
//
// 被積分関数
//
double func(double x)
{
  return 1.0 / (1.0 + x * x);
}
//
int main(void)
{
  int nmax0 = 50; // Gauss-Legendre公式の点数の最大値
  GLquad GLquad1(nmax0); // Gauss-Legendre数値積分公式のクラス
  //
  double xmin, xmax, s, se, error;
  xmin = 0.2; // (xmin, xmax) : 積分区間
  xmax = 2.0;
  cout << "xmin = " << xmin << ", xmax = " << xmax << endl;
  se = atan(xmax) - atan(xmin); // 積分の厳密値
  std::cout << std::setiosflags(std::ios::scientific);
  cout << setprecision(15);
  cout << "I(exact) = " << se << endl;
  cout << "  n I(Gauss-Legendre rule) relative error" << endl;
  cout << "-----------------------------------------" << endl;
  //
  int nmax = 16;
  for (int n=4; n<=nmax; n+=2)
    {
      s = GLquad1.integral(func, xmin, xmax, n); // n点Gauss-Legendre公式による数値積分値
      error = fabs(s - se) / fabs(se);           // 数値積分値の相対誤差
      cout << setw(3) << right << n << " ";
      std::cout << std::setiosflags(std::ios::scientific);
      cout << setprecision(15);
      cout << s << " ";
      cout << setprecision(3);
      cout << error << endl;
    }
  //
  return 0;
}