' pips2.bas BCX 4.75 bitsieve of Eratosthenes routine
' Cino Hilliard Rev 05/02/05
'^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
'This algorithm is a modification and improvement of sieve.c by Frank Pilhofer
'<fp@fpx.de> Sieve.c is included in the zip and compiles with lcc. I sped up
'the org routine by stopping the outer loop of the masking routine at the
'sqr(maxn) vs maxn and starting the mask at l*l below. this cuts the time in
'half.
'^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
'$nowin
const bitset(a,b) = (*(a+(b)/16) |= *(ptr_mask2+imod(b,16)/2)) 'define bit set.
const bitread(a,b) = (*(a+(b)/16) & *(ptr_mask2+imod(b,16)/2)) 'define bit read
set mask2[16] as ULONG
128,64,32,16,8,4,2,1 'force forward bits
end set
'set mask[16] as ULONG
'1,2,4,8,16,32,64,128 'required for nope
'end set
set ptr_mask2 as PULONG
mask2
end set
dim t1!,t2!,c$,n1$,l2,tmp$
'declare some variables c style
! unsigned long j,j1,primes,tprimes,index,pips,k,l,alloc;
! unsigned long maxn,n;
! unsigned char *ptr_a=NULL, *ptr_m;
cls
c$ = command$ 'get command tail
if val(c$) > 10000 then
maxn = val(c$)+10000
else
maxn = alloc = 100000000L ' Can be up to to 4.29 billion+
end if
ptr_a = ptr_m = malloc(maxn)
print " Bit Sieve of Eratosthenes "
print " An analysis of Primes,PIPS2 and Twin Primes"
print
!xprintf(" Sieving %lu numbers start address = %p\n",maxn,ptr_a);
print
t1! = timer
for j=0 to alloc-1 'clear allocated memory
*ptr_m++ = 0
next j
'*************** The incredible masking routine *************************;
for l = 3 to sqr(maxn) step 2 'sqr(maxn) speeds it up
l2=l+l
if not bitread(ptr_a,l) then 'skip if already set
for j = l*l to maxn step l2 'start at square of l vs 3*l
bitset(ptr_a,j)
next j
end if
next l
for j = 0 to 50 'Print some of the binary string comma delimited
tmp$ = bin$(*(ptr_a+j))
tmp$ = repeat$(8-len(tmp$),"0") & tmp$
for k = 1 to 8
print mid$(tmp$,k,1);",";
next k
next
print "first 32 forms " 'be careful! These bit patterns
for j = 0 to 32 'are backwards. I fixed this!
tmp$ = bin$(*(ptr_a+j))
tmp$ = repeat$(8-len(tmp$),"0") & tmp$
print *(ptr_a+j);",";
next j
print
t2!=timer
!xprintf("It took %f sec to mask out multiples of primes < %lu \n",t2-t1,maxn);
print " Press <ENTER>on blank to exit"
print
' ******************** c like print out *********************
do
input "number ",n1$
n = val(n1$)
if n=0 then break
primes=1
tprimes=0
pips=1
t1=timer
for k = 3 to n-1 step 2
if not bitread(ptr_a,k) then 'test for primes
primes++ 'count of regular prime
if not bitread(ptr_a,(k+1)) then 'test for twin primes
incr tprimes 'count twin primes
end if
if not bitread(ptr_a,primes) and mod(primes,2) then 'test for pips
pips++ 'track no of pips
end if
end if
next k
t2=timer
print
print "Primes,twins,pips in a range statistics"
print
!xprintf("primes < %lu = %lu\n",n,primes);
!xprintf("twin p < %lu = %lu\n",n,tprimes);
print "twins/primes ";(float) tprimes/primes
!xprintf("pips < %lu = %lu\n",n,pips);
print "pips/primes ";(float) pips/primes
print "pips/twins ";(float) pips/tprimes
print "twins - pips ";(float) tprimes-pips
print "took ";t2!-t1!;" sec to print"
loop until n = 0
free(ptr_m)
end
