/*************************    Program 2.2    ****************************/
/*                                                                      */
/************************************************************************/
/* Please Note:                                                         */
/*                                                                      */
/* (1) This computer program is written by Tao Pang in conjunction with */
/*     his book, "An Introduction to Computational Physics," published  */
/*     by Cambridge University Press in 1997.                           */
/*                                                                      */
/* (2) No warranties, express or implied, are made for this program.    */
/*                                                                      */
/************************************************************************/

#define NMAX 21

void pfit (n,m,x,f,a,u)
/* Subroutine generating orthonormal polynomials U(M,N) up to
   (M-1)th order and coefficients A(M), for the least squares
   approximation of the function F(N) at X(N). Other variables
   used: G(K) for g_k, H(K) for h_k, S(K) for <u_k|u_k>.
   Copyright (c) Tao Pang 1997. */

int n,m;
double x[],f[];
double a[],u[NMAX][NMAX];
{
int i,j;
double stmp;
double s[NMAX],g[NMAX],h[NMAX];

if (n > NMAX)
  {
  printf ("Too many points\n");
  exit (1);
  }
if (m > n)
  {
  printf ("Order too high\n");
  exit (2);
  }

/* Set up the zeroth order polynomial u_0 */

s[0] = 0;  g[0] = 0;  a[0] = 0;

for (i = 0; i < n; ++i)
  {
  u[0][i] = 1;
  stmp   = u[0][i]*u[0][i];
  s[0]   = s[0]+stmp;
  g[0]   = g[0]+x[i]*stmp;
  a[0]   = a[0]+u[0][i]*f[i];
  }
g[0] = g[0]/s[0];
a[0] = a[0]/s[0];
h[0] = 0;

/* Set up the first order polynomial u_1 */

s[1] = 0;  g[1] = 0;  h[1] = 0; a[1] = 0;

for (i = 0; i < n; ++i)
  {
  u[1][i] = x[i]*u[0][i]-g[0]*u[0][i];
  s[1]   = s[1]+u[1][i]*u[1][i];
  g[1]   = g[1]+x[i]*u[1][i]*u[1][i];
  h[1]   = h[1]+x[i]*u[1][i]*u[0][i];
  a[1]   = a[1]+u[1][i]*f[i];
  }
g[1] = g[1]/s[1];
a[1] = a[1]/s[1];
h[1] = h[1]/s[0];

/* Higher order polynomials u_k from the recursive relation */

if (m >= 3)
  {
  for (i = 1; i < m-1; ++i)
    {
    s[i+1] = 0;  g[i+1] = 0;  h[i+1] = 0; a[i+1] = 0;
    for (j = 0; j < n; ++j)
      {
      u[i+1][j] = x[j]*u[i][j]-g[i]*u[i][j]-h[i]*u[i-1][j];
      s[i+1]   = s[i+1]+u[i+1][j]*u[i+1][j];
      g[i+1]   = g[i+1]+x[j]*u[i+1][j]*u[i+1][j];
      h[i+1]   = h[i+1]+x[j]*u[i+1][j]*u[i][j];
      a[i+1]   = a[i+1]+u[i+1][j]*f[j];
      }
    g[i+1] = g[i+1]/s[i+1];
    a[i+1] = a[i+1]/s[i+1];
    h[i+1] = h[i+1]/s[i];
    }
  }
}