! File: scaled2h.f90
! Public domain 2004 James Van Buskirk

! PURPOSE: Generates full complex FFT algorithms of order 2**n.
! Prompts the user for the exponent of the problem, and if the
! user wants to generate FFT code for problems of order 2**n, he
! enters n at the prompt.  Two files are then created:
! 1) fft###.i90 contains a subroutine fft###(h) that converts a
!    complex vector h of order ### = 2**n into its DFT:
!    h(0:2**n-1) = &
!    (/(sum((/(exp(-2*pi*(0,1)*i*j/2**n),j=0,2**n-1)/)),i=0,2**n-1)/)
! 2) fft###.f90 contains a test program that shows how one has
!    to INCLUDE the above file in a module that defines the
!    parameter wp to be a kind type parameter suitable for real
!    types, and how to call fft### to convert a vector into its
!    DFT.
! If everything works OK, the fft###.f90 file should compile
! to a program that prints out 2*2**n lines of output that are
! near EPSILON(1.0d0) and a final line of output that is the
! maximum of the previous 2*2**n lines.

! ***************************************************************
! module mykinds
! Sets up kind type parameter for double precision and
! (if available) quadruple precision.  If there is no
! quad precision type, double precision will be used.
! ***************************************************************
module mykinds
   implicit none
   integer, parameter :: dp = selected_real_kind(15,300)
   integer, parameter :: qp_preferred = selected_real_kind(30,1000)
   integer, parameter :: qp = (1+sign(1,qp_preferred))/2*qp_preferred+ &
                              (1-sign(1,qp_preferred))/2*dp
end module mykinds
! ***************************************************************

! ***************************************************************
! module avl_mod
! Contains tools and derived types to store randomly generated
! data in an AVL tree and convert the tree to an array.
! NOTE ABOUT RE-USE: This is not a general-purpose AVL tree
! module: it was designed to work on only positive data, and
! there is no deletion procedure nor would it be simple to
! implement one for this module without adding an extra data
! structure to handle the indices.
! ***************************************************************
module avl_mod
   use mykinds
   implicit none

!................................................................
! type avl_node
! A node in the AVL tree like any avl node, but also has
! an index that tells us where the key will end up in
! the eventual array.
!................................................................
   type avl_node
      real(qp) key                               ! Actual (positive) data
      integer index                              ! Index into ultimate array
      integer :: balance = 0                     ! -1 = heavy left, +1 = heavey right
      type(avl_node), pointer :: left => NULL()  ! Pointer to left child
      type(avl_node), pointer :: right => NULL() ! Pointer to right child
   end type avl_node
!................................................................

!................................................................
! type avl_bag
! Contains an AVL tree root and a counter of its number of nodes
!................................................................
   type avl_bag
      type(avl_node), pointer :: root => NULL()   ! Root of AVL tree
      integer :: num_nodes = 0                    ! Node counter
   end type avl_bag
!................................................................

!................................................................
! Interface operator(<)
! Allows us to fuzzy compare a real(qp) with an avl_node
!................................................................
   interface operator(<)
      module procedure lessn
   end interface operator(<)
!................................................................

!................................................................
! Interface operator(>)
! Allows us to fuzzy compare a real(qp) with an avl_node
!................................................................
   interface operator(>)
      module procedure greatern
   end interface operator(>)
!................................................................

! Fuzzy compares above will be up to 100 ulps
   integer, parameter :: ulps = 100

! Module entities private by default
   private
! Need the below
   public avl_node
! Need the above
! Only public entities:
   public avl_bag, key2index, bag2array

   contains

!................................................................
! function lessn(key,root)
! Performs fuzzy compare of real(qp) key with type(avl_tree) root
!................................................................
      pure function lessn(key,root)
         real(qp), intent(in) :: key        ! Left operand of comparison
         type(avl_node), intent(in) :: root ! Right operand of comparison
         logical lessn                      ! Result of comparison

         lessn = abs(key) < abs(root%key)-ulps*spacing(root%key)
      end function lessn
!................................................................

!................................................................
! function greatern(key,root)
! Performs fuzzy compare of real(qp) key with type(avl_tree) root
!................................................................
      pure function greatern(key,root)
         real(qp), intent(in) :: key        ! Left operand of comparison
         type(avl_node), intent(in) :: root ! Right operand of comparison
         logical greatern                   ! Result of comparison

         greatern = abs(key) > abs(root%key)+ulps*spacing(root%key)
      end function greatern
!................................................................

!................................................................
! function key2index(key,bag,index)
! Stores key as avl_node in bag and returns index of key
! in eventual array.  Does not store duplicate copies.
!................................................................
      subroutine key2index(key,bag,index)
         real(qp), intent(in) :: key         ! Key to be compared and stored
         type(avl_bag), intent(inout) :: bag ! Where the key will be stored
         integer, intent(out) :: index       ! Index of incoming key
         logical changed                     ! Required for avl_insert; ignored

         index = bag%num_nodes
         if(associated(bag%root)) then  ! We have an AVL tree to search
            call avl_insert(key,bag%root,index,changed)
            bag%num_nodes = max(bag%num_nodes, index)
         else                           ! Start new AVL tree
            bag%num_nodes = 1
            allocate(bag%root)
            bag%root%key = key
            index = 1
            bag%root%index = index
         end if
      end subroutine key2index
!................................................................

!................................................................
! subroutine avl_insert(key,root,index,changed)
! Searched AVL tree starting at root for key.  If found, returns
! index of matching key and clears changed.  If not found, creates
! new node containing key and next index.  Sets change accordingly.
!................................................................
      recursive subroutine avl_insert(key,root,index,changed)
         real(qp), intent(in) :: key     ! Key to search for in AVL tree
         type(avl_node), pointer :: root ! Root of AVL tree (must be associated)
         integer, intent(inout) :: index ! On input, last valid index
                                         ! On output, index of key
         logical, intent(out) :: changed ! Set if tree grew longer

         if(key < root) then
            if(associated(root%left)) then
               call avl_insert(key,root%left,index,changed)
! After insertion, if tree grew longer we must recalculate
! balance factor and rotate to maintain AVL condition if
! necessary.
               if(changed) then
                  if(root%balance == -1) then
                     if(root%left%balance == 1) then
                        call rot2right(root)
                     else ! root%left%balance <= 0
                        call rot1right(root)
                     end if
                     changed = .FALSE.
                  else if(root%balance == 0) then
                     root%balance = -1
                  else ! root%balance == 1
                     root%balance = 0
                     changed = .FALSE.
                  end if
               end if
            else ! .NOT.associated(root%left)
! If there is an empty subtree to search, we know the key
! goes here.
               index = index+1
               allocate(root%left)
               root%left%key = key
               root%left%index = index
               changed = root%balance == 0
               root%balance = root%balance-1
            end if
         else if(key > root) then
            if(associated(root%right)) then
               call avl_insert(key,root%right,index,changed)
! After insertion, if tree grew longer we must recalculate
! balance factor and rotate to maintain AVL condition if
! necessary.
               if(changed) then
                  if(root%balance == -1) then
                     root%balance = 0
                     changed = .FALSE.
                  else if(root%balance == 0) then
                     root%balance = 1
                  else ! root%balance == 1
                     if(root%right%balance == -1) then
                        call rot2left(root)
                     else ! root%right%balance <= 0
                        call rot1left(root)
                     end if
                     changed = .FALSE.
                  end if
               end if
            else ! .NOT.associated(root%right)
! If there is an empty subtree to search, we know the key
! goes here.
               index = index+1
               allocate(root%right)
               root%right%key = key
               root%right%index = index
               changed = root%balance == 0
               root%balance = root%balance+1
            end if
         else ! key == root
! If we found the key, return the matching index.
            changed = .FALSE.
            index = root%index
         end if
      end subroutine avl_insert
!................................................................

!................................................................
! Subroutine rot2right(root)
! Performs a double right rotation of the AVL tree at root
!................................................................
      subroutine rot2right(root)
         type(avl_node), pointer :: root        ! Root of AVL tree
         type(avl_node), pointer :: rotate_temp ! Temporary pointer to new root

! Create new connections
         rotate_temp => root%left%right
         root%left%right => rotate_temp%left
         rotate_temp%left => root%left
         root%left => rotate_temp%right
         rotate_temp%right => root
         root => rotate_temp
! Recalculate balance factors
         root%right%balance = max(-root%balance, 0)
         root%left%balance = -(root%balance+root%right%balance)
         root%balance = 0
      end subroutine rot2right
!................................................................

!................................................................
! Subroutine rot2left(root)
! Performs a double left rotation of the AVL tree at root
!................................................................
      subroutine rot2left(root)
         type(avl_node), pointer :: root        ! Root of AVL tree
         type(avl_node), pointer :: rotate_temp ! Temporary pointer to new root

! Create new connections
         rotate_temp => root%right%left
         root%right%left => rotate_temp%right
         rotate_temp%right => root%right
         root%right => rotate_temp%left
         rotate_temp%left => root
         root => rotate_temp
! Recalculate balance factors
         root%left%balance = min(-root%balance, 0)
         root%right%balance = -(root%balance+root%left%balance)
         root%balance = 0
      end subroutine rot2left
!................................................................

!................................................................
! Subroutine rot1right(root)
! Performs a single right rotation of the AVL tree at root
!................................................................
      subroutine rot1right(root)
         type(avl_node), pointer :: root        ! Root of AVL tree
         type(avl_node), pointer :: rotate_temp ! Temporary pointer to new root

! Create new connections
         rotate_temp => root%left
         root%left => rotate_temp%right
         rotate_temp%right => root
         root => rotate_temp
! Recalculate balance factors
         root%balance = root%balance+1
         root%right%balance = -root%balance
      end subroutine rot1right
!................................................................

!................................................................
! Subroutine rot1left(root)
! Performs a single left rotation of the AVL tree at root
!................................................................
      subroutine rot1left(root)
         type(avl_node), pointer :: root        ! Root of AVL tree
         type(avl_node), pointer :: rotate_temp ! Temporary pointer to new root

! Create new connections
         rotate_temp => root%right
         root%right => rotate_temp%left
         rotate_temp%left => root
         root => rotate_temp
! Recalculate balance factors
         root%balance = root%balance-1
         root%left%balance = -root%balance
      end subroutine rot1left
!................................................................

!................................................................
! subroutine bag2array(bag,array)
! Destroys AVL tree rooted in bag and stores the keys at
! the appropriate indices in array
!................................................................
      subroutine bag2array(bag, array)
         type(avl_bag), intent(inout) :: bag           ! Has tree to be converted
         real(qp), intent(out) :: array(bag%num_nodes) ! Resultant array of keys

! Set up call to first level of avl_detroy
         if(associated(bag%root)) then
            call avl_destroy(bag%root,array)
            bag%num_nodes = 0
         end if
      end subroutine bag2array
!................................................................

!................................................................
! subroutine avl_destroy(root,array)
! Destroys AVL tree rooted at root and stores keys at the
! appropriate indices in array.
!................................................................
      recursive subroutine avl_destroy(root,array)
         type(avl_node), pointer :: root     ! Root of AVL tree to destroy
         real(qp), intent(inout) :: array(*) ! Array to store contents of tree

! Set up call to next level on left
         if(associated(root%left)) then
            call avl_destroy(root%left,array)
         end if
! Set up call to next level on right
         if(associated(root%right)) then
            call avl_destroy(root%right,array)
         end if
! Store current key in array in destroy current node
         array(root%index) = root%key
         deallocate(root)
      end subroutine avl_destroy
!................................................................

end module avl_mod
! ***************************************************************

! ***************************************************************
! module code_gen
! Contains tools and derived types to generate full complex DFT code
! ***************************************************************
module code_gen
   use mykinds
   use avl_mod
   implicit none

!................................................................
! type get_factor_stuff
! Contains a pointer to a quad-precision array
!................................................................
   type get_factor_stuff
      real(qp), pointer :: factor(:) ! Pointer to quad array for current order
   end type get_factor_stuff
!................................................................

!................................................................
! type odd_tail_stuff
! Contains pointers to 2-d arrays for sine and cosine array entries
! Entries will contain the indices of the constants
!................................................................
   type odd_tail_stuff
      integer, pointer :: matrixc(:,:) ! Matrix of cosine coefficient indices
      integer, pointer :: matrixs(:,:) ! Matrix of sine coefficient indices
   end type odd_tail_stuff
!................................................................

!................................................................
! type scaled_odd_tail_stuff
! Contains pointers to arrays of tangents, sine and cosine entries,
! and the last diagonal.  Entries will contain the indices of the constants
!................................................................
   type scaled_odd_tail_stuff
      integer, pointer :: matrixc(:)    ! Array of cosine tangent indices
      integer, pointer :: matrixs(:)    ! Array of sine tangent indices
      integer, pointer :: matrix1c(:,:) ! Matrix of cosine coefficient indices
      integer, pointer :: matrix1s(:,:) ! Matrix of sin coefficient indices
      integer, pointer :: diag(:)       ! Array of last diagonal indices
   end type scaled_odd_tail_stuff
!................................................................

!................................................................
! type correct_stuff
! Contains data to compute odd_tail(8).  Entries
! will contain the indices of the constants.
!................................................................
   type correct_stuff
      integer rsq2  ! Index of 1/sqrt(2)
   end type correct_stuff
!................................................................

!................................................................
! type var_stuff
! Contains latitude, longtude, and label of current variable.
! Note LF95 pads at the end to fit into arrays
!................................................................
   type var_stuff
      integer depth      ! Number of LOC executed to reach variable
      integer index      ! Index of variable among current set
      character(2) label ! xr = real even, yr = real odd
                         ! xi = imaginary even, yi = imaginary odd
                         ! 'y' labels not yet implemented
   end type var_stuff
!................................................................

   contains
!................................................................
! Subroutine get_factor(factor,n,old_factor)
! Computes scale factor thrown out by scaled_odd_tail(n)
!................................................................
      recursive subroutine get_factor(factor, n, old_factor)
         integer, intent(in) :: n                     ! Order of the problem
         real(qp), intent(out) :: factor(2**(n-4))    ! Scale factor
         real(qp), intent(in) :: old_factor(2**(n-6)) ! Previous scale factor
         integer i                                    ! Multiple of fundamental angle
         integer k                                    ! Current block in factor
         real(qp) pi

         pi = 4*atan(1.0_qp)
         if(n <= 3) then
            factor = 1
         else
            factor = (/(cos(2*pi*i/2**n),i=1,2**(n-3),2)/)
            if(n > 5) then
               k = 2**(n-6)
               factor(0*k+1:0*k+k) = factor(0*k+1:0*k+k)*old_factor
               factor(1*k+1:1*k+k) = factor(1*k+1:1*k+k)*old_factor(k:1:-1)
               factor(2*k+1:2*k+k) = factor(2*k+1:2*k+k)*old_factor
               factor(3*k+1:3*k+k) = factor(3*k+1:3*k+k)*old_factor(k:1:-1)
            end if
         end if
      end subroutine get_factor
!................................................................

!................................................................
! Subroutine prepare_correct(n,body,correct_data)
! Prepares constants for odd_tail algorithms leading to degenerate
! recursion
!................................................................
      subroutine prepare_correct(n, body, correct_data)
         type(avl_bag), intent(inout) :: body ! Holds any keys created
         integer, intent(in) :: n             ! Order of overall problem
         type(correct_stuff) correct_data     ! Holds indices of keys created
         real(qp) constant                    ! Current constant created
         real(qp) pi                          ! Pi

         if(n > 2) then
            pi = 4*atan(1.0_qp)
            constant = cos(pi/4)
            call key2index(constant, body, correct_data%rsq2)
         end if
      end subroutine prepare_correct
!................................................................

!................................................................
! Subroutine odd_tail_correct(h,n,correct_data,iunit)
! Computes odd_tail(n)*h using the theoretical formula for
! orders too degenerate for recursion and immotalizes the
! result in code
!................................................................
      subroutine odd_tail_correct(h,n,correct_data,iunit)
         integer, intent(in) :: n                      ! Order of the problem
         type(var_stuff), intent(inout) :: h(2**(n-1)) ! Array of variable mnemonics
         type(correct_stuff) correct_data              ! Data for degenerate problem
         integer, intent(in) :: iunit                  ! Write code to this unit
         type(var_stuff) h1(2**(n-1)), h2(2**(n-1))    ! Temporary variable mnemonics
         integer max_depth                             ! Previous maximum depth
         integer i                                     ! General index

         select case(n)
            case(1)
            case(2)
            case(3)
               max_depth = maxval(h%depth)
               h1(1) = var_stuff(max_depth+1, h(1)%index, h(1)%label)
               call WX('   '//write_var(h1(1))//' = '//write_const(correct_data%rsq2)// &
                              '*('//write_var(h(2))//'-'//write_var(h(4))//')',iunit)
!               xc = const_array(correct_data%rsq2)*(h(2)-h(4))
               h1(2) = var_stuff(max_depth+1, h(2)%index, h(2)%label)
               call WX('   '//write_var(h1(2))//' = '//write_const(correct_data%rsq2)// &
                              '*('//write_var(h(2))//'+'//write_var(h(4))//')',iunit)
!               xs = const_array(correct_data%rsq2)*(h(2)+h(4))
               do i = 1, 2**(n-1)
                  h2(i) = var_stuff(max_depth+2,h(i)%index,h(i)%label)
               end do
               call WX('   '//write_var(h2(1))//' = '// &
                              write_var(h(1))//'+'//write_var(h1(1)),iunit)
               call WX('   '//write_var(h2(2))//' = '// &
                              write_var(h(1))//'-'//write_var(h1(1)),iunit)
               call WX('   '//write_var(h2(3))//' = '// &
                              write_var(h1(2))//'+'//write_var(h(3)),iunit)
               call WX('   '//write_var(h2(4))//' = '// &
                              write_var(h1(2))//'-'//write_var(h(3)),iunit)
!               h = (/h(1)+xc,h(1)-xc,xs+h(3),xs-h(3)/)
               do i = 1, 2**(n-1)
                  h(i) = var_stuff(max_depth+2,h2(i)%index,h2(i)%label)
               end do
            case default
               stop ' bigger n than 4 in odd_tail_correct'
         end select
      end subroutine odd_tail_correct
!................................................................

!................................................................
! Subroutine odd_tail_gen(stuff,n,body,get_factor_data)
! Compute odd_tail constants and places their indices in arrays
!................................................................
      subroutine odd_tail_gen(stuff,n,body,get_factor_data)
         integer, intent(in) :: n                   ! Order of the problem
         type(avl_bag) body                         ! AVL tree to hold constants
         type(odd_tail_stuff), intent(out) :: stuff ! Indices of data for
                                                    ! current order
         type(get_factor_stuff), intent(in) :: get_factor_data(n-2) ! Factors
                                                    ! to absorb
         integer i                                  ! Current row
         integer j                                  ! Column select
         integer k                                  ! Current blocksize
         real(qp) pi                                ! Pi
         real(qp) matrix(2,2**(n-3))                ! Array of actual data

         pi = 4*atan(1.0_qp)
         allocate(stuff%matrixc(2,2**(n-3)))
         allocate(stuff%matrixs(2,2**(n-3)))
         if(n >= 4) then
            k = 2**(n-4)
            do i = 1, k
               matrix(1,i) = cos(2*pi*(2*i-1)/2**n)
               matrix(1,2*k+1-i) = sin(2*pi*(2*i-1)/2**n)
               matrix(2,i) = -sin(2*pi*(2*i-1)/2**n)
               matrix(2,2*k+1-i) = cos(2*pi*(2*i-1)/2**n)
            end do
            k = 2**(n-6)
            do i = 1, k
               matrix(:,0*k+i) = matrix(:,0*k+i)* &
                                get_factor_data(n-2)%factor(i)
               matrix(:,8*k+1-i) = matrix(:,8*k+1-i)* &
                                get_factor_data(n-2)%factor(i)
               matrix(:,1*k+i) = matrix(:,1*k+i)* &
                                get_factor_data(n-2)%factor(k+1-i)
               matrix(:,7*k+1-i) = matrix(:,7*k+1-i)* &
                                get_factor_data(n-2)%factor(k+1-i)
               matrix(:,2*k+i) = matrix(:,2*k+i)* &
                                get_factor_data(n-2)%factor(i)
               matrix(:,6*k+1-i) = matrix(:,6*k+1-i)* &
                                get_factor_data(n-2)%factor(i)
               matrix(:,3*k+i) = matrix(:,3*k+i)* &
                                get_factor_data(n-2)%factor(k+1-i)
               matrix(:,5*k+1-i) = matrix(:,5*k+1-i)* &
                                get_factor_data(n-2)%factor(k+1-i)
            end do
            do i = 1, 2**(n-3)
               do j = 1, 2
                  call key2index(abs(matrix(j,i)), body, stuff%matrixc(j,i))
               end do
            end do
            k = 2**(n-4)
            do i = 1, k
               matrix(1,i) = sin(2*pi*(2*i-1)/2**n)
               matrix(1,2*k+1-i) = cos(2*pi*(2*i-1)/2**n)
               matrix(2,i) = cos(2*pi*(2*i-1)/2**n)
               matrix(2,2*k+1-i) = -sin(2*pi*(2*i-1)/2**n)
            end do
            k = 2**(n-6)
            do i = 1, k
               matrix(:,0*k+i) = matrix(:,0*k+i)* &
                                get_factor_data(n-2)%factor(i)
               matrix(:,8*k+1-i) = matrix(:,8*k+1-i)* &
                                get_factor_data(n-2)%factor(i)
               matrix(:,1*k+i) = matrix(:,1*k+i)* &
                                get_factor_data(n-2)%factor(k+1-i)
               matrix(:,7*k+1-i) = matrix(:,7*k+1-i)* &
                                get_factor_data(n-2)%factor(k+1-i)
               matrix(:,2*k+i) = matrix(:,2*k+i)* &
                                get_factor_data(n-2)%factor(i)
               matrix(:,6*k+1-i) = matrix(:,6*k+1-i)* &
                                get_factor_data(n-2)%factor(i)
               matrix(:,3*k+i) = matrix(:,3*k+i)* &
                                get_factor_data(n-2)%factor(k+1-i)
               matrix(:,5*k+1-i) = matrix(:,5*k+1-i)* &
                                get_factor_data(n-2)%factor(k+1-i)
            end do
            do i = 1, 2**(n-3)
               do j = 1, 2
                  call key2index(abs(matrix(j,i)), body, stuff%matrixs(j,i))
               end do
            end do
         end if
      end subroutine odd_tail_gen
!................................................................

!................................................................
! Subroutine odd_tail(h,n,correct_data, &
!      odd_tail_data,scaled_odd_tail_data,iunit)
! Compute action of odd tail matrix of order 2**(n-1)
! and immortalize in code
!................................................................
      recursive subroutine odd_tail(h,n,correct_data, &
                     odd_tail_data,scaled_odd_tail_data,iunit)
         integer, intent(in) :: n                       ! Order of the problem
         type(var_stuff), intent(inout) :: h(2**(n-1))  ! Array of variable mnemonics
         type(correct_stuff) correct_data               ! Indices of data for
                                                        ! degenerate orders
         type(odd_tail_stuff), intent(in) :: odd_tail_data(n) ! Indices of data
                                                        ! For the problem
         type(scaled_odd_tail_stuff) scaled_odd_tail_data(n-2) ! Indices of data for
                                                        ! scaled version of problem
         integer, intent(in) :: iunit                   ! Write code to this unit
         type(var_stuff) hc(2**(n-3)), hc1(2**(n-3))    ! Temporary varaible mnemonics
         type(var_stuff) hs(2**(n-3)), hs1(2**(n-3))    ! Temporary varaible mnemonics
         type(var_stuff) h2(2**(n-2))                   ! Temporary varaible mnemonics
         integer i                                      ! Row number
         integer k                                      ! Blocksize
         integer M                                      ! Blocksize
         integer max_depth                              ! Previous maximum depth

         if(n < 4) then
            call odd_tail_correct(h,n,correct_data,iunit)
         else
            max_depth = maxval(h%depth)
            M = 2**(n-3)
            do i = 1, M
               h2(i) = h(2*i-1)
               hc(i) = var_stuff(max_depth+1,h(2*i)%index,h(2*i)%label)
               hs(i) = var_stuff(max_depth+1,h(2*(i+M))%index,h(2*(i+M))%label)
               h2(i+M) = h(2*(i+M)-1)
            end do
            do i = 1, M
               call WX('   '//write_var(hc(i))//' = '// &
                  write_var(h(4*i-2))//'-'//write_var(h(4*(M+1-i))),iunit)
            end do
!            hc = h(2::4)-h(2**(n-1):1:-4)
            call scaled_odd_tail(hc,n-2,correct_data,scaled_odd_tail_data,iunit)
            max_depth = maxval(hc%depth)
            do i = 1, M
               hc1(i) = var_stuff(max_depth+1,hc(i)%index,hc(i)%label)
            end do
            M = 2**(n-4)
            do i = 1, M
               call WX('   '//write_var(hc1(i))//' = '// &
                  write_const(odd_tail_data(n)%matrixc(1,i))//'*'// &
                  write_var(hc(i))//'-'//write_const(odd_tail_data(n)%matrixc(2,i)) &
                  //'*'//write_var(hc(M+i)),iunit)
            end do
            do i = M+1,2*M
               call WX('   '//write_var(hc1(i))//' = '// &
                  write_const(odd_tail_data(n)%matrixc(1,i))//'*'// &
                  write_var(hc(2*M+1-i))//'+'// &
                  write_const(odd_tail_data(n)%matrixc(2,i))//'*'// &
                  write_var(hc(3*M+1-i)),iunit)
            end do
!            hc = (/(/(odd_tail_data(n)%matrixc(1,i)*hc(i)- &
!                      odd_tal_data(n)%matrixc(2,i)*hc(M+i),i=1,M)/), &
!                   (/(odd_tail_data(n)%matrixc(1,i)*hc(2*M+1-i)+ &
!                      odd_tail_data(n)%matrixc(2,i)*hc(3*M+1-i),i=M+1,2*M)/)/)
            M = 2**(n-3)
            do i = 1, M
               call WX('   '//write_var(hs(i))//' = '// &
                  write_var(h(4*i-2))//'+'//write_var(h(4*(M+1-i))),iunit)
            end do
!            hs = h(2::4)+h(2**(n-1):1:-4)
            call scaled_odd_tail(hs,n-2,correct_data,scaled_odd_tail_data,iunit)
            max_depth = maxval(hs%depth)
            do i = 1, M
               hs1(i) = var_stuff(max_depth+1,hs(i)%index,hs(i)%label)
            end do
            M = 2**(n-4)
            do i = 1, M
               call WX('   '//write_var(hs1(i))//' = '// &
                  write_const(odd_tail_data(n)%matrixs(1,i))//'*'// &
                  write_var(hs(i))//'+'// &
                  write_const(odd_tail_data(n)%matrixs(2,i))//'*'// &
                  write_var(hs(M+i)),iunit)
            end do
            do i = M+1,2*M
               call WX('   '//write_var(hs1(i))//' = '// &
                  write_const(odd_tail_data(n)%matrixs(1,i))//'*'// &
                  write_var(hs(2*M+1-i))//'-'// &
                  write_const(odd_tail_data(n)%matrixs(2,i))//'*'// &
                  write_var(hs(3*M+1-i)),iunit)
            end do
!            hs = (/(/(odd_tail_data(n)%matrixs(1,i)*hs(i)+ &
!                      odd_tail_data(n)%matrixs(2,i)*hs(M+i),i=1,M)/), &
!                   (/(odd_tail_data(n)%matrixs(1,i)*hs(2*M+1-i)- &
!                      odd_tail_data(n)%matrixs(2,i)*hs(3*M+1-i),i=M+1,2*M)/)/)
!            h2 = h(::2)
            call odd_tail(h2,n-1,correct_data,odd_tail_data,scaled_odd_tail_data,iunit)
            k = 2**(n-3)
            max_depth = max(maxval(hc1%depth),maxval(hs1%depth),maxval(h2%depth))
            do i = 1, k
               h(2*i-1) = var_stuff(max_depth+1,h2(i)%index,h2(i)%label)
               h(2*i) = var_stuff(max_depth+1,hc1(i)%index,hc1(i)%label)
               h(2*(i+k)) = var_stuff(max_depth+1,hs1(i)%index,hs1(i)%label)
               h(2*(i+k)-1) = var_stuff(max_depth+1,h2(i+k)%index,h2(i+k)%label)
            end do
            do i = 1, k
               call WX('   '//write_var(h(i))//' = '// &
                  write_var(h2(i))//'+'//write_var(hc1(i)),iunit)
            end do
!            h(0*k+1:1*k) = h2(:k)+hc
            do i = 1, k
               call WX('   '//write_var(h(k+i))//' = '// &
                  write_var(h2(k+1-i))//'-'//write_var(hc1(k+1-i)),iunit)
            end do
!            h(1*k+1:2*k) = h2(k:1:-1)-hc(k:1:-1)
            do i = 1, k
               call WX('   '//write_var(h(2*k+i))//' = '// &
                  write_var(hs1(i))//'+'//write_var(h2(k+i)),iunit)
            end do
!            h(2*k+1:3*k) = hs+h2(k+1:)
            do i = 1, k
               call WX('   '//write_var(h(3*k+i))//' = '// &
                  write_var(hs1(k+1-i))//'-'//write_var(h2(2*k+1-i)),iunit)
            end do
!            h(3*k+1:4*k) = hs(k:1:-1)-h2(2*k:k+1:-1)
         end if
      end subroutine odd_tail
!................................................................

!................................................................
! Subroutine scaled_odd_tail_gen(stuff,n,body,get_factor_data)
! Generates constants for scaled_odd_tail and store their indices
!................................................................
      subroutine scaled_odd_tail_gen(stuff,n,body,get_factor_data)
         integer, intent(in) :: n          ! Order of the problem
         type(avl_bag) body                ! Holds actual data
         type(scaled_odd_tail_stuff), intent(out) :: stuff ! Holds
                                           ! Indices of coefficients
         type(get_factor_stuff), intent(in) :: get_factor_data(n) ! Factors
                                           ! thrown out by scaled_odd_tail
         real(qp) factor(2**(n-4))         ! Factors passed through
         integer i                         ! Row number
         integer j                         ! Column select
         integer k                         ! Blocksize
         real(qp) pi                       ! Pi
         real(qp) matrix(2**(n-3))         ! Actual array of coefficients
         real(qp) matrix1(min(2,2**(n-4)),2**(n-4)) ! Actual matrix of coefficients
         real(qp) diag(2**(n-3))           ! Actual array of coefficients

         pi = 4*atan(1.0_qp)
         allocate(stuff%matrixc(2**(n-3)))
         allocate(stuff%matrixs(2**(n-3)))
         allocate(stuff%matrix1c(min(2,2**(n-4)),2**(n-4)))
         allocate(stuff%matrix1s(min(2,2**(n-4)),2**(n-4)))
         allocate(stuff%diag(2**(n-3)))
         if(n == 3) then
! Insert garbage values: should never be used
            stuff%matrixc = 0
            stuff%matrixs = 0
            stuff%diag = 0
         else if(n > 3) then
            factor = get_factor_data(n)%factor
            k = 2**(n-4)
            do i = 1, k
               matrix(2*k+1-i) = tan(2*pi*(2*i-1)/2**n)
               matrix(i) = -tan(2*pi*(2*i-1)/2**n)
            end do
            do i = 1, 2**(n-3)
               call key2index(abs(matrix(i)), body, stuff%matrixc(i))
            end do
            k = 2**(n-4)
            do i = 1, k
               matrix(i) = tan(2*pi*(2*i-1)/2**n)
               matrix(2*k+1-i) = -tan(2*pi*(2*i-1)/2**n)
            end do
            do i = 1, 2**(n-3)
               call key2index(abs(matrix(i)), body, stuff%matrixs(i))
            end do
            k = 2**(n-5)
            matrix1(1,1) = cos(pi/4)
            do i = 1, k
               matrix1(1,i) = cos(2*pi*(2*i-1)/2**(n-1))
               matrix1(1,2*k+1-i) = sin(2*pi*(2*i-1)/2**(n-1))
               matrix1(2,i) = -sin(2*pi*(2*i-1)/2**(n-1))
               matrix1(2,2*k+1-i) = cos(2*pi*(2*i-1)/2**(n-1))
            end do
            k = 2**(n-7)
            do i = 1, k
               matrix1(:,0*k+i) = matrix1(:,0*k+i)* &
                               get_factor_data(n-3)%factor(i)
               matrix1(:,8*k+1-i) = matrix1(:,8*k+1-i)* &
                               get_factor_data(n-3)%factor(i)
               matrix1(:,1*k+i) = matrix1(:,1*k+i)* &
                               get_factor_data(n-3)%factor(k+1-i)
               matrix1(:,7*k+1-i) = matrix1(:,7*k+1-i)* &
                               get_factor_data(n-3)%factor(k+1-i)
               matrix1(:,2*k+i) = matrix1(:,2*k+i)* &
                               get_factor_data(n-3)%factor(i)
               matrix1(:,6*k+1-i) = matrix1(:,6*k+1-i)* &
                               get_factor_data(n-3)%factor(i)
               matrix1(:,3*k+i) = matrix1(:,3*k+i)* &
                               get_factor_data(n-3)%factor(k+1-i)
               matrix1(:,5*k+1-i) = matrix1(:,5*k+1-i)* &
                               get_factor_data(n-3)%factor(k+1-i)
            end do
            k = 2**(n-4)
            do i = 1, k
               matrix1(:,0*k+i) = matrix1(:,0*k+i)/factor(i)
            end do
            do i = 1, 2**(n-4)
               do j = 1, min(2,2**(n-4))
                  call key2index(abs(matrix1(j,i)), body, stuff%matrix1c(j,i))
               end do
            end do
            k = 2**(n-5)
            matrix1(1,1) = cos(pi/4)
            do i = 1, k
               matrix1(1,i) = sin(2*pi*(2*i-1)/2**(n-1))
               matrix1(1,2*k+1-i) = cos(2*pi*(2*i-1)/2**(n-1))
               matrix1(2,i) = cos(2*pi*(2*i-1)/2**(n-1))
               matrix1(2,2*k+1-i) = -sin(2*pi*(2*i-1)/2**(n-1))
            end do
            k = 2**(n-7)
            do i = 1, k
               matrix1(:,0*k+i) = matrix1(:,0*k+i)* &
                               get_factor_data(n-3)%factor(i)
               matrix1(:,8*k+1-i) = matrix1(:,8*k+1-i)* &
                               get_factor_data(n-3)%factor(i)
               matrix1(:,1*k+i) = matrix1(:,1*k+i)* &
                               get_factor_data(n-3)%factor(k+1-i)
               matrix1(:,7*k+1-i) = matrix1(:,7*k+1-i)* &
                               get_factor_data(n-3)%factor(k+1-i)
               matrix1(:,2*k+i) = matrix1(:,2*k+i)* &
                               get_factor_data(n-3)%factor(i)
               matrix1(:,6*k+1-i) = matrix1(:,6*k+1-i)* &
                               get_factor_data(n-3)%factor(i)
               matrix1(:,3*k+i) = matrix1(:,3*k+i)* &
                               get_factor_data(n-3)%factor(k+1-i)
               matrix1(:,5*k+1-i) = matrix1(:,5*k+1-i)* &
                               get_factor_data(n-3)%factor(k+1-i)
            end do
            k = 2**(n-4)
            do i = 1, k
               matrix1(:,0*k+i) = matrix1(:,0*k+i)/factor(i)
            end do
            do i = 1, 2**(n-4)
               do j = 1, min(2,2**(n-4))
                  call key2index(abs(matrix1(j,i)), body, stuff%matrix1s(j,i))
               end do
            end do
            k = 2**(n-6)
            diag = 1
            do i = 1, k
               diag(0*k+i) = get_factor_data(n-2)%factor(i)
               diag(1*k+i) = get_factor_data(n-2)%factor(k+1-i)
               diag(2*k+i) = get_factor_data(n-2)%factor(i)
               diag(3*k+i) = get_factor_data(n-2)%factor(k+1-i)
               diag(4*k+i) = get_factor_data(n-2)%factor(i)
               diag(5*k+i) = get_factor_data(n-2)%factor(k+1-i)
               diag(6*k+i) = get_factor_data(n-2)%factor(i)
               diag(7*k+i) = get_factor_data(n-2)%factor(k+1-i)
            end do
            k = 2**(n-4)
            do i = 1, k
               diag(0*k+i) = diag(0*k+i)/factor(i)
               diag(1*k+i) = diag(1*k+i)/factor(i)
            end do
            do i = 1, 2**(n-3)
               call key2index(abs(diag(i)), body, stuff%diag(i))
            end do
         end if
      end subroutine scaled_odd_tail_gen

!................................................................
! Subroutine scaled_odd_tail(h,n,correct_data,scaled_odd_tail_data,iunit)
! Computes action of scaled odd tail matrix of order 2**(n-1) and
! immortalizes in code.
!................................................................
      recursive subroutine scaled_odd_tail(h,n,correct_data,scaled_odd_tail_data,iunit)
         integer, intent(in) :: n                      ! Order of the problem
         type(var_stuff), intent(inout) :: h(2**(n-1)) ! Array of variable mnemonics
         type(scaled_odd_tail_stuff), intent(in) :: scaled_odd_tail_data(n) ! Indices
                                                       ! of data for scaled_odd_tail
         type(correct_stuff) correct_data              ! Indices for degenerate problems
         integer, intent(in) :: iunit                  ! Write code to this unit
         type(var_stuff) hc(2**(n-3)), hcx(2**(n-3))   ! Temporary variable mnemonics
         type(var_stuff) hs(2**(n-3)), hsx(2**(n-3))   ! Temporary variable mnemonics
         type(var_stuff) h2(2**(n-2))                  ! Temporary variable mnemonics
         type(var_stuff) hc1(2**(n-4)), hc1x(2**(n-4)) ! Temporary variable mnemonics
         type(var_stuff) hs1(2**(n-4)), hs1x(2**(n-4)) ! Temporary variable mnemonics
         type(var_stuff) h21(2**(n-3)), h21x(2**(n-3)) ! Temporary variable mnemonics
         integer i                                     ! Row number
         integer k                                     ! Blocksize
         integer M                                     ! Blocksize
         integer max_depth                             ! Previous maximum depth

         if(n <= 3) then
            call odd_tail_correct(h,n,correct_data,iunit)
         else
            max_depth = maxval(h%depth)
            M = 2**(n-3)
            do i = 1, M
               h2(i) = h(2*i-1)
               hc(i) = var_stuff(max_depth+1,h(2*i)%index,h(2*i)%label)
               hs(i) = var_stuff(max_depth+1,h(2*(i+M))%index,h(2*(i+M))%label)
               h2(i+M) = h(2*(i+M)-1)
            end do
            do i = 1, M
               call WX('   '//write_var(hc(i))//' = '// &
                  write_var(h(4*i-2))//'-'//write_var(h(4*(M+1-i))),iunit)
            end do
!            hc = h(2::4)-h(2**(n-1):1:-4)
            call scaled_odd_tail(hc,n-2,correct_data,scaled_odd_tail_data,iunit)
            max_depth = maxval(hc%depth)
            do i = 1, M
               hcx(i) = var_stuff(max_depth+1,hc(i)%index,hc(i)%label)
            end do
            M = 2**(n-4)
            do i = 1, M
               call WX('   '//write_var(hcx(i))//' = '// &
                  write_var(hc(i))//'-'// &
                  write_const(scaled_odd_tail_data(n)%matrixc(i))//'*'// &
                  write_var(hc(M+i)),iunit)
            end do
            do i = M+1,2*M
               call WX('   '//write_var(hcx(i))//' = '// &
                  write_const(scaled_odd_tail_data(n)%matrixc(i))//'*'// &
                  write_var(hc(2*M+1-i))//'+'//write_var(hc(3*M+1-i)),iunit)
            end do
!            hc = (/(/(hc(i)-scaled_odd_tail_data(n)%matrixc(i)*hc(M+i),i=1,M)/), &
!                   (/(scaled_odd_tail_data(n)%matrixc(i)*hc(2*M+1-i)+hc(3*M+1-i),i=M+1,2*M)/)/)
            M = 2**(n-3)
            do i = 1, M
               call WX('   '//write_var(hs(i))//' = '// &
                  write_var(h(4*i-2))//'+'//write_var(h(4*(M+1-i))),iunit)
            end do
!            hs = h(2::4)+h(2**(n-1):1:-4)
            call scaled_odd_tail(hs,n-2,correct_data,scaled_odd_tail_data,iunit)
            max_depth = maxval(hs%depth)
            do i = 1, M
               hsx(i) = var_stuff(max_depth+1,hs(i)%index,hs(i)%label)
            end do
            M = 2**(n-4)
            do i = 1, M
               call WX('   '//write_var(hsx(i))//' = '// &
                  write_const(scaled_odd_tail_data(n)%matrixs(i))//'*'// &
                  write_var(hs(i))//'+'//write_var(hs(M+i)),iunit)
            end do
            do i = M+1,2*M
               call WX('   '//write_var(hsx(i))//' = '// &
                  write_var(hs(2*M+1-i))//'-'// &
                  write_const(scaled_odd_tail_data(n)%matrixs(i))//'*'// &
                  write_var(hs(3*M+1-i)),iunit)
            end do
!            hs = (/(/(scaled_odd_tail_data(n)%matrixs(i)*hs(i)+hs(M+i),i=1,M)/), &
!                   (/(hs(2*M+1-i)-scaled_odd_tail_data(n)%matrixs(i)*hs(3*M+1-i),i=M+1,2*M)/)/)
!            h2 = h(::2)
            max_depth = maxval(h2%depth)
            M = 2**(n-4)
            do i = 1, M
               h21(i) = h2(2*i-1)
               hc1(i) = var_stuff(max_depth+1,h2(2*i)%index,h2(2*i)%label)
               hs1(i) = var_stuff(max_depth+1,h2(2*(i+M))%index,h2(2*(i+M))%label)
               h21(i+M) = h2(2*(i+M)-1)
            end do
            do i = 1, M
               call WX('   '//write_var(hc1(i))//' = '// &
                  write_var(h2(4*i-2))//'-'//write_var(h2(4*(M+1-i))),iunit)
            end do
!            hc1 = h2(2::4)-h2(2**(n-2):1:-4)
            call scaled_odd_tail(hc1,n-3,correct_data,scaled_odd_tail_data,iunit)
            max_depth = maxval(hc1%depth)
            do i = 1, M
               hc1x(i) = var_stuff(max_depth+1,hc1(i)%index,hc1(i)%label)
            end do
            M = 2**(n-5)
            if(n == 4) then
               call WX('   '//write_var(hc1x(1))//' = '// &
                  write_const(scaled_odd_tail_data(n)%matrix1c(1,1))//'*'// &
                  write_var(hc1(1)),iunit)
            else
               do i = 1, M
                  call WX('   '//write_var(hc1x(i))//' = '// &
                     write_const(scaled_odd_tail_data(n)%matrix1c(1,i))//'*'// &
                     write_var(hc1(i))//'-'// &
                     write_const(scaled_odd_tail_data(n)%matrix1c(2,i))// &
                     '*'//write_var(hc1(M+i)),iunit)
               end do
               do i = M+1,2*M
                  call WX('   '//write_var(hc1x(i))//' = '// &
                     write_const(scaled_odd_tail_data(n)%matrix1c(1,i))//'*'// &
                     write_var(hc1(2*M+1-i))//'+'// &
                     write_const(scaled_odd_tail_data(n)%matrix1c(2,i))//'*'// &
                     write_var(hc1(3*M+1-i)),iunit)
               end do
            end if
!            if(n == 4) then
!               hc1(1) = scaled_odd_tail_data(n)%matrix1c(1,1)*hc1(1)
!            else
!               hc1 = (/(/(scaled_odd_tail_data(n)%matrix1c(1,i)*hc1(i)- &
!                          scaled_odd_tail_data(n)%matrix1c(2,i)*hc1(M+i),i=1,M)/), &
!                       (/(scaled_odd_tail_data(n)%matrix1c(1,i)*hc1(2*M+1-i)+ &
!                          scaled_odd_tail_data(n)%matrix1c(2,i)*hc1(3*M+1-i),i=M+1,2*M)/)/)
!            end if
            M = 2**(n-4)
            do i = 1, M
               call WX('   '//write_var(hs1(i))//' = '// &
                  write_var(h2(4*i-2))//'+'//write_var(h2(4*(M+1-i))),iunit)
            end do
!            hs1 = h2(2::4)+h2(2**(n-2):1:-4)
            call scaled_odd_tail(hs1,n-3,correct_data,scaled_odd_tail_data,iunit)
            max_depth = maxval(hs1%depth)
            do i = 1, M
               hs1x(i) = var_stuff(max_depth+1,hs1(i)%index,hs1(i)%label)
            end do
            M = 2**(n-5)
            if(n == 4) then
               call WX('   '//write_var(hs1x(1))//' = '// &
                  write_const(scaled_odd_tail_data(n)%matrix1s(1,1))//'*'// &
                  write_var(hs1(1)),iunit)
            else
               do i = 1, M
                  call WX('   '//write_var(hs1x(i))//' = '// &
                     write_const(scaled_odd_tail_data(n)%matrix1s(1,i))//'*'// &
                     write_var(hs1(i))//'+'// &
                     write_const(scaled_odd_tail_data(n)%matrix1s(2,i))//'*'// &
                     write_var(hs1(M+i)),iunit)
               end do
               do i = M+1,2*M
                  call WX('   '//write_var(hs1x(i))//' = '// &
                     write_const(scaled_odd_tail_data(n)%matrix1s(1,i))//'*'// &
                     write_var(hs1(2*M+1-i))//'-'// &
                     write_const(scaled_odd_tail_data(n)%matrix1s(2,i))//'*'// &
                     write_var(hs1(3*M+1-i)),iunit)
               end do
            end if
!            if(n == 4) then
!               hs1(1) = scaled_odd_tail_data(n)%matrix1s(1,1)*hs1(1)
!            else
!            hs1 = (/(/(scaled_odd_tail_data(n)%matrix1s(1,i)*hs1(i)+ &
!                       scaled_odd_tail_data(n)%matrix1s(2,i)*hs1(M+i),i=1,M)/), &
!                    (/(scaled_odd_tail_data(n)%matrix1s(1,i)*hs1(2*M+1-i)- &
!                       scaled_odd_tail_data(n)%matrix1s(2,i)*hs1(3*M+1-i),i=M+1,2*M)/)/)
!            end if
!            h21 = h2(::2)
            call scaled_odd_tail(h21,n-2,correct_data,scaled_odd_tail_data,iunit)
            max_depth = maxval(h21%depth)
            M = 2**(n-3)
            do i = 1, M
               h21x(i) = var_stuff(max_depth+1,h21(i)%index,h21(i)%label)
            end do
            do i = 1, M
               call WX('   '//write_var(h21x(i))//' = '// &
                  write_const(scaled_odd_tail_data(n)%diag(i))//'*'// &
                  write_var(h21(i)),iunit)
            end do
!            h21 = scaled_odd_tail_data(n)%diag*h21
            k = 2**(n-4)
            max_depth = max(maxval(hc1x%depth),maxval(hs1x%depth),maxval(h21x%depth))
            do i = 1, k
               h2(2*i-1) = var_stuff(max_depth+1,h21x(i)%index,h21x(i)%label)
               h2(2*i) = var_stuff(max_depth+1,hc1x(i)%index,hc1x(i)%label)
               h2(2*(i+k)) = var_stuff(max_depth+1,hs1x(i)%index,hs1x(i)%label)
               h2(2*(i+k)-1) = var_stuff(max_depth+1,h21x(i+k)%index,h21x(i+k)%label)
            end do
            do i = 1, k
               call WX('   '//write_var(h2(i))//' = '// &
                  write_var(h21x(i))//'+'//write_var(hc1x(i)),iunit)
            end do
!            h2(0*k+1:1*k) = h21(:k)+hc1
            do i = 1, k
               call WX('   '//write_var(h2(k+i))//' = '// &
                  write_var(h21x(k+1-i))//'-'//write_var(hc1x(k+1-i)),iunit)
            end do
!            h2(1*k+1:2*k) = h21(k:1:-1)-hc1(k:1:-1)
            do i = 1, k
               call WX('   '//write_var(h2(2*k+i))//' = '// &
                  write_var(hs1x(i))//'+'//write_var(h21x(k+i)),iunit)
            end do
!            h2(2*k+1:3*k) = hs1+h21(k+1:)
            do i = 1, k
               call WX('   '//write_var(h2(3*k+i))//' = '// &
                  write_var(hs1x(k+1-i))//'-'//write_var(h21x(2*k+1-i)),iunit)
            end do
!            h2(3*k+1:4*k) = hs1(k:1:-1)-h21(2*k:k+1:-1)
            k = 2**(n-3)
            max_depth = max(maxval(hcx%depth),maxval(hsx%depth),maxval(h2%depth))
            do i = 1, k
               h(2*i-1) = var_stuff(max_depth+1,h2(i)%index,h2(i)%label)
               h(2*i) = var_stuff(max_depth+1,hcx(i)%index,hcx(i)%label)
               h(2*(i+k)) = var_stuff(max_depth+1,hsx(i)%index,hsx(i)%label)
               h(2*(i+k)-1) = var_stuff(max_depth+1,h2(i+k)%index,h2(i+k)%label)
            end do
            do i = 1, k
               call WX('   '//write_var(h(i))//' = '// &
                  write_var(h2(i))//'+'//write_var(hcx(i)),iunit)
            end do
!            h(0*k+1:1*k) = h2(:k)+hc
            do i = 1, k
               call WX('   '//write_var(h(k+i))//' = '// &
                  write_var(h2(k+1-i))//'-'//write_var(hcx(k+1-i)),iunit)
            end do
!            h(1*k+1:2*k) = h2(k:1:-1)-hc(k:1:-1)
            do i = 1, k
               call WX('   '//write_var(h(2*k+i))//' = '// &
                  write_var(hsx(i))//'+'//write_var(h2(k+i)),iunit)
            end do
!            h(2*k+1:3*k) = hs+h2(k+1:)
            do i = 1, k
               call WX('   '//write_var(h(3*k+i))//' = '// &
                  write_var(hsx(k+1-i))//'-'//write_var(h2(2*k+1-i)),iunit)
            end do
!            h(3*k+1:4*k) = hs(k:1:-1)-h2(2*k:k+1:-1)
         end if
      end subroutine scaled_odd_tail
!................................................................

!................................................................
! Subroutine real_fft_correct(h,n,iunit)
! Calculates real-half complex dft of order 2**n from definition
! for degenerate cases and immortalizes the result in code
!................................................................
      subroutine real_fft_correct(h,n,iunit)
         integer, intent(in) :: n                      ! Order of the problem
         type(var_stuff), intent(inout) :: h(0:2**n-1) ! Array of variable mnemonics
         integer, intent(in) :: iunit                  ! Write code to this unit
         type(var_stuff) h1(0:2**n-1)                  ! Temporary variable mnemonics
         type(var_stuff) h2(0:2**n-1)                  ! Temporary variable mnemonics
         integer i                                     ! Row index
         integer max_depth                             ! Depth of previous problem

!         M = 2**n
!         h = (/(/(sum((/(h(j)*cos(2*pi*mod(i*j,M)/M),j=0,M-1)/)),i=0,M/2)/), &
!               (/(sum((/(h(j)*sin(2*pi*mod(i*j,M)/M),j=0,M-1)/)),i=1,(M-1)/2)/)/)
         select case(n)
            case(1)
               max_depth = maxval(h%depth)
               do i = 0, 2**n-1
                  h1(i) = var_stuff(max_depth+1,h(i)%index,h(i)%label)
               end do
               if(max_depth == 0) then
                  call WX('   '//write_var(h1(0))//' = '// &
                     write_input(h(0))//'+'//write_input(h(1)),iunit)
                  call WX('   '//write_var(h1(1))//' = '// &
                     write_input(h(0))//'-'//write_input(h(1)),iunit)
               else
                  call WX('   '//write_var(h1(0))//' = '// &
                     write_var(h(0))//'+'//write_var(h(1)),iunit)
                  call WX('   '//write_var(h1(1))//' = '// &
                     write_var(h(0))//'-'//write_var(h(1)),iunit)
               end if
               do i = 0, 2**n-1
                  h(i) = h1(i)
               end do
            case(2)
               max_depth = maxval(h%depth)
               do i = 0, 2**n-1
                  h1(i) = var_stuff(max_depth+1,h(i)%index,h(i)%label)
                  h2(i) = var_stuff(max_depth+2,h(i)%index,h(i)%label)
               end do
               if(max_depth == 0) then
                  call WX('   '//write_var(h1(0))//' = '// &
                     write_input(h(0))//'+'//write_input(h(2)),iunit)
                  call WX('   '//write_var(h1(2))//' = '// &
                     write_input(h(1))//'+'//write_input(h(3)),iunit)
                  call WX('   '//write_var(h2(0))//' = '// &
                     write_var(h1(0))//'+'//write_var(h1(2)),iunit)
                  call WX('   '//write_var(h2(1))//' = '// &
                     write_input(h(0))//'-'//write_input(h(2)),iunit)
                  call WX('   '//write_var(h2(2))//' = '// &
                     write_var(h1(0))//'-'//write_var(h1(2)),iunit)
                  call WX('   '//write_var(h2(3))//' = '// &
                     write_input(h(1))//'-'//write_input(h(3)),iunit)
               else
                  call WX('   '//write_var(h1(0))//' = '// &
                     write_var(h(0))//'+'//write_var(h(2)),iunit)
                  call WX('   '//write_var(h1(2))//' = '// &
                     write_var(h(1))//'+'//write_var(h(3)),iunit)
                  call WX('   '//write_var(h2(0))//' = '// &
                     write_var(h1(0))//'+'//write_var(h1(2)),iunit)
                  call WX('   '//write_var(h2(1))//' = '// &
                     write_var(h(0))//'-'//write_var(h(2)),iunit)
                  call WX('   '//write_var(h2(2))//' = '// &
                     write_var(h1(0))//'-'//write_var(h1(2)),iunit)
                  call WX('   '//write_var(h2(3))//' = '// &
                     write_var(h(1))//'-'//write_var(h(3)),iunit)
               end if
               do i = 0, 2**n-1
                  h(i) = h2(i)
               end do
            case default
               stop 'bigger n than 2 in real_fft_correct'
         end select
      end subroutine real_fft_correct
!................................................................

!................................................................
! Subroutine real_fft(h,n,correct_data,odd_tail_data, &
!                                   scaled_odd_tail_data,iunit)
! Calculates real-half complex fft of order 2**n and
! immortalizes the result in code
!................................................................
      recursive subroutine real_fft(h,n,correct_data,odd_tail_data, &
                                    scaled_odd_tail_data,iunit)
         integer, intent(in) :: n                              ! Order of the problem
         type(var_stuff), intent(inout) :: h(0:2**n-1)         ! Array of variable mnemonics
         type(correct_stuff) correct_data                      ! Indices for degenerate problems
         type(odd_tail_stuff) odd_tail_data(n)                 ! Indices for odd_tail
         type(scaled_odd_tail_stuff) scaled_odd_tail_data(n-2) ! Indices for scaled_odd_tail
         integer, intent(in) :: iunit                          ! Write code to this unit
         type(var_stuff) h0(0:2**(n-1)-1)                      ! even freq. mnemonics
         type(var_stuff) h1(0:2**(n-1)-1)                      ! odd freq. mnemonics
         integer max_depth                                     ! Max depth of previous
         integer M                                             ! Blocksize
         integer i                                             ! Current index

         if(n < 3) then
            call real_fft_correct(h,n,iunit)
         else
            max_depth = maxval(h%depth)
            if(max_depth == 0) then
               M = 2**(n-1)
               do i = 0, 2**(n-1)-1
                  h0(i) = var_stuff(max_depth+1, h(i)%index, h(i)%label)
                  call WX('   '//write_var(h0(i))//' = '// &
                                 write_input(h(i))//'+'//write_input(h(i+M)),iunit)
               end do
!               h0 = h(0:2**(n-1)-1)+h(2**(n-1):)
               call real_fft(h0,n-1,correct_data,odd_tail_data,scaled_odd_tail_data,iunit)
               do i = 0, 2**(n-1)-1
                  h1(i) = var_stuff(max_depth+1, h(i+M)%index, h(i+M)%label)
                  call WX('   '//write_var(h1(i))//' = '// &
                                 write_input(h(i))//'-'//write_input(h(i+M)),iunit)
               end do
!               h1 = h(0:2**(n-1)-1)-h(2**(n-1):)
            else
               M = 2**(n-1)
               do i = 0, 2**(n-1)-1
                  h0(i) = var_stuff(max_depth+1, h(i)%index, h(i)%label)
                  call WX('   '//write_var(h0(i))//' = '// &
                                 write_var(h(i))//'+'//write_var(h(i+M)),iunit)
               end do
!               h0 = h(0:2**(n-1)-1)+h(2**(n-1):)
               call real_fft(h0,n-1,correct_data,odd_tail_data,scaled_odd_tail_data,iunit)
               do i = 0, 2**(n-1)-1
                  h1(i) = var_stuff(max_depth+1, h(i+M)%index, h(i+M)%label)
                  call WX('   '//write_var(h1(i))//' = '// &
                                 write_var(h(i))//'-'//write_var(h(i+M)),iunit)
               end do
!               h1 = h(0:2**(n-1)-1)-h(2**(n-1):)
            end if
               call odd_tail(h1,n,correct_data,odd_tail_data,scaled_odd_tail_data,iunit)
            do i = 0, M-1
               h(2*i) = h0(i)
               h(2*i+1) = h1(i)
            end do
!            h(0::2) = h0
!            h(1::2) = h1
         end if
      end subroutine real_fft
!................................................................

!................................................................
! Subroutine complex_fft(n,correct_data,odd_tail_data,scaled_odd_tail_data,iunit)
! Computes full complex fft of order 2**n and immortalizes results in code
!      subroutine complex_fft(h,n)
!................................................................
      subroutine complex_fft(n,correct_data,odd_tail_data,scaled_odd_tail_data,iunit)
         integer, intent(in) :: n               ! Order of the problem
         type(correct_stuff) correct_data       ! Indices for degenerate problems
         type(odd_tail_stuff) odd_tail_data(n)  ! Indices for odd tail
         type(scaled_odd_tail_stuff) scaled_odd_tail_data(n-2) ! Indices for
                                                ! scaled form of odd tail
         integer, intent(in) :: iunit           ! Write code to this unit
         type(var_stuff) hreal(0:2**n-1)        ! Array of real variable mnemonics
         type(var_stuff) himag(0:2**n-1)        ! Array of imaginary variable mnemonics
         integer M                              ! Blocksize
         integer i                              ! Row number

         M = 2**n
!         hreal = real(h)
         do i = 0, 2**n-1
            hreal(i) = var_stuff(0,i,'xr')
         end do
         call real_fft(hreal,n,correct_data,odd_tail_data,scaled_odd_tail_data,iunit)
!         himag = aimag(h)
         do i = 0, 2**n-1
            himag(i) = var_stuff(0,i,'xi')
         end do
         call real_fft(himag,n,correct_data,odd_tail_data,scaled_odd_tail_data,iunit)
!         h(0) = cmplx(hreal(0),himag(0),dp)
         call WX('   h(0) = cmplx('// &
                        write_var(hreal(0))//','//write_var(himag(0))//',wp)',iunit)
         do i = 1, (M-1)/2
            call WX('   '//write_output(i)//' = cmplx('// &
                        write_var(hreal(i))//'+'//write_var(himag(M/2+i))//','// &
                        write_var(himag(i))//'-'//write_var(hreal(M/2+i))//',wp)',iunit)
!            h(i) = cmplx(hreal(i)-himag(M/2+i),hreal(M/2+i)+himag(i),dp)
            call WX('   '//write_output(M-i)//' = cmplx('// &
                        write_var(hreal(i))//'-'//write_var(himag(M/2+i))//','// &
                        write_var(hreal(M/2+i))//'+'//write_var(himag(i))//',wp)',iunit)
!            h(M-i) = cmplx(hreal(i)+himag(M/2+i),-hreal(M/2+i)+himag(i),dp)
         end do
         call WX('   '//write_output(M/2)//' = cmplx('// &
                        write_var(hreal(M/2))//','//write_var(himag(M/2))//',wp)',iunit)
!         h(M/2) = cmplx(hreal(M/2),himag(M/2),dp)
      end subroutine complex_fft
!................................................................

!................................................................
! function write_var(var)
! Writes out variable name corresponding to mnemonic
!................................................................
      function write_var(var)
         type(var_stuff), intent(in) :: var      ! Variable mnemonic
         character(write_var_len(var)) write_var ! Printed variable name

         write(write_var,'(a,i0,a,i0)') trim(adjustl(var%label)), &
                                        var%depth, '_', var%index
      end function write_var
!................................................................

!................................................................
! function write_var_len(var)
! Specification function to compute len(write_var)
!................................................................
      pure function write_var_len(var)
         integer write_var_len              ! Len(write_var(var))
         type(var_stuff), intent(in) :: var ! Variable mnemonic
         character(80) test_string          ! Holds name temporarily

         write(test_string,'(a,i0,a,i0)') trim(adjustl(var%label)), &
                                        var%depth, '_', var%index
         write_var_len = len_trim(test_string)
      end function write_var_len
!................................................................

!................................................................
! function write_input(var)
! Write initial variable name in h array
!................................................................
      function write_input(var)
         type(var_stuff), intent(in) :: var          ! Variable mnemonic
         character(write_input_len(var)) write_input ! Variable name

         if(scan(var%label,'iI') == 0) then
            write(write_input,'(a,i0,a)') 'real(h(',var%index,'))'
         else
            write(write_input,'(a,i0,a)') 'aimag(h(',var%index,'))'
         end if
      end function write_input
!................................................................

!................................................................
! function write_input_len(var)
! Specification function to compute len(write_input(var))
!................................................................
      pure function write_input_len(var)
         type(var_stuff), intent(in) :: var ! Variable mnemonic
         integer write_input_len            ! Len(write_input(var))
         character(80) test_string          ! Holds name temporarily

         if(scan(var%label,'iI') == 0) then
            write(test_string,'(a,i0,a)') 'real(h(',var%index,'))'
         else
            write(test_string,'(a,i0,a)') 'aimag(h(',var%index,'))'
         end if
         write_input_len = len_trim(test_string)
      end function write_input_len
!................................................................

!................................................................
! function write_const(var)
! Write constant name
!................................................................
      function write_const(var)
         integer, intent(in) :: var                  ! Constant index
         character(write_const_len(var)) write_const ! Constant name

         write(write_const,'(a,i0)') 'd',var
      end function write_const
!................................................................

!................................................................
! function write_const_len(var)
! Specification function to compute len(write_const(var))
!................................................................
      pure function write_const_len(var)
         integer, intent(in) :: var         ! Constant index
         integer write_const_len            ! Len(write_const(var))
         character(80) test_string          ! Holds name temporarily

         write(test_string,'(a,i0)') 'd',var
         write_const_len = len_trim(test_string)
      end function write_const_len
!................................................................

!................................................................
! function write_output(var)
! Write array element name
!................................................................
      function write_output(var)
         integer, intent(in) :: var                    ! Array element index
! The following line causes CVF 6.6C to crash at runtime.  The
! workaround, commented out below, was found by James Giles.
         character(write_output_len(var)) write_output ! Array element name
!         character(write_output_len((var))) write_output ! Array element name

         write(write_output,'(a,i0,a)') 'h(',var,')'
      end function write_output
!................................................................

!................................................................
! function write_output_len(var)
! Specification function to compute len(write_output(var))
!................................................................
      pure function write_output_len(var)
         integer, intent(in) :: var         ! Array element index
         integer write_output_len           ! Len(write_output(var))
         character(80) test_string          ! Holds name temporarily

         write(test_string,'(a,i0,a)') 'h(',var,')'
         write_output_len = len_trim(test_string)
      end function write_output_len
!................................................................

!................................................................
! subroutine WX(string,iunit)
! Writes string to iunit
!................................................................
      subroutine WX(string,iunit)
         integer, intent(in) :: iunit
         character(*), intent(in) :: string

         write(iunit,'(a)') string
      end subroutine WX
!................................................................

end module code_gen
! ***************************************************************

! ***************************************************************
! program test
! Prompts user for size of the problem and creates a dft subroutine
! of order 2**n along with a test program.
! ***************************************************************
program test
   use mykinds
   use code_gen
   use avl_mod
   implicit none
   integer n                                ! Order of the problem
   integer i                                ! Subproblem order
   real(dp) t0, t1                          ! Timer start and stop times
   real(qp), allocatable :: qconst_array(:) ! Array of required constants
   character(80) filename1                  ! Name of output *.f90 file
   character(80) filename2                  ! Name of output *.i90 file
   character(8) date                        ! Date for DATE_AND_TIME
   character(10) time                       ! Time for DATE_AND_TIME
   integer An, Mn                           ! Theoretical opcounts
   type(get_factor_stuff), allocatable :: get_factor_data(:) ! Scale
                                            ! factors from scaled_odd_tail
   type(correct_stuff) correct_data         ! Indices for degenerate problems
   type(odd_tail_stuff), allocatable :: odd_tail_data(:) ! Indices for
                                            ! unscaled odd_tail problem
   type(scaled_odd_tail_stuff), allocatable :: scaled_odd_tail_data(:) ! Indices
                                            ! for scaled odd tail problem
   type(avl_bag) body                       ! Holds actual data
   integer iunit                            ! write code to this unit
   integer junit                            ! test code unit
   character(80) temp_out                   ! Staging area for numeric output

   write(*,'(a)',advance='no') ' Enter the exponent of the problem:> '
   read(*,*) n
   if(n <= 0) then
      write(*,'(a)') ' Sorry, the exponent must be positive.'
      stop
   else if(n > 24) then
      write(*,'(a)') ' Sorry, the maximum value is 2**24.'
      write(*,'(a,i0,a)') ' Your input corresponds to 2**',n,'.'
      stop
   end if
   write(*,*) ' Creating DFT code:'
   call cpu_time(t0)
! Write test program
   write(filename1,'(a,i0,a)') 'fft', 2**n, '.f90'
   junit = 11
   open(junit,file=filename1,status='replace')
   write(junit,'(a,a)') '! File: ', trim(filename1)
   call date_and_time(date, time)
   write(junit,'(12a)') '! Generated by scaled2h.f90 ', date(5:6), &
                    '/', date(7:8), '/', date(1:4), ' ', time(1:2), &
                    ':', time(3:4), ':', time(5:10)
   write(junit,'(a)') 'module mykinds'
   write(junit,'(a)') '   implicit none'
   write(junit,'(a)') '   integer, parameter :: sp = selected_real_kind(6,30)'
   write(junit,'(a)') '   integer, parameter :: dp = selected_real_kind(15,300)'
   write(junit,'(a)') 'end module mykinds'
   write(junit,'(a)') ''
   write(junit,'(a)') 'module dp_mod'
   write(junit,'(a)') '   use mykinds, only: wp => dp'
   write(junit,'(a)') '   implicit none'
   write(junit,'(a)') '   private wp'
   write(junit,'(a)') '   contains'
   write(junit,'(a,i0,a)') 'include ''fft',2**n,'.i90'''
   write(junit,'(a)') 'end module dp_mod'
   write(junit,'(a)') ''
   write(junit,'(a)') 'program test'
   write(junit,'(a)') '   use mykinds'
   write(junit,'(a)') '   use dp_mod'
   write(junit,'(a)') '   implicit none'
   write(junit,'(a,i0)') '   integer, parameter :: N = ',2**n
   write(junit,'(a)') '   complex(dp) h1(0:N-1)'
   write(junit,'(a)') '   complex(dp) h2(0:N-1)'
   write(junit,'(a)') '   complex(dp), parameter :: ri(2) = (/(1,0),(0,1)/)'
   write(junit,'(a)') '   integer i, j, k, L'
   write(junit,'(a)') '   real(dp) pi'
   write(junit,'(a)') '   real(dp) err, maxerr'
   write(junit,'(a)') ''
   write(junit,'(a)') '   pi = 4*atan(1.0_dp)'
   write(junit,'(a)') '   maxerr = 0'
   write(junit,'(a)') '   do i = 0, N-1'
   write(junit,'(a)') '      do j = 1, 2'
   write(junit,'(a)') '         h2 = 0'
   write(junit,'(a)') '         h2(i) = ri(j)'
   write(junit,'(a)') '         h1 = (/(h2(i)*exp(-2*pi*(0,1)*mod(i*k,N)/N),k=0,N-1)/)'
   write(junit,'(a,i0,a)') '         call fft',2**n,'(h2)'
   write(junit,'(a)') '         err = max(maxval(abs(real(h1-h2))),maxval(abs(aimag(h1-h2))))'
   write(junit,'(a)') '         write(*,*) i, j, err'
   write(junit,'(a)') '         maxerr = max(err, maxerr)'
   write(junit,'(a)') '      end do'
   write(junit,'(a)') '   end do'
   write(junit,'(a)') '   write(*,*)'
   write(junit,'(a)') "   write(*,*) 'maxerr = ',maxerr"
   write(junit,'(a)') 'end program test'
   close(junit)
! Write DFT subroutine file
   write(filename2,'(a,i0,a)') 'fft', 2**n, '.i90'
   iunit = 10
   open(iunit,file=filename2,status='replace')
   write(iunit,'(a,a)') '! File: ', trim(filename2)
   call date_and_time(date, time)
   write(iunit,'(12a)') '! Generated by scaled2h.f90 ', date(5:6), &
                    '/', date(7:8), '/', date(1:4), ' ', time(1:2), &
                    ':', time(3:4), ':', time(5:10)
   write(iunit,'(a,i0,a)') 'subroutine fft', 2**n, '(h)'
   An = (8*3*n*2**n-16*2**n-2*(-1)**n)/9+2
   Mn = (3*10*n*2**n-76*2**n-2*3*n*(-1)**n+22*(-1)**n)/27-2*n+6
   write(iunit,'(a,i0,a,i0,a)') '! ',An,' adds, ',Mn,' muls'
   write(iunit,'(a)') '   implicit real(wp) (x,y)'
   write(iunit,'(a,i0,a)') '   complex(wp), intent(inout) :: h(0:', 2**n-1, ')'
   allocate(get_factor_data(n-2))
   do i = 1, n-2
      allocate(get_factor_data(i)%factor(2**(i-4)))
      if(i >= 4) then
         call get_factor(get_factor_data(i)%factor,i,get_factor_data(i-2)%factor)
      end if
   end do
   call prepare_correct(n, body, correct_data)
   allocate(odd_tail_data(n))
   do i = 1, n
      call odd_tail_gen(odd_tail_data(i),i,body,get_factor_data)
   end do
   allocate(scaled_odd_tail_data(n-2))
   do i = 1, n-2
      call scaled_odd_tail_gen(scaled_odd_tail_data(i),i,body,get_factor_data)
   end do
   allocate(qconst_array(body%num_nodes))
   call bag2array(body, qconst_array)
   write(*,'()')
   call cpu_time(t1)
   do i = 1, size(qconst_array)
      write(temp_out,'(g41.34)') qconst_array(i)
      write(iunit,'(a,i0,a,a,a)') '   real(wp), parameter :: d', &
                                  i, ' = ', trim(adjustl(temp_out)), '_wp'
   end do
   deallocate(qconst_array)
   write(*,'(a,f0.3)') ' Prep time = ', t1-t0
   call cpu_time(t0)
   call complex_fft(n,correct_data,odd_tail_data,scaled_odd_tail_data,iunit)
   call cpu_time(t1)
   write(*,'(a,f0.3)') ' Execution time = ', t1-t0
   write(iunit,'(a,i0)') 'end subroutine fft', 2**n
   write(*,'(a)') ' Files created: '//trim(filename1)//', '//trim(filename2)
end program test
! ***************************************************************

