'****************************************************
'*  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].*
'*                                                  *
'* ------------------------------------------------ *
'* SAMPLE RUN:                                      *
'*                                                  *
'*  Lagrange interpolation of SIN(X):               *
'*                                                  *
'*    Input X: 0.5                                  *
'*    Order of interpolation: 4                     *
'*                                                  *
'*    SIN(X) = 0.47942552 (exact value: 0.47942554) *
'*    Error check:  4                               *
'*                                                  *
'****************************************************
defint i-n
defdbl a-h,o-z
dim XL(9),Y(15),X(15)
  cls
  print
  print " Lagrange interpolation of SIN(X):"
  iv=14
  'Input table
  for i=1 to iv
    read X(i), Y(i)
  next i
  '------------------------------------------------------------
  'Sine table values from  Handbook of mathematical
  'functions by M. Abramowitz and I.A. Stegun,
  'NBS, june 1964
  data 0.000#,0.00000000#,0.125#,0.12467473#
  data 0.217#,0.21530095#,0.299#,0.29456472#,0.376#,0.36720285#
  data 0.450#,0.43496553#,0.520#,0.49688014#,0.589#,0.55552980#
  data 0.656#,0.60995199#,0.721#,0.66013615#
  data 0.7853981634#,0.7071067812#
  data 0.849#,0.75062005#,0.911#,0.79011709#
  data 0.972#,0.82601466#
  '------------------------------------------------------------
  'Input interpolation point
100 print

  input "   Input X: ",xx
  input "   Order of interpolation: ",n

  gosub 1000

  if n<>0 then
    print
    print using "   SIN(X) = #.########"; yy;
    print using " (exact value: #.########)"; SIN(xx)
    print "   Error check: ";n
    goto 100
  end if
  print
end


'********************************************************
'*          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.   *
'********************************************************
1000 'Check to see if interpolation point is correct
  if xx<X(1) then goto 1100
  if xx<=X(iv-n) then goto 1200
  'An error has been encountered
1100 n=0 : return
  'Find the relevant table interval
1200 i=0
1300 i=i+1
  if xx>X(i) then goto 1300
  i=i-1
  'Begin interpolation
  for j=0 to n
    XL(j)=1#
  next j
  yy=0
  for k=0 to n
    for j=0 to n
      if j=k then goto 1400
      XL(k)=XL(k)*(xx-X(j+i))/(X(i+k)-X(j+i))
1400 next j
    yy=yy+XL(k)*Y(i+k)
  next k
return

'End of file Lagrange.bas