'********************************************************************
'*       BEST APPROXIMATION OF A DISCRETE REAL FUNCTION F(X)        *
'*             BY STIEFEL-REMES'S POLYNOMIAL METHOD                 *
'* ---------------------------------------------------------------- *
'* SAMPLE RUN:                                                      *
'* (Approximate Exp(x) from 0.0 to 1.0).                            *
'*                                                                  *
'* AFTER 7  ITERATIONS                                              *
'* POLYNOMIAL COEFFICIENTS ARE:                                     *
'*   0   0.99999891                                                 *
'*   1   1.00008075                                                 *
'*   2   0.49908690                                                 *
'*   3   0.17042178                                                 *
'*   4   0.03478373                                                 *
'*   5   0.01390866                                                 *
'*           X          Y         YAPP       DIFF                   *
'*   1   0.000000   1.000000   0.999999  -0.00000109                *  
'*   2   0.100000   1.105171   1.105172   0.00000097                *
'*   3   0.200000   1.221403   1.221402  -0.00000074                *
'*   4   0.300000   1.349859   1.349858  -0.00000092                *
'*   5   0.400000   1.491825   1.491825   0.00000030                *
'*   6   0.500000   1.648721   1.648722   0.00000109                *
'*   7   0.600000   1.822119   1.822119   0.00000046                *
'*   8   0.700000   2.013753   2.013752  -0.00000082                *
'*   9   0.800000   2.225541   2.225540  -0.00000084                *
'*  10   0.900000   2.459603   2.459604   0.00000088                *
'*  11   1.000000   2.718282   2.718281  -0.00000109                *
'*                                                                  *
'* ---------------------------------------------------------------- *
'* Reference: From Numath Library By Tuan Dang Trong in Fortran 77  *
'*            [BIBLI 18].                                           *
'*                                                                  *
'*                             Basic Release By J-P Moreau, Paris   *
'*                                      (www.jpmoreau.fr)           *
'********************************************************************
'PROGRAM TEST_CHEAPP

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

  NMAX = 30

  DIM X(NMAX), Y(NMAX), A(NMAX)
  DIM T(NMAX), AX(NMAX), AY(NMAX), AH(NMAX)
  DIM BY(NMAX), BH(NMAX), IN0(NMAX)

  M = 21        'NUMBER OF POINTS OF THE FUNCTION
  N = 5         'DEGREE OF APPROXIMATION POLYNOMIAL
  NZ = 11       'NUMBER OF POINTS FOR VERIFICATION
  ITMAX = 20    'MAXIMUM NUMBER OF ITERATIONS

  ' X: TABLE OF SIZE (M) STORING THE ABSCISSAS
  ' Y: TABLE OF SIZE (M) STORING THE ORDINATES
  ' A: TABLE OF SIZE (0:N) STORING THE COEFFICIENTS
  '    OF APPROXIMATION POLYNOMIAL
  
  ' DEFINE CURVE (M PTS)
  FOR I = 1 TO M
    X(I) = 1! * (I - 1) / (1! * (M - 1))
    Y(I) = EXP(X(I))
  NEXT I

  GOSUB 1000 'call CHEFIT(M,N,X,Y,A,T,AX,AY,AH,BY,BH,IN0,IT,ITMAX)

  F$ = "####"
  F1$ = "####.######"
  F2$ = "####.########"

' PRINT RESULTS
  CLS
  PRINT
  PRINT " AFTER "; IT; " ITERATIONS"
  PRINT " POLYNOMIAL COEFFICIENTS ARE:"
  FOR I = 0 TO N
    PRINT USING F$; I;
    PRINT USING F2$; A(I)
  NEXT I

' VERIFICATION (NZ PTS)
  PRINT "           X          Y         YAPP       DIFF"
  Z0 = 0!
  DZ = .1

  FOR I = 1 TO NZ
    Z = Z0 + (I - 1) * DZ
    V = A(N)
    FOR J = N - 1 TO 0 STEP -1
      V = A(J) + V * Z
    NEXT J
    PRINT USING F$; I;
    PRINT USING F1$; Z;
    PRINT USING F1$; EXP(Z);
    PRINT USING F1$; V;
    PRINT USING F2$; V - EXP(Z)
  NEXT I

  PRINT

END 'of main program



1000  'Subroutine CHEFIT(M,N,X,Y,A,T,AX,AY,AH,BY,BH,IN0,IT,ITMAX)
'===================================================================
'     BEST APPROXIMATION OF A CONTINUOUS FUNCTION WITH REAL VARIABLE
'     IN THE SENSE OF MAXIMUM NORM.
'     THE STIEFEL-REMES"S METHOD IS USED.
'----- CALLING MODE:
'          CHEAPP(M,N,X,Y,ITMAX,IT,A,W,IW);
'----- INPUTS:
'     M:  NUMBER OF POINTS OF THE FUNCTION
'     N:  DEGREE OF APPROXIMATION POLYNOMIAL
'     X:  TABLE OF SIZE (M) STORING THE ABSCISSAS
'     Y:  TABLE OF SIZE (M) STORING THE ORDINATES
'     ITMAX: MAXIMUM NUMBER OF ITERATIONS
'     IT:    ACTUAL NUMBER OF ITERATIONS MADE
'----- OUTPUTS:
'     A:  TABLE OF SIZE (0:N) STORING THE COEFFICIENTS
'         OF APPROXIMATION POLYNOMIAL
'----- WORKING ZONES:
'     W : TABLE OF SIZE (M+5N+10)
'     (here divided into T,AX,AY,AH,BY,BH).
'     IW: TABLE OF INTEGERS OF SIZE (N+2)
'----- NOTE:
'     ABSCISSAS X(I) MUST BE IN INCREASING ORDER
'----- REFERENCE:
'     PROCEDURES ALGOL EN ANALYSE NUMERIQUE
'     EDITIONS DU CNRS, 1967
'===================================================================
      K = (M - 1) / (N + 1)
      FOR I = 1 TO N + 1
        IN0(I) = (I - 1) * K + 1
      NEXT I
      IN0(N + 2) = M
      IT = 0
1     IT = IT + 1
      FOR I = 1 TO N + 2
        AX(I) = X(IN0(I))
        AY(I) = Y(IN0(I))
        AH(I) = (-1) ^ (I - 1)
      NEXT I
      FOR I = 2 TO N + 2
        FOR J = I - 1 TO N + 2
          BY(J) = AY(J)
          BH(J) = AH(J)
        NEXT J
        FOR J = I TO N + 2
          DX = AX(J) - AX(J - I + 1)
          AY(J) = (BY(J) - BY(J - 1)) / DX
          AH(J) = (BH(J) - BH(J - 1)) / DX
        NEXT J
      NEXT I
      H = -AY(N + 2) / AH(N + 2)
      FOR I = 0 TO N
        A(I) = AY(I + 1) + AH(I + 1) * H
        BY(I + 1) = 0!
      NEXT I
      BY(1) = 1!
      TMAX = ABS(H)
      IMAX = IN0(1)
      FOR I = 1 TO N
        FOR J = 0 TO I - 1
          W = BY(I + 1 - J) - BY(I - J) * X(IN0(I))
          A(J) = A(J) + A(I) * W
          BY(I + 1 - J) = W
        NEXT J
      NEXT I
      FOR I = 1 TO M
        T(I) = A(N)
        FOR J = 1 TO N
          T(I) = T(I) * X(I) + A(N - J)
        NEXT J
        T(I) = T(I) - Y(I)
        IF ABS(T(I)) > TMAX THEN
          TMAX = ABS(T(I))
          IMAX = I
        END IF
      NEXT I
      FOR I = 1 TO N + 2
        L = I
        IF IMAX < IN0(L) THEN GOTO 2
        'Normal return
        IF (IMAX = IN0(L)) OR (IT = ITMAX) THEN RETURN
      NEXT I
      L = N + 2
2     P = T(IMAX) * T(IN0(L))
      IF P < 0! THEN GOTO 3
      IN0(L) = IMAX
      GOTO 1

3     IF IN0(1) < IMAX THEN GOTO 4
      FOR I = 1 TO N + 1
        IN0(N + 3 - I) = IN0(N + 2 - I)
      NEXT I
      IN0(1) = IMAX
      GOTO 1

4     IF IN0(N + 2) <= IMAX THEN GOTO 5
      IN0(L - 1) = IMAX
      GOTO 1

5     FOR I = 1 TO N + 1
        IN0(I) = IN0(I + 1)
      NEXT I
      IN0(N + 2) = IMAX
      GOTO 1
RETURN

'end of file tcheapp.bas