Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
Prime Contest
#11
Now seriously:
Provisional scores for FB:

my (Rich's) entry Execution time: 2.585 s
yetifoot's second entry Execution time: 4.377 s
yetifoot's second entry Execution time: 2.131 s

Rich Geldreich's method is not that slow after all....
Antoni
Reply
#12
Indeed, I found rich geldreichs method fast, even considering the fact its aimed at qb 1.1 and has to account for that.

I spotted another optimization for my second one, i've just edited the post, as it wasn't enough to warrent posting the whole code again. I wasn't checking n1 <= sqr(maxval), even though i knew about that trick, i hadn't been able to find where to add it until know.

I think my entry is now as fast as i can get it, but i already thought that about 5 times already, so who knows? Smile

Is anyone else thinking of entering some code?

Maybe its just me and you Antoni, i fear i'll be in last place if that is the case Smile
EVEN MEN OF STEEL RUST.
[Image: chav.gif]
Reply
#13
Your modified last entry scores 1.27 seconds, it nearly doubles its previous speed...

Im' surprised other people has'nt showed up. Is someone coding something?
Antoni
Reply
#14
I think i have killed the competition by myself...

First I forgot QB1.1 can't do huge arrays so a sieve up to 1.6M to find the 100,000th prime is out of limits.

Then I posted Rich Geldeich's solution for QB4.5, thinking a sieve was faster. In fact it is..in FreeBASIC. With all that huge arrays, index offseting and bit twiddling a sieve in QB4.5 is actually slower than Rich's solution.

So perhaps this is the reason because no one else enters.

Here is a sieve for QB 4.5


Code:
'Prime sieve for QB4.5  by Antoni Gual 10/2006
'Run QB4.5 with /AH
  
DEFLNG A-H
DEFINT I-M
amaxp = 1300000
'index offset of -4000 so it does not reach the limit of 16737...
REDIM p(-4000 TO amaxp \ 64 - 4000) AS LONG

'to avoid bit rotations
DIM powers2(31) AS LONG
b = 1
FOR i = 0 TO 31
  powers2(i) = b
  IF i = 30 THEN EXIT FOR
  b = b + b
NEXT
powers2(31) = &H80000000

ctarget = 1000
cnt = 1
asqrt = INT(SQR(amaxp))

FOR b = 3 TO amaxp STEP 2
IF (p((b \ 64) - 4000) AND powers2((b AND 63) \ 2)) = 0 THEN
   cnt = cnt + 1
   IF cnt = ctarget THEN
     PRINT ctarget; b
     IF ctarget = 100000 THEN SYSTEM
     ctarget = ctarget * 10
   END IF
   IF b <= asqrt THEN
    FOR bb = b * b TO amaxp STEP 2 * b
      p((bb \ 64) - 4000) = p((bb \ 64) - 4000) OR powers2((bb AND 63) \ 2)
    NEXT
   END IF
END IF
NEXT
Antoni
Reply
#15
There are a few optimizations in there that I hadn't thought of.

1. Using a lookup table for the powers, to make the bit-twiddling faster.

2. Starting the inner loop at n * n, i had discovered that i only needed to start an 3n, but n ^ 2 is even better!

3. Makeing the bitarray \ 64, because theres no need to store data for even numbers, I had thought of this, but I never got it working in my program, that halves the memory overheads.

I added these to my program, but it was only a very very small speed increase (less than 5%), so I haven't updated my code.
EVEN MEN OF STEEL RUST.
[Image: chav.gif]
Reply
#16
There is a "killer" optimization possible yet.
If you make the program to print what prime candidate is sieving, you will find it spends a lot of time in sieving the small primes as 3,5,7 because of the high number of bits to set.
As the bit pattern is periodical for a few primes sieved, you can sieve for those small integers in a multiple of the period that fills compeletely a number of integers, then memcopy that small part to the rest of the array.

Another optimization I have found is to separate the count from the sieve.

This way you have to sieve only up to sqr(n). The primes above it do not have multiples in the table (remember we start sieving at p*p), so no check is needed, the separate count finds them.
The count can be very fast if you use a bit counting trick as the one I posted at The Code Post.

EDITED:
I have a source ...and it's fast
Code:
'n-sieve-bits in FreeBASIC by Antoni Gual

#include "crt.bi"
option explicit
option escape

dim shared p as uinteger ptr, m as integer
#define subset 3*5*7*11
#define precount 5
#define lastpre  11
#define firstsieve 13
'
'-----------------------------------------------
sub mymemcopy
dim as uinteger ptr p1,p2
dim as integer i,j
p1=p+subset
  for i=subset to m\64+1-subset step subset
    p2=p1-subset
    for j=1 to subset
      *p1=*p2
      p1+=1
      p2+=1
    next
  next
  p2=p1-subset
  while p1<m\64+1
    *p1=*p2
    p1+=1
    p2+=1
  wend
end sub
'
'----------------------------------------------------
sub presieve
'sieve a small part of the sieve for 3,5,7,11, marking
'p*1 as composites

dim as integer i,nn,cs
i=3
  for nn=i to subset*64 step 2*i
    p[nn\64] or= 1u shl((nn and 63)\2)
  next
i=5
cs=2
while i<=lastpre
  for nn=i to (subset*64) step 2*i
     p[nn\64] or= 1u shl((nn and 63)\2)
  next
  i=i+cs
  cs=6-cs
wend
end sub
'
'---------------------------------------------------
sub dosieve
'does not sieve even nrs or multiples of 3
'starts at 11, previous sieving has been made by preload and mymemcopy
'starts checking multiples at p*p
'stops trying at sqrt(maxprime), as p*p multiples of bigger
'  candidates  would fall past the end of the sieve.


dim as integer sm,i,nn,cs,ccs,ccs1
  sm=int(sqr(m))
  i=firstsieve
  cs=iif((i mod 3)=2,2,4)
  'sieve  
  do
    if (p[i\64] and (1u shl((i and 63)\2)))=0u then
      nn=i*i
      ccs=I* iif ((I mod 3)=2,2,4)
      ccs1=6*I
      do
        p[nn\64] or= 1u shl((nn and 63)\2)
        nn+=ccs
        ccs=ccs1-ccs
      loop until nn>m
    end if
    i+=cs
    cs=6-cs
  loop until i>sm
end sub
'
'------------------------------------------------------
sub docount
'count unset bytes
'
dim as integer cnt,cnt1,i,i1,i2,n=1000
dim as unsigned integer u,uc
cnt=precount  
i=0
do
'count unset bytes by 4 bytes groups
do
   u=p[i]
   '
   'bit count snippet (See it at the Code post)

   uc= u-((u shr 1) and &o33333333333)-((u shr 2) and &o11111111111)
   cnt+=32-( ((uc+(uc shr 3)) and &o30707070707) mod 63 )

   i+=1
loop until cnt>=n

'if count too high, count down single unset bytes in the last
'sieve 4 bytes group counted
if cnt>n then
   cnt1=cnt
   i2=32
   do
     i2-=1
     if (u and (1u shl i2))=0 then cnt1-=1
   loop until cnt1=n  
  
   printf("the %7dth prime is %9d\n",cnt1,i*64-(31-i2)*2-1)
else
   printf("the %7dth prime is %9d\n", cnt, i*64+1)

endif

n*=10
loop until n>1000000
end sub
'
'-------------------------------------------------------
dim as integer m1
m=16e6
m1=(m\64+1) *sizeof(uinteger)
p =callocate(m1)
presieve
mymemcopy
dosieve
docount
deallocate p
Antoni
Reply
#17
Ok so here's a QB entry - you have to use /AH on the command line to stop the overflow.
Code:
TT = TIMER
DEFLNG I-N
N = 21000
DIM A(0 TO N) AS LONG
DIM B(0 TO 30) AS LONG
B(0) = 1
FOR I = 1 TO 30
B(I) = B(I - 1) + B(I - 1)
NEXT I
FOR J = 2 TO 571
J1 = J \ 31: J2 = J MOD 31
IF (A(J1) AND B(J2)) = 0 THEN
K1 = J + J - 1
K2 = K1 * J - K1 + J
FOR I = K2 TO 650000 STEP K1
I1 = I \ 31: I2 = I MOD 31
A(I1) = (A(I1) OR B(I2))
NEXT I
END IF
NEXT J
FOR I = 2 TO 650000
I1 = I \ 31: I2 = I MOD 31
L = A(I1) AND B(I2)
IF L = 0 THEN
L1 = L1 + 1
IF L1 + 1 = 1000 THEN PRINT I + I - 1
IF L1 + 1 = 10000 THEN PRINT I + I - 1
IF L1 + 1 = 100000 THEN PRINT I + I - 1
END IF
NEXT I
PRINT TIMER - TT
Reply
#18
Quibbler:
Your QB sieve is faster than mine!
Yours uses the idea of separing the count and sieving only to sqrt(n). That had to give good results!
Antoni
Reply


Forum Jump:


Users browsing this thread: 1 Guest(s)