/* legendre.c  Legendre polynomials and quadrature                   */
/*                                                                   */
/*  P0(x)=1  P1(x)=x  Pn(x)=x*((2n-1)/n)*Pn-1(x) - ((n-1)/n)*Pn-2(x) */
/*              or  Pn+1(x)=x*)((2n+1)/(n+1))Pn(x)- (n/(n+1))Pn-1(x) */
/*  P2(x)=3/2 x^2 - 1/2   P3(x)=5/2 x^3 - 3/2 x                      */
/*  P4(x)=35/8 x^4 -30/8 x^2 + 3/8                                   */
/*  int -1 to 1 of f(x)dx = sum i=1,n  w[i] F(x[i])                  */
/*  w[i]= -2/((n+1)P'n(x[i])Pn+1(x[i]))                              */
/*  x[i]= roots of Pn(x)                                             */
/*                                                                   */
/*  n=2 x= +/- .577350269189625 w=1.000000000000000                  */
/*  n=3 x= +/- .774596669241483 w= .555555555555556                  */
/*            0.0                  .888888888888889                  */
/*  n=4 x= +/- .861136311594053 w= .347854845137454                  */
/*         +/- .339981043584856    .652145154862546                  */
/*  n=5 x= +/- .906179845938664    .236926885056189                  */
/*         +/- .538469310105683    .478628670499355                  */
/*            0.0                  .568888888888889                  */
/*                                                                   */
/* conversion/scaling  integral a to b of f(x)dx =                   */
/*             (b-a)/2 integral -1 to 1 of f((a+b+x*(b-a))/2)dx      */
/*                                                                   */
/* Gauss-Legendre quadrature, integrate f(x) from -1 to 1            */
/*   integral -1 to 1 f(x) dx = sum i=0,n-1 w[i]f(x[i])              */

#include <stdio.h>
#include <math.h>

#undef  abs
#define abs(x) ((x)<0.0?(-(x)):(x))

void LRoot(int n, double P[17][17], double R[17][17]);
double Evaluate(int n, double A[], double x);

int main(int argc, char *argv[])
{
  int N=16;          /* change to max polynomial size */
  double P[17][17];  /* coefficients of Pn(x)         */
  double DP[17][17]; /* derivatives of Pn(x)          */
  double R[17][17];  /* roots of Pn(x) the x[i]       */
  double W[17][17];  /* weights of Pn(x) the w[i]     */
  int i, j, n;
  double sumw=0.0;
  double a=1.0;
  double b=5.0;
  double exact=0.0;
  double approx=0.0;
  double err=0.0;
  
  printf("Legendre\n");
  /* generate family of Legendre polynomials Pn(x)    */
  /* as P[degree][coefficients]                       */
  /* build coefficients of Pn(x)                      */
  for(n=0; n<N+1; n++) for(i=0; i<N; i++) P[n][i]=0.0;
  P[0][0]=1.0;
  P[1][1]=1.0; /* start for recursion */
  for(n=2; n<N+1; n++)
  {
    for(i=0; i<=n; i++)
    {
      if(i<n-1) P[n][i]=-((double)(n-1)/(double)n)*P[n-2][i];
      if(i>0)   P[n][i]=P[n][i]+((double)(2*n-1)/(double)n)*P[n-1][i-1];
      printf("P[%d][%d]=%23.16E\n", n, i, P[n][i]);
    }
    printf("\n");
  }

  /* compute derivatives of Pn(x) */
  for(n=0; n<N; n++) for(i=0; i<N; i++) DP[n][i]=0.0;
  for(n=1; n<N; n++)
  {
    for(i=0; i<n; i++)
    {
      DP[n][i]=(double)(i+1)*P[n][i+1];
      printf("DP[%d][%d]=%23.16E", n, i, DP[n][i]);
    }
    printf("\n");
  }

  /* find roots of Pn(x), the x[i] */
  for(n=0; n<N; n++) for(i=0; i<N; i++) R[n][i]=0.0;
  R[1][0]=1.0;
  for(n=2; n<N; n++)
  {
    LRoot(n, P, R);
    for(i=0; i<n; i++)
    {
      printf("R[%d][%d]=%23.16E\n", n, i, R[n][i]);
    }
    printf("\n");
  }

  /* compute weights of Pn(x) */
  for(n=0; n<N; n++) for(i=0; i<N; i++) W[n][i]=0.0;
  W[1][0]=1.0;
  for(n=2; n<N; n++)
  {
    sumw = 0.0;
    for(i=0; i<n; i++)
    {
      W[n][i]=-2.0/((double)(n+1)*Evaluate(n, DP[n], R[n][i])*
                                  Evaluate(n+1, P[n+1], R[n][i]));
      sumw = sumw + W[n][i];
      printf("R[%d][%d]=%23.16E,   W[%d][%d]=%23.16E\n",
             n, i, R[n][i], n, i, W[n][i]);
    }
    printf("                            sum W=%23.16E\n",sumw);
  }
  printf("\n");

  /* integrate e^x from a to b */

  exact=exp(1.0)-exp(-1.0);
  printf("Gauss-Legendre quadrature of e^x from -1 to 1 = %23.16E\n",
         exact);
  for(n=2; n<N; n++)
  {
    approx=0.0;
    for(i=0; i<n; i++)
    {
      approx = approx + W[n][i]*exp(R[n][i]);
    }
    err = approx-exact;
    printf("n=%d,  approx=%23.16E, err=%23.16E\n", n, approx, err);
  }
  printf("\n");
  exact=exp(b)-exp(a);
  printf("Gauss-Legendre quadrature of e^x from %g to %g =%23.16E\n",
         a, b, exact);
  for(n=2; n<N; n++)
  {
    approx=0.0;
    for(i=0; i<n; i++)
    {
      approx = approx + W[n][i]*exp((a+b+R[n][i]*(b-a))/2.0);
    }
    approx = approx*(b-a)/2.0;
    err = approx-exact;
    printf("n=%d,  approx=%23.16E, err=%23.16E\n", n, approx, err);
  }
  printf("\n");
  return 0;
} /* end main of legendre.c */
  
void LRoot(int n, double P[17][17], double R[17][17])
{
  /* specialized for roots of Legendre polynomials */
  /* check first for zero root, reduce into T.     */
  /* T=initial polynomial, store root at R[n][i]   */
  double T[17];
  double DT[17];
  double rold, r;
  int i, j, k;

  i = 0;
  if(P[n][0]==0.0)
  { 
    R[n][i]=0.0;
    i++;
    for(j=1; j<=n; j++) T[j-1] = P[n][j]; /* reduce */
  }
  else
  {
    for(j=0; j<=n; j++) T[j] = P[n][j];
  }
  while(i<n)
  {
    for(j=0; j<n-i; j++)
    {
      DT[j]=(double)(j+1)*T[j+1];
    }
    /* do root Newton */
    rold=0.0;
    r=1.0;
    
    for(k=0; k<50; k++)
    {
      r = r - Evaluate(n-i, T, r)/Evaluate(n-i-1, DT, r);
      printf("k=%d, r=%23.16E, y=%23.16E, dy=%23.16E\n",
             k, r, Evaluate(n-i, T, r), Evaluate(n-i-1, DT, r));
      if(abs(r-rold)<1.0E-15) break;
      rold = r;
    }
    for(j=n-i; j>0; j--)
    {
      T[j-1]=T[j-1]+r*T[j];
    }
    for(j=1; j<=n-i; j++) T[j-1] = T[j]; /* reduce */
    R[n][i] = r;
    i++;
    r = -r;
    for(j=n-i; j>0; j--)
    {
      T[j-1]=T[j-1]+r*T[j];
    }
    for(j=1; j<=n-i; j++) T[j-1] = T[j]; /* reduce */
    R[n][i] = r;
    i++;
  }
} /* end LRoot */

double Evaluate(int n, double A[], double x)
{
  int i;
  double val=A[n]*x;
  for(i=n-1; i>0; i--)
  {
    val = (A[i] + val)*x;
  }
  return val+A[0];
} /* end Evaluate */

/* end legendre.c */