C/ FIGURE 5.3.4
      SUBROUTINE NRFFT(F,M)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMPLEX*16 F(2**M),FODD(2**M/2),FEVEN(2**M/2) 
      COMPLEX*16 D,H,U
C
C  SUBROUTINE NRFFT PERFORMS A FAST FOURIER TRANSFORM ON THE COMPLEX
C     VECTOR F, OF LENGTH N=2**M.  THE FOURIER TRANSFORM IS DEFINED BY
C       Y(K) = SUM FROM J=1 TO N OF: EXP[I*2*PI*(K-1)*(J-1)/N]*F(J)
C     WHERE I = SQRT(-1).
C
C  ARGUMENTS
C
C             ON INPUT                          ON OUTPUT
C             --------                          ---------
C
C    F      - THE COMPLEX VECTOR OF LENGTH      THE TRANSFORMED VECTOR
C             2**M TO BE TRANSFORMED.           Y.
C
C    M      - THE LENGTH OF THE VECTOR F
C             IS ASSUMED TO BE 2**M.
C
C-----------------------------------------------------------------------
      N = 2**M
      PI = 3.141592653589793 D0
C                              REORDER ELEMENTS OF F AS THEY WILL BE 
C                              IN TRANSFORMED VECTOR
      DO 40 I=M,2,-1
         NI = 2**I
         N2 = NI/2
         DO 30 J=0,N-1,NI
            DO 10 K=1,N2
               FODD(K) = F(J+2*K-1)
               FEVEN(K) = F(J+2*K)
   10       CONTINUE
            DO 20 K=1,N2
               F(J+K) = FODD(K)
               F(J+N2+K) = FEVEN(K)
   20       CONTINUE
   30    CONTINUE
   40 CONTINUE
C                              PROCESS ELEMENTS BY GROUPS OF 2, 
C                              THEN BY GROUPS OF 4,8,...,2**M 
      DO 70 I=1,M
         NI = 2**I
         N2 = NI/2
         H = COS(2*PI/NI) + SIN(2*PI/NI)*CMPLX(0.0,1.0)
         DO 60 J=0,N-1,NI
            D = 1.0
            DO 50 K=1,N2
               FODD(K) = F(J+K)
               FEVEN(K) = F(J+N2+K)
               U = D*FEVEN(K)
               F(J+K) = FODD(K) + U
               F(J+N2+K) = FODD(K) - U
               D = D*H
   50       CONTINUE
   60    CONTINUE
   70 CONTINUE
      RETURN
      END
