!********************************************************************
!*     BEST APPROXIMATION OF A CONTINUOUS REAL FUNCTION F(X)        *
!*     BY STIEFEL-REMES'S METHOD                                    *
!* ---------------------------------------------------------------- *
!* SAMPLE RUN:                                                      *
!* (Approximate Exp(x) from 0.0 to 1.0).                            *
!*                                                                  *
!* AFTER   7 ITERATIONS                                             *
!* POLYNOMIAL COEFFICIENTS ARE:                                     *
!*   0   0.999999E+00                                               *
!*   1   0.100008E+01                                               *
!*   2   0.499087E+00                                               *
!*   3   0.170422E+00                                               *
!*   4   0.347837E-01                                               *
!*   5   0.139087E-01                                               *
!*            X              Y             YAPP           DIFF      *
!*   1   0.000000E+00   0.100000E+01   0.999999E+00  -0.109150E-05  *   
!*   2   0.100000E+00   0.110517E+01   0.110517E+01   0.974011E-06  *
!*   3   0.200000E+00   0.122140E+01   0.122140E+01  -0.743961E-06  * 
!*   4   0.300000E+00   0.134986E+01   0.134986E+01  -0.917507E-06  *
!*   5   0.400000E+00   0.149182E+01   0.149182E+01   0.298756E-06  *
!*   6   0.500000E+00   0.164872E+01   0.164872E+01   0.109150E-05  *
!*   7   0.600000E+00   0.182212E+01   0.182212E+01   0.458606E-06  *
!*   8   0.700000E+00   0.201375E+01   0.201375E+01  -0.816261E-06  *
!*   9   0.800000E+00   0.222554E+01   0.222554E+01  -0.841892E-06  * 
!*  10   0.900000E+00   0.245960E+01   0.245960E+01   0.875494E-06  *
!*  11   0.100000E+01   0.271828E+01   0.271828E+01  -0.109150E-05  *
!*                                                                  *
!* ---------------------------------------------------------------- *
!* Reference: From Numath Library By Tuan Dang Trong in Fortran 77  *
!*            [BIBLI 18].                                           *
!*                                                                  *
!*                               F90 Release By J-P Moreau, Paris   *
!*                                       (www.jpmoreau.fr)          *
!******************************************************************** 
PROGRAM TEST_CHEAPP
  PARAMETER (M=21,N=5,N2=N+2,NZ=11,ITMAX=20,NW=M+5*N2)
  REAL*8 X(M),Y(M),A(0:N),W(NW),V,Z,Z0,DZ
  INTEGER IW(N2)

! DEFINE CURVE (21 PTS)
  DO I=1,M
    X(I)=FLOAT(I-1)/FLOAT(M-1)
    Y(I)=EXP(X(I))
  ENDDO

  CALL CHEAPP(M,N,X,Y,ITMAX,IT,A,W,IW)

! RESULTS
  WRITE(*,'(1X,A5,I4,A11)') 'AFTER',IT,' ITERATIONS'
  WRITE(*,'(1X,A)') 'POLYNOMIAL COEFFICIENTS ARE:'
  WRITE(*,'(I4,E15.6)') (I,A(I),I=0,N)

! VERIFICATION (11 PTS)
  WRITE(*,'(1X,A)') '           X              Y             YAPP           DIFF'
  Z0=0.D0
  DZ=0.1D0
  DO I=1,NZ
    Z=Z0+(I-1)*DZ
    V=A(N)
    DO J=N-1,0,-1
      V=A(J)+V*Z
    ENDDO
    WRITE(*,'(I4,4E15.6)') I,Z,EXP(Z),V,V-EXP(Z)
  ENDDO
END

      SUBROUTINE CHEAPP(M,N,X,Y,ITMAX,IT,A,W,IW)
!===================================================================
!     BEST APPROXIMATION OF A CONTINUOUS FUNCTION WITH REAL VARIABLE
!     IN THE SENSE OF MAXIMUM NORM.
!     THE STIEFEL-REMES'S METHOD IS USED.
!----- CALLING MODE:
!     CALL CHEAPP(M,N,X,Y,ITMAX,IT,A,W,IW)
!----- INPUTS:
!     M:  NUMBER OF POINTS OF THE FUNCTION
!     N:  DEGREE OF CHEBYCHEV 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)
!     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
!===================================================================
      REAL*8 X(M),Y(M),A(0:N),W(M+5*N+10)
      INTEGER IW(N+2)
      N1=M+1
      N2=N1+N+2
      N3=N2+N+2
      N4=N3+N+2
      N5=N4+N+2
      CALL CHEFIT(M,N,X,Y,A,W(1),W(N1),W(N2),W(N3),W(N4),W(N5),IW,IT,  &
      ITMAX)
      RETURN
      END
!
      SUBROUTINE CHEFIT(M,N,X,Y,A,T,AX,AY,AH,BY,BH,IN,IT,ITMAX)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION X(M),Y(M),A(0:N),T(M),AX(1:N+2),AY(1:N+2),AH(1:N+2),   &
      BY(1:N+2),BH(1:N+2),IN(1:N+2)
      K=(M-1)/(N+1)
      DO I=1,N+1
      IN(I)=(I-1)*K+1
      ENDDO
      IN(N+2)=M
      IT=0
    1 IT=IT+1
      DO I=1,N+2
      AX(I)=X(IN(I))
      AY(I)=Y(IN(I))
      AH(I)=(-1)**(I-1)
      ENDDO
      DO I=2,N+2
      DO J=I-1,N+2
      BY(J)=AY(J)
      BH(J)=AH(J)
      ENDDO
      DO J=I,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
      ENDDO
      ENDDO
      H=-AY(N+2)/AH(N+2)
      DO I=0,N
      A(I)=AY(I+1)+AH(I+1)*H
      BY(I+1)=0.D0
      ENDDO
      BY(1)=1.D0
      TMAX=ABS(H)
      IMAX=IN(1)
      DO I=1,N
      DO J=0,I-1
      W=BY(I+1-J)-BY(I-J)*X(IN(I))
      A(J)=A(J)+A(I)*W
      BY(I+1-J)=W
      ENDDO
      ENDDO
      DO I=1,M
      T(I)=A(N)
      DO J=1,N
      T(I)=T(I)*X(I)+A(N-J)
      ENDDO
      T(I)=T(I)-Y(I)
      IF(ABS(T(I)).GT.TMAX) THEN
      TMAX=ABS(T(I))
      IMAX=I
      ENDIF
      ENDDO
      DO I=1,N+2
      L=I
      IF(IMAX.LT.IN(L)) GO TO 2
      IF(IMAX.EQ.IN(L).OR.IT.EQ.ITMAX) RETURN
      ENDDO
      L=N+2
    2 P=T(IMAX)*T(IN(L))
      IF(P.LT.0.D0) GO TO 3
      IN(L)=IMAX
      GO TO 1

    3 IF(IN(1).LT.IMAX) GO TO 4
      DO I=1,N+1
      IN(N+3-I)=IN(N+2-I)
      ENDDO
      IN(1)=IMAX
      GO TO 1

    4 IF(IN(N+2).LE.IMAX) GO TO 5
      IN(L-1)=IMAX
      GO TO 1

    5 DO I=1,N+1
      IN(I)=IN(I+1)
      ENDDO
      IN(N+2)=IMAX
      GO TO 1
      END

!end of file tcheapp.f90