{█▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀█}
{█                                                       █}
{█      Virtual Pascal Examples. Version 1.10            █}
{█      MathPak Least Squares example                    █}
{█      ─────────────────────────────────────────────────█}
{█      Copyright (C) 1995-96 fPrint UK Ltd              █}
{█                                                       █}
{▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀}

program LeastSq;
{$M 32768}

uses
  MP87Lite;

{$PMTYPE VIO}
{$IFDEF VPDEMO}
  {&Dynamic VP11Demo.Lib}
{$ENDIF}

  { This is an example program to show how MATHPAK 87 procedures can
    be used to program a polynomial least-squares data fitting program.
    See the Appendix of your manual for a description of the program.  }

Const
  MaxDeg    = 10;    { Maximum degree of fitting polynomial }
  MaxDegx2  = 20;
  MaxNumPoints = 500;

Type
  Vector      = Array[1..MaxNumPoints] of Double;
  CoeffV      = Array[0..MaxDeg]       of Double;

Var
  x,y         : Vector;     { Data Vectors }
  NumPoints,
  Degree,i    : Integer;
  a           : CoeffV;
  f           : Double;


procedure Lsq(x,y : Vector; N,m: Integer; Var a: CoeffV);
 Var
  i,j,k : Integer;
  B   : Array [0..MaxDeg,0..MaxDeg] of Double;
  ip  : Array [0..MaxDeg] of LongInt;
  r   : CoeffV;
  Det : Double;    { Matrix Determinant }
  sx,sxy : Double; { Temporaries        }
  tx,txy : Array [1..MaxNumPoints] of Double;
begin
  { Fill Matrix A for Least Squares Calculation }
  KFillV(tx[1],1.0,N);
  CopyV(txy[1],y[1],N);
  for i := 0 to 2*m do begin
    sx  := SumV(tx[1],N);  VxV(tx[1],tx[1],x[1],N);
    sxy := SumV(txy[1],N); VxV(txy[1],txy[1],x[1],N);
    for j := 0 to i do begin
      k := i - j;
      if (k <= m) and (j <= m) then
        B[j,k] := sx;
    end;
    if i <= m then begin
      r[i] := sxy;
    end;
  end;

  InvMatr(B[0][0],m+1,MaxDeg+1,ip[0],Det);    { Calculate Inverse of B in place }
  if det = 0.0 then WriteLn('Singular matrix');
  MxV(a[0],B[0][0],r[0],m+1,m+1,MaxDeg+1);    { a := Inverse(B)*r               }
end;

begin
  WriteLn('LEASTSQ - Example Program For Least Squares Data Fitting With MATHPAK 87');
  WriteLn;
  Write('How many points              ? '); ReadLn(NumPoints);
  Repeat
    Write('Degree of Fitting Polynomial ? '); ReadLn(Degree);
    if NumPoints < Degree + 2 then begin
      WriteLn;
      WriteLn('Error  --  Degree must be less than ',NumPoints-1);
      WriteLn;
    end;
  Until NumPoints >= Degree + 2;
  WriteLn;
  WriteLn('Enter (x,y) data pairs separated by space');
  for i := 1 to NumPoints do begin
    Write('Point ',i:2,' ==> ');
    ReadLn(x[i],y[i]);
    a[i-1] := i;
  end;
  Lsq(x,y,NumPoints,Degree,a);
  WriteLn;
  WriteLn('Fitted Coefficients');
  for i := 0 to Degree do WriteLn('a[',i,'] = ',a[i]:14);
  WriteLn;
  WriteLn(' i  x[i]        y[i]        f[i]        Residual');
  WriteLn('--  ----        ----        ----        --------');
  for i := 1 to NumPoints do begin
    f := poly87(a[0],x[i],Degree);
    WriteLn(i:2,'  ',x[i]:10,'  ',y[i]:10,'  ',f:10,'  ',(f-y[i]):10);
  end;
end.