'*************************************************************
'*    Interpolate a function F(x) by continuous fractions    *
'* --------------------------------------------------------- *
'* SAMPLE RUN:                                               *
'* (Interpolate function e(x) between x=0 and x=2)           *
'*                                                           *
'* Number of points: 3                                       *
'* X, Y: 0,1                                                 *
'* X, Y: 1,2.71828                                           *
'* X, Y: 2,7.38906                                           *
'*                                                           *
'* Coefficients D(K):                                        *
'* D(0) =      1.000000                                      *
'* D(1) =      0.581977                                      *
'* D(2) =     -3.718271                                      *
'*                                                           *
'* X = 1.5                                                   *
'*                                                           *
'* For X = 1.5    Y =    4.351909                            *
'*                                                           *
'* --------------------------------------------------------- *
'* Ref.: "Methodes de calcul numerique, Tome 2 By Claude     *
'*        Nowakowski, PSI Edition, 1984" [BIBLI 04].         *
'*                                                           *
'*                       Basic Release By J-P Moreau, Paris. *
'*                                (www.jpmoreau.fr)          *
'*************************************************************
DEFINT I-N
DEFDBL A-H, O-Z

CLS
PRINT
INPUT " Number of points: ", N1

N = N1 - 1
M = N

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

'Read data from screen
FOR K = 0 TO N
  INPUT " X, Y: ", X(K), Y(K)
NEXT K

'Calculate coefficients D(K)
FOR K = 0 TO N
  D(K) = Y(K)
NEXT K

FOR L = 1 TO M
  FOR K = L TO N
    DD = (X(K) - X(L - 1)) / (D(K) - D(L - 1))
    IF K <> L THEN
      D(K) = DD
    ELSE
      DL = DD
    END IF
  NEXT K
  D(L) = DL
NEXT L

F$ = "####.######"

'print coefficients
PRINT
PRINT " Coefficients D(K):"
FOR K = 0 TO N
  PRINT " D("; K; ") = ";
  PRINT USING F$; D(K)
NEXT K

'Interpolate for X=XX
PRINT
INPUT " X=", XX

'Evaluate continuous fraction
S = (XX - X(N - 1)) / D(N)
FOR K = N - 1 TO 1 STEP -1
  S = (XX - X(K - 1)) / (D(K) + S)
NEXT K

S = S + D(0)
PRINT
PRINT " For X= "; XX; "  Y=";
PRINT USING F$; S

END