/* lagrange.c  simple definition and build 
 *             L_n,j(x) = product i=0,n i!=j (x-x_i)/(x_j-x_i)  
 *             let 0 < x < 1 for points x_0, x_1, ... , x_n
 *             let h = x_i - x_i-1,  x_i = i * h,  h = 1/n
 *             L_n(x) = sum j=0,n f(x_j) L_n,j(x)
 *             This L_n(x) fits f(x) with a n th order polynomial
 *
 *  generate L_n,j(x) for n=1,6  j=0,n
 */
 
#include <stdio.h>
#include <math.h>
#undef  max
#define max(a,b) ((a)>(b)?(a):(b))
#undef  min
#define min(a,b) ((a)<(b)?(a):(b))
#undef  abs
#define abs(a)   ((a)<0.0?(-(a)):(a))

int main(int argc, char *argv[])
{
  int N = 6;
  int i, j, k, n, m;
  double h, termx, termc;
  double p[20]; /* coefficients of polynomial */
  double q[20]; /* temp */
  double L[20]; /* accumulate fit of e^x */
  double x, fit, err, maxerr, rmserr;

  printf("lagrange.c running on 0<x<1 for N=%d \n", N);
  printf(" L_order,point  order is highest power of x, point 0..n \n");
  printf("      L_n,j(x) = product i=0,n i!=j (x-x_i)/(x_j-x_i)  \n");
  printf("      let 0 < x < 1 for points x_0, x_1, ... , x_n \n");
  printf("      let h = x_i - x_i-1,  x_i = i * h,  h = 1/n \n");
  printf("      L_n(x) = sum j=0,n f(x_j) L_n,j(x) \n");
  printf("      This L_n(x) fits f(x) with a n th order polynomial \n");

  for(n=1; n<=N; n++)
  {
    h = 1.0/n;
    for(j=0; j<=n; j++)
    {
      /* printf("computing terms of L_%d,%d \n", n, j); */
      p[0] = 1.0; /* identity polynomial to start product */
      for(i=1; i<=n; i++) p[i] = 0.0;
      for(i=0; i<=n; i++)
      {
        if(i!=j)
        {
          termx = 1.0/((j-i)*h); /* (x    )/(x_j-x_i) */
          termc = -h*i*termx;    /* ( -x_i)/(x_j-x_i) */
          /* do product p[] * (termx * x + termc) */
          for(k=0; k<=n; k++) q[k] = p[k] * termc;
          for(k=1; k<=n; k++) q[k] = q[k] + p[k-1]*termx;
          for(k=0; k<=n; k++) p[k] = q[k];
        }
      }
      /* print L_n,j(x) = p[...] */
      printf("L_%d,%d(x) = %g x^0 +\n", n, j, p[0]);
      for(i=1; i<=n; i++)
         printf("           %g x^%d + \n", p[i], i);
      printf("\n");
    } /* end loop on j */
    printf("\n");
  } /* end loop on n */

  /* repeat generation of L_5,* and fit e^x  0<x<1 */
  printf("\n");
  n = 5;
  h = 1.0/n;
  for(i=0; i<=n; i++) L[i] = 0.0; /* accumulate sum */
  for(j=0; j<=n; j++)
  {
    p[0] = 1.0; /* identity polynomial to start product */
    for(i=1; i<=n; i++) p[i] = 0.0;
    for(i=0; i<=n; i++)
    {
      if(i!=j)
      {
        termx = 1.0/((j-i)*h); /* (x    )/(x_j-x_i) */
        termc = -h*i*termx;    /* ( -x_i)/(x_j-x_i) */
        /* do product p[] * (termx * x + termc) */
        for(k=0; k<=n; k++) q[k] = p[k] * termc;
        for(k=1; k<=n; k++) q[k] = q[k] + p[k-1]*termx;
        for(k=0; k<=n; k++) p[k] = q[k];
      }
    }
    for(k=0; k<=n; k++) L[k] = L[k]+exp(j*h)*p[k];
  } /* end loop on j */
  /* check results at computed points and intermediate points */
  printf("check Lagrange 5 point fit of exp(x) in 0 <= x <=1 \n");
  printf("L_%d(x) = %g x^0 +\n", n, L[0]);
  for(i=1; i<=n; i++)
      printf("        %g x^%d + \n", L[i], i);
  printf("\n");
  m = 2*n;
  h = 1.0/m;
  maxerr = 0.0;
  rmserr = 0.0;
  for(k=0; k<=m; k++)
  {
    x = k*h;
    fit = L[n]*x; /* evaluate L at x */
    for(i=n-1; i>0; i--) fit = (L[i]+fit)*x;
    fit = fit + L[0];
    err = fit - exp(x);
    maxerr = max(maxerr, abs(err));
    rmserr = rmserr + err*err;
    printf("check exp(%g)=%g, fit=%g, err=%g \n", x, exp(x), fit, err);
  }
  rmserr = sqrt(rmserr/(double)(m+1));
  printf("maxerr=%g, rmserr=%g \n", maxerr, rmserr);

  /* repeat generation of L_6,* and fit e^x  0<x<1 */
  printf("\n");
  n = 6;
  h = 1.0/n;
  for(i=0; i<=n; i++) L[i] = 0.0; /* accumulate sum */
  for(j=0; j<=n; j++)
  {
    p[0] = 1.0; /* identity polynomial to start product */
    for(i=1; i<=n; i++) p[i] = 0.0;
    for(i=0; i<=n; i++)
    {
      if(i!=j)
      {
        termx = 1.0/((j-i)*h); /* (x    )/(x_j-x_i) */
        termc = -h*i*termx;    /* ( -x_i)/(x_j-x_i) */
        /* do product p[] * (termx * x + termc) */
        for(k=0; k<=n; k++) q[k] = p[k] * termc;
        for(k=1; k<=n; k++) q[k] = q[k] + p[k-1]*termx;
        for(k=0; k<=n; k++) p[k] = q[k];
      }
    }
    for(k=0; k<=n; k++) L[k] = L[k]+exp(j*h)*p[k];
  } /* end loop on j */
  /* check results at computed points and intermediate points */
  printf("check Lagrange 6 point fit of exp(x) in 0 <= x <=1 \n");
  printf("L_%d(x) = %g x^0 +\n", n, L[0]);
  for(i=1; i<=n; i++)
      printf("        %g x^%d + \n", L[i], i);
  printf("\n");
  m = 2*n;
  h = 1.0/m;
  maxerr = 0.0;
  rmserr = 0.0;
  for(k=0; k<=m; k++)
  {
    x = k*h;
    fit = L[n]*x; /* evaluate L at x */
    for(i=n-1; i>0; i--) fit = (L[i]+fit)*x;
    fit = fit + L[0];
    err = fit - exp(x);
    maxerr = max(maxerr, abs(err));
    rmserr = rmserr + err*err;
    printf("check exp(%g)=%g, fit=%g, err=%g \n", x, exp(x), fit, err);
  }
  rmserr = sqrt(rmserr/(double)(m+1));
  printf("maxerr=%g, rmserr=%g \n", maxerr, rmserr);

  printf("\n end lagrange.c \n");
  return 0;
} /* end lagrange.c */