c
c     Numerical Analysis:
c     The Mathematics of Scientific Computing
c     D.R. Kincaid & E.W. Cheney
c     Brooks/Cole Publ., 1990
c
c     Section 6.12
c
c     Example of Fast Fourier Transform method
c
c
c     file: fft.f
c
      parameter (m=4,n3=15)
      complex z(0:n3),c(0:n3),d(0:n3),w,f,u,v
c
      print *
      print *,' Fast Fourier Transform example'
      print *,' Section 6.12, Kincaid-Cheney'
      print *
c
      nn=2**m
      m2=nn/2
      pi=4.0*atan(1.0)
      y=2.0*pi/float(nn)
      w=cmplx(0.0 , -y)
      w=cexp(w)
      z(0)=cmplx(1.0,0.0)
      c(0)=f(0.0)
      do 2 k=1,nn-1
         z(k)=w*z(k-1)
 	 c(k)=f(y*float(k))
 2    continue
      print *,' Initial c-vector'
      print 5,c
      print *
      do 4 n=0,m-1
         print *,' n = ',n
	 print 5,c
         print *
	 mn=2**(m-n-1)
	 n2=2**(n-1)
	    do 3 k=0,mn-1
	       do 3 j=0,n2
		  u=c(2*n2*k+j)
		  v=z(j*mn)*c(2*n2*k+j+m2)
		  d(4*n2*k+j)=(u+v)*cmplx(0.5,0.0)
 		  d(4*n2*k+j+2*n2)=(u-v)*cmplx(0.5,0.0)
 3             continue
       	    do 4 j=0,nn-1
 	       c(j)=d(j)
 4          continue
            print *,' Final c-vector'
	    print 5,c
c
 5    format (2x,16( '(',e13.6,',',e13.6,')',4x))
      stop
      end
c
      function f(x)
      complex f,z,d
      u=cos(x)
      v=sin(x)
      z=cmplx(u,v)
      d=cmplx(16.0,0.0)
      do 2 k=1,15
         y=float(k)
         d=d*z + cmplx(16.0-y,0.0)
 2    continue
      f=d
c
      return
      end