C Using the D.P. discrete Fourier transform, find the approximate 
C Fourier coefficients to Runge's function on [-1,1] with N=16 and 
C N=17. 
C 
       PARAMETER (MCOEF=17)
       DOUBLE PRECISION  A(MCOEF/2),B(MCOEF/2),R(MCOEF),
     *WSAVE(3*MCOEF+15),TN,ER,PI,DEL,XJ,DFTA(MCOEF/2),
     *DFTB(MCOEF/2),F,C(MCOEF/2),S(MCOEF/2),X,RUNGE,X0,AZERO 
C
C Arithmetic statement function for Runge's function. 
       RUNGE(X) = 1.0D0/(1.0D0+25.0D0*X*X) 
C 
       X0 = -1.0D0 
       PI = ASIN(1.0D0)*2.0D0
C
       DO 10 N = MCOEF-1,MCOEF 
          CALL DEZFTI (N,WSAVE) 
C Function assumed to be periodic on [-1,1], of length 2.
          DEL = 2.0D0/N 
          F = 2.0D0*PI/(N*DEL) 
          DO 1 J = 1,N 
C First sample point at -1, last at 1-DEL
              XJ = (-1.0D0) + (J-1)*DEL 
              R(J) = RUNGE(XJ) 
C Compute sines and cosines to adjust output of DEZFTF to give
C approximate Fourier coefficients.
              IF (J .LE. N/2) THEN 
                  C(J) = COS(J*F*X0) 
                  S(J) = SIN(J*F*X0) 
              END IF 
    1     CONTINUE 
          CALL DEZFTF (N,R,AZERO,A,B,WSAVE) 
C 
C As a convenience this loop can go to N/2. If N is even last B is 
C zero.
          DO 11 J = 1,N/2 
              DFTA(J) = A(J)*C(J) - B(J)*S(J) 
              DFTB(J) = A(J)*S(J) + B(J)*C(J)  
   11     CONTINUE 
          WRITE (*,'(A,I3,A,D18.10)') ' DEZFTF RESULTS FOR N= ' ,
     *N, ' AZERO = ',AZERO
          WRITE (*,*) '    J          DFTA(J)              DFTB(J) '
          DO 12 J = 1,N/2
              WRITE(*,'(1X,I5,3X,D18.10,3X,D18.10)') J,
     *DFTA(J),DFTB(J)
   12     CONTINUE

          M = 101 
C
          WRITE (*,*) 'FOR BREVITY 101 EVALUATION POINTS OMITTED'
          IF (M.GT.0)GOTO 10
C Evaluate interpolant at 101 points on [-1,1]     
          WRITE (*,*) ' RESULTS FOR N= ',N 
          DO 20 K = 1,M 
              X = -1.0D0 + 2.0D0*(K-1.0D0)/(M-1.0D0) 
              TN = AZERO 
              DO 19 J = 1,N/2 
                  TN = TN + DFTA(J)*COS(J*F*X) + DFTB(J)*SIN(J*F*X) 
   19         CONTINUE 
              ER = TN - RUNGE(X) 
              WRITE (*,'(2X,D15.8,2X,D15.8,2X,D15.8)') X,TN,ER 
   20     CONTINUE 
C
          WRITE (*,*)
C
   10  CONTINUE 
       WRITE (*,*)
       WRITE (*,*) 'REFERENCE RESULTS FROM IBM PC/AT'
       WRITE (*,*)
     *  ' DEZFTF RESULTS FOR N=  17    AZERO =    0.274581'    
       WRITE (*,*)'    J         DFTA(J)              DFTB(J) '
       WRITE (*,*)'    1     0.3442428079E+00     0.5030951647E-16'
       WRITE (*,*)'    2     0.1755197395E+00     0.4668315002E-16'
       WRITE (*,*)'    3     0.9699027603E-01     0.3740740688E-16'
       WRITE (*,*)'    4     0.4964416088E-01     0.2403131256E-16'
       WRITE (*,*)'    5     0.2759681609E-01     0.8202373493E-17'
       WRITE (*,*)'    6     0.1323313221E-01     0.4148226049E-17'
       WRITE (*,*)'    7     0.7099464912E-02     0.1053997863E-16'
       WRITE (*,*)'    8     0.1413251200E-02     0.5512226224E-17'
C
       STOP 
       END