/* poly_fit.c  given n data points (xi,yi)        */
/*  exactly fit a n-1 order polynomial            */
/*  y = c0 + c1*x + c2*x^2 + ... + cn-1 * x^(n-1) */
/*                                                */
/*  create system of equations:                   */
/*                                                */
/*  | 1  x1  x1^2  ... x1^n |   | c0   |   | y1 | */
/*  | 1  x2  x2^2  ... x2^n |   | c2   |   | y2 | */
/*  |              ...      | * | ...  | = |... | */
/*  |              ...      |   | ...  |   |... | */
/*  | 1  xn  xn^2  ... xn^n |   | cn-1 |   | yn | */
/*                                                */
/*  solve for c0, c1, ... , cn-1                  */

#include <stdio.h>
#include "simeq.h"
double peval(int n, double x, double c[]);

int main(int argc, char *argv[])
{
  int n = 4;    /* number of equations */
  double X[4] = { 0.0, 0.25, 0.5, 1.0 }; /* (x,y) points to fit */
  double Y[4] = { 1.0, 0.0,  0.0, 1.0 };
  double A[16]; /* matrix */
  double C[4];  /* coefficients */
  double xx, yy, err;
  int i, j, k;

  printf("poly_fit.c running \n");
  for(i=0; i<n; i++) printf("data %d, (x=%g,y=%g) \n",i, X[i], Y[i]);
  printf("\n");

  for(i=0; i<n; i++) /* build matrix */
  {
    /* Y[i] is already Y[i] */
    xx = 1.0;
    for(j=0; j<n; j++)
    {
      A[i*n+j] = xx;
      xx = xx * X[i];
    }
  }  
  /* solve system of equations */
  simeq(n, A, Y, C);
  printf("fit equation,  y = %g + \n", C[0]);
  for(i=1; i<n; i++) printf("                   %g * x^%d + \n", C[i], i);

  printf("\n check polynomial against given points \n");
  for(i=0; i<n; i++)
  {
    yy = peval(n-1, X[i], C);
    err = Y[i]-yy;
    printf("evaluate at x=%g, yy=%g, Y[%d]=%g, error=%g \n",
           X[i], yy, i, Y[i], err);
  }
  printf("\n poly_fit.c ends \n");

  return 0;
} /* end main */

/* peval.c Horners method for evaluating a polynomial */
double peval(int n, double x, double c[])
{ /* Horners method for evaluating a polynomial   */
  /* an nth order polynomial has n+1 coefficients */
  /* stored in c[0] through c[n]                  */
  /* y = c[0] + c[1]*x + c[2]*x^2 +...+ c[n]*x^n  */
  int i;
  double y = c[n]*x;
  if(n<=0) return c[0];
  for(i=n-1; i>0; i--) y = (c[i]+y)*x;
  return y+c[0];
} /* end peval */

/* end poly_fit.c */