common /time/ f
      complex f(0:10000),ai,q,g(0:8),c,c2
      pi=acos(-1.0)
      ai = (0,1.0)
      n=8
      delta=0.4
      q = cos(2.0 * pi/n) + ai*sin(2.0 * pi/n)
      w = (2.0 * pi/n)/delta

      do j=1,8
       t = (j-4) * delta
       th1 = j * pi * (7.0/8.0)
       c = cos(th1) + ai * sin(th1)
       f(j) = c * (exp(-4.0 * t**2.0) + 0.1 * cos(5.0 * t))
      enddo
      call eight(q,g,1,0)
      do j=1,7
       th2 = j * pi * (7.0/8.0)
       c2 = cos(th2) + ai*sin(th2)
       g(j) = g(j) * c2
       print *,j,w*(j-4),abs(g(j))
      enddo

      stop
      end


      subroutine four (t,x,k,i)
       common /time/ f
       complex t,f(0:10000),x(0:3),s1,s2,s3,s4
       s1=f(i)+f(i+2*k)
       s2=f(i+k)+f(i+3*k)
       s3=f(i)-f(i+2*k)
       s4=f(i+k)-f(i+3*k)
       x(0)=s1+s2
       x(1)=s3+t*s4
       x(2)=s1-s2
       x(3)=s3-t*s4
       return
      end

      subroutine eight (q,g,k,i)
       common /time/ f
       complex q,g(0:7),x(0:3),y(0:3),t,f(0:10000)
       call four(q**2,x,2*k,i)
       call four(q**2,y,2*k,i+k)
       t=1.0
       do j=0,3
        g(j)=x(j)+t*y(j)
        g(j+4)=x(j)-t*y(j)
        t=t*q
       enddo
       return
      end

      subroutine sixteen(q,g,k,i)
       common /time/ f
       complex q,g(0:15),x(0:7),y(0:7),t,f(0:10000)
       call eight(q**2,x,2*k,i)
       call eight(q**2,y,2*k,i+k)
       t = 1.0
       do j=0,7
        g(j)=x(j)+t*y(j)
        g(j+8)=x(j)-t*y(j)
        t=t*q
       enddo
       return
      end