C Count all complete rulers of a given length (OEIS A103295)
C and all complete rulers with a given number of segments (A103296)
C
C Author: Hugo Pfoertner http://www.pfoertner.org/
C URL of this file:
C http://www.randomwalk.de/sequences/a103295.txt
C
C Version history
C Mar 4 2005 Difference counting by shifts
C Mar 3 2005 Initial version
C
C This program depends on the availability of 64-bit integers:
implicit integer*8 (a-z)
logical btest
C Length of lookup table for bit counting
parameter(tbits=21, mask = 2**tbits-1)
C mmax has to be set to (range of ruler lengths to be investigated)+1
parameter (mmax=37)
dimension bitset(0:mask), rulcnt(0:mmax)
C Prepare lookup table for bit counting
do 20 i = 0, mask
bc = 0
do 25 k = 0,tbits
25 if ( btest(i,k) ) bc = bc + 1
bitset(i) = bc
20 continue
C array for counting rulers
do 15 i = 0, mmax
rulcnt(i) = 0
15 continue
C
C Loop over number of interior positions
do 100 m = 1, mmax
last = 2**m-1
count = 0
C Loop over all possible rulers with length m+1
do 10 r = 0, last
C set all bits (no difference covered)
difmis = last
C
C Check differences of set bits
C Distances of set bits from left boundary
C This clears the bits at all positions of the marks
difmis = iand ( difmis, not(-not(r)-1) )
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
do 30 j = 1, m
C Clear bit positions measured from right boundary
if ( btest(r,m-j) ) difmis = ibclr ( difmis, m-j )
C Check only positions not cleared so far to avoid the expensive shift
if (btest(difmis,j)) then
C If there is at least one common position in the original bit field and
C the bit field shifted by j positions, then difference j is cleared
if ( iand( r, ishft(r,j)) .ne. 0 ) difmis=ibclr ( difmis, j )
endif
30 continue
C check if all bits have been cleared
if ( difmis .eq. 0 ) then
count = count + 1
C number of set bits in this ruler
bits = bitset(iand(r,mask))+ bitset(ishft(r,-tbits))
rulcnt(bits) = rulcnt(bits) + 1
endif
C End of loop over all rulers of length m
10 continue
write (*,*) m, count
100 continue
do 200 i = 1, mmax
if ( rulcnt(i) .gt. 0 ) write(*,*) i, rulcnt(i)
200 continue
end