! Program to print the longest consecutive sequence of integers in
! A330904 up to 3**l_max.  This is to verify the conjecture that
! there are at most four of these.

! If you want to run this efficiently, make sure to switch on the
! highest level of optimization and to select your architecture
! (using gfortran, that would be done with "-O3 -march=native") so
! the popcnt intrinsic is directly expressed as a machine instruction
! instead of a call to a library routine.

! Public domain, no warranty of any kind, etc.
! Thomas König, 2020

program main
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  integer, parameter :: ip = selected_int_kind(15)
  integer, parameter :: l_max = 30
  integer, parameter :: n_seq = 10  ! Let's be optimistic...
  integer :: seq_max, i_seq, i
  integer(kind=ip), dimension(n_seq) :: seq
  real (kind=dp) :: t0, t1, t2
  seq_max = 1
  i_seq = 1
  seq(1) = 0
  seq (2:) = -1
  call cpu_time (t0)
  t1 = t0
  do i=0, l_max
     write (*,'(A,I3)', advance="no") "Testing level ", i
     call btsum (1_ip, 1_ip, i)
     call cpu_time (t2)
     write (*,'(A,G12.5,A,ES12.5,A)') " cpu_time = ", t2-t1, ", ", &
          3_dp**i / (t2-t0), " total iterations per second"
     t1 = t2
  end do
contains
  ! This is the fastest way I know of to generate A065363.
  recursive subroutine btsum (n, s, level)
    integer(kind=ip), value :: n, s
    integer, value :: level
    if (level /= 0) then
       call btsum (3*n-1, s-1, level-1)
       call btsum (3*n  , s  , level-1)
       call btsum (3*n+1, s+1, level-1)
       return
    end if
    if (s > 0 .and. popcnt (n) == s) then
       call do_something (n)
    end if
  end subroutine btsum
  subroutine do_something (n)
    integer (kind=ip), value :: n
    if (i_seq >= n_seq) then
       stop "Overreaching success!"
    end if
    if (n - 1 == seq(i_seq)) then
       i_seq = i_seq + 1
       if (i_seq > seq_max) then
          write (*,'()')
          print *, seq(1:i_seq-1),n
          seq_max = i_seq
       end if
    else
       i_seq = 1
    end if
    seq(i_seq) = n
  end subroutine do_something
end program main