'*********************************************************
'*       Polynomial Interpolation or Extrapolation       *
'*              of a Discrete Function F(x)              *
'* ----------------------------------------------------- *
'* SAMPLE RUN:                                           *
'* (Example: Function sin(x) - 2*cos(x) is given by 12   *
'*          points from x=0 to x=1.1.                    *
'*          Extrapolate for x=1.255).                    *
'*                                                       *
'*  For X             =  1.255                           *
'*  Estimated Y value =  .3294023272245815               *
'*  Estimated Error   = -8.273064603451457D-11           *
'*  Exact Y value     =  .3294023272200048               *
'*                                                       *
'* ----------------------------------------------------- *
'* REFERENCE: "Numerical Recipes, The Art of Scientific  *
'*             Computing By W.H. Press, B.P. Flannery,   *
'*             S.A. Teukolsky and W.T. Vetterling,       *
'*             Cambridge University Press, 1986"         *
'*             [BIBLI 08].                               *
'*                                                       *
'*                  Basic Release By J-P Moreau, Paris.  *
'*                           (www.jpmoreau.fr)           *
'*********************************************************  
'PROGRAM TEST_POLINT

DEFDBL A-H, O-Z
DEFINT I-N

N = 12  'Number of points

DIM X(N), Y(N), C(N), D(N)

'define tables X and Y
  X(1) = 0#
  X(2) = .1
  X(3) = .2
  X(4) = .3
  X(5) = .4
  X(6) = .5
  X(7) = .6
  X(8) = .7
  X(9) = .8
  X(10) = .9
  X(11) = 1#
  X(12) = 1.1
  FOR I = 1 TO N
    X1 = X(I): GOSUB 500: Y(I) = FCT
  NEXT I

' define interpolation abscissa
  XX = 1.255#

  CLS
' call interpolation procedure
  GOSUB 1000  'call POLINT(X,Y,N,XX,YY,DY)
  
' print results
  PRINT
  PRINT "  For X             = "; XX
  PRINT "  Estimated Y value = "; YY
  PRINT "  Estimated Error   = "; DY
  X1 = XX: GOSUB 500
  PRINT "  Exact Y value     = "; FCT
  PRINT

END  'of main program


500 'FUNCTION FCT(X1)
  FCT = SIN(X1) - 2# * COS(X1)
RETURN

1000 'Subroutine POLINT(X,Y,N,XX,YY,DY)
'*****************************************************
'*     Polynomial Interpolation or Extrapolation     *
'*            of a Discreet Function                 *
'* ------------------------------------------------- *
'* INPUTS:                                           *
'*    X:    Table of abscissas (N)                   *
'*    Y:    Table of ordinates (N)                   *
'*    N:    Number of points                         *
'*   XX:    Interpolation abscissa value             *
'* OUTPUTS:                                          *
'*   YY:    Returned estimation of function for X    *
'*   DY:    Estimated error for YY                   *
'*****************************************************
  NS = 1
  DIF = ABS(XX - X(1))
  FOR I = 1 TO N
    DIFT = ABS(XX - X(I))
    IF DIFT < DIF THEN
      NS = I               'index of closest table entry
      DIF = DIFT
    END IF
    C(I) = Y(I)            'Initialize the C"s and D"s
    D(I) = Y(I)
  NEXT I
  YY = Y(NS)               'Initial approximation of Y
  NS = NS - 1
  FOR M = 1 TO N - 1
    FOR I = 1 TO N - M
      HO = X(I) - XX
      HP = X(I + M) - XX
      W = C(I + 1) - D(I)
      DEN = HO - HP
      IF DEN = 0# THEN
        PRINT
        PRINT " Error: two identical abscissas."
        STOP               'stop program
      END IF
      DEN = W / DEN
      D(I) = HP * DEN      'Update the C's and D's
      C(I) = HO * DEN
    NEXT I
    IF 2 * NS < N - M THEN 'After each column in the tableau XA is completed,
      DY = C(NS + 1)       'we decide which correction, C or D, we want to
    ELSE                   'add to our accumulating value of Y, i.e. which
      DY = D(NS)           'path to take through the tableau, forking up or
      NS = NS - 1          'down. We do this in such a way as to take the
    END IF                 'most "straight line" route through the tableau to
                           'its apex, updating NS accordingly to keep track
    YY = YY + DY           'of where we are. This route keeps the partial
  NEXT M                   'approximations centered (insofar as possible) on
                           'the target X.The last DY added is thus the error
                           'indication.
RETURN

'end of file tpolint.bas