anchor.diff:  module sort
   
   use const, only: int_32, int_64, p
  +
  +implicit none
   
   private
   public :: qsort
  
positive.code: interface qsort
      module procedure qsort_int_32_list
      module procedure qsort_int_64_list
  end interface qsort
  
  ! Insertion sort...
  
negative.code: 
      end subroutine qsort_int_64_list
  
      pure subroutine qsort_psip_info(nstates, states, pops, dat)
  
          ! Sort a set of psip information (states, populations and data) in order according
          ! to the state labels.
  
          ! states(:,i) is regarded as greater than states(:,j) if the first
          ! non-identical element between states(:,i) and states(:,j) is greater in
          ! states(:,i).
  
          ! In:
          !    nstates: number of occupied states.
          ! In/Out:
          !    states: 2D array of i0 integers containing the state label for each occupied state.
          !        Sorted on output.
          !    pops, dat: population and data arrays for each state.  Sorted by states on output.
  
          use const, only: int_p, i0, p
          use bit_utils, only: operator(.bitstrge.), operator(.bitstrgt.)
  
          integer, intent(in) :: nstates
          integer(i0), intent(inout) :: states(:,:)
          integer(int_p), intent(inout) :: pops(:,:)
          real(p), intent(inout) :: dat(:,:)
  
          ! Threshold.  When a substates gets to this length, switch to using
          ! insertion sort to sort the substates.
          integer, parameter :: switch_threshold = 7
  
          ! sort needs auxiliary storage of length 2*log_2(n).
          integer, parameter :: stack_max = 50
  
          integer :: pivot, lo, hi, i, j
          integer(i0) :: tmp_state(ubound(states,dim=1))
          integer(int_p) :: tmp_pop(ubound(pops,dim=1))
          real(p) :: tmp_dat(ubound(dat,dim=1))
  
          ! Stack.  This is the auxilliary memory required by quicksort.
          integer :: stack(2,stack_max), nstack
  
          hi = nstates
  
          nstack = 0
          lo = 1
  
          do
              ! If the section/partition we are looking at is smaller than
              ! switch_threshold then perform an insertion sort.
              if (hi - lo < switch_threshold) then
                  do j = lo + 1, hi
                      tmp_state = states(:,j)
                      tmp_pop = pops(:,j)
                      tmp_dat = dat(:,j)
                      do i = j - 1, 1, -1
                          if (tmp_state .bitstrge. states(:,i)) exit
                          states(:,i+1) = states(:,i)
                          pops(:,i+1) = pops(:,i)
                          dat(:,i+1) = dat(:,i)
                      end do
                      states(:,i+1) = tmp_state
                      pops(:,i+1) = tmp_pop
                      dat(:,i+1) = tmp_dat
                  end do
  
                  if (nstack == 0) exit
                  hi = stack(2,nstack)
                  lo = stack(1,nstack)
                  nstack = nstack - 1
  
              else
                  ! Otherwise start partitioning with quicksort.
  
                  ! Pick the pivot element to be the median of states(:,lo), states(:,hi)
                  ! and states(:,(lo+hi)/2).
                  ! This largely overcomes a major problem with quicksort, where it
                  ! degrades if the pivot is always the smallest element.
                  pivot = (lo + hi)/2
                  call swap_states(states(:,pivot), pops(:,pivot), dat(:,pivot), states(:,lo+1), pops(:,lo+1), dat(:,lo+1))
                  if (states(:,lo) .bitstrgt. states(:,hi)) then
                      call swap_states(states(:,lo), pops(:,lo), dat(:,lo), states(:,hi), pops(:,hi), dat(:,hi))
                  end if
                  if (states(:,lo+1) .bitstrgt. states(:,hi)) then
                      call swap_states(states(:,lo+1), pops(:,lo+1), dat(:,lo+1), states(:,hi), pops(:,hi), dat(:,hi))
                  end if
                  if (states(:,lo) .bitstrgt. states(:,lo+1)) then
                      call swap_states(states(:,lo), pops(:,lo), dat(:,lo), states(:,lo+1), pops(:,lo+1), dat(:,lo+1))
                  end if
  
                  i = lo + 1
                  j = hi
                  tmp_state = states(:,lo+1) ! a is the pivot value
                  tmp_pop = pops(:,lo+1)
                  tmp_dat = dat(:,lo+1)
                  do while (.true.)
                      ! Scan down states to find element > a.
                      i = i + 1
                      do while (tmp_state .bitstrgt. states(:,i))
                          i = i + 1
                      end do
  
                      ! Scan down states to find element < a.
                      j = j - 1
                      do while (states(:,j) .bitstrgt. tmp_state)
                          j = j - 1
                      end do
  
                      ! When the pointers crossed, partitioning is complete.
                      if (j < i) exit
  
                      ! Swap the elements, so that all elements < a end up
                      ! in lower indexed variables.
                      call swap_states(states(:,i), pops(:,i), dat(:,i), states(:,j), pops(:,j), dat(:,j))
                  end do
  
                  ! Insert partitioning element
                  states(:,lo + 1) = states(:,j)
                  pops(:,lo + 1) = pops(:,j)
                  dat(:,lo + 1) = dat(:,j)
                  states(:,j) = tmp_state
                  pops(:,j) = tmp_pop
                  dat(:,j) = tmp_dat
  
                  ! Push the larger of the partitioned sections onto the stack
                  ! of sections to look at later.
                  ! --> need fewest stack elements.
                  nstack = nstack + 1
  
                  ! With a stack_max of 50, we can sort arrays of length
                  ! 1125899906842624.  It is safe to say this will never be
                  ! exceeded, and so this test can be skipped.
  !                if (nstack > stack_max) call stop_all('qsort_int_64_states', "parameter stack_max too small")
  
                  if (hi - i + 1 >= j - lo) then
                      stack(2,nstack) = hi
                      stack(1,nstack) = i
                      hi = j - 1
                  else
                      stack(2,nstack) = j - 1
                      stack(1,nstack) = lo
                      lo = i
                  end if
  
              end if
          end do
  
      contains
  
          pure subroutine swap_states(s1,p1,d1,s2,p2,d2)
  
              integer(i0), intent(inout) :: s1(:), s2(:)
              integer(int_p), intent(inout) :: p1(:), p2(:)
              real(p), intent(inout) :: d1(:), d2(:)
              integer(i0) :: tmp_state(ubound(s1,dim=1))
              integer(int_p) :: tmp_pop(ubound(p1,dim=1))
              real(p) :: tmp_dat(ubound(d1,dim=1))
  
              tmp_state = s1
              s1 = s2
              s2 = tmp_state
  
              tmp_pop = p1
              p1 = p2
              p2 = tmp_pop
  
              tmp_dat = d1
              d1 = d2
              d2 = tmp_dat
  
          end subroutine swap_states
  
      end subroutine qsort_psip_info
  
  !--- In-place insertion sort.  ---
  
      pure subroutine insert_sort_int_32(list)
  