module OouraFFT
! Fast Fourier/Cosine/Sine Transform
!     dimension   :one
!     data length :power of 2
!     decimation  :frequency
!     radix       :split-radix
!     data        :inplace
!     table       :use
! subroutines
!     cdft: Complex Discrete Fourier Transform
!     rdft: Real Discrete Fourier Transform
!     ddct: Discrete Cosine Transform
!     ddst: Discrete Sine Transform
!     dfct: Cosine Transform of RDFT (Real Symmetric DFT)
!     dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
!
!
! -------- Complex DFT (Discrete Fourier Transform) --------
!     [definition]
!         <case1>
!             X(k) = sum_j=0^n-1 x(j)*exp(2*pi*i*j*k/n), 0<=k<n
!         <case2>
!             X(k) = sum_j=0^n-1 x(j)*exp(-2*pi*i*j*k/n), 0<=k<n
!         (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
!     [usage]
!         <case1>
!             ip(0) = 0  ! first time only
!             call cdft(2*n, 1, a, ip, w)
!         <case2>
!             ip(0) = 0  ! first time only
!             call cdft(2*n, -1, a, ip, w)
!     [parameters]
!         2*n          :data length (integer)
!                       n >= 1, n = power of 2
!         a(0:2*n-1)   :input/output data (real*8)
!                       input data
!                           a(2*j) = Re(x(j)),
!                           a(2*j+1) = Im(x(j)), 0<=j<n
!                       output data
!                           a(2*k) = Re(X(k)),
!                           a(2*k+1) = Im(X(k)), 0<=k<n
!         ip(0:*)      :work area for bit reversal (integer)
!                       length of ip >= 2+sqrt(n)
!                       strictly,
!                       length of ip >=
!                           2+2**(int(log(n+0.5)/log(2.0))/2).
!                       ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:n/2-1)   :cos/sin table (real*8)
!                       w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of
!             call cdft(2*n, -1, a, ip, w)
!         is
!             call cdft(2*n, 1, a, ip, w)
!             do j = 0, 2 * n - 1
!                 a(j) = a(j) / n
!             end do
!         .
!
!
! -------- Real DFT / Inverse of Real DFT --------
!     [definition]
!         <case1> RDFT
!             R(k) = sum_j=0^n-1 a(j)*cos(2*pi*j*k/n), 0<=k<=n/2
!             I(k) = sum_j=0^n-1 a(j)*sin(2*pi*j*k/n), 0<k<n/2
!         <case2> IRDFT (excluding scale)
!             a(k) = (R(0) + R(n/2)*cos(pi*k))/2 +
!                    sum_j=1^n/2-1 R(j)*cos(2*pi*j*k/n) +
!                    sum_j=1^n/2-1 I(j)*sin(2*pi*j*k/n), 0<=k<n
!     [usage]
!         <case1>
!             ip(0) = 0  ! first time only
!             call rdft(n, 1, a, ip, w)
!         <case2>
!             ip(0) = 0  ! first time only
!             call rdft(n, -1, a, ip, w)
!     [parameters]
!         n            :data length (integer)
!                       n >= 2, n = power of 2
!         a(0:n-1)     :input/output data (real*8)
!                       <case1>
!                           output data
!                               a(2*k) = R(k), 0<=k<n/2
!                               a(2*k+1) = I(k), 0<k<n/2
!                               a(1) = R(n/2)
!                       <case2>
!                           input data
!                               a(2*j) = R(j), 0<=j<n/2
!                               a(2*j+1) = I(j), 0<j<n/2
!                               a(1) = R(n/2)
!         ip(0:*)      :work area for bit reversal (integer)
!                       length of ip >= 2+sqrt(n/2)
!                       strictly,
!                       length of ip >=
!                           2+2**(int(log(n/2+0.5)/log(2.0))/2).
!                       ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:n/2-1)   :cos/sin table (real*8)
!                       w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of
!             call rdft(n, 1, a, ip, w)
!         is
!             call rdft(n, -1, a, ip, w)
!             do j = 0, n - 1
!                 a(j) = a(j) * 2 / n
!             end do
!         .
!
!
! -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
!     [definition]
!         <case1> IDCT (excluding scale)
!             C(k) = sum_j=0^n-1 a(j)*cos(pi*j*(k+1/2)/n), 0<=k<n
!         <case2> DCT
!             C(k) = sum_j=0^n-1 a(j)*cos(pi*(j+1/2)*k/n), 0<=k<n
!     [usage]
!         <case1>
!             ip(0) = 0  ! first time only
!             call ddct(n, 1, a, ip, w)
!         <case2>
!             ip(0) = 0  ! first time only
!             call ddct(n, -1, a, ip, w)
!     [parameters]
!         n            :data length (integer)
!                       n >= 2, n = power of 2
!         a(0:n-1)     :input/output data (real*8)
!                       output data
!                           a(k) = C(k), 0<=k<n
!         ip(0:*)      :work area for bit reversal (integer)
!                       length of ip >= 2+sqrt(n/2)
!                       strictly,
!                       length of ip >=
!                           2+2**(int(log(n/2+0.5)/log(2.0))/2).
!                       ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:n*5/4-1) :cos/sin table (real*8)
!                       w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of
!             call ddct(n, -1, a, ip, w)
!         is
!             a(0) = a(0) / 2
!             call ddct(n, 1, a, ip, w)
!             do j = 0, n - 1
!                 a(j) = a(j) * 2 / n
!             end do
!         .
!
!
! -------- DST (Discrete Sine Transform) / Inverse of DST --------
!     [definition]
!         <case1> IDST (excluding scale)
!             S(k) = sum_j=1^n A(j)*sin(pi*j*(k+1/2)/n), 0<=k<n
!         <case2> DST
!             S(k) = sum_j=0^n-1 a(j)*sin(pi*(j+1/2)*k/n), 0<k<=n
!     [usage]
!         <case1>
!             ip(0) = 0  ! first time only
!             call ddst(n, 1, a, ip, w)
!         <case2>
!             ip(0) = 0  ! first time only
!             call ddst(n, -1, a, ip, w)
!     [parameters]
!         n            :data length (integer)
!                       n >= 2, n = power of 2
!         a(0:n-1)     :input/output data (real*8)
!                       <case1>
!                           input data
!                               a(j) = A(j), 0<j<n
!                               a(0) = A(n)
!                           output data
!                               a(k) = S(k), 0<=k<n
!                       <case2>
!                           output data
!                               a(k) = S(k), 0<k<n
!                               a(0) = S(n)
!         ip(0:*)      :work area for bit reversal (integer)
!                       length of ip >= 2+sqrt(n/2)
!                       strictly,
!                       length of ip >=
!                           2+2**(int(log(n/2+0.5)/log(2.0))/2).
!                       ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:n*5/4-1) :cos/sin table (real*8)
!                       w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of
!             call ddst(n, -1, a, ip, w)
!         is
!             a(0) = a(0) / 2
!             call ddst(n, 1, a, ip, w)
!             do j = 0, n - 1
!                 a(j) = a(j) * 2 / n
!             end do
!         .
!
!
! -------- Cosine Transform of RDFT (Real Symmetric DFT) --------
!     [definition]
!         C(k) = sum_j=0^n a(j)*cos(pi*j*k/n), 0<=k<=n
!     [usage]
!         ip(0) = 0  ! first time only
!         call dfct(n, a, t, ip, w)
!     [parameters]
!         n            :data length - 1 (integer)
!                       n >= 2, n = power of 2
!         a(0:n)       :input/output data (real*8)
!                       output data
!                           a(k) = C(k), 0<=k<=n
!         t(0:n/2)     :work area (real*8)
!         ip(0:*)      :work area for bit reversal (integer)
!                       length of ip >= 2+sqrt(n/4)
!                       strictly,
!                       length of ip >=
!                           2+2**(int(log(n/4+0.5)/log(2.0))/2).
!                       ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:n*5/8-1) :cos/sin table (real*8)
!                       w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of
!             a(0) = a(0) / 2
!             a(n) = a(n) / 2
!             call dfct(n, a, t, ip, w)
!         is
!             a(0) = a(0) / 2
!             a(n) = a(n) / 2
!             call dfct(n, a, t, ip, w)
!             do j = 0, n
!                 a(j) = a(j) * 2 / n
!             end do
!         .
!
!
! -------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
!     [definition]
!         S(k) = sum_j=1^n-1 a(j)*sin(pi*j*k/n), 0<k<n
!     [usage]
!         ip(0) = 0  ! first time only
!         call dfst(n, a, t, ip, w)
!     [parameters]
!         n            :data length + 1 (integer)
!                       n >= 2, n = power of 2
!         a(0:n-1)     :input/output data (real*8)
!                       output data
!                           a(k) = S(k), 0<k<n
!                       (a(0) is used for work area)
!         t(0:n/2-1)   :work area (real*8)
!         ip(0:*)      :work area for bit reversal (integer)
!                       length of ip >= 2+sqrt(n/4)
!                       strictly,
!                       length of ip >=
!                           2+2**(int(log(n/4+0.5)/log(2.0))/2).
!                       ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:n*5/8-1) :cos/sin table (real*8)
!                       w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of
!             call dfst(n, a, t, ip, w)
!         is
!             call dfst(n, a, t, ip, w)
!             do j = 1, n - 1
!                 a(j) = a(j) * 2 / n
!             end do
!         .
!
!
! Appendix :
!     The cos/sin table is recalculated when the larger table required.
!     w() and ip() are compatible with all routines.
!
!
contains
   subroutine cdft(n, isgn, a, ip, w)
      integer n, isgn, ip(0:*), nw
      real*8 a(0:n - 1), w(0:*)
      nw = ip(0)
      if (n .gt. 4*nw) then
         nw = n/4
         call makewt(nw, ip, w)
      end if
      if (isgn .ge. 0) then
         call cftfsub(n, a, ip, nw, w)
      else
         call cftbsub(n, a, ip, nw, w)
      end if
   end
!
   subroutine rdft(n, isgn, a, ip, w)
      integer n, isgn, ip(0:*), nw, nc
      real*8 a(0:n - 1), w(0:*), xi
      nw = ip(0)
      if (n .gt. 4*nw) then
         nw = n/4
         call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n .gt. 4*nc) then
         nc = n/4
         call makect(nc, ip, w(nw))
      end if
      if (isgn .ge. 0) then
         if (n .gt. 4) then
            call cftfsub(n, a, ip, nw, w)
            call rftfsub(n, a, nc, w(nw))
         else if (n .eq. 4) then
            call cftfsub(n, a, ip, nw, w)
         end if
         xi = a(0) - a(1)
         a(0) = a(0) + a(1)
         a(1) = xi
      else
         a(1) = 0.5d0*(a(0) - a(1))
         a(0) = a(0) - a(1)
         if (n .gt. 4) then
            call rftbsub(n, a, nc, w(nw))
            call cftbsub(n, a, ip, nw, w)
         else if (n .eq. 4) then
            call cftbsub(n, a, ip, nw, w)
         end if
      end if
   end
!
   subroutine ddct(n, isgn, a, ip, w)
      integer n, isgn, ip(0:*), j, nw, nc
      real*8 a(0:n - 1), w(0:*), xr
      nw = ip(0)
      if (n .gt. 4*nw) then
         nw = n/4
         call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n .gt. nc) then
         nc = n
         call makect(nc, ip, w(nw))
      end if
      if (isgn .lt. 0) then
         xr = a(n - 1)
         do j = n - 2, 2, -2
            a(j + 1) = a(j) - a(j - 1)
            a(j) = a(j) + a(j - 1)
         end do
         a(1) = a(0) - xr
         a(0) = a(0) + xr
         if (n .gt. 4) then
            call rftbsub(n, a, nc, w(nw))
            call cftbsub(n, a, ip, nw, w)
         else if (n .eq. 4) then
            call cftbsub(n, a, ip, nw, w)
         end if
      end if
      call dctsub(n, a, nc, w(nw))
      if (isgn .ge. 0) then
         if (n .gt. 4) then
            call cftfsub(n, a, ip, nw, w)
            call rftfsub(n, a, nc, w(nw))
         else if (n .eq. 4) then
            call cftfsub(n, a, ip, nw, w)
         end if
         xr = a(0) - a(1)
         a(0) = a(0) + a(1)
         do j = 2, n - 2, 2
            a(j - 1) = a(j) - a(j + 1)
            a(j) = a(j) + a(j + 1)
         end do
         a(n - 1) = xr
      end if
   end
!
   subroutine ddst(n, isgn, a, ip, w)
      integer n, isgn, ip(0:*), j, nw, nc
      real*8 a(0:n - 1), w(0:*), xr
      nw = ip(0)
      if (n .gt. 4*nw) then
         nw = n/4
         call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n .gt. nc) then
         nc = n
         call makect(nc, ip, w(nw))
      end if
      if (isgn .lt. 0) then
         xr = a(n - 1)
         do j = n - 2, 2, -2
            a(j + 1) = -a(j) - a(j - 1)
            a(j) = a(j) - a(j - 1)
         end do
         a(1) = a(0) + xr
         a(0) = a(0) - xr
         if (n .gt. 4) then
            call rftbsub(n, a, nc, w(nw))
            call cftbsub(n, a, ip, nw, w)
         else if (n .eq. 4) then
            call cftbsub(n, a, ip, nw, w)
         end if
      end if
      call dstsub(n, a, nc, w(nw))
      if (isgn .ge. 0) then
         if (n .gt. 4) then
            call cftfsub(n, a, ip, nw, w)
            call rftfsub(n, a, nc, w(nw))
         else if (n .eq. 4) then
            call cftfsub(n, a, ip, nw, w)
         end if
         xr = a(0) - a(1)
         a(0) = a(0) + a(1)
         do j = 2, n - 2, 2
            a(j - 1) = -a(j) - a(j + 1)
            a(j) = a(j) - a(j + 1)
         end do
         a(n - 1) = -xr
      end if
   end
!
   subroutine dfct(n, a, t, ip, w)
      integer n, ip(0:*), j, k, l, m, mh, nw, nc
      real*8 a(0:n), t(0:n/2), w(0:*), xr, xi, yr, yi
      nw = ip(0)
      if (n .gt. 8*nw) then
         nw = n/8
         call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n .gt. 2*nc) then
         nc = n/2
         call makect(nc, ip, w(nw))
      end if
      m = n/2
      yi = a(m)
      xi = a(0) + a(n)
      a(0) = a(0) - a(n)
      t(0) = xi - yi
      t(m) = xi + yi
      if (n .gt. 2) then
         mh = m/2
         do j = 1, mh - 1
            k = m - j
            xr = a(j) - a(n - j)
            xi = a(j) + a(n - j)
            yr = a(k) - a(n - k)
            yi = a(k) + a(n - k)
            a(j) = xr
            a(k) = yr
            t(j) = xi - yi
            t(k) = xi + yi
         end do
         t(mh) = a(mh) + a(n - mh)
         a(mh) = a(mh) - a(n - mh)
         call dctsub(m, a, nc, w(nw))
         if (m .gt. 4) then
            call cftfsub(m, a, ip, nw, w)
            call rftfsub(m, a, nc, w(nw))
         else if (m .eq. 4) then
            call cftfsub(m, a, ip, nw, w)
         end if
         a(n - 1) = a(0) - a(1)
         a(1) = a(0) + a(1)
         do j = m - 2, 2, -2
            a(2*j + 1) = a(j) + a(j + 1)
            a(2*j - 1) = a(j) - a(j + 1)
         end do
         l = 2
         m = mh
         do while (m .ge. 2)
            call dctsub(m, t, nc, w(nw))
            if (m .gt. 4) then
               call cftfsub(m, t, ip, nw, w)
               call rftfsub(m, t, nc, w(nw))
            else if (m .eq. 4) then
               call cftfsub(m, t, ip, nw, w)
            end if
            a(n - l) = t(0) - t(1)
            a(l) = t(0) + t(1)
            k = 0
            do j = 2, m - 2, 2
               k = k + 4*l
               a(k - l) = t(j) - t(j + 1)
               a(k + l) = t(j) + t(j + 1)
            end do
            l = 2*l
            mh = m/2
            do j = 0, mh - 1
               k = m - j
               t(j) = t(m + k) - t(m + j)
               t(k) = t(m + k) + t(m + j)
            end do
            t(mh) = t(m + mh)
            m = mh
         end do
         a(l) = t(0)
         a(n) = t(2) - t(1)
         a(0) = t(2) + t(1)
      else
         a(1) = a(0)
         a(2) = t(0)
         a(0) = t(1)
      end if
   end
!
   subroutine dfst(n, a, t, ip, w)
      integer n, ip(0:*), j, k, l, m, mh, nw, nc
      real*8 a(0:n - 1), t(0:n/2 - 1), w(0:*), xr, xi, yr, yi
      nw = ip(0)
      if (n .gt. 8*nw) then
         nw = n/8
         call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n .gt. 2*nc) then
         nc = n/2
         call makect(nc, ip, w(nw))
      end if
      if (n .gt. 2) then
         m = n/2
         mh = m/2
         do j = 1, mh - 1
            k = m - j
            xr = a(j) + a(n - j)
            xi = a(j) - a(n - j)
            yr = a(k) + a(n - k)
            yi = a(k) - a(n - k)
            a(j) = xr
            a(k) = yr
            t(j) = xi + yi
            t(k) = xi - yi
         end do
         t(0) = a(mh) - a(n - mh)
         a(mh) = a(mh) + a(n - mh)
         a(0) = a(m)
         call dstsub(m, a, nc, w(nw))
         if (m .gt. 4) then
            call cftfsub(m, a, ip, nw, w)
            call rftfsub(m, a, nc, w(nw))
         else if (m .eq. 4) then
            call cftfsub(m, a, ip, nw, w)
         end if
         a(n - 1) = a(1) - a(0)
         a(1) = a(0) + a(1)
         do j = m - 2, 2, -2
            a(2*j + 1) = a(j) - a(j + 1)
            a(2*j - 1) = -a(j) - a(j + 1)
         end do
         l = 2
         m = mh
         do while (m .ge. 2)
            call dstsub(m, t, nc, w(nw))
            if (m .gt. 4) then
               call cftfsub(m, t, ip, nw, w)
               call rftfsub(m, t, nc, w(nw))
            else if (m .eq. 4) then
               call cftfsub(m, t, ip, nw, w)
            end if
            a(n - l) = t(1) - t(0)
            a(l) = t(0) + t(1)
            k = 0
            do j = 2, m - 2, 2
               k = k + 4*l
               a(k - l) = -t(j) - t(j + 1)
               a(k + l) = t(j) - t(j + 1)
            end do
            l = 2*l
            mh = m/2
            do j = 1, mh - 1
               k = m - j
               t(j) = t(m + k) + t(m + j)
               t(k) = t(m + k) - t(m + j)
            end do
            t(0) = t(m + mh)
            m = mh
         end do
         a(l) = t(0)
      end if
      a(0) = 0
   end
!
! -------- initializing routines --------
!
   subroutine makewt(nw, ip, w)
      integer nw, ip(0:*), j, nwh, nw0, nw1
      real*8 w(0:nw - 1), delta, wn4r, wk1r, wk1i, wk3r, wk3i
      ip(0) = nw
      ip(1) = 1
      if (nw .gt. 2) then
         nwh = nw/2
         delta = atan(1.0d0)/nwh
         wn4r = cos(delta*nwh)
         w(0) = 1
         w(1) = wn4r
         if (nwh .eq. 4) then
            w(2) = cos(delta*2)
            w(3) = sin(delta*2)
         else if (nwh .gt. 4) then
            call makeipt(nw, ip)
            w(2) = 0.5d0/cos(delta*2)
            w(3) = 0.5d0/cos(delta*6)
            do j = 4, nwh - 4, 4
               w(j) = cos(delta*j)
               w(j + 1) = sin(delta*j)
               w(j + 2) = cos(3*delta*j)
               w(j + 3) = -sin(3*delta*j)
            end do
         end if
         nw0 = 0
         do while (nwh .gt. 2)
            nw1 = nw0 + nwh
            nwh = nwh/2
            w(nw1) = 1
            w(nw1 + 1) = wn4r
            if (nwh .eq. 4) then
               wk1r = w(nw0 + 4)
               wk1i = w(nw0 + 5)
               w(nw1 + 2) = wk1r
               w(nw1 + 3) = wk1i
            else if (nwh .gt. 4) then
               wk1r = w(nw0 + 4)
               wk3r = w(nw0 + 6)
               w(nw1 + 2) = 0.5d0/wk1r
               w(nw1 + 3) = 0.5d0/wk3r
               do j = 4, nwh - 4, 4
                  wk1r = w(nw0 + 2*j)
                  wk1i = w(nw0 + 2*j + 1)
                  wk3r = w(nw0 + 2*j + 2)
                  wk3i = w(nw0 + 2*j + 3)
                  w(nw1 + j) = wk1r
                  w(nw1 + j + 1) = wk1i
                  w(nw1 + j + 2) = wk3r
                  w(nw1 + j + 3) = wk3i
               end do
            end if
            nw0 = nw1
         end do
      end if
   end
!
   subroutine makeipt(nw, ip)
      integer nw, ip(0:*), j, l, m, m2, p, q
      ip(2) = 0
      ip(3) = 16
      m = 2
      l = nw
      do while (l .gt. 32)
         m2 = 2*m
         q = 8*m2
         do j = m, m2 - 1
            p = 4*ip(j)
            ip(m + j) = p
            ip(m2 + j) = p + q
         end do
         m = m2
         l = l/4
      end do
   end
!
   subroutine makect(nc, ip, c)
      integer nc, ip(0:*), j, nch
      real*8 c(0:nc - 1), delta
      ip(1) = nc
      if (nc .gt. 1) then
         nch = nc/2
         delta = atan(1.0d0)/nch
         c(0) = cos(delta*nch)
         c(nch) = 0.5d0*c(0)
         do j = 1, nch - 1
            c(j) = 0.5d0*cos(delta*j)
            c(nc - j) = 0.5d0*sin(delta*j)
         end do
      end if
   end
!
! -------- child routines --------
!
   subroutine cftfsub(n, a, ip, nw, w)
      integer n, ip(0:*), nw
      real*8 a(0:n - 1), w(0:nw - 1)
      if (n .gt. 8) then
         if (n .gt. 32) then
            call cftf1st(n, a, w(nw - n/4))
            if (n .gt. 512) then
               call cftrec4(n, a, nw, w)
            else if (n .gt. 128) then
               call cftleaf(n, 1, a, nw, w)
            else
               call cftfx41(n, a, nw, w)
            end if
            call bitrv2(n, ip, a)
         else if (n .eq. 32) then
            call cftf161(a, w(nw - 8))
            call bitrv216(a)
         else
            call cftf081(a, w)
            call bitrv208(a)
         end if
      else if (n .eq. 8) then
         call cftf040(a)
      else if (n .eq. 4) then
         call cftx020(a)
      end if
   end
!
   subroutine cftbsub(n, a, ip, nw, w)
      integer n, ip(0:*), nw
      real*8 a(0:n - 1), w(0:nw - 1)
      if (n .gt. 8) then
         if (n .gt. 32) then
            call cftb1st(n, a, w(nw - n/4))
            if (n .gt. 512) then
               call cftrec4(n, a, nw, w)
            else if (n .gt. 128) then
               call cftleaf(n, 1, a, nw, w)
            else
               call cftfx41(n, a, nw, w)
            end if
            call bitrv2conj(n, ip, a)
         else if (n .eq. 32) then
            call cftf161(a, w(nw - 8))
            call bitrv216neg(a)
         else
            call cftf081(a, w)
            call bitrv208neg(a)
         end if
      else if (n .eq. 8) then
         call cftb040(a)
      else if (n .eq. 4) then
         call cftx020(a)
      end if
   end
!
   subroutine bitrv2(n, ip, a)
      integer n, ip(0:*), j, j1, k, k1, l, m, nh, nm
      real*8 a(0:n - 1), xr, xi, yr, yi
      m = 1
      l = n/4
      do while (l .gt. 8)
         m = m*2
         l = l/4
      end do
      nh = n/2
      nm = 4*m
      if (l .eq. 8) then
         do k = 0, m - 1
            do j = 0, k - 1
               j1 = 4*j + 2*ip(m + k)
               k1 = 4*k + 2*ip(m + j)
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nm
               k1 = k1 + 2*nm
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nm
               k1 = k1 - nm
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nm
               k1 = k1 + 2*nm
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nh
               k1 = k1 + 2
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nm
               k1 = k1 - 2*nm
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nm
               k1 = k1 + nm
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nm
               k1 = k1 - 2*nm
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + 2
               k1 = k1 + nh
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nm
               k1 = k1 + 2*nm
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nm
               k1 = k1 - nm
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nm
               k1 = k1 + 2*nm
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nh
               k1 = k1 - 2
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nm
               k1 = k1 - 2*nm
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nm
               k1 = k1 + nm
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nm
               k1 = k1 - 2*nm
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
            end do
            k1 = 4*k + 2*ip(m + k)
            j1 = k1 + 2
            k1 = k1 + nh
            xr = a(j1)
            xi = a(j1 + 1)
            yr = a(k1)
            yi = a(k1 + 1)
            a(j1) = yr
            a(j1 + 1) = yi
            a(k1) = xr
            a(k1 + 1) = xi
            j1 = j1 + nm
            k1 = k1 + 2*nm
            xr = a(j1)
            xi = a(j1 + 1)
            yr = a(k1)
            yi = a(k1 + 1)
            a(j1) = yr
            a(j1 + 1) = yi
            a(k1) = xr
            a(k1 + 1) = xi
            j1 = j1 + nm
            k1 = k1 - nm
            xr = a(j1)
            xi = a(j1 + 1)
            yr = a(k1)
            yi = a(k1 + 1)
            a(j1) = yr
            a(j1 + 1) = yi
            a(k1) = xr
            a(k1 + 1) = xi
            j1 = j1 - 2
            k1 = k1 - nh
            xr = a(j1)
            xi = a(j1 + 1)
            yr = a(k1)
            yi = a(k1 + 1)
            a(j1) = yr
            a(j1 + 1) = yi
            a(k1) = xr
            a(k1 + 1) = xi
            j1 = j1 + nh + 2
            k1 = k1 + nh + 2
            xr = a(j1)
            xi = a(j1 + 1)
            yr = a(k1)
            yi = a(k1 + 1)
            a(j1) = yr
            a(j1 + 1) = yi
            a(k1) = xr
            a(k1 + 1) = xi
            j1 = j1 - nh + nm
            k1 = k1 + 2*nm - 2
            xr = a(j1)
            xi = a(j1 + 1)
            yr = a(k1)
            yi = a(k1 + 1)
            a(j1) = yr
            a(j1 + 1) = yi
            a(k1) = xr
            a(k1 + 1) = xi
         end do
      else
         do k = 0, m - 1
            do j = 0, k - 1
               j1 = 4*j + ip(m + k)
               k1 = 4*k + ip(m + j)
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nm
               k1 = k1 + nm
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nh
               k1 = k1 + 2
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nm
               k1 = k1 - nm
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + 2
               k1 = k1 + nh
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nm
               k1 = k1 + nm
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nh
               k1 = k1 - 2
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nm
               k1 = k1 - nm
               xr = a(j1)
               xi = a(j1 + 1)
               yr = a(k1)
               yi = a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
            end do
            k1 = 4*k + ip(m + k)
            j1 = k1 + 2
            k1 = k1 + nh
            xr = a(j1)
            xi = a(j1 + 1)
            yr = a(k1)
            yi = a(k1 + 1)
            a(j1) = yr
            a(j1 + 1) = yi
            a(k1) = xr
            a(k1 + 1) = xi
            j1 = j1 + nm
            k1 = k1 + nm
            xr = a(j1)
            xi = a(j1 + 1)
            yr = a(k1)
            yi = a(k1 + 1)
            a(j1) = yr
            a(j1 + 1) = yi
            a(k1) = xr
            a(k1 + 1) = xi
         end do
      end if
   end
!
   subroutine bitrv2conj(n, ip, a)
      integer n, ip(0:*), j, j1, k, k1, l, m, nh, nm
      real*8 a(0:n - 1), xr, xi, yr, yi
      m = 1
      l = n/4
      do while (l .gt. 8)
         m = m*2
         l = l/4
      end do
      nh = n/2
      nm = 4*m
      if (l .eq. 8) then
         do k = 0, m - 1
            do j = 0, k - 1
               j1 = 4*j + 2*ip(m + k)
               k1 = 4*k + 2*ip(m + j)
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nm
               k1 = k1 + 2*nm
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nm
               k1 = k1 - nm
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nm
               k1 = k1 + 2*nm
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nh
               k1 = k1 + 2
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nm
               k1 = k1 - 2*nm
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nm
               k1 = k1 + nm
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nm
               k1 = k1 - 2*nm
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + 2
               k1 = k1 + nh
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nm
               k1 = k1 + 2*nm
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nm
               k1 = k1 - nm
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nm
               k1 = k1 + 2*nm
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nh
               k1 = k1 - 2
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nm
               k1 = k1 - 2*nm
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nm
               k1 = k1 + nm
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nm
               k1 = k1 - 2*nm
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
            end do
            k1 = 4*k + 2*ip(m + k)
            j1 = k1 + 2
            k1 = k1 + nh
            a(j1 - 1) = -a(j1 - 1)
            xr = a(j1)
            xi = -a(j1 + 1)
            yr = a(k1)
            yi = -a(k1 + 1)
            a(j1) = yr
            a(j1 + 1) = yi
            a(k1) = xr
            a(k1 + 1) = xi
            a(k1 + 3) = -a(k1 + 3)
            j1 = j1 + nm
            k1 = k1 + 2*nm
            xr = a(j1)
            xi = -a(j1 + 1)
            yr = a(k1)
            yi = -a(k1 + 1)
            a(j1) = yr
            a(j1 + 1) = yi
            a(k1) = xr
            a(k1 + 1) = xi
            j1 = j1 + nm
            k1 = k1 - nm
            xr = a(j1)
            xi = -a(j1 + 1)
            yr = a(k1)
            yi = -a(k1 + 1)
            a(j1) = yr
            a(j1 + 1) = yi
            a(k1) = xr
            a(k1 + 1) = xi
            j1 = j1 - 2
            k1 = k1 - nh
            xr = a(j1)
            xi = -a(j1 + 1)
            yr = a(k1)
            yi = -a(k1 + 1)
            a(j1) = yr
            a(j1 + 1) = yi
            a(k1) = xr
            a(k1 + 1) = xi
            j1 = j1 + nh + 2
            k1 = k1 + nh + 2
            xr = a(j1)
            xi = -a(j1 + 1)
            yr = a(k1)
            yi = -a(k1 + 1)
            a(j1) = yr
            a(j1 + 1) = yi
            a(k1) = xr
            a(k1 + 1) = xi
            j1 = j1 - nh + nm
            k1 = k1 + 2*nm - 2
            a(j1 - 1) = -a(j1 - 1)
            xr = a(j1)
            xi = -a(j1 + 1)
            yr = a(k1)
            yi = -a(k1 + 1)
            a(j1) = yr
            a(j1 + 1) = yi
            a(k1) = xr
            a(k1 + 1) = xi
            a(k1 + 3) = -a(k1 + 3)
         end do
      else
         do k = 0, m - 1
            do j = 0, k - 1
               j1 = 4*j + ip(m + k)
               k1 = 4*k + ip(m + j)
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nm
               k1 = k1 + nm
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nh
               k1 = k1 + 2
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nm
               k1 = k1 - nm
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + 2
               k1 = k1 + nh
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 + nm
               k1 = k1 + nm
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nh
               k1 = k1 - 2
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
               j1 = j1 - nm
               k1 = k1 - nm
               xr = a(j1)
               xi = -a(j1 + 1)
               yr = a(k1)
               yi = -a(k1 + 1)
               a(j1) = yr
               a(j1 + 1) = yi
               a(k1) = xr
               a(k1 + 1) = xi
            end do
            k1 = 4*k + ip(m + k)
            j1 = k1 + 2
            k1 = k1 + nh
            a(j1 - 1) = -a(j1 - 1)
            xr = a(j1)
            xi = -a(j1 + 1)
            yr = a(k1)
            yi = -a(k1 + 1)
            a(j1) = yr
            a(j1 + 1) = yi
            a(k1) = xr
            a(k1 + 1) = xi
            a(k1 + 3) = -a(k1 + 3)
            j1 = j1 + nm
            k1 = k1 + nm
            a(j1 - 1) = -a(j1 - 1)
            xr = a(j1)
            xi = -a(j1 + 1)
            yr = a(k1)
            yi = -a(k1 + 1)
            a(j1) = yr
            a(j1 + 1) = yi
            a(k1) = xr
            a(k1 + 1) = xi
            a(k1 + 3) = -a(k1 + 3)
         end do
      end if
   end
!
   subroutine bitrv216(a)
      real*8 a(0:31), x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i
      real*8 x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i
      real*8 x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i
      x1r = a(2)
      x1i = a(3)
      x2r = a(4)
      x2i = a(5)
      x3r = a(6)
      x3i = a(7)
      x4r = a(8)
      x4i = a(9)
      x5r = a(10)
      x5i = a(11)
      x7r = a(14)
      x7i = a(15)
      x8r = a(16)
      x8i = a(17)
      x10r = a(20)
      x10i = a(21)
      x11r = a(22)
      x11i = a(23)
      x12r = a(24)
      x12i = a(25)
      x13r = a(26)
      x13i = a(27)
      x14r = a(28)
      x14i = a(29)
      a(2) = x8r
      a(3) = x8i
      a(4) = x4r
      a(5) = x4i
      a(6) = x12r
      a(7) = x12i
      a(8) = x2r
      a(9) = x2i
      a(10) = x10r
      a(11) = x10i
      a(14) = x14r
      a(15) = x14i
      a(16) = x1r
      a(17) = x1i
      a(20) = x5r
      a(21) = x5i
      a(22) = x13r
      a(23) = x13i
      a(24) = x3r
      a(25) = x3i
      a(26) = x11r
      a(27) = x11i
      a(28) = x7r
      a(29) = x7i
   end
!
   subroutine bitrv216neg(a)
      real*8 a(0:31), x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i
      real*8 x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i
      real*8 x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i
      real*8 x13r, x13i, x14r, x14i, x15r, x15i
      x1r = a(2)
      x1i = a(3)
      x2r = a(4)
      x2i = a(5)
      x3r = a(6)
      x3i = a(7)
      x4r = a(8)
      x4i = a(9)
      x5r = a(10)
      x5i = a(11)
      x6r = a(12)
      x6i = a(13)
      x7r = a(14)
      x7i = a(15)
      x8r = a(16)
      x8i = a(17)
      x9r = a(18)
      x9i = a(19)
      x10r = a(20)
      x10i = a(21)
      x11r = a(22)
      x11i = a(23)
      x12r = a(24)
      x12i = a(25)
      x13r = a(26)
      x13i = a(27)
      x14r = a(28)
      x14i = a(29)
      x15r = a(30)
      x15i = a(31)
      a(2) = x15r
      a(3) = x15i
      a(4) = x7r
      a(5) = x7i
      a(6) = x11r
      a(7) = x11i
      a(8) = x3r
      a(9) = x3i
      a(10) = x13r
      a(11) = x13i
      a(12) = x5r
      a(13) = x5i
      a(14) = x9r
      a(15) = x9i
      a(16) = x1r
      a(17) = x1i
      a(18) = x14r
      a(19) = x14i
      a(20) = x6r
      a(21) = x6i
      a(22) = x10r
      a(23) = x10i
      a(24) = x2r
      a(25) = x2i
      a(26) = x12r
      a(27) = x12i
      a(28) = x4r
      a(29) = x4i
      a(30) = x8r
      a(31) = x8i
   end
!
   subroutine bitrv208(a)
      real*8 a(0:15), x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i
      x1r = a(2)
      x1i = a(3)
      x3r = a(6)
      x3i = a(7)
      x4r = a(8)
      x4i = a(9)
      x6r = a(12)
      x6i = a(13)
      a(2) = x4r
      a(3) = x4i
      a(6) = x6r
      a(7) = x6i
      a(8) = x1r
      a(9) = x1i
      a(12) = x3r
      a(13) = x3i
   end
!
   subroutine bitrv208neg(a)
      real*8 a(0:15), x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i
      real*8 x5r, x5i, x6r, x6i, x7r, x7i
      x1r = a(2)
      x1i = a(3)
      x2r = a(4)
      x2i = a(5)
      x3r = a(6)
      x3i = a(7)
      x4r = a(8)
      x4i = a(9)
      x5r = a(10)
      x5i = a(11)
      x6r = a(12)
      x6i = a(13)
      x7r = a(14)
      x7i = a(15)
      a(2) = x7r
      a(3) = x7i
      a(4) = x3r
      a(5) = x3i
      a(6) = x5r
      a(7) = x5i
      a(8) = x1r
      a(9) = x1i
      a(10) = x6r
      a(11) = x6i
      a(12) = x2r
      a(13) = x2i
      a(14) = x4r
      a(15) = x4i
   end
!
   subroutine cftf1st(n, a, w)
      integer n, j, j0, j1, j2, j3, k, m, mh
      real*8 a(0:n - 1), w(0:*)
      real*8 wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i
      real*8 wd1r, wd1i, wd3r, wd3i
      real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
      real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i
      mh = n/8
      m = 2*mh
      j1 = m
      j2 = j1 + m
      j3 = j2 + m
      x0r = a(0) + a(j2)
      x0i = a(1) + a(j2 + 1)
      x1r = a(0) - a(j2)
      x1i = a(1) - a(j2 + 1)
      x2r = a(j1) + a(j3)
      x2i = a(j1 + 1) + a(j3 + 1)
      x3r = a(j1) - a(j3)
      x3i = a(j1 + 1) - a(j3 + 1)
      a(0) = x0r + x2r
      a(1) = x0i + x2i
      a(j1) = x0r - x2r
      a(j1 + 1) = x0i - x2i
      a(j2) = x1r - x3i
      a(j2 + 1) = x1i + x3r
      a(j3) = x1r + x3i
      a(j3 + 1) = x1i - x3r
      wn4r = w(1)
      csc1 = w(2)
      csc3 = w(3)
      wd1r = 1
      wd1i = 0
      wd3r = 1
      wd3i = 0
      k = 0
      do j = 2, mh - 6, 4
         k = k + 4
         wk1r = csc1*(wd1r + w(k))
         wk1i = csc1*(wd1i + w(k + 1))
         wk3r = csc3*(wd3r + w(k + 2))
         wk3i = csc3*(wd3i + w(k + 3))
         wd1r = w(k)
         wd1i = w(k + 1)
         wd3r = w(k + 2)
         wd3i = w(k + 3)
         j1 = j + m
         j2 = j1 + m
         j3 = j2 + m
         x0r = a(j) + a(j2)
         x0i = a(j + 1) + a(j2 + 1)
         x1r = a(j) - a(j2)
         x1i = a(j + 1) - a(j2 + 1)
         y0r = a(j + 2) + a(j2 + 2)
         y0i = a(j + 3) + a(j2 + 3)
         y1r = a(j + 2) - a(j2 + 2)
         y1i = a(j + 3) - a(j2 + 3)
         x2r = a(j1) + a(j3)
         x2i = a(j1 + 1) + a(j3 + 1)
         x3r = a(j1) - a(j3)
         x3i = a(j1 + 1) - a(j3 + 1)
         y2r = a(j1 + 2) + a(j3 + 2)
         y2i = a(j1 + 3) + a(j3 + 3)
         y3r = a(j1 + 2) - a(j3 + 2)
         y3i = a(j1 + 3) - a(j3 + 3)
         a(j) = x0r + x2r
         a(j + 1) = x0i + x2i
         a(j + 2) = y0r + y2r
         a(j + 3) = y0i + y2i
         a(j1) = x0r - x2r
         a(j1 + 1) = x0i - x2i
         a(j1 + 2) = y0r - y2r
         a(j1 + 3) = y0i - y2i
         x0r = x1r - x3i
         x0i = x1i + x3r
         a(j2) = wk1r*x0r - wk1i*x0i
         a(j2 + 1) = wk1r*x0i + wk1i*x0r
         x0r = y1r - y3i
         x0i = y1i + y3r
         a(j2 + 2) = wd1r*x0r - wd1i*x0i
         a(j2 + 3) = wd1r*x0i + wd1i*x0r
         x0r = x1r + x3i
         x0i = x1i - x3r
         a(j3) = wk3r*x0r + wk3i*x0i
         a(j3 + 1) = wk3r*x0i - wk3i*x0r
         x0r = y1r + y3i
         x0i = y1i - y3r
         a(j3 + 2) = wd3r*x0r + wd3i*x0i
         a(j3 + 3) = wd3r*x0i - wd3i*x0r
         j0 = m - j
         j1 = j0 + m
         j2 = j1 + m
         j3 = j2 + m
         x0r = a(j0) + a(j2)
         x0i = a(j0 + 1) + a(j2 + 1)
         x1r = a(j0) - a(j2)
         x1i = a(j0 + 1) - a(j2 + 1)
         y0r = a(j0 - 2) + a(j2 - 2)
         y0i = a(j0 - 1) + a(j2 - 1)
         y1r = a(j0 - 2) - a(j2 - 2)
         y1i = a(j0 - 1) - a(j2 - 1)
         x2r = a(j1) + a(j3)
         x2i = a(j1 + 1) + a(j3 + 1)
         x3r = a(j1) - a(j3)
         x3i = a(j1 + 1) - a(j3 + 1)
         y2r = a(j1 - 2) + a(j3 - 2)
         y2i = a(j1 - 1) + a(j3 - 1)
         y3r = a(j1 - 2) - a(j3 - 2)
         y3i = a(j1 - 1) - a(j3 - 1)
         a(j0) = x0r + x2r
         a(j0 + 1) = x0i + x2i
         a(j0 - 2) = y0r + y2r
         a(j0 - 1) = y0i + y2i
         a(j1) = x0r - x2r
         a(j1 + 1) = x0i - x2i
         a(j1 - 2) = y0r - y2r
         a(j1 - 1) = y0i - y2i
         x0r = x1r - x3i
         x0i = x1i + x3r
         a(j2) = wk1i*x0r - wk1r*x0i
         a(j2 + 1) = wk1i*x0i + wk1r*x0r
         x0r = y1r - y3i
         x0i = y1i + y3r
         a(j2 - 2) = wd1i*x0r - wd1r*x0i
         a(j2 - 1) = wd1i*x0i + wd1r*x0r
         x0r = x1r + x3i
         x0i = x1i - x3r
         a(j3) = wk3i*x0r + wk3r*x0i
         a(j3 + 1) = wk3i*x0i - wk3r*x0r
         x0r = y1r + y3i
         x0i = y1i - y3r
         a(j3 - 2) = wd3i*x0r + wd3r*x0i
         a(j3 - 1) = wd3i*x0i - wd3r*x0r
      end do
      wk1r = csc1*(wd1r + wn4r)
      wk1i = csc1*(wd1i + wn4r)
      wk3r = csc3*(wd3r - wn4r)
      wk3i = csc3*(wd3i - wn4r)
      j0 = mh
      j1 = j0 + m
      j2 = j1 + m
      j3 = j2 + m
      x0r = a(j0 - 2) + a(j2 - 2)
      x0i = a(j0 - 1) + a(j2 - 1)
      x1r = a(j0 - 2) - a(j2 - 2)
      x1i = a(j0 - 1) - a(j2 - 1)
      x2r = a(j1 - 2) + a(j3 - 2)
      x2i = a(j1 - 1) + a(j3 - 1)
      x3r = a(j1 - 2) - a(j3 - 2)
      x3i = a(j1 - 1) - a(j3 - 1)
      a(j0 - 2) = x0r + x2r
      a(j0 - 1) = x0i + x2i
      a(j1 - 2) = x0r - x2r
      a(j1 - 1) = x0i - x2i
      x0r = x1r - x3i
      x0i = x1i + x3r
      a(j2 - 2) = wk1r*x0r - wk1i*x0i
      a(j2 - 1) = wk1r*x0i + wk1i*x0r
      x0r = x1r + x3i
      x0i = x1i - x3r
      a(j3 - 2) = wk3r*x0r + wk3i*x0i
      a(j3 - 1) = wk3r*x0i - wk3i*x0r
      x0r = a(j0) + a(j2)
      x0i = a(j0 + 1) + a(j2 + 1)
      x1r = a(j0) - a(j2)
      x1i = a(j0 + 1) - a(j2 + 1)
      x2r = a(j1) + a(j3)
      x2i = a(j1 + 1) + a(j3 + 1)
      x3r = a(j1) - a(j3)
      x3i = a(j1 + 1) - a(j3 + 1)
      a(j0) = x0r + x2r
      a(j0 + 1) = x0i + x2i
      a(j1) = x0r - x2r
      a(j1 + 1) = x0i - x2i
      x0r = x1r - x3i
      x0i = x1i + x3r
      a(j2) = wn4r*(x0r - x0i)
      a(j2 + 1) = wn4r*(x0i + x0r)
      x0r = x1r + x3i
      x0i = x1i - x3r
      a(j3) = -wn4r*(x0r + x0i)
      a(j3 + 1) = -wn4r*(x0i - x0r)
      x0r = a(j0 + 2) + a(j2 + 2)
      x0i = a(j0 + 3) + a(j2 + 3)
      x1r = a(j0 + 2) - a(j2 + 2)
      x1i = a(j0 + 3) - a(j2 + 3)
      x2r = a(j1 + 2) + a(j3 + 2)
      x2i = a(j1 + 3) + a(j3 + 3)
      x3r = a(j1 + 2) - a(j3 + 2)
      x3i = a(j1 + 3) - a(j3 + 3)
      a(j0 + 2) = x0r + x2r
      a(j0 + 3) = x0i + x2i
      a(j1 + 2) = x0r - x2r
      a(j1 + 3) = x0i - x2i
      x0r = x1r - x3i
      x0i = x1i + x3r
      a(j2 + 2) = wk1i*x0r - wk1r*x0i
      a(j2 + 3) = wk1i*x0i + wk1r*x0r
      x0r = x1r + x3i
      x0i = x1i - x3r
      a(j3 + 2) = wk3i*x0r + wk3r*x0i
      a(j3 + 3) = wk3i*x0i - wk3r*x0r
   end
!
   subroutine cftb1st(n, a, w)
      integer n, j, j0, j1, j2, j3, k, m, mh
      real*8 a(0:n - 1), w(0:*)
      real*8 wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i
      real*8 wd1r, wd1i, wd3r, wd3i
      real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
      real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i
      mh = n/8
      m = 2*mh
      j1 = m
      j2 = j1 + m
      j3 = j2 + m
      x0r = a(0) + a(j2)
      x0i = -a(1) - a(j2 + 1)
      x1r = a(0) - a(j2)
      x1i = -a(1) + a(j2 + 1)
      x2r = a(j1) + a(j3)
      x2i = a(j1 + 1) + a(j3 + 1)
      x3r = a(j1) - a(j3)
      x3i = a(j1 + 1) - a(j3 + 1)
      a(0) = x0r + x2r
      a(1) = x0i - x2i
      a(j1) = x0r - x2r
      a(j1 + 1) = x0i + x2i
      a(j2) = x1r + x3i
      a(j2 + 1) = x1i + x3r
      a(j3) = x1r - x3i
      a(j3 + 1) = x1i - x3r
      wn4r = w(1)
      csc1 = w(2)
      csc3 = w(3)
      wd1r = 1
      wd1i = 0
      wd3r = 1
      wd3i = 0
      k = 0
      do j = 2, mh - 6, 4
         k = k + 4
         wk1r = csc1*(wd1r + w(k))
         wk1i = csc1*(wd1i + w(k + 1))
         wk3r = csc3*(wd3r + w(k + 2))
         wk3i = csc3*(wd3i + w(k + 3))
         wd1r = w(k)
         wd1i = w(k + 1)
         wd3r = w(k + 2)
         wd3i = w(k + 3)
         j1 = j + m
         j2 = j1 + m
         j3 = j2 + m
         x0r = a(j) + a(j2)
         x0i = -a(j + 1) - a(j2 + 1)
         x1r = a(j) - a(j2)
         x1i = -a(j + 1) + a(j2 + 1)
         y0r = a(j + 2) + a(j2 + 2)
         y0i = -a(j + 3) - a(j2 + 3)
         y1r = a(j + 2) - a(j2 + 2)
         y1i = -a(j + 3) + a(j2 + 3)
         x2r = a(j1) + a(j3)
         x2i = a(j1 + 1) + a(j3 + 1)
         x3r = a(j1) - a(j3)
         x3i = a(j1 + 1) - a(j3 + 1)
         y2r = a(j1 + 2) + a(j3 + 2)
         y2i = a(j1 + 3) + a(j3 + 3)
         y3r = a(j1 + 2) - a(j3 + 2)
         y3i = a(j1 + 3) - a(j3 + 3)
         a(j) = x0r + x2r
         a(j + 1) = x0i - x2i
         a(j + 2) = y0r + y2r
         a(j + 3) = y0i - y2i
         a(j1) = x0r - x2r
         a(j1 + 1) = x0i + x2i
         a(j1 + 2) = y0r - y2r
         a(j1 + 3) = y0i + y2i
         x0r = x1r + x3i
         x0i = x1i + x3r
         a(j2) = wk1r*x0r - wk1i*x0i
         a(j2 + 1) = wk1r*x0i + wk1i*x0r
         x0r = y1r + y3i
         x0i = y1i + y3r
         a(j2 + 2) = wd1r*x0r - wd1i*x0i
         a(j2 + 3) = wd1r*x0i + wd1i*x0r
         x0r = x1r - x3i
         x0i = x1i - x3r
         a(j3) = wk3r*x0r + wk3i*x0i
         a(j3 + 1) = wk3r*x0i - wk3i*x0r
         x0r = y1r - y3i
         x0i = y1i - y3r
         a(j3 + 2) = wd3r*x0r + wd3i*x0i
         a(j3 + 3) = wd3r*x0i - wd3i*x0r
         j0 = m - j
         j1 = j0 + m
         j2 = j1 + m
         j3 = j2 + m
         x0r = a(j0) + a(j2)
         x0i = -a(j0 + 1) - a(j2 + 1)
         x1r = a(j0) - a(j2)
         x1i = -a(j0 + 1) + a(j2 + 1)
         y0r = a(j0 - 2) + a(j2 - 2)
         y0i = -a(j0 - 1) - a(j2 - 1)
         y1r = a(j0 - 2) - a(j2 - 2)
         y1i = -a(j0 - 1) + a(j2 - 1)
         x2r = a(j1) + a(j3)
         x2i = a(j1 + 1) + a(j3 + 1)
         x3r = a(j1) - a(j3)
         x3i = a(j1 + 1) - a(j3 + 1)
         y2r = a(j1 - 2) + a(j3 - 2)
         y2i = a(j1 - 1) + a(j3 - 1)
         y3r = a(j1 - 2) - a(j3 - 2)
         y3i = a(j1 - 1) - a(j3 - 1)
         a(j0) = x0r + x2r
         a(j0 + 1) = x0i - x2i
         a(j0 - 2) = y0r + y2r
         a(j0 - 1) = y0i - y2i
         a(j1) = x0r - x2r
         a(j1 + 1) = x0i + x2i
         a(j1 - 2) = y0r - y2r
         a(j1 - 1) = y0i + y2i
         x0r = x1r + x3i
         x0i = x1i + x3r
         a(j2) = wk1i*x0r - wk1r*x0i
         a(j2 + 1) = wk1i*x0i + wk1r*x0r
         x0r = y1r + y3i
         x0i = y1i + y3r
         a(j2 - 2) = wd1i*x0r - wd1r*x0i
         a(j2 - 1) = wd1i*x0i + wd1r*x0r
         x0r = x1r - x3i
         x0i = x1i - x3r
         a(j3) = wk3i*x0r + wk3r*x0i
         a(j3 + 1) = wk3i*x0i - wk3r*x0r
         x0r = y1r - y3i
         x0i = y1i - y3r
         a(j3 - 2) = wd3i*x0r + wd3r*x0i
         a(j3 - 1) = wd3i*x0i - wd3r*x0r
      end do
      wk1r = csc1*(wd1r + wn4r)
      wk1i = csc1*(wd1i + wn4r)
      wk3r = csc3*(wd3r - wn4r)
      wk3i = csc3*(wd3i - wn4r)
      j0 = mh
      j1 = j0 + m
      j2 = j1 + m
      j3 = j2 + m
      x0r = a(j0 - 2) + a(j2 - 2)
      x0i = -a(j0 - 1) - a(j2 - 1)
      x1r = a(j0 - 2) - a(j2 - 2)
      x1i = -a(j0 - 1) + a(j2 - 1)
      x2r = a(j1 - 2) + a(j3 - 2)
      x2i = a(j1 - 1) + a(j3 - 1)
      x3r = a(j1 - 2) - a(j3 - 2)
      x3i = a(j1 - 1) - a(j3 - 1)
      a(j0 - 2) = x0r + x2r
      a(j0 - 1) = x0i - x2i
      a(j1 - 2) = x0r - x2r
      a(j1 - 1) = x0i + x2i
      x0r = x1r + x3i
      x0i = x1i + x3r
      a(j2 - 2) = wk1r*x0r - wk1i*x0i
      a(j2 - 1) = wk1r*x0i + wk1i*x0r
      x0r = x1r - x3i
      x0i = x1i - x3r
      a(j3 - 2) = wk3r*x0r + wk3i*x0i
      a(j3 - 1) = wk3r*x0i - wk3i*x0r
      x0r = a(j0) + a(j2)
      x0i = -a(j0 + 1) - a(j2 + 1)
      x1r = a(j0) - a(j2)
      x1i = -a(j0 + 1) + a(j2 + 1)
      x2r = a(j1) + a(j3)
      x2i = a(j1 + 1) + a(j3 + 1)
      x3r = a(j1) - a(j3)
      x3i = a(j1 + 1) - a(j3 + 1)
      a(j0) = x0r + x2r
      a(j0 + 1) = x0i - x2i
      a(j1) = x0r - x2r
      a(j1 + 1) = x0i + x2i
      x0r = x1r + x3i
      x0i = x1i + x3r
      a(j2) = wn4r*(x0r - x0i)
      a(j2 + 1) = wn4r*(x0i + x0r)
      x0r = x1r - x3i
      x0i = x1i - x3r
      a(j3) = -wn4r*(x0r + x0i)
      a(j3 + 1) = -wn4r*(x0i - x0r)
      x0r = a(j0 + 2) + a(j2 + 2)
      x0i = -a(j0 + 3) - a(j2 + 3)
      x1r = a(j0 + 2) - a(j2 + 2)
      x1i = -a(j0 + 3) + a(j2 + 3)
      x2r = a(j1 + 2) + a(j3 + 2)
      x2i = a(j1 + 3) + a(j3 + 3)
      x3r = a(j1 + 2) - a(j3 + 2)
      x3i = a(j1 + 3) - a(j3 + 3)
      a(j0 + 2) = x0r + x2r
      a(j0 + 3) = x0i - x2i
      a(j1 + 2) = x0r - x2r
      a(j1 + 3) = x0i + x2i
      x0r = x1r + x3i
      x0i = x1i + x3r
      a(j2 + 2) = wk1i*x0r - wk1r*x0i
      a(j2 + 3) = wk1i*x0i + wk1r*x0r
      x0r = x1r - x3i
      x0i = x1i - x3r
      a(j3 + 2) = wk3i*x0r + wk3r*x0i
      a(j3 + 3) = wk3i*x0i - wk3r*x0r
   end
!
   subroutine cftrec4(n, a, nw, w)
      integer n, nw, cfttree, isplt, j, k, m
      real*8 a(0:n - 1), w(0:nw - 1)
      m = n
      do while (m .gt. 512)
         m = m/4
         call cftmdl1(m, a(n - m), w(nw - m/2))
      end do
      call cftleaf(m, 1, a(n - m), nw, w)
      k = 0
      do j = n - m, m, -m
         k = k + 1
         isplt = cfttree(m, j, k, a, nw, w)
         call cftleaf(m, isplt, a(j - m), nw, w)
      end do
   end
!
   integer function cfttree(n, j, k, a, nw, w)
      integer n, j, k, nw, i, isplt, m
      real*8 a(0:j - 1), w(0:nw - 1)
      if (mod(k, 4) .ne. 0) then
         isplt = mod(k, 2)
         if (isplt .ne. 0) then
            call cftmdl1(n, a(j - n), w(nw - n/2))
         else
            call cftmdl2(n, a(j - n), w(nw - n))
         end if
      else
         m = n
         i = k
         do while (mod(i, 4) .eq. 0)
            m = m*4
            i = i/4
         end do
         isplt = mod(i, 2)
         if (isplt .ne. 0) then
            do while (m .gt. 128)
               call cftmdl1(m, a(j - m), w(nw - m/2))
               m = m/4
            end do
         else
            do while (m .gt. 128)
               call cftmdl2(m, a(j - m), w(nw - m))
               m = m/4
            end do
         end if
      end if
      cfttree = isplt
   end
!
   subroutine cftleaf(n, isplt, a, nw, w)
      integer n, isplt, nw
      real*8 a(0:n - 1), w(0:nw - 1)
      if (n .eq. 512) then
         call cftmdl1(128, a, w(nw - 64))
         call cftf161(a, w(nw - 8))
         call cftf162(a(32), w(nw - 32))
         call cftf161(a(64), w(nw - 8))
         call cftf161(a(96), w(nw - 8))
         call cftmdl2(128, a(128), w(nw - 128))
         call cftf161(a(128), w(nw - 8))
         call cftf162(a(160), w(nw - 32))
         call cftf161(a(192), w(nw - 8))
         call cftf162(a(224), w(nw - 32))
         call cftmdl1(128, a(256), w(nw - 64))
         call cftf161(a(256), w(nw - 8))
         call cftf162(a(288), w(nw - 32))
         call cftf161(a(320), w(nw - 8))
         call cftf161(a(352), w(nw - 8))
         if (isplt .ne. 0) then
            call cftmdl1(128, a(384), w(nw - 64))
            call cftf161(a(480), w(nw - 8))
         else
            call cftmdl2(128, a(384), w(nw - 128))
            call cftf162(a(480), w(nw - 32))
         end if
         call cftf161(a(384), w(nw - 8))
         call cftf162(a(416), w(nw - 32))
         call cftf161(a(448), w(nw - 8))
      else
         call cftmdl1(64, a, w(nw - 32))
         call cftf081(a, w(nw - 8))
         call cftf082(a(16), w(nw - 8))
         call cftf081(a(32), w(nw - 8))
         call cftf081(a(48), w(nw - 8))
         call cftmdl2(64, a(64), w(nw - 64))
         call cftf081(a(64), w(nw - 8))
         call cftf082(a(80), w(nw - 8))
         call cftf081(a(96), w(nw - 8))
         call cftf082(a(112), w(nw - 8))
         call cftmdl1(64, a(128), w(nw - 32))
         call cftf081(a(128), w(nw - 8))
         call cftf082(a(144), w(nw - 8))
         call cftf081(a(160), w(nw - 8))
         call cftf081(a(176), w(nw - 8))
         if (isplt .ne. 0) then
            call cftmdl1(64, a(192), w(nw - 32))
            call cftf081(a(240), w(nw - 8))
         else
            call cftmdl2(64, a(192), w(nw - 64))
            call cftf082(a(240), w(nw - 8))
         end if
         call cftf081(a(192), w(nw - 8))
         call cftf082(a(208), w(nw - 8))
         call cftf081(a(224), w(nw - 8))
      end if
   end
!
   subroutine cftmdl1(n, a, w)
      integer n, j, j0, j1, j2, j3, k, m, mh
      real*8 a(0:n - 1), w(0:*)
      real*8 wn4r, wk1r, wk1i, wk3r, wk3i
      real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
      mh = n/8
      m = 2*mh
      j1 = m
      j2 = j1 + m
      j3 = j2 + m
      x0r = a(0) + a(j2)
      x0i = a(1) + a(j2 + 1)
      x1r = a(0) - a(j2)
      x1i = a(1) - a(j2 + 1)
      x2r = a(j1) + a(j3)
      x2i = a(j1 + 1) + a(j3 + 1)
      x3r = a(j1) - a(j3)
      x3i = a(j1 + 1) - a(j3 + 1)
      a(0) = x0r + x2r
      a(1) = x0i + x2i
      a(j1) = x0r - x2r
      a(j1 + 1) = x0i - x2i
      a(j2) = x1r - x3i
      a(j2 + 1) = x1i + x3r
      a(j3) = x1r + x3i
      a(j3 + 1) = x1i - x3r
      wn4r = w(1)
      k = 0
      do j = 2, mh - 2, 2
         k = k + 4
         wk1r = w(k)
         wk1i = w(k + 1)
         wk3r = w(k + 2)
         wk3i = w(k + 3)
         j1 = j + m
         j2 = j1 + m
         j3 = j2 + m
         x0r = a(j) + a(j2)
         x0i = a(j + 1) + a(j2 + 1)
         x1r = a(j) - a(j2)
         x1i = a(j + 1) - a(j2 + 1)
         x2r = a(j1) + a(j3)
         x2i = a(j1 + 1) + a(j3 + 1)
         x3r = a(j1) - a(j3)
         x3i = a(j1 + 1) - a(j3 + 1)
         a(j) = x0r + x2r
         a(j + 1) = x0i + x2i
         a(j1) = x0r - x2r
         a(j1 + 1) = x0i - x2i
         x0r = x1r - x3i
         x0i = x1i + x3r
         a(j2) = wk1r*x0r - wk1i*x0i
         a(j2 + 1) = wk1r*x0i + wk1i*x0r
         x0r = x1r + x3i
         x0i = x1i - x3r
         a(j3) = wk3r*x0r + wk3i*x0i
         a(j3 + 1) = wk3r*x0i - wk3i*x0r
         j0 = m - j
         j1 = j0 + m
         j2 = j1 + m
         j3 = j2 + m
         x0r = a(j0) + a(j2)
         x0i = a(j0 + 1) + a(j2 + 1)
         x1r = a(j0) - a(j2)
         x1i = a(j0 + 1) - a(j2 + 1)
         x2r = a(j1) + a(j3)
         x2i = a(j1 + 1) + a(j3 + 1)
         x3r = a(j1) - a(j3)
         x3i = a(j1 + 1) - a(j3 + 1)
         a(j0) = x0r + x2r
         a(j0 + 1) = x0i + x2i
         a(j1) = x0r - x2r
         a(j1 + 1) = x0i - x2i
         x0r = x1r - x3i
         x0i = x1i + x3r
         a(j2) = wk1i*x0r - wk1r*x0i
         a(j2 + 1) = wk1i*x0i + wk1r*x0r
         x0r = x1r + x3i
         x0i = x1i - x3r
         a(j3) = wk3i*x0r + wk3r*x0i
         a(j3 + 1) = wk3i*x0i - wk3r*x0r
      end do
      j0 = mh
      j1 = j0 + m
      j2 = j1 + m
      j3 = j2 + m
      x0r = a(j0) + a(j2)
      x0i = a(j0 + 1) + a(j2 + 1)
      x1r = a(j0) - a(j2)
      x1i = a(j0 + 1) - a(j2 + 1)
      x2r = a(j1) + a(j3)
      x2i = a(j1 + 1) + a(j3 + 1)
      x3r = a(j1) - a(j3)
      x3i = a(j1 + 1) - a(j3 + 1)
      a(j0) = x0r + x2r
      a(j0 + 1) = x0i + x2i
      a(j1) = x0r - x2r
      a(j1 + 1) = x0i - x2i
      x0r = x1r - x3i
      x0i = x1i + x3r
      a(j2) = wn4r*(x0r - x0i)
      a(j2 + 1) = wn4r*(x0i + x0r)
      x0r = x1r + x3i
      x0i = x1i - x3r
      a(j3) = -wn4r*(x0r + x0i)
      a(j3 + 1) = -wn4r*(x0i - x0r)
   end
!
   subroutine cftmdl2(n, a, w)
      integer n, j, j0, j1, j2, j3, k, kr, m, mh
      real*8 a(0:n - 1), w(0:*)
      real*8 wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i
      real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
      real*8 y0r, y0i, y2r, y2i
      mh = n/8
      m = 2*mh
      wn4r = w(1)
      j1 = m
      j2 = j1 + m
      j3 = j2 + m
      x0r = a(0) - a(j2 + 1)
      x0i = a(1) + a(j2)
      x1r = a(0) + a(j2 + 1)
      x1i = a(1) - a(j2)
      x2r = a(j1) - a(j3 + 1)
      x2i = a(j1 + 1) + a(j3)
      x3r = a(j1) + a(j3 + 1)
      x3i = a(j1 + 1) - a(j3)
      y0r = wn4r*(x2r - x2i)
      y0i = wn4r*(x2i + x2r)
      a(0) = x0r + y0r
      a(1) = x0i + y0i
      a(j1) = x0r - y0r
      a(j1 + 1) = x0i - y0i
      y0r = wn4r*(x3r - x3i)
      y0i = wn4r*(x3i + x3r)
      a(j2) = x1r - y0i
      a(j2 + 1) = x1i + y0r
      a(j3) = x1r + y0i
      a(j3 + 1) = x1i - y0r
      k = 0
      kr = 2*m
      do j = 2, mh - 2, 2
         k = k + 4
         wk1r = w(k)
         wk1i = w(k + 1)
         wk3r = w(k + 2)
         wk3i = w(k + 3)
         kr = kr - 4
         wd1i = w(kr)
         wd1r = w(kr + 1)
         wd3i = w(kr + 2)
         wd3r = w(kr + 3)
         j1 = j + m
         j2 = j1 + m
         j3 = j2 + m
         x0r = a(j) - a(j2 + 1)
         x0i = a(j + 1) + a(j2)
         x1r = a(j) + a(j2 + 1)
         x1i = a(j + 1) - a(j2)
         x2r = a(j1) - a(j3 + 1)
         x2i = a(j1 + 1) + a(j3)
         x3r = a(j1) + a(j3 + 1)
         x3i = a(j1 + 1) - a(j3)
         y0r = wk1r*x0r - wk1i*x0i
         y0i = wk1r*x0i + wk1i*x0r
         y2r = wd1r*x2r - wd1i*x2i
         y2i = wd1r*x2i + wd1i*x2r
         a(j) = y0r + y2r
         a(j + 1) = y0i + y2i
         a(j1) = y0r - y2r
         a(j1 + 1) = y0i - y2i
         y0r = wk3r*x1r + wk3i*x1i
         y0i = wk3r*x1i - wk3i*x1r
         y2r = wd3r*x3r + wd3i*x3i
         y2i = wd3r*x3i - wd3i*x3r
         a(j2) = y0r + y2r
         a(j2 + 1) = y0i + y2i
         a(j3) = y0r - y2r
         a(j3 + 1) = y0i - y2i
         j0 = m - j
         j1 = j0 + m
         j2 = j1 + m
         j3 = j2 + m
         x0r = a(j0) - a(j2 + 1)
         x0i = a(j0 + 1) + a(j2)
         x1r = a(j0) + a(j2 + 1)
         x1i = a(j0 + 1) - a(j2)
         x2r = a(j1) - a(j3 + 1)
         x2i = a(j1 + 1) + a(j3)
         x3r = a(j1) + a(j3 + 1)
         x3i = a(j1 + 1) - a(j3)
         y0r = wd1i*x0r - wd1r*x0i
         y0i = wd1i*x0i + wd1r*x0r
         y2r = wk1i*x2r - wk1r*x2i
         y2i = wk1i*x2i + wk1r*x2r
         a(j0) = y0r + y2r
         a(j0 + 1) = y0i + y2i
         a(j1) = y0r - y2r
         a(j1 + 1) = y0i - y2i
         y0r = wd3i*x1r + wd3r*x1i
         y0i = wd3i*x1i - wd3r*x1r
         y2r = wk3i*x3r + wk3r*x3i
         y2i = wk3i*x3i - wk3r*x3r
         a(j2) = y0r + y2r
         a(j2 + 1) = y0i + y2i
         a(j3) = y0r - y2r
         a(j3 + 1) = y0i - y2i
      end do
      wk1r = w(m)
      wk1i = w(m + 1)
      j0 = mh
      j1 = j0 + m
      j2 = j1 + m
      j3 = j2 + m
      x0r = a(j0) - a(j2 + 1)
      x0i = a(j0 + 1) + a(j2)
      x1r = a(j0) + a(j2 + 1)
      x1i = a(j0 + 1) - a(j2)
      x2r = a(j1) - a(j3 + 1)
      x2i = a(j1 + 1) + a(j3)
      x3r = a(j1) + a(j3 + 1)
      x3i = a(j1 + 1) - a(j3)
      y0r = wk1r*x0r - wk1i*x0i
      y0i = wk1r*x0i + wk1i*x0r
      y2r = wk1i*x2r - wk1r*x2i
      y2i = wk1i*x2i + wk1r*x2r
      a(j0) = y0r + y2r
      a(j0 + 1) = y0i + y2i
      a(j1) = y0r - y2r
      a(j1 + 1) = y0i - y2i
      y0r = wk1i*x1r - wk1r*x1i
      y0i = wk1i*x1i + wk1r*x1r
      y2r = wk1r*x3r - wk1i*x3i
      y2i = wk1r*x3i + wk1i*x3r
      a(j2) = y0r - y2r
      a(j2 + 1) = y0i - y2i
      a(j3) = y0r + y2r
      a(j3 + 1) = y0i + y2i
   end
!
   subroutine cftfx41(n, a, nw, w)
      integer n, nw
      real*8 a(0:n - 1), w(0:nw - 1)
      if (n .eq. 128) then
         call cftf161(a, w(nw - 8))
         call cftf162(a(32), w(nw - 32))
         call cftf161(a(64), w(nw - 8))
         call cftf161(a(96), w(nw - 8))
      else
         call cftf081(a, w(nw - 8))
         call cftf082(a(16), w(nw - 8))
         call cftf081(a(32), w(nw - 8))
         call cftf081(a(48), w(nw - 8))
      end if
   end
!
   subroutine cftf161(a, w)
      real*8 a(0:31), w(0:*), wn4r, wk1r, wk1i
      real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
      real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i
      real*8 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i
      real*8 y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i
      real*8 y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i
      wn4r = w(1)
      wk1r = w(2)
      wk1i = w(3)
      x0r = a(0) + a(16)
      x0i = a(1) + a(17)
      x1r = a(0) - a(16)
      x1i = a(1) - a(17)
      x2r = a(8) + a(24)
      x2i = a(9) + a(25)
      x3r = a(8) - a(24)
      x3i = a(9) - a(25)
      y0r = x0r + x2r
      y0i = x0i + x2i
      y4r = x0r - x2r
      y4i = x0i - x2i
      y8r = x1r - x3i
      y8i = x1i + x3r
      y12r = x1r + x3i
      y12i = x1i - x3r
      x0r = a(2) + a(18)
      x0i = a(3) + a(19)
      x1r = a(2) - a(18)
      x1i = a(3) - a(19)
      x2r = a(10) + a(26)
      x2i = a(11) + a(27)
      x3r = a(10) - a(26)
      x3i = a(11) - a(27)
      y1r = x0r + x2r
      y1i = x0i + x2i
      y5r = x0r - x2r
      y5i = x0i - x2i
      x0r = x1r - x3i
      x0i = x1i + x3r
      y9r = wk1r*x0r - wk1i*x0i
      y9i = wk1r*x0i + wk1i*x0r
      x0r = x1r + x3i
      x0i = x1i - x3r
      y13r = wk1i*x0r - wk1r*x0i
      y13i = wk1i*x0i + wk1r*x0r
      x0r = a(4) + a(20)
      x0i = a(5) + a(21)
      x1r = a(4) - a(20)
      x1i = a(5) - a(21)
      x2r = a(12) + a(28)
      x2i = a(13) + a(29)
      x3r = a(12) - a(28)
      x3i = a(13) - a(29)
      y2r = x0r + x2r
      y2i = x0i + x2i
      y6r = x0r - x2r
      y6i = x0i - x2i
      x0r = x1r - x3i
      x0i = x1i + x3r
      y10r = wn4r*(x0r - x0i)
      y10i = wn4r*(x0i + x0r)
      x0r = x1r + x3i
      x0i = x1i - x3r
      y14r = wn4r*(x0r + x0i)
      y14i = wn4r*(x0i - x0r)
      x0r = a(6) + a(22)
      x0i = a(7) + a(23)
      x1r = a(6) - a(22)
      x1i = a(7) - a(23)
      x2r = a(14) + a(30)
      x2i = a(15) + a(31)
      x3r = a(14) - a(30)
      x3i = a(15) - a(31)
      y3r = x0r + x2r
      y3i = x0i + x2i
      y7r = x0r - x2r
      y7i = x0i - x2i
      x0r = x1r - x3i
      x0i = x1i + x3r
      y11r = wk1i*x0r - wk1r*x0i
      y11i = wk1i*x0i + wk1r*x0r
      x0r = x1r + x3i
      x0i = x1i - x3r
      y15r = wk1r*x0r - wk1i*x0i
      y15i = wk1r*x0i + wk1i*x0r
      x0r = y12r - y14r
      x0i = y12i - y14i
      x1r = y12r + y14r
      x1i = y12i + y14i
      x2r = y13r - y15r
      x2i = y13i - y15i
      x3r = y13r + y15r
      x3i = y13i + y15i
      a(24) = x0r + x2r
      a(25) = x0i + x2i
      a(26) = x0r - x2r
      a(27) = x0i - x2i
      a(28) = x1r - x3i
      a(29) = x1i + x3r
      a(30) = x1r + x3i
      a(31) = x1i - x3r
      x0r = y8r + y10r
      x0i = y8i + y10i
      x1r = y8r - y10r
      x1i = y8i - y10i
      x2r = y9r + y11r
      x2i = y9i + y11i
      x3r = y9r - y11r
      x3i = y9i - y11i
      a(16) = x0r + x2r
      a(17) = x0i + x2i
      a(18) = x0r - x2r
      a(19) = x0i - x2i
      a(20) = x1r - x3i
      a(21) = x1i + x3r
      a(22) = x1r + x3i
      a(23) = x1i - x3r
      x0r = y5r - y7i
      x0i = y5i + y7r
      x2r = wn4r*(x0r - x0i)
      x2i = wn4r*(x0i + x0r)
      x0r = y5r + y7i
      x0i = y5i - y7r
      x3r = wn4r*(x0r - x0i)
      x3i = wn4r*(x0i + x0r)
      x0r = y4r - y6i
      x0i = y4i + y6r
      x1r = y4r + y6i
      x1i = y4i - y6r
      a(8) = x0r + x2r
      a(9) = x0i + x2i
      a(10) = x0r - x2r
      a(11) = x0i - x2i
      a(12) = x1r - x3i
      a(13) = x1i + x3r
      a(14) = x1r + x3i
      a(15) = x1i - x3r
      x0r = y0r + y2r
      x0i = y0i + y2i
      x1r = y0r - y2r
      x1i = y0i - y2i
      x2r = y1r + y3r
      x2i = y1i + y3i
      x3r = y1r - y3r
      x3i = y1i - y3i
      a(0) = x0r + x2r
      a(1) = x0i + x2i
      a(2) = x0r - x2r
      a(3) = x0i - x2i
      a(4) = x1r - x3i
      a(5) = x1i + x3r
      a(6) = x1r + x3i
      a(7) = x1i - x3r
   end
!
   subroutine cftf162(a, w)
      real*8 a(0:31), w(0:*)
      real*8 wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i
      real*8 x0r, x0i, x1r, x1i, x2r, x2i
      real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i
      real*8 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i
      real*8 y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i
      real*8 y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i
      wn4r = w(1)
      wk1r = w(4)
      wk1i = w(5)
      wk3r = w(6)
      wk3i = -w(7)
      wk2r = w(8)
      wk2i = w(9)
      x1r = a(0) - a(17)
      x1i = a(1) + a(16)
      x0r = a(8) - a(25)
      x0i = a(9) + a(24)
      x2r = wn4r*(x0r - x0i)
      x2i = wn4r*(x0i + x0r)
      y0r = x1r + x2r
      y0i = x1i + x2i
      y4r = x1r - x2r
      y4i = x1i - x2i
      x1r = a(0) + a(17)
      x1i = a(1) - a(16)
      x0r = a(8) + a(25)
      x0i = a(9) - a(24)
      x2r = wn4r*(x0r - x0i)
      x2i = wn4r*(x0i + x0r)
      y8r = x1r - x2i
      y8i = x1i + x2r
      y12r = x1r + x2i
      y12i = x1i - x2r
      x0r = a(2) - a(19)
      x0i = a(3) + a(18)
      x1r = wk1r*x0r - wk1i*x0i
      x1i = wk1r*x0i + wk1i*x0r
      x0r = a(10) - a(27)
      x0i = a(11) + a(26)
      x2r = wk3i*x0r - wk3r*x0i
      x2i = wk3i*x0i + wk3r*x0r
      y1r = x1r + x2r
      y1i = x1i + x2i
      y5r = x1r - x2r
      y5i = x1i - x2i
      x0r = a(2) + a(19)
      x0i = a(3) - a(18)
      x1r = wk3r*x0r - wk3i*x0i
      x1i = wk3r*x0i + wk3i*x0r
      x0r = a(10) + a(27)
      x0i = a(11) - a(26)
      x2r = wk1r*x0r + wk1i*x0i
      x2i = wk1r*x0i - wk1i*x0r
      y9r = x1r - x2r
      y9i = x1i - x2i
      y13r = x1r + x2r
      y13i = x1i + x2i
      x0r = a(4) - a(21)
      x0i = a(5) + a(20)
      x1r = wk2r*x0r - wk2i*x0i
      x1i = wk2r*x0i + wk2i*x0r
      x0r = a(12) - a(29)
      x0i = a(13) + a(28)
      x2r = wk2i*x0r - wk2r*x0i
      x2i = wk2i*x0i + wk2r*x0r
      y2r = x1r + x2r
      y2i = x1i + x2i
      y6r = x1r - x2r
      y6i = x1i - x2i
      x0r = a(4) + a(21)
      x0i = a(5) - a(20)
      x1r = wk2i*x0r - wk2r*x0i
      x1i = wk2i*x0i + wk2r*x0r
      x0r = a(12) + a(29)
      x0i = a(13) - a(28)
      x2r = wk2r*x0r - wk2i*x0i
      x2i = wk2r*x0i + wk2i*x0r
      y10r = x1r - x2r
      y10i = x1i - x2i
      y14r = x1r + x2r
      y14i = x1i + x2i
      x0r = a(6) - a(23)
      x0i = a(7) + a(22)
      x1r = wk3r*x0r - wk3i*x0i
      x1i = wk3r*x0i + wk3i*x0r
      x0r = a(14) - a(31)
      x0i = a(15) + a(30)
      x2r = wk1i*x0r - wk1r*x0i
      x2i = wk1i*x0i + wk1r*x0r
      y3r = x1r + x2r
      y3i = x1i + x2i
      y7r = x1r - x2r
      y7i = x1i - x2i
      x0r = a(6) + a(23)
      x0i = a(7) - a(22)
      x1r = wk1i*x0r + wk1r*x0i
      x1i = wk1i*x0i - wk1r*x0r
      x0r = a(14) + a(31)
      x0i = a(15) - a(30)
      x2r = wk3i*x0r - wk3r*x0i
      x2i = wk3i*x0i + wk3r*x0r
      y11r = x1r + x2r
      y11i = x1i + x2i
      y15r = x1r - x2r
      y15i = x1i - x2i
      x1r = y0r + y2r
      x1i = y0i + y2i
      x2r = y1r + y3r
      x2i = y1i + y3i
      a(0) = x1r + x2r
      a(1) = x1i + x2i
      a(2) = x1r - x2r
      a(3) = x1i - x2i
      x1r = y0r - y2r
      x1i = y0i - y2i
      x2r = y1r - y3r
      x2i = y1i - y3i
      a(4) = x1r - x2i
      a(5) = x1i + x2r
      a(6) = x1r + x2i
      a(7) = x1i - x2r
      x1r = y4r - y6i
      x1i = y4i + y6r
      x0r = y5r - y7i
      x0i = y5i + y7r
      x2r = wn4r*(x0r - x0i)
      x2i = wn4r*(x0i + x0r)
      a(8) = x1r + x2r
      a(9) = x1i + x2i
      a(10) = x1r - x2r
      a(11) = x1i - x2i
      x1r = y4r + y6i
      x1i = y4i - y6r
      x0r = y5r + y7i
      x0i = y5i - y7r
      x2r = wn4r*(x0r - x0i)
      x2i = wn4r*(x0i + x0r)
      a(12) = x1r - x2i
      a(13) = x1i + x2r
      a(14) = x1r + x2i
      a(15) = x1i - x2r
      x1r = y8r + y10r
      x1i = y8i + y10i
      x2r = y9r - y11r
      x2i = y9i - y11i
      a(16) = x1r + x2r
      a(17) = x1i + x2i
      a(18) = x1r - x2r
      a(19) = x1i - x2i
      x1r = y8r - y10r
      x1i = y8i - y10i
      x2r = y9r + y11r
      x2i = y9i + y11i
      a(20) = x1r - x2i
      a(21) = x1i + x2r
      a(22) = x1r + x2i
      a(23) = x1i - x2r
      x1r = y12r - y14i
      x1i = y12i + y14r
      x0r = y13r + y15i
      x0i = y13i - y15r
      x2r = wn4r*(x0r - x0i)
      x2i = wn4r*(x0i + x0r)
      a(24) = x1r + x2r
      a(25) = x1i + x2i
      a(26) = x1r - x2r
      a(27) = x1i - x2i
      x1r = y12r + y14i
      x1i = y12i - y14r
      x0r = y13r - y15i
      x0i = y13i + y15r
      x2r = wn4r*(x0r - x0i)
      x2i = wn4r*(x0i + x0r)
      a(28) = x1r - x2i
      a(29) = x1i + x2r
      a(30) = x1r + x2i
      a(31) = x1i - x2r
   end
!
   subroutine cftf081(a, w)
      real*8 a(0:15), w(0:*)
      real*8 wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
      real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i
      real*8 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i
      wn4r = w(1)
      x0r = a(0) + a(8)
      x0i = a(1) + a(9)
      x1r = a(0) - a(8)
      x1i = a(1) - a(9)
      x2r = a(4) + a(12)
      x2i = a(5) + a(13)
      x3r = a(4) - a(12)
      x3i = a(5) - a(13)
      y0r = x0r + x2r
      y0i = x0i + x2i
      y2r = x0r - x2r
      y2i = x0i - x2i
      y1r = x1r - x3i
      y1i = x1i + x3r
      y3r = x1r + x3i
      y3i = x1i - x3r
      x0r = a(2) + a(10)
      x0i = a(3) + a(11)
      x1r = a(2) - a(10)
      x1i = a(3) - a(11)
      x2r = a(6) + a(14)
      x2i = a(7) + a(15)
      x3r = a(6) - a(14)
      x3i = a(7) - a(15)
      y4r = x0r + x2r
      y4i = x0i + x2i
      y6r = x0r - x2r
      y6i = x0i - x2i
      x0r = x1r - x3i
      x0i = x1i + x3r
      x2r = x1r + x3i
      x2i = x1i - x3r
      y5r = wn4r*(x0r - x0i)
      y5i = wn4r*(x0r + x0i)
      y7r = wn4r*(x2r - x2i)
      y7i = wn4r*(x2r + x2i)
      a(8) = y1r + y5r
      a(9) = y1i + y5i
      a(10) = y1r - y5r
      a(11) = y1i - y5i
      a(12) = y3r - y7i
      a(13) = y3i + y7r
      a(14) = y3r + y7i
      a(15) = y3i - y7r
      a(0) = y0r + y4r
      a(1) = y0i + y4i
      a(2) = y0r - y4r
      a(3) = y0i - y4i
      a(4) = y2r - y6i
      a(5) = y2i + y6r
      a(6) = y2r + y6i
      a(7) = y2i - y6r
   end
!
   subroutine cftf082(a, w)
      real*8 a(0:15), w(0:*)
      real*8 wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i
      real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i
      real*8 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i
      wn4r = w(1)
      wk1r = w(2)
      wk1i = w(3)
      y0r = a(0) - a(9)
      y0i = a(1) + a(8)
      y1r = a(0) + a(9)
      y1i = a(1) - a(8)
      x0r = a(4) - a(13)
      x0i = a(5) + a(12)
      y2r = wn4r*(x0r - x0i)
      y2i = wn4r*(x0i + x0r)
      x0r = a(4) + a(13)
      x0i = a(5) - a(12)
      y3r = wn4r*(x0r - x0i)
      y3i = wn4r*(x0i + x0r)
      x0r = a(2) - a(11)
      x0i = a(3) + a(10)
      y4r = wk1r*x0r - wk1i*x0i
      y4i = wk1r*x0i + wk1i*x0r
      x0r = a(2) + a(11)
      x0i = a(3) - a(10)
      y5r = wk1i*x0r - wk1r*x0i
      y5i = wk1i*x0i + wk1r*x0r
      x0r = a(6) - a(15)
      x0i = a(7) + a(14)
      y6r = wk1i*x0r - wk1r*x0i
      y6i = wk1i*x0i + wk1r*x0r
      x0r = a(6) + a(15)
      x0i = a(7) - a(14)
      y7r = wk1r*x0r - wk1i*x0i
      y7i = wk1r*x0i + wk1i*x0r
      x0r = y0r + y2r
      x0i = y0i + y2i
      x1r = y4r + y6r
      x1i = y4i + y6i
      a(0) = x0r + x1r
      a(1) = x0i + x1i
      a(2) = x0r - x1r
      a(3) = x0i - x1i
      x0r = y0r - y2r
      x0i = y0i - y2i
      x1r = y4r - y6r
      x1i = y4i - y6i
      a(4) = x0r - x1i
      a(5) = x0i + x1r
      a(6) = x0r + x1i
      a(7) = x0i - x1r
      x0r = y1r - y3i
      x0i = y1i + y3r
      x1r = y5r - y7r
      x1i = y5i - y7i
      a(8) = x0r + x1r
      a(9) = x0i + x1i
      a(10) = x0r - x1r
      a(11) = x0i - x1i
      x0r = y1r + y3i
      x0i = y1i - y3r
      x1r = y5r + y7r
      x1i = y5i + y7i
      a(12) = x0r - x1i
      a(13) = x0i + x1r
      a(14) = x0r + x1i
      a(15) = x0i - x1r
   end
!
   subroutine cftf040(a)
      real*8 a(0:7), x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
      x0r = a(0) + a(4)
      x0i = a(1) + a(5)
      x1r = a(0) - a(4)
      x1i = a(1) - a(5)
      x2r = a(2) + a(6)
      x2i = a(3) + a(7)
      x3r = a(2) - a(6)
      x3i = a(3) - a(7)
      a(0) = x0r + x2r
      a(1) = x0i + x2i
      a(2) = x1r - x3i
      a(3) = x1i + x3r
      a(4) = x0r - x2r
      a(5) = x0i - x2i
      a(6) = x1r + x3i
      a(7) = x1i - x3r
   end
!
   subroutine cftb040(a)
      real*8 a(0:7), x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i
      x0r = a(0) + a(4)
      x0i = a(1) + a(5)
      x1r = a(0) - a(4)
      x1i = a(1) - a(5)
      x2r = a(2) + a(6)
      x2i = a(3) + a(7)
      x3r = a(2) - a(6)
      x3i = a(3) - a(7)
      a(0) = x0r + x2r
      a(1) = x0i + x2i
      a(2) = x1r + x3i
      a(3) = x1i - x3r
      a(4) = x0r - x2r
      a(5) = x0i - x2i
      a(6) = x1r - x3i
      a(7) = x1i + x3r
   end
!
   subroutine cftx020(a)
      real*8 a(0:3), x0r, x0i
      x0r = a(0) - a(2)
      x0i = a(1) - a(3)
      a(0) = a(0) + a(2)
      a(1) = a(1) + a(3)
      a(2) = x0r
      a(3) = x0i
   end
!
   subroutine rftfsub(n, a, nc, c)
      integer n, nc, j, k, kk, ks, m
      real*8 a(0:n - 1), c(0:nc - 1), wkr, wki, xr, xi, yr, yi
      m = n/2
      ks = 2*nc/m
      kk = 0
      do j = 2, m - 2, 2
         k = n - j
         kk = kk + ks
         wkr = 0.5d0 - c(nc - kk)
         wki = c(kk)
         xr = a(j) - a(k)
         xi = a(j + 1) + a(k + 1)
         yr = wkr*xr - wki*xi
         yi = wkr*xi + wki*xr
         a(j) = a(j) - yr
         a(j + 1) = a(j + 1) - yi
         a(k) = a(k) + yr
         a(k + 1) = a(k + 1) - yi
      end do
   end
!
   subroutine rftbsub(n, a, nc, c)
      integer n, nc, j, k, kk, ks, m
      real*8 a(0:n - 1), c(0:nc - 1), wkr, wki, xr, xi, yr, yi
      m = n/2
      ks = 2*nc/m
      kk = 0
      do j = 2, m - 2, 2
         k = n - j
         kk = kk + ks
         wkr = 0.5d0 - c(nc - kk)
         wki = c(kk)
         xr = a(j) - a(k)
         xi = a(j + 1) + a(k + 1)
         yr = wkr*xr + wki*xi
         yi = wkr*xi - wki*xr
         a(j) = a(j) - yr
         a(j + 1) = a(j + 1) - yi
         a(k) = a(k) + yr
         a(k + 1) = a(k + 1) - yi
      end do
   end
!
   subroutine dctsub(n, a, nc, c)
      integer n, nc, j, k, kk, ks, m
      real*8 a(0:n - 1), c(0:nc - 1), wkr, wki, xr
      m = n/2
      ks = nc/n
      kk = 0
      do j = 1, m - 1
         k = n - j
         kk = kk + ks
         wkr = c(kk) - c(nc - kk)
         wki = c(kk) + c(nc - kk)
         xr = wki*a(j) - wkr*a(k)
         a(j) = wkr*a(j) + wki*a(k)
         a(k) = xr
      end do
      a(m) = c(0)*a(m)
   end
!
   subroutine dstsub(n, a, nc, c)
      integer n, nc, j, k, kk, ks, m
      real*8 a(0:n - 1), c(0:nc - 1), wkr, wki, xr
      m = n/2
      ks = nc/n
      kk = 0
      do j = 1, m - 1
         k = n - j
         kk = kk + ks
         wkr = c(kk) - c(nc - kk)
         wki = c(kk) + c(nc - kk)
         xr = wki*a(k) - wkr*a(j)
         a(k) = wkr*a(k) + wki*a(j)
         a(j) = xr
      end do
      a(m) = c(0)*a(m)
   end
!
end module OouraFFT
