! Program to generate terms for A212222. Input: The number l_max. ! Output: A b-file for A212222 up to 3**l_max - 1. ! This program will work, in principle, up to l_max = 40 ! (around 2^63), but calculation time grows exponentially with l_max. ! Currently (2020), l_max = 30 will take around 200 CPU hours on ! a fast system (8-10 CPU cycles per number). If your CPU has ! the POPCNT instruction, be sure to use -march=native or equivalent ! to gain speed. ! This program is in the public domain, no warranty of any ! kind etc. ! Thomas Koenig, 2020-08-17 program main implicit none integer, parameter :: ip = selected_int_kind(15) integer :: l_max integer(kind=ip) :: n_level integer :: num num = 0 read (*,*) l_max write (0,'(A,I0)') "Checking up to ", l_max call btsum (0_ip, 0_ip, l_max) write (0,'(A)') "Stop!" 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, s , level-1) call btsum (3*n+1,s+1, level-1) call btsum (3*n+2,s+2, level-1) return end if if (popcnt (n) == s) then call do_something (n, s) end if end subroutine btsum subroutine do_something (n,s) integer (kind=ip), value :: n, s if (qs5(n) /= s) return if (qs7(n) /= s) return if (qs11(n) /=s) return num = num + 1 print '(I0," ",I0)', num, n flush (6) end subroutine do_something function qs5(n) result(qs) integer(kind=ip), parameter :: div = 5 integer(kind=ip), value :: n integer(kind=ip) :: qs qs = 0 do while (n > 0) qs = qs + mod(n, div) n = n / div end do end function qs5 function qs7(n) result(qs) integer(kind=ip), parameter :: div = 7 integer(kind=ip), value :: n integer(kind=ip) :: qs qs = 0 do while (n > 0) qs = qs + mod(n, div) n = n / div end do end function qs7 function qs11(n) result(qs) integer(kind=ip), parameter :: div = 11 integer(kind=ip), value :: n integer(kind=ip) :: qs qs = 0 do while (n > 0) qs = qs + mod(n, div) n = n / div end do end function qs11 end program main