'****************************************************
'*   Program to demonstrate Newton 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:                                      *
'*                                                  *
'*  Newton interpolation of SIN(X):                 *
'*                                                  *
'*    Input X: 0.5                                  *
'*                                                  *
'*    SIN(X) = 0.47961461 (exact value: 0.47942554) *
'*    Error estimate: -0.00000028                   *
'*                                                  *
'****************************************************
defint i-n
defdbl a-h,o-z
  iv=14
  dim X(iv+1),Y(iv+1),Y1(iv),Y2(iv-1),Y3(iv-2)
  cls
  print
  print " Newton interpolation of SIN(X):"
  '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

  gosub 1000

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


'********************************************************
'* Newton divided differences interpolation subroutine  *
'* ---------------------------------------------------- *
'* Calculates cubic interpolation for a given table.    *
'* iv 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 arbitrarily spa- *
'* ced. xx is the interpolation point which is assumed  *
'* to be in the interval  with at least one table value *
'* to the left, and 3 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. e is the error   *
'* estimate. The returned value is yy.                  *
'********************************************************
1000 'Check to see if interpolation point is correct
  n=1
  if xx>=X(1) then goto 1100
  n=0 : return
1100 if xx<=X(iv-3) then goto 1200
  n=0 : return
  'Generate divided differences
1200 for i=1 to iv-1
    Y1(i)=(Y(i+1)-Y(i))/(X(i+1)-X(i))
  next i
  for i=1 to iv-2
    Y2(i)=(Y1(i+1)-Y1(i))/(X(i+1)-X(i))
  next i
  for i=1 to iv-3
    Y3(i)=(Y2(i+1)-Y2(i))/(X(i+1)-X(i))
  next i
  'Find relevant table interval
  i=0
1300 i=i+1
  if xx>X(i) then goto 1300
  i=i-1
  'Begin interpolation
  a=xx-X(i)
  b=a*(xx-X(i+1))
  c=b*(xx-X(i+2))
  yy=Y(i)+a*Y1(i)+b*Y2(i)+c*Y3(i)
  'Calculate next term in the expansion for an error estimate
  e=c*(xx-X(i+3))*yy/24#
return

'End of file Newton.bas