c---------------------------------------------------------------------
c compute the roots-of-unity array that will be used for subsequent FFTs. 
c---------------------------------------------------------------------
      subroutine CompExp (n, exponent)

      implicit none
      integer n
      double complex exponent(n) 
      integer ilog2
      external ilog2      
     
      integer m,nu,ku,i,j,ln
      double precision t, ti, pi 
      data pi /3.141592653589793238d0/

      nu = n
      m = ilog2(n)
      exponent(1) = m
      ku = 2
      ln = 1
      do j = 1, m
         t = pi / ln
         do i = 0, ln - 1
            ti = i * t
            exponent(i+ku) = dcmplx(cos(ti),sin(ti))
         enddo        
         ku = ku + ln
         ln = 2 * ln
      enddo
            
      return
      end

c---------------------------------------------------------------------
c---------------------------------------------------------------------
c---------------------------------------------------------------------
      integer function ilog2(n)
      implicit none
      integer n
c---------------------------------------------------------------------
c---------------------------------------------------------------------  
      integer nn, lg
      if (n .eq. 1) then
         ilog2=0
         return
      endif
      lg = 1
      nn = 2
      do while (nn .lt. n)
         nn = nn*2
         lg = lg+1
      end do
      ilog2 = lg
      return
      end
c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine ipow46(a, exponent, result)
c---------------------------------------------------------------------
c compute a^exponent mod 2^46
c---------------------------------------------------------------------

      implicit none
      double precision a, result, dummy, q, r
      integer exponent, n, n2
      external randlc
      double precision randlc
c---------------------------------------------------------------------
c Use
c   a^n = a^(n/2)*a^(n/2) if n even else
c   a^n = a*a^(n-1)       if n odd
c---------------------------------------------------------------------
      result = 1
      if (exponent .eq. 0) return
      q = a
      r = 1
      n = exponent

      do while (n .gt. 1)
         n2 = n/2
         if (n2 * 2 .eq. n) then
            dummy = randlc(q, q) 
            n = n2
         else
            dummy = randlc(r, q)
            n = n-1
         endif
      end do
      dummy = randlc(r, q)
      result = r
      return
      end
c---------------------------------------------------------------------
      subroutine CalculateChecksum(csum,iterN,u,d1,d2,d3)
        implicit none
        integer iterN
        integer d1,d2,d3
        double complex csum
        double complex u(d1+1,d2,d3)
        integer i, i1, ii, ji, ki
        csum = dcmplx (0.0, 0.0)
        do i = 1, 1024
          i1 = i
          ii = mod (i1, d1) + 1
          ji = mod (3 * i1, d2) + 1
          ki = mod (5 * i1, d3) + 1
          csum = csum + u(ii,ji,ki)
        end do
        csum = csum/dble(d1*d2*d3)
        write(*,30) iterN, csum
 30     format (' T =',I5,5X,'Checksum =',1P2D22.12)
      return
      end
c---------------------------------------------------------------------
c---------------------------------------------------------------------
      subroutine compute_initial_conditions(u0,d1,d2,d3)

      implicit none
      include 'npbparams.h'
      integer d1,d2,d3
      double complex u0(d1+1,d2,d3), tmp(maxdim)
      double precision x0, start, an, dummy
      double precision RanStarts(maxdim)

      integer i,j,k
      double precision seed, a
      parameter (seed = 314159265.d0, a = 1220703125.d0)
      external randlc
      double precision randlc
      
      start = seed                                    
c---------------------------------------------------------------------
c Jump to the starting element for our first plane.
c---------------------------------------------------------------------
      call ipow46(a, 0, an)
      dummy = randlc(start, an)
      call ipow46(a, 2*d1*d2, an)
c---------------------------------------------------------------------
c Go through by z planes filling in one square at a time.
c---------------------------------------------------------------------
      RanStarts(1) = start
      do k = 2, d3 
         dummy = randlc(start, an)
         RanStarts(k) = start
      end do
      
      do k = 1, d3 
         x0 = RanStarts(k)
         do j = 1, d2 
           call vranlc(2*d1, x0, a, tmp)
           do i = 1, d1 
             u0(i,j,k)=tmp(i)
           end do
         end do
      end do

      return
      end
c---------------------------------------------------------------------
c---------------------------------------------------------------------
      subroutine evolve(x,y,twiddle,nx,ny,nz)
      implicit none
      integer nx,ny,nz
      double complex x(nx+1,ny,nz),y(nx+1,ny,nz)
      real*8 twiddle(nx+1,ny,nz)
      integer i,j,k
           do i = 1, nz
             do k = 1, ny
               do j = 1, nx
                   y(j,k,i)=y(j,k,i)*twiddle(j,k,i)
                   x(j,k,i)=y(j,k,i)
                 end do
              end do
           end do
      
      return
      end
c---------------------------------------------------------------------
c---------------------------------------------------------------------