C Using complex discrete Fourier transform, find the approximate Fourier
C coefficients to Runge's function on [-1,1] with N=16 and N=17. 
C 
       DOUBLE PRECISION  WSAVE(150),RUNGE,X,X0,PI,DEL,F,XJ
       COMPLEX*16 COEFF(0:16), SQTM1
C
C Note: Complex Double Precision is not Fortran Standard and many
C       compiler vendors may fail to implement it or implement it
C       in some other manner.  If the COMPLEX*16 declaration does
C       not work on your compiler we recommend that you check
C       with the vendor's instructions and modify the above accordingly.
C
C Note 0 subscript, which makes indexing easier, allowed in Fortran 77.
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 
       SQTM1 = CMPLX(0.0D0,-1.0D0)

C
       DO 10 N = 16,17 
          CALL DCFFTI (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 = 0,N-1
C First sample point at -1, last at 1-DEL
              XJ = (-1.0D0) + J*DEL 
              COEFF(J) = CMPLX(RUNGE(XJ),0.0D0) 
    1     CONTINUE 
          CALL DCFFTF (N,COEFF,WSAVE) 
C Returned coefficients must be divided by N for correct normaliziation.
C 
C Note repetition after N/2 in original coefficients. Scaling because
C X0 not at origin destroys this to some extent. 
C
          WRITE (*,*) ' DCFFTF RESULTS FOR N = ' ,N 
          WRITE (*,'(A,1X,2D17.8)') ' CZERO = ',COEFF(0D0)/N*2.0D0
          WRITE (*,*)  
     *' J          OUTPUT FROM DCFFTF,         SCALED COEFFICIENTS'   
          DO 11 J = 1,N-1
              WRITE (*,'(I5,2D17.8,5X,2D17.8)') 
     *              J, COEFF(J), EXP(-SQTM1*J*F*X0) * COEFF(J)/N *2
   11     CONTINUE 
          WRITE (*,*)

   10  CONTINUE 
       WRITE (*,*) 'REFERENCE RESULTS (PARTIAL) FROM IBM PC/AT'
       WRITE (*,*) ' DCFFTF RESULTS FOR N =           17'
       WRITE (*,*) ' CZERO =  (0.54916124E+00, 0.00000000E+00)'
       WRITE (*,*) '    .  '
       WRITE (*,*) '    .  '
       WRITE (*,*) '  10  -0.60345452E-01  -0.37837322E-16 '  
       WRITE (*,*) '  11   0.11248162E+00  -0.47424091E-16 '  
       WRITE (*,*) '  12  -0.23457294E+00   0.73973161E-16 '  
       WRITE (*,*) '  13   0.42197537E+00  -0.25268310E-17 '  
       WRITE (*,*) '  14  -0.82441735E+00  -0.14952867E-16 '
       WRITE (*,*) '  15   0.14919178E+01   0.31242621E-16 '
       WRITE (*,*) '  16  -0.29260639E+01  -0.69144627E-16 '
       STOP 
       END