C     TITLE: FIND A NEAR-MINIMAX APPROXIMATION
C
C     THIS PROGRAM CALCULATES A NEAR-MINIMAX APPROXIMATION
C     TO A GIVEN FUNCTION 'FCN' ON THE INTERVAL [-1,1].
C     THE APPROXIMATION IS BASED ON INTERPOLATION AT THE
C     CHEBYSHEV ZEROS ON [-1,1].  THE PROGRAM THEN
C     ESTIMATES THE MAXIMUM ERROR BY CALCULATING THE ERROR
C     AT A LARGE SAMPLE OF EVENLY SPACED POINTS ON [-1,1].
C     THE PROGRAM CAN BE EASILY MODIFIED TO PRINT THE
C     DIVIDED DIFFERENCES AND THE ERRORS.
C
      PARAMETER (MAXN=20)
      PARAMETER (ZERO=0.0, ONE=1.0, P02=.02)
      DIMENSION FX(0:MAXN), X(0:MAXN)
C
C     THE FOLLOWING STATEMENT DEFINES THE FUNCTION TO BE
C     APPROXIMATED. IT CAN ALSO BE SPECIFIED USING A
C     FUNCTION SUBPROGRAM.
      FCN(X) = EXP(X)
C
C     INITIALIZE VARIABLES.
      PI = 4*ATAN(ONE)
      PRINT *, ' WHAT IS THE DEGREE N?' 
      PRINT *, ' THE DEGREE SHOULD BE .LE. ', MAXN
      READ *, N
      H = PI/(2*(N+1))
C
C     SET UP INTERPOLATION NODES AND CORRESPONDING
C     FUNCTION VALUES.
      DO 10 I=0,N
        X(I) = COS((2*I+1)*H)
        FX(I) = FCN(X(I))
10      CONTINUE
C
C     CALCULATE THE DIVIDED DIFFERENCES FOR THE
C     INTERPOLATION POLYNOMIAL.
      CALL DIVDIF(N,X,FX)
C
C     CALCULATE THE MAXIMUM INTERPOLATION ERROR AT 101
C     EVENLY SPACED POINTS ON [-1,1].
      ERRMAX = ZERO
      DO 20 I=0,100
        Z = -ONE + P02*I
C       USE NESTED MULTIPLICATION TO EVALUATE THE
C       INTERPOLATION POLYNOMIAL AT Z.
        P = FX(N)
        DO 30 J=N-1,0,-1
30        P = FX(J) + (Z - X(J))*P
        ERROR = FCN(Z) - P
        ERRMAX = MAX(ABS(ERROR),ERRMAX)
20      CONTINUE
C
C     PRINT MAXIMUM ERROR.
      PRINT *, ' FOR N =', N
      PRINT *, ' MAXIMUM ERROR =', ERRMAX
      STOP
      END

      SUBROUTINE DIVDIF(N,X,F)
C     -----------------
C
C     INPUT: (1) 'N' DENOTES A POSITIVE INTEGER.
C            (2) 'X' DENOTES A VECTOR OF POINTS X(0),...,X(N).
C            (3) 'F' DENOTES A VECTOR OF VALUES OF SOME FUNCTION
C                EVALUATED AT THE NODES IN X.
C     OUTPUT: THE VECTOR F IS CONVERTED TO A VECTOR OF DIVIDED
C                DIFFERENCES: F(I)=F[X(0),...,X(I)], I=0,...,N
C
C     WARNING: NOTE THE FORM OF THE DIMENSION STATEMENT.
C
      DIMENSION X(0:*), F(0:*)
C
      DO 1 I=1,N
      DO 1 J=N,I,-1
1       F(J) = (F(J) - F(J-1))/(X(J) - X(J-I))
      RETURN
      END