      
      program fft1
      
      integer nmax
      real*8 pi
      parameter (nmax = 8192)
      integer n, z, v
      real*8 a(0 : nmax)
      pi = 4.0 * ATAN(1.0)

      write(*,*) 'length of data n = ?  '
      read(*,*) n
      
      
      write(*,*) 'input values are :  '
      do z = 0, n - 1
         a(z) = (z + 1)
         write(*,*) a(z)
      end do

      call cdft(n, cos(2 * pi / n), sin(2 * pi / n), a)
      call cdft(n, cos(2 * pi / n), -sin(2 * pi / n), a)
      

      write(*,*) 'output values are :  '
      do v = 0, n - 1
         write(*,*) a(v)
         
      end do
      
      return
      end


**********************************************************************
      subroutine bitrv(n, a)
      integer n, j, j1, k1, l, m, m2, n2
      real*8 a(0 : n-1), xr, xi
      m = n / 4
      m2 = 2 * m
      n2 = n - 2
      k = 0
      do j = 0, m2 - 4, 4
         if (j .lt. k) then
            xr = a(j)
            xi = a(j + 1)
            a(j) = a(k)
            a(j + 1) = a(k + 1)
            a(k) = xr
            a(k + 1) = xi
         else if (j .gt. k) then
            j1 = n2 - j
            k1 = n2 - k
            xr = a(j1)
            xi = a(j1 + 1)
            a(j1) = a(k1)
            a(j1 + 1) = a(k1 + 1)
            a(k1) = xr
            a(k1 + 1) = xi
         end if
         k1 = m2 + k
         xr = a(j + 2)
         xi = a(j + 3)
         a(j + 2) = a(k1)
         a(j + 3) = a(k1 + 1)
         a(k1) = xr
         a(k1 + 1) = xi
         l = m
         do while (k .ge. l)
            k = k - l
            l = l / 2
         end do
         k = k + 1
      end do
      end
**********************************************************************
      subroutine cdft(n, wr, wi, a)
      integer n, i, j, k, l, m
      real*8 wr, wi, a(0 : n -1), wmr, wmi, wkr, wki 
      real*8 wdr, wdi, ss, xr, xi
      wmr = wr
      wmi = wi
      m = n
      do while (m .gt. 4)
         l = m / 2
         wkr = 1
         wki = 0
         wdr = 1 - 2 * wmi * wmi
         wdi = 2 * wmi * wmr
         ss = 2 * wdi
         wmr = wdr
         wmi = wdi
         do j = 0, n -m, m
            i = j + l
            xr = a(j) - a(i)
            xi = a(j + 1) - a(i + 1)
            a(j) = a(j) + a(i)
            a(j + 1) = a(j + 1) + a(i + 1)
            a(i) = xr
            a(i + 1) = xi
            xr = a(j + 2) - a(i + 2)
            xi = a(j + 3) - a(i + 3)
            a(j + 2) = a(j + 2) + a(i + 2)
            a(j + 3) = a(j + 3) + a(i + 3)
            a(i + 2) = wdr * xr - wdi * xi
            a(i + 3) = wdr * xi + wdi * xr
         end do
         do k = 4, l - 4, 4
            wkr = wkr - ss * wdi
            wki = wki + ss * wdr
            wdr = wdr - ss * wki
            wdi = wdi + ss * wkr
            do j = k, n - m + k, m
               i = j + l
               xr = a(j) - a(i)
               xi = a(j + 1) - a(i + 1)
               a(j) = a(j) + a(i)
               a(j + 1) = a(j + 1) + a(i + 1)
               a(i) = wkr * xr - wki * xi
               a(i + 1) = wkr * xi + wki * xr
               xr = a(j + 2) - a(i + 2)
               xi = a(j + 3) - a(i + 3)
               a(j + 2) = a(j + 2) + a(i + 2)
               a(j + 3) = a(j + 3) + a(i + 3)
               a(i + 2) = wdr * xr - wdi *xi
               a(i + 3) = wdr * xi + wdi *xr
            end do
         end do
         m = l
      end do
      if(m .gt. 2) then
         do j = 0, n - 4, 4
            xr = a(j) - a(j + 2)
            xi = a(j + 1) - a(j + 3)
            a(j) = a(j + 1) + a(j + 3)
            a(j + 2) = xr
            a(j + 3) = xi
         end do
      end if
      if (n .gt. 4) call bitrv(n, a)
      end
      
