// polyfit.c    polyfit(pwr, n, x, y, p)   p computed
// i=0, i<n     y[i] = p[0] + p[1]*x[i] + p[2]*x[i]^2 + ... x[i]^pwr
// The problem is stated as follows :
// Given measured data for values of Y based on values of X. e.g.
//
//  i    y_actual         x
// ---   --------       -----
//  0     32.5           1.0
//  1      7.2           2.0
//  2      6.9           3.2
//  3      8.4           2.7
//  4      9.6           1.4
//
// Find p[] such that  
// y_approximate =  p[0] + p[1]*x[i] + p[2]*x[i]^2 ... +p[pwr]*x[i]^pwr
// and such that the sum of (y_actual - y_approximate) squared is minimized.
// To find the coefficients  p  solve the linear system of equations
//
// | SUM(1*1)   SUM(1*x)   SUM(1*x^2)   SUM(1*x^3)   | | p[0] | | SUM(1*y)   |
// | SUM(x*1)   SUM(x*x)   SUM(x*x^2)   SUM(x*x^3)   |x| p[1] |=| SUM(x*y)   |
// | SUM(x^2*1) SUM(x^2*x) SUM(x^2*x^2) SUM(x^2*x^2) | | p[2] | | SUM(x^2*y) |
// | SUM(x^3*1) SUM(x^3*x) SUM(x^3*x^2) SUM(x^3*X^2) | | p[3] | | SUM(x^3*y) |
//   solve  |A|*|p|=|Y| using simeq
  
#include  "simeq.h" // compile with  simeq.c

void polyfit(int pwr, int n, double x[], double y[], double p[])
{
  int pn = pwr+1;
  double A[pn*pn];
  double Y[pn];
  double xpwr[pn];
  int i, j, k;

  for(i=0; i<pn; i++) // zero A and Y
  {
    for(j=0; j<pn; j++)
    {
      A[i*pn+j] = 0.0;
    }
    Y[i] = 0.0;
  }
  for(k=0; k<n; k++)
  {
    xpwr[0] = 1.0;
    for(i=1; i<pn; i++) xpwr[i] = xpwr[i-1]*x[k];
    for(i=0; i<pn; i++)
    {
      for(j=0; j<pn; j++)
      {
        A[i*pn+j] = A[i*pn+j] + xpwr[i]*xpwr[j];
      }
      Y[i] = Y[i] + y[k]*xpwr[i];
    }
  }

  simeq(pn, A, Y, p);
  return;
} // end polyfit.c