/* legendre_fit.c                                                   */
/*    fit  f(x) = a_0 g0 P0(x) + a_1 g1 P1(x) + ... a_n gn Pn(x)    */
/*         where Pi(x) as P(i,x) is a Legendre polynomial           */
/*         the weight function is 1.0                               */
/*         the integral range is -1 to 1                            */
/*         a_i = integral -1 to 1 of f(x) Pi(x) dx                 */
/*         gi = (2*i+1)/2                                           */
/* test on e^x = 1 + x + 1/2 x^2 + 1/6 x^3 + 1/24 x^4 + 1/120 x^5   */
/* (just using Pi(x) does not work. (2i+1)/2 could be the weight)   */
/* using Phi(i,x) = sqrt((2i+1)/2) P(i,x) for integration and       */
/* evaluation works also.                                           */

/*    first few Legendre polynomials:                               */
/*    P(0,x) = 1.0                                                  */
/*    P(1,x) = x                                                    */
/*    P(2,x) = 1/2 (3x^2 - 1)                                       */
/*    P(3,x) = 1/2 (5x^3-3.0*x)                                     */
/*    P(4,x) = 1/8 (35x^4-30x^2+3)                                  */

#include <stdio.h>
#include <math.h>
#undef  max
#define max(a,b) ((a)>(b)?(a):(b))
#undef  abs
#define abs(a)   ((a)<0.0?(-(a)):(a))

static double P(int i, double x);
static double Phi(int i, double x);

int main(int argc, char * argv[])
{
  double a[10];
  double exact[10];
  int n = 4;
  int nn = 2048; /* for trapezoid integration */
  int m = 10; /* for check */
  double err, maxerr, rmserr;
  double xmin = -1.0;
  double xmax = 1.0;
  double h, hh, x, val;
  int i, j;
  
  printf("legendre_fit.c running \n");
  printf("fit  f(x) = a_0 g0 P0(x) + a_1 g1 P1(x) + ... a_n gn Pn(x)     \n");
  printf("where Pi(x) as P(i,x) is a Legendre polynomial                 \n");
  printf("the weight function is 1.0                                     \n");
  printf("the integral range is -1 to 1                                  \n");
  printf("a_i = integral -1 to 1 of f(x) Pi(x) dx                       \n");
  printf("gi = (2*i+1)/2                                                 \n");
  printf("test on e^x = 1 + x + 1/2 x^2 + 1/6 x^3 + 1/24 x^4 + 1/120 x^5 \n");
  printf("using Phi(i,x) = sqrt((2i+1)/2 P(i,x) both places works also   \n");
  printf("\n");
  
  exact[0] = 1.0;
  for(i=1; i<10; i++)
  {
    exact[i] = exact[i-1]/(double)i;
  }
  
  h = (xmax-xmin)/(double)nn;
  hh = (double)nn;
  printf("using P0 thru P4, trapezoidal nn=%d \n", nn);
  for(i=0; i<=n; i++)
  {
    a[i] = exp(xmin)*P(i,xmin)/2.0; /* trapezoidal */
    a[i] = a[i]+exp(xmax)*P(i,xmax)/2.0;
    for(j=1; j<nn; j++)
    {
      x = xmin + (double)j*h;
      a[i] = a[i] + exp(x)*P(i,x);
      /* printf("x=%8.5f, exp(x)=%8.5f, P(%d,x)=%8.5f \n",
             x, exp(x), i, P(i,x)); */
    }
    a[i] = a[i]/hh;
    a[i] = a[i]*(xmax-xmin);
    printf("a[%d]=%8.5f, exp coefficient=%8.5f \n", i, a[i], exact[i]);
  }

  printf("P(i,x) with no correction does not work \n");
  h = (xmax-xmin)/(double)(m-1);
  maxerr = 0.0;
  rmserr = 0.0;
  for(j=0; j<m; j++)
  {
    x = xmin + (double)j*h;
    val = 0.0;
    for(i=0; i<=n; i++)
    {
      val = val + a[i]*P(i,x);
    }
    err = val-exp(x);
    maxerr = max(maxerr, abs(err));
    rmserr = rmserr + err*err;
    printf("x=%8.5f, val=%8.5f, exp(x)=%8.5f, err=%8.5f \n",
           x, val, exp(x), err);
  }
  printf("                      maxerr=%8.5f, rmserr=%8.5f \n\n",
         maxerr, sqrt(rmserr/(double)m));

  printf("using evaluation factor (2i+1)/2 a_i P(i,x) works \n");
  h = (xmax-xmin)/(double)(m-1);
  maxerr = 0.0;
  rmserr = 0.0;
  for(j=0; j<m; j++)
  {
    x = xmin + (double)j*h;
    val = 0.0;
    for(i=0; i<=n; i++)
    {
      val = val + ((double)(2*i+1)/2.0)*a[i]*P(i,x);
    }
    err = val-exp(x);
    maxerr = max(maxerr, abs(err));
    rmserr = rmserr + err*err;
    printf("x=%8.5f, val=%8.5f, exp(x)=%8.5f, err=%8.5f \n",
           x, val, exp(x), err);
  }
  printf("                      maxerr=%8.5f, rmserr=%8.5f \n\n",
         maxerr, sqrt(rmserr/(double)m));
  
  h = (xmax-xmin)/(double)nn;
  hh = (double)nn;
  printf("using Phi(i,x)=sqrt((2i+1)/2) P(i,x) for integration \n");
  for(i=0; i<=n; i++)
  {
    a[i] = exp(xmin)*Phi(i,xmin)/2.0;
    a[i] = a[i]+exp(xmax)*Phi(i,xmax)/2.0;
    for(j=1; j<nn; j++)
    {
      x = xmin + (double)j*h;
      a[i] = a[i] + exp(x)*Phi(i,x);
      /* printf("x=%8.5f, exp(x)=%8.5f, Phi(%d,x)=%8.5f \n",
             x, exp(x), i, Phi(i,x)); */
    }
    a[i] = a[i]/hh;
    a[i] = a[i]*(xmax-xmin);
    printf("a[%d]=%8.5f, exp coefficient=%8.5f \n", i, a[i], exact[i]);
  }

  printf("using Phi(i,x)=sqrt((2i+1)/2) P(i,x) for evaluation \n");
  h = (xmax-xmin)/(double)(m-1);
  maxerr = 0.0;
  rmserr = 0.0;
  for(j=0; j<m; j++)
  {
    x = xmin + (double)j*h;
    val = 0.0;
    for(i=0; i<=n; i++)
    {
      val = val + a[i]*Phi(i,x);
    }
    err = val-exp(x);
    maxerr = max(maxerr, abs(err));
    rmserr = rmserr + err*err;
    printf("x=%8.5f, val=%8.5f, exp(x)=%8.5f, err=%8.5f \n",
           x, val, exp(x), err);
  }
  printf("                      maxerr=%8.5f, rmserr=%8.5f \n\n",
         maxerr, sqrt(rmserr/(double)m));

  nn = 8192;
  n = 6;
  printf("using P0 thru P6, trapezoidal nn=%d \n", nn);
  h = (xmax-xmin)/(double)nn;
  hh = (double)nn;
  printf("using Phi(i,x)=sqrt((2i+1)/2) P(i,x) for integration \n");
  for(i=0; i<=n; i++)
  {
    a[i] = exp(xmin)*P(i,xmin)/2.0;
    a[i] = a[i]+exp(xmax)*P(i,xmax)/2.0;
    for(j=1; j<nn; j++)
    {
      x = xmin + (double)j*h;
      a[i] = a[i] + exp(x)*P(i,x);
      /* printf("x=%9.6f, exp(x)=%9.6f, P(%d,x)=%9.6f \n",
             x, exp(x), i, P(i,x)); */
    }
    a[i] = a[i]/hh;
    a[i] = a[i]*(xmax-xmin);
    printf("a[%d]=%9.6f, exp coefficient=%9.6f \n", i, a[i], exact[i]);
  }

  printf("using (2i+1)/2 P(i,x) for evaluation \n");
  h = (xmax-xmin)/(double)(m-1);
  maxerr = 0.0;
  rmserr = 0.0;
  for(j=0; j<m; j++)
  {
    x = xmin + (double)j*h;
    val = 0.0;
    for(i=0; i<=n; i++)
    {
      val = val + ((double)(2*i+1)/2.0)*a[i]*P(i,x);
    }
    err = val-exp(x);
    maxerr = max(maxerr, abs(err));
    rmserr = rmserr + err*err;
    printf("x=%9.6f, val=%9.6f, exp(x)=%9.6f, err=%9.7f \n",
           x, val, exp(x), err);
  }
  printf("                      maxerr=%10.7f, rmserr=%10.7f \n\n",
         maxerr, sqrt(rmserr/(double)m));

  printf("\nlegendre_fit.c finished \n");  
  return 0;
}

static double P(int i, double x)
{
  switch(i)
  {
    case 0: return 1.0;
    case 1: return x;
    case 2: return (3.0*x*x-1.0)/2.0;
    case 3: return (5.0*x*x*x-3.0*x)/2.0;
    case 4: return (35.0*x*x*x*x-30.0*x*x+3.0)/8.0;
    case 5: return 1.875*x-8.75*x*x*x+7.875*x*x*x*x*x;
    case 6: return -0.3125+6.5625*x*x-19.6875*x*x*x*x+14.4375*x*x*x*x*x*x;
    default: return 0.0;
  }
}

static double Phi(int i, double x)
{
  switch(i)
  {
    case 0: return sqrt(1.0/2.0)*1.0;
    case 1: return sqrt(3.0/2.0)*x;
    case 2: return sqrt(5.0/2.0)*(3.0*x*x-1.0)/2.0;
    case 3: return sqrt(7.0/2.0)*(5.0*x*x*x-3.0*x)/2.0;
    case 4: return sqrt(9.0/2.0)*(35.0*x*x*x*x-30.0*x*x+3.0)/8.0;
    default: return 0.0;
  }
}