'*********************************************************
'*       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 =  .3294023272088324               *
'*  Estimated Error   = -1.327291621149679D-10           *
'*  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 RATINT(X,Y,N,XX,YY,DY)
'*****************************************************
'*    Interpolation or Extrapolation of a Discrete   *
'*       Function By a Quotient of Polynomials       *
'* ------------------------------------------------- *
'* 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                   *
'*****************************************************
  TINY = .000000000000001#
  NS = 1
  HH = ABS(XX - X(1))
  FOR I = 1 TO N
    H = ABS(XX - X(I))
    IF H = 0# THEN
      Y = Y(I)
      DY = 0#
      STOP
    ELSEIF H < HH THEN
      NS = I               'index of closest table entry
      HH = H
    END IF
    C(I) = Y(I)
    D(I) = Y(I) + TINY     '{TINY is to prevent a zero-over-zero
  NEXT I                   'condition.
  YY = Y(NS)
  NS = NS - 1
  FOR M = 1 TO N - 1
    FOR I = 1 TO N - M
      W = C(I + 1) - D(I)
      H = X(I + M) - XX      'H<>0 (tested before)
      T = (X(I) - XX) * D(I) / H
      DD = T - C(I + 1)
      IF DD = 0# THEN
        PRINT
        PRINT " Error: Pole at requested value of X."
        STOP                 'stop program
      END IF
      DD = W / DD
      D(I) = C(I + 1) * DD
      C(I) = T * DD
    NEXT I
    IF 2 * NS < N - M THEN
      DY = C(NS + 1)
    ELSE
      DY = D(NS)
      NS = NS - 1
    END IF
    YY = YY + DY
  NEXT M

RETURN

'end of file tratint.bas