ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c  A real-valued, in place, split-radix FFT program
c  Real input and output in data array X
c  Length is N = 2 ** M
c  Decimation-in-time, cos/sin in second loop
c  Output in order:
c         [ Re(0), Re(1), ..., Re(N/2), Im(N/2-1), ..., Im(1) ]
c
c  This FFT computes
c     X(k) = sum_{j=0}^{N-1} x(j)*exp(-2ijk*pi/N)
c
c
c  H.V. Sorensen, Rice University, Oct. 1985
c
c  Reference:  H.V. Sorensen, D.L. Jones, M.T. Heideman, & C.S. Burrus;
c              Real-Valued Fast Fourier Transform Algorithms; IEEE
c              Trans. Acoust., Speech, Signal Process.; Vol ASSP-35,
c              June 1987, pp. 849-863.
c
c  This code was originally named RVFFT.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      SUBROUTINE SFFTCF( X, N, M )
c ... Scalar arguments ...
      INTEGER N, M
c ... Array arguments ...
      REAL  X(*)
c ... Local scalars ...
      INTEGER  J, I, K, IS, ID, I0, I1, I2, I3, I4, I5, I6, I7, I8
      INTEGER  N1, N2, N4, N8
      REAL  SQRT2, TWOPI, XT, R1, T1, T2, T3, T4, T5, T6
      REAL  A, A3, E, CC1, SS1, CC3, SS3
c ... Intrinsic functions ...
      INTRINSIC  SIN, COS
c ... Parameters ...
      PARAMETER ( SQRT2 = 1.4142135623730950488,
     +            TWOPI = 6.2831853071795864769 )
c
c ... Exe. statements ...
c
      IF ( N .EQ. 1 ) RETURN
c
 100  J = 1
      N1 = N - 1
      DO 104, I = 1, N1
         IF ( I .GE. J ) GOTO 101
         XT = X(J)
         X(J) = X(I)
         X(I) = XT
 101     K = N / 2
 102     IF ( K .GE. J ) GOTO 103
            J = J - K
            K = K / 2
            GOTO 102
 103     J = J + K
 104  CONTINUE
c 
      IS = 1
      ID = 4
 70   DO 60, I0 = IS, N, ID
         I1 = I0 + 1
         R1 = X(I0)
         X(I0) = R1 + X(I1)
         X(I1) = R1 - X(I1)
 60   CONTINUE
      IS = 2 * ID - 1
      ID = 4 * ID
      IF ( IS .LT. N ) GOTO 70
c
      N2 = 2
      DO 10, K = 2, M
         N2 = N2 * 2
         N4 = N2 / 4
         N8 = N2 / 8
         E = TWOPI / N2
         IS = 0
         ID = N2 * 2
 40      DO 38, I = IS, N-1, ID
            I1 = I + 1
            I2 = I1 + N4
            I3 = I2 + N4
            I4 = I3 + N4
            T1 = X(I4) + X(I3)
            X(I4) = X(I4) - X(I3)
            X(I3) = X(I1) - T1
            X(I1) = X(I1) + T1
            IF ( N4 .EQ. 1 ) GOTO 38
            I1 = I1 + N8
            I2 = I2 + N8
            I3 = I3 + N8
            I4 = I4 + N8
            T1 = ( X(I3) + X(I4) ) / SQRT2
            T2 = ( X(I3) - X(I4) ) / SQRT2
            X(I4) = X(I2) - T1
            X(I3) = - X(I2) - T1
            X(I2) = X(I1) - T2
            X(I1) = X(I1) + T2
 38      CONTINUE
         IS = 2 * ID - N2
         ID = 4 * ID
         IF ( IS .LT. N ) GOTO 40
         A = E
         DO 32, J = 2, N8
            A3 = 3 * A
            CC1 = COS(A)
            SS1 = SIN(A)
            CC3 = COS(A3)
            SS3 = SIN(A3)
            A = J * E
            IS = 0
            ID = 2 * N2
 36         DO 30, I = IS, N-1, ID
               I1 = I + J
               I2 = I1 + N4
               I3 = I2 + N4
               I4 = I3 + N4
               I5 = I + N4 - J + 2
               I6 = I5 + N4
               I7 = I6 + N4
               I8 = I7 + N4
               T1 = X(I3) * CC1 + X(I7) * SS1
               T2 = X(I7) * CC1 - X(I3) * SS1
               T3 = X(I4) * CC3 + X(I8) * SS3
               T4 = X(I8) * CC3 - X(I4) * SS3
               T5 = T1 + T3
               T6 = T2 + T4
               T3 = T1 - T3
               T4 = T2 - T4
               T2 = X(I6) + T6
               X(I3) = T6 - X(I6)
               X(I8) = T2
               T2 = X(I2) - T3
               X(I7) = - X(I2) - T3
               X(I4) = T2
               T1 = X(I1) + T5
               X(I6) = X(I1) - T5
               X(I1) = T1
               T1 = X(I5) + T4
               X(I5) = X(I5) - T4
               X(I2) = T1
 30         CONTINUE
            IS = 2 * ID - N2
            ID = 4 * ID
            IF ( IS .LT. N ) GOTO 36
 32      CONTINUE
 10   CONTINUE
      RETURN
c
c ... End of subroutine SFFTCF ...
c
      END