// test simple  polyfit.c  only  y = f(x)  given x,y pairs
#include "polyfit.h"
#include "polyval.h"
#include "simeq.h"
#include <stdio.h>
#include <stdlib.h>

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

int main(int argc, char argv[])
{
  #define n 10
  #define pwr 3
  double x[n];
  double y[n];
  double yc[n]; // for checking
  double p[n+1];
  double err;
  int i, j;

  printf("Testing polyval and polyfit.c \n");
  printf("given x,y   i=0 to n-1  y[i] = f(x[i]) \n");
  printf("find p[] polynomial coefficients that fit  f \n");
  printf("y[i] = p[0] + p[1]*x[i] + p[2]*x[i]^2 + p[3]*x[i]^3 ... +p[pwr]*x[i]^pwr \n");
  for(j=0; j<=pwr; j++)
  {
    p[j] = 1.2*(double)(j+1); // for computing y
    printf("j=%d  p[j]=%f \n", j, p[j]);
  }
  for(i=0; i<n; i++)
  {
     x[i] = 0.25*(double)(i+1);
     printf("i=%d x[i]=%f \n", i, x[i]);
  }
  for(i=0; i<n; i++)
  {
    y[i] = polyval(pwr, p, x[i]);
    yc[i] = y[i];
    printf("i=%d  x[i]=%f  y[i]=%f \n", i, x[i], y[i]);
  }
  for(j=0; j<=pwr; j++) p[j] = -1.0;
  
  polyfit(pwr, n, x, y, p);
  printf("computed values \n");
  for(j=0; j<=pwr; j++)
  {
    printf("j=%d  p[j]=%f \n", j, p[j]);
  }
  for(i=0; i<n; i++)
  {
    y[i] = polyval(pwr, p, x[i]);
    err = abs(y[i]-yc[i]);
    printf("i=%d  x[i]=%f  y[i]=%f err=%e \n", i, x[i], y[i], err);
  }

  return 0;
} // end test polyfit.c