{***************************************************
*      Program to demonstrate least squares        *
*         polynomial fitting subroutine            *
* ------------------------------------------------ *
* Reference: BASIC Scientific Subroutines, Vol. II *
* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].*
*                                                  *
*           Pascal version by J-P Moreau, Paris.   *
*                    (www.jpmoreau.fr)             *
* ------------------------------------------------ *
* SAMPLE RUN:                                      *
*                                                  *
* LEAST SQUARES POLYNOMIAL FITTING                 *
*                                                  *
* What is the order of the fit      : 4            *
* What is the error reduction factor: 0            *
* How many data pooints are there   : 11           *
*                                                  *
* Input the data points as prompted:               *
*                                                  *
*  1   X, Y = 1,1                                  *
*  2   X, Y = 2,2                                  *
*  3   X, Y = 3,3                                  *
*  4   X, Y = 4,4                                  *
*  5   X, Y = 5,5                                  *
*  6   X, Y = 6,6                                  *
*  7   X, Y = 7,7                                  *
*  8   X, Y = 8,8                                  *
*  9   X, Y = 9,9                                  *
*  10   X, Y = 0,0                                 *
*  11   X, Y = 1,1                                 *
*                                                  *
* Coefficients are:                                *
*  0    -0.000000                                  *
*  1    1                                          *
*  2    -0.000000                                  *
*  3    0                                          *
*  4    -0.000000                                  *
*                                                  *
* Standard deviation =  0.00000000                 *
*                                                  *
***************************************************}
PROGRAM LSQRPOLY;
Uses WinCrt;

TYPE    TAB = Array[0..25] of DOUBLE;

VAR
        i,l,m,n : integer;
        dd,e1,vv : DOUBLE;

        x,y,v,a,b,c,d,c2,e,f : TAB;



{****************************************************************
*        LEAST SQUARES POLYNOMIAL FITTING PROCEDURE             *
* ------------------------------------------------------------- *
* This program least squares fits a polynomial to input data.   *
* Forsythe orthogonal polynomials are used in the fitting.      *
* The number of data points is n.                               *
* The data is input to the subroutine in x[i], y[i] pairs.      *
* The coefficients are returned in c[i],                        *
* the smoothed data is returned in v[i],                        *
* the order of the fit is specified by m.                       *
* The standard deviation of the fit is returned in d.           *
* There are two options available by use of the parameter e:    *
*  1. if e = 0, the fit is to order m,                          *
*  2. if e > 0, the order of fit increases towards m, but will  *
*     stop if the relative standard deviation does not decrease *
*     by more than e between successive fits.                   *
* The order of the fit then obtained is l.                      *
****************************************************************}
PROCEDURE LS_POLY;
Label 10,15,20,30,50,fin;
Var i,l2,n1 : integer;
    a1,a2,b1,b2,c1,d1,f1,f2,v1,v2,w : DOUBLE;
Begin
  n1 := m + 1;
  v1 := 1e7;
  {Initialize the arrays}
  FOR i := 1 TO n1 DO
  begin
    a[i] := 0; b[i] := 0; f[i] := 0
  end;
  FOR i := 1 TO n DO
  begin
    v[i] := 0; d[i] := 0
  end;
  d1 := SQRT(n); w := d1;
  FOR i := 1 TO n DO
  begin
    e[i] := 1 / w
  end;
  f1 := d1; a1 := 0;
  FOR i := 1 TO n DO
  begin
    a1 := a1 + x[i] * e[i] * e[i]
  end;
  c1 := 0;
  FOR i := 1 TO n DO
  begin
    c1 := c1 + y[i] * e[i]
  end;
  b[1] := 1 / f1; f[1] := b[1] * c1;
  FOR i := 1 TO n DO
  begin
    v[i] := v[i] + e[i] * c1
  end;
  m := 1;
10: {Save latest results}
  FOR i := 1 TO l DO c2[i] := c[i];
  l2 := l; v2 := v1; f2 := f1; a2 := a1; f1 := 0;
  FOR i := 1 TO n DO
  begin
    b1 := e[i];
    e[i] := (x[i] - a2) * e[i] - f2 * d[i];
    d[i] := b1;
    f1 := f1 + e[i] * e[i]
  end;
  f1 := SQRT(f1);
  FOR i := 1 TO n DO e[i] := e[i] / f1;
  a1 := 0;
  FOR i := 1 TO n DO a1 := a1 + x[i] * e[i] * e[i];
  c1 := 0;
  FOR i := 1 TO n DO c1 := c1 + e[i] * y[i]; 
  m := m + 1; i := 0;
15: l := m - i; b2 := b[l]; d1 := 0;
  IF l > 1 THEN d1 := b[l - 1];
  d1 := d1 - a2 * b[l] - f2 * a[l];
  b[l] := d1 / f1; a[l] := b2; i := i + 1;
  IF i <> m THEN GOTO 15;
  FOR i := 1 TO n DO v[i] := v[i] + e[i] * c1; 
  FOR i := 1 TO n1 DO
  begin
    f[i] := f[i] + b[i] * c1;
    c[i] := f[i]
  end;
  vv := 0;
  FOR i := 1 TO n DO vv := vv + (v[i] - y[i]) * (v[i] - y[i]); 
  {Note the division is by the number of degrees of freedom}
  vv := SQRT(vv / (n - l - 1)); l := m;
  IF e1 = 0 THEN GOTO 20;
  {Test for minimal improvement}
  IF ABS(v1 - vv) / vv < e1 THEN GOTO 50;
  {If error is larger, quit}
  IF e1 * vv > e1 * v1 THEN GOTO 50;
  v1 := vv;
20: IF m = n1 THEN GOTO 30;
  GOTO 10;
30: {Shift the c[i] down, so c(0) is the constant term}
  FOR i := 1 TO l DO c[i - 1] := c[i];
  c[l] := 0;
  {l is the order of the polynomial fitted}
  l := l - 1; dd := vv;
  goto fin;
50: {Aborted sequence, recover last values}
  l := l2; vv := v2;
  FOR i := 1 TO l DO c[i] := c2[i];
  GOTO 30;
fin: End;


BEGIN
  Clrscr;
  Writeln(' LEAST SQUARES POLYNOMIAL FITTING');
  Writeln;
  write(' What is the order of the fit      : '); Read(m);
  write(' What is the error reduction factor: '); Read(e1);
  write(' How many data pooints are there   : '); Read(n);
  Writeln;
  Writeln(' Input the data points as prompted:');
  Writeln;
  FOR i := 1 TO n DO
  begin
    Write(' ',i,'  X  Y = '); Read(x[i], y[i])
  end;
  Writeln;

  LS_POLY;

  Writeln(' Coefficients are:');
  FOR i := 0 TO l DO Writeln('  ',i,'  ',c[i]:9:6);
  Writeln;
  Writeln(' Standard deviation: ',dd:9:8);
  writeln;
  ReadKey; DoneWinCrt
END.

{End lsqrpoly.pas}