I would like to generate in an efficient way a list of integers (preferably ordered) with the following defining properties:
All integers have the same number of bit set
N.All integers have the same sum of bit indices
K.
To be definite, for an integer I
its binary representation is:
$I=\sum_{j=0}^M c_j 2^j$ where $c_j=0$ or $1$
The number of bit sets is:
$N(I)=\sum_{j=0}^M c_j$
The sum of bit indices is:
$K(I)=\sum_{j=0}^M j c_j$
I have an inefficient way to generate the list as follows: make a do/for loop over integers incrementing by use of a "snoob" function - smallest next integer with same number of bit set and at each increment checking if it has the correct value of K
this is grossly inefficient because in general starting from an integer
with the correct N and K value the snoob integer from I does not have the correct K and one has to make many snoob calculations to get the next integer
with both N and K equal to the chosen values.
Using snoob gives an ordered list which is handy for dichotomic search but
not absolutely compulsory.
Counting the number of elements in this list is easily done by recursion when viewed as a partition numner counting. here is a recursive function in fortran 90 doing that job:
=======================================================================
recursive function BoundedPartitionNumberQ(N, M, D) result (res)
implicit none
! number of partitions of N into M distinct integers, bounded by D
! appropriate for Fermi counting rules
integer(8) :: N, M, D, Nmin
integer(8) :: res
Nmin = M*(M+1)/2 ! the Fermi sea
if(N < Nmin) then
res = 0
else if((N == Nmin) .and. (D >= M)) then
res = 1
else if(D < M) then
res = 0
else if(D == M) then
if(N == Nmin) then
res = 1
else
res = 0
endif
else if(M == 0) then
res = 0
else
res = BoundedPartitionNumberQ(N-M,M-1,D-1)+BoundedPartitionNumberQ(N-M,M,D-1)
endif
end function BoundedPartitionNumberQ
========================================================================================
My present solution is inefficient when I want to generate lists with several $10^7$
elements. Ultimately I want to stay within the realm of C/C++/Fortran and reach lists of lengths
up to a few $10^9$
my present f90 code is the following:
program test
implicit none
integer(8) :: Nparticles
integer(8) :: Nmax, TmpL, CheckL, Nphi
integer(8) :: i, k, counter
integer(8) :: NextOne
Nphi = 31 ! word size is Nphi+1
Nparticles = 16 ! number of bit set
print*,Nparticles,Nphi
Nmax = ishft(1_8, Nphi + 1) - ishft(1_8, Nphi + 1 - Nparticles)
i = ishft(1, Nparticles) - 1
counter = 0
! integer CheckL is the sum of bit indices
CheckL = Nparticles*Nphi/2 ! the value of the sum giving the largest list
do while(i .le. Nmax) ! we increment the integer
TmpL = 0
do k=0,Nphi
if (btest(i,k)) TmpL = TmpL + k
end do
if (TmpL == CheckL) then ! we check whether the sum of bit indices is OK
counter = counter + 1
end if
i = NextOne(i) ! a version of "snoob" described below
end do
print*,counter
end program
!==========================================================================
function NextOne (state)
implicit none
integer(8) :: bit
integer(8) :: counter
integer(8) :: NextOne,state,pstate
bit = 1
counter = -1
! find first one bit
do while (iand(bit,state) == 0)
bit = ishft(bit,1)
end do
! find next zero bit
do while (iand(bit,state) /= 0)
counter = counter + 1
bit = ishft(bit,1)
end do
if (bit == 0) then
print*,'overflow in NextOne'
NextOne = not(0)
else
state = iand(state,not(bit-1)) ! clear lower bits i &= (~(bit-1));
pstate = ishft(1_8,counter)-1 ! needed by IBM/Zahir compiler
! state = ior(state,ior(bit,ishft(1,counter)-1)) ! short version OK with gcc
state = ior(state,ior(bit,pstate))
NextOne = state
end if
end function NextOne
