!  Demonstration program for the 3x+1 calculation.  This version is
!  typically used for very small-scale calculations.  In particular,
!  the INTEGER type is not large enough for serious exploration.
!  Incorporates test for sequences which cycle through X[0].
!  Ernest Sibert - September 1996

!  Perform n steps of the 3x+1 calculation.  The current state of the
!  calculation is
!
!    X    - Current value in sequence at each position
!    X0   - Starting value for sequence at each position
!    C    - Number of steps in sequence at each position
!    Cmax - Maximum steps till stopping at each position
!    Cy   - Total number of iterations performed so far

      subroutine stepper(n)
      integer :: n       ! Number of steps to perform
   

      INCLUDE 'threex.h' ! Global search state (common)
      logical :: cycle_found
      integer :: k           ! Local loop counter
      integer :: Y(0:NC-1)   ! Working array
      
      cycle_found = .false.
      k = 1
      do while (cycle_found = .false. .and. k <= n)
      
        Y = X/2                      ! Integer division
        WHERE  ( btest(X, 0) )       ! Odd X
          X = X + Y + 1              ! (3*X + 1)/2
         ELSEWHERE
          X = Y                      ! X/2
        END WHERE
        C = C + 1                    ! One more step at each position
        if ( ANY(X == X0) ) then     ! Sequence has a cycle (unlikely)
          call report_cycles()        ! Report the cycle(s)
          cycle_found = .true.
        end if
        WHERE ( X < X0 )             ! Sequence finished
          X0 = X0 + step             ! Advance to new sequence
          X = X0                     ! Begin exploration
          Cmax = max(Cmax, C)        ! Record longest sequence
          C = 0                      ! Reset to count next sequence
        END WHERE
        k = k + 1
      end do
      Cy = Cy + n                    ! Total number of steps
      end subroutine stepper


!  Report cycles in the current sequences.  It is extremely unlikely that
!  this routine will ever be used except in an artificial test.

      subroutine report_cycles()

      INCLUDE 'threex.h'

      integer I(0:NC-1), p     ! Coordinates, "processor" number
      logical cyc(0:NC-1)      ! Flags locations of cycles

      print *, "CYCLE(S) DETECTED AT Processor      X0          C"
      I = (/ (k, k=0,NC-1) /);  cyc = (X == X0)   ! Marks positions with cycles
      do while ( ANY(cyc) )
        p = MINVAL(i, mask=cyc)         ! Lowest remaining (unreported) cycle
        print 21, p, X0(p), C(p)        ! Report the cycle
        cyc(p) = .false.                ! This cycle has been reported
      end do
  21  FORMAT (22X, I7, 3X, I10, 2X, I10)
      print *
      end subroutine report_cycles


!  Artificial test of report_cycles.  This routine is careful to avoid
!  any net chage of the computation state.

      subroutine testrc()

      INCLUDE 'threex.h'

      integer I(0:NC-1), Xsv(0:NC-1)     ! To save X

      print *, "Test of cycle report -- odd positions should report cycles."
      Xsv = X;  I = (/ (k, k=0,NC-1) /)
      where ( btest(I, 0) )              ! Odd-numbered entries
        X = X0
      elsewhere
        X = X0 + 25
      end where
      call report_cycles()
      X = Xsv                            ! So X isn't ultimately changed
      end subroutine testrc


!  Initialize the search

      subroutine start()

      INCLUDE 'threex.h'

      X0 = (/ (4*k + 3, k=0,NC-1) /)
      X = X0;  C = 0;  Cmax = 0;  Cy = 0
      step = 4*NC
      print *, "Initialized with ", NC, " classes, step = ", step
      end subroutine start


!  Run n steps of the calculation, report results

      subroutine run(n)

      INCLUDE 'threex.h'

      integer n

      call stepper(n);
      print *, "Confirmed to ", minval(X0) - 1, &
               ",  longest sequence ", maxval(Cmax), &
               ",  ", Cy, " cycles."
      end subroutine run


!  Report complete state of computation.  Works in any situation, but
!  intended for use with small-scale test calculations.

      subroutine report_state()

      INCLUDE 'threex.h'

      integer p

      print *, "Processor      X0           X         Cmax          C"
      do p = 0, NC-1
        print 21, p, X0(p), X(p), Cmax(p), C(p)
  21    FORMAT (I8, 4(2X, I10))
      end do
      end subroutine report_state


!  Program to manage interactive experiments
!  Commands are:
!                r     - run
!                  <n>   -   n steps (separate prompt)
!                i     - (re)initialize calculation
!                s     - report complete state
!                t     - test cycle report
!                q     - quit

      program threex

      INCLUDE 'threex.h'

      integer   ns         ! Number of steps
      character ibuff*80   ! Input buffer
      logical   more

      print 21
  21  FORMAT ("Commands are:"/  &
              "   r - run (prompts for steps)      s - print state"/  &
              "   i - (re)initialize calculation   q - quit       "/  &
              "   t - test cycle report"                            )
      ibuff(1:1) = 'i';  more = 'true'       ! To force initialization
      print *, "Initializing ..."
      do while (more)
        if      ( ibuff(1:1) == 'i' ) then
          call start()
        else if ( ibuff(1:1) == 's' ) then
          call report_state()
        else if ( ibuff(1:1) == 't' ) then
          call testrc()
        else if ( ibuff(1:1) == 'r' ) then
          print '("steps > "$)';  read *, ns
          call run(ns)                       ! Run the steps
        else if ( ibuff(1:1) == 'q' ) then
          more = false
        else
          print *, "Inappropriate command - ignored."
        end if
        if ( more ) then
          print '("i r s t q > "$)';  read *, ibuff
        end if
      end do
      end program threex
