{***************************************************
*  Program to demonstrate Lagrange interpolation   *
*     of Function SIN(X) in double precision       *
* ------------------------------------------------ *
* 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:                                      *
*                                                  *
*  Lagrange interpolation of SIN(X):               *
*                                                  *
*    Input X: 0.5                                  *
*    Order of interpolation: 4                     *
*                                                  *
*    SIN(X) = 0.47942552 (exact value: 0.47942554) *
*                                                  *
***************************************************}
PROGRAM Lagrange;
Uses WinCrt;

Label   100;

Var
        
        X,Y : Array[0..14] of DOUBLE;
        xx,yy : DOUBLE;
        iv,n  : integer;


{*******************************************************
*          Lagrange interpolation subroutine           *
* n is the level of the interpolation ( Ex. n=2 is     *
* quadratic ). v is the total number of table values.  *
* X(i), Y(i) are the coordinate table values, Y(i)     *
* being the dependant variable. The X(i) may be arbi-  *
* trarily spaced.  x is the interpolation point which  *
* is assumed to be in the interval  with at least one  *
* table value to the left, and n to the right. If this *
* is violated, n will be set to zero. It is assumed    *
* that the table values are in ascending X(i) order.   *
*******************************************************}
PROCEDURE Interpol_Lagrange;
Label 100,200,300,400,fin;
Var i,j,k : INTEGER;
    XL    : Array[0..9] of DOUBLE;
Begin
  {Check to see if interpolation point is correct}
  if xx < X[1] then goto 100; 
  if xx <= X[iv-n] then goto 200;
  {An error has been encountered}
100: n:=0 ; goto fin;
  {Find the relevant table interval}
200: i:=0;
300: i:=i+1;
  if xx > X[i] then goto 300;
  i:=i-1;
  {Begin interpolation}
  for j:=0 to n do  XL[j]:=1.0;
  yy:=0.0;
  for k:=0 to n do
  begin
    for j:=0 to n do
    begin
      if j=k then goto 400;
      XL[k]:=XL[k]*(xx-X[j+i])/(X[i+k]-X[j+i]);
400: end;
    yy:=yy+XL[k]*Y[i+k]
  end;
fin: End;


{main program}
BEGIN
  clrscr;
  writeln;
  writeln(' Lagrange interpolation of SIN(X):');
  iv:=14;
  {Input sine table
  -----------------------------------------------------------------
   Sine table values from  Handbook of mathematical functions
   by M. Abramowitz and I.A. Stegun, NBS, june 1964
  ----------------------------------------------------------------}
  X[1]:=0.000; Y[1]:=0.00000000; X[2]:=0.125; Y[2]:=0.12467473;
  X[3]:=0.217; Y[3]:=0.21530095; X[4]:=0.299; Y[4]:=0.29456472;
  X[5]:=0.376; Y[5]:=0.36720285; X[6]:=0.450; Y[6]:=0.43496553;
  X[7]:=0.520; Y[7]:=0.49688014; X[8]:=0.589; Y[8]:=0.55552980;
  X[9]:=0.656; Y[9]:=0.60995199; X[10]:=0.721; Y[10]:=0.66013615;
  X[11]:=0.7853981634; Y[11]:=0.7071067812;
  X[12]:=0.849; Y[12]:=0.75062005; X[13]:=0.911; Y[13]:=0.79011709;
  X[14]:=0.972; Y[14]:=0.82601466;
  {----------------------------------------------------------------

   Input interpolation point }
100: writeln;

  write('   Input X: '); read(xx);
  write('   Order of interpolation: '); read(n);

  Interpol_Lagrange;

  if n<>0 then
  begin
    writeln;
    writeln('   SIN(X) = ',yy:10:8,' (exact value: ',sin(xx):10:8,')');
    goto 100
  end;
  writeln;
  DoneWinCrt
END.

{End of file lagrange.pas}