program main
    implicit none
    integer,parameter :: FTSIZE=8
    real*8,parameter :: M_PI=3.141592653589793238460D0
    integer :: i
    complex*16,dimension(0:FTSIZE-1):: X,B,C
    do i=0,FTSIZE-1
        X(i)=cmplx(1D0/(i+1),1D0/(FTSIZE-i),kind(X(i)));
    end do
    print '("N=",I0)',FTSIZE
    call fft(X,B)
    print '("fft_B=")'
    call cvecprint(B)
    return

contains 

recursive subroutine fftwork(N,s,o,x,p,b)
    integer :: k,N2,N,s,o,p
    complex*16,dimension(0:) :: x,b
    complex*16 :: even,odd,w
    if(N.eq.1) then
        b(p)=x(o)
        return
    end if
    if(mod(N,2).ne.0) then
        print '("N not divisible by 2!")'
        stop
    end if
    N2=N/2
    call fftwork(N2,2*s,o,x,p,b)
    call fftwork(N2,2*s,o+s,x,p+N2,b)
    do k=0,N2-1
        even=b(p+k)
        odd=b(p+k+N2)
        w=exp(cmplx(0.0D0,-2*M_PI*k/N,kind(w)))
        b(p+k)=even+w*odd
        b(p+k+N2)=even-w*odd
    end do
end subroutine

subroutine fft(x,b)
    complex*16,dimension(0:) :: x,b
    integer :: N
    N=size(x)
    call fftwork(N,1,0,x,0,b)
end subroutine

subroutine cvecprint(x)
    complex*16,dimension(0:) :: x
    integer :: N,i
    N=size(x)
    do i=0,N-1
        print '("(",G0," ",G0,")")',x(i)
    end do
end subroutine

end program
