1,000,000 primes in 6.9 seconds!

BlitzMax Forums/BlitzMax Programming/1,000,000 primes in 6.9 seconds!

Nate the Great

(Posted 2010) [#1]

so I was reading some wiki articles on number theory and prime numbers and I had an idea for a faster prime number program. This one finds the first 1,000,000 prime numbers in ~~6.9 2.3~~ 0.4 seconds on my computer which is not much faster than average. Don't be fooled, it finishes calculating the numbers before it displays them. My IDE freezes when bmax is printing the numbers and then unfreezes several seconds later after it has printed all the numbers so I have left the thing that displays the numbers out but you can just uncomment it if you want to see them.

edit: new faster method

Global num:Int = 15485864
Global sqrnum:Int = Sqr(num)+1
Global PrimeFlags:Byte[num]
Global Primes:Int[num]
Local tim:Int
Local mill:Int = MilliSecs()
For Local x:Int = 2 Until sqrnum
	If primeflags[x] = 0 Then
		tim = x
		Repeat
			tim:+x
			primeflags[tim] = True
		Until x+tim => num
	EndIf
Next 
Local count:Int
For Local y:Int = 2 Until num
	If primeflags[y] = False Then
		primes[count] = y
		count:+1
		
	EndIf
Next
mill = MilliSecs()-mill
Print mill/1000.0 +" seconds"
Rem
For Local i:Int = 0 Until count
	Print primes[i]
Next
EndRem
Print count + " primes found!"

heres the old code:

Local PR:Int = 3
Local primelist:Int[1000001]
Local prcnt:Int = 1
primelist[0] = 2

Local n:Int = 0
Local flag:Byte = False

Local c:Int = MilliSecs()
Repeat
	n = 0
	flag = False
	While primelist[n]*primelist[n] <= PR
		If ((pr Mod primelist[n])=0) Then
			flag = True
			Exit
		EndIf
		n:+1
	Wend
	If Not flag
		primelist[prcnt] = pr
		prcnt:+1
	EndIf
	pr:+2
Until prcnt = 1000000
c = MilliSecs()-c

'For i:Int = 0 Until 1000000
'	Print primelist[i]
'Next

Print "First 1,000,000 primes in "+c/1000.0+" seconds!"

UPDATED

if you want to see the numbers without risking freezing your computer, you could just display the first 1000 or the last 1000

be sure NOT to run in DEBUG MODE

Last edited 2010

Last edited 2010

Last edited 2010

Last edited 2010

Last edited 2010

xlsior

(Posted 2010) [#2]

Just out of curiosity, what processor do you have in your PC?
It takes 12.91 seconds on my computer (Core 2 Duo E6600 @ 2.4GHz, windows 7 x64)

Nate the Great

(Posted 2010) [#3]

core 2 duo t8300 @ 2.4 ghz & 2.4 ghz
this is what my control panel says... although I thought it was 2.6 ghz until now...

and its win vista 32 bit

make sure you are not in debug and you have no other programs running... i consistently get 6.9 something, never more or less than that for the first 2 digits

ok just did it in debug and got 14.5 seconds.. maybe thats youre problem?

youre processor http://ark.intel.com/Product.aspx?id=27250
mine http://ark.intel.com/Product.aspx?id=33099

they look almost identicle to my untrained eye

edit: looks like the time it takes to do this is more affected by read/write speeds for your ram considering it reads a ton of prime numbers off of your ram for every number it tests.

Last edited 2010

Dabhand

(Posted 2010) [#4]

First 1,000,000 primes in 5.90299988 seconds!

On bootcamped XP!

Intel i3 CPU 550@...

Dabz

PowerPC603

(Posted 2010) [#5]

I got 15.5 seconds by running this code on my laptop. (non-debug mode)
With debug mode on, I got 23.97 seconds.

Win Vista 32 bit, processor: Core 2 Duo T7200 @ 2GHz.

Last edited 2010

Pengwin

(Posted 2010) [#6]

Just for comparison, I ran this on my iMac (2.66 GHz Intel Core 2 Duo) with OSX 10.5.8

I got 6.339 seconds for non-debug and 15.796 seconds for debug mode.

Czar Flavius

(Posted 2010) [#7]

If you want a challenge make a multithreaded one!

zzz	(Posted 2010) [#8]

The oiginal code ran in 8.71700001 seconds on a Core 2 Duo T6400. Changing the line

If (pr/(primelist[n]*1.0)) = Int(pr/primelist[n]) Then

If ((pr Mod primelist[n])=0) Then

made it run in 2.86400008 seconds instead :p

Last edited 2010

Last edited 2010

Arabia

(Posted 2010) [#9]

It took 2hrs and 23minutes on my Vic 20 :(

Nate the Great

(Posted 2010) [#10]

@zzz thanks! I cant believe I overlooked that! thanks

edit: updated!

Last edited 2010

Hummelpups

(Posted 2010) [#11]

:D

First 1,000,000 primes in 1.79700005 seconds!

Nate the Great

(Posted 2010) [#12]

hmmm translated it to ti-basic and it takes at least several hours if not days on my calculator. I didnt want to run my batteries down!

Shortwind

(Posted 2010) [#13]

:D

Last edited 2010

Me.32

(Posted 2010) [#14]

translated to C:
<5.809s> on my Slow 2GHz Core 2 Duo iMac :D

BMax:
<5.918s> on the same machine

iMac, Intel Core 2 Dou 2.0 GHZ

<update>
C:
<5.726s> with LLVM 1.5 on same Machine :D

Last edited 2010

Shortwind

(Posted 2010) [#15]

I guess my computer sucks. I've got a (dual) 2.5ghz cpu, and in both win 7 (32bit) and ubuntu (32bit) it runs in 10 seconds, over and over again. hmm, time to upgrade I guess.

;(

PowerPC603

(Posted 2010) [#16]

After the code-change (suggested by zzz), it now runs in 5.97 seconds instead of 15.5 seconds on my laptop.

Shortwind

(Posted 2010) [#17]

Nate, or someone else, here is the sieve of eratosthene. Super fast.

Done in 0.962000012 seconds.
Total primes generated: 1000000

This is the time on my computer.

Can someone help me clean up the code. It's needs to be better able to find a specific number of primes.

But I think you'll agree this is much faster.
After the prime arrary is built, it could be turned into a bit table and take up far less ram.

The end for loop is to print sections of the prime array to verify the numbers are actually primes.

Have fun.

P.S. Don't program when your drunk! ha ha

SuperStrict

Const stopme:Int=100000000
Local flags:Int[stopme+1]
Local i:Int
Local theprimes:Int[1000000]

For i=0 To stopme
	flags[i]=-1
Next

Local z:Int=0
Local endprime:Int=15485863
Local test:Int=2
Local x:Int=0
Local count:Int=0

z=MilliSecs()

For Local test:Int=2 To endprime

	If flags[test]<>-1 Then
		Continue
	EndIf
	
	x=test
	If flags[test]=-1 Then
		theprimes[count]=test
		count=count+1
		flags[test]=test
	EndIf

	Repeat
		flags[test+x]=test
		x=x+test
	Until x>=endprime

Next

z=MilliSecs()-z

Print "Done in "+z/1000.0+" seconds."
Print "Total primes generated: "+count


For i=25 To 50
	Print theprimes[i]
Next

Last edited 2010

Last edited 2010

Last edited 2010

Last edited 2010

Last edited 2010

Shortwind

(Posted 2010) [#18]

Hmm. my bit array for the first 1,000,000 primes is still 125,000 bytes long.... :)

Nate the Great

(Posted 2010) [#19]

thats interesting, I dont see how using that method would be twice as fast as mine... but i get 1.17 seconds with yours which is twice as fast as mine and 7 times as fast as the original! :p I wonder what the fastest way ever though of is to find primes... I cant seem to find many methods for finding prime numbers.

Matty

(Posted 2010) [#20]

Isn't there something called the "Sieve of Eratosthenes " or something - very quick to find primes....

http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes

Shortwind

(Posted 2010) [#21]

Um, LOL, Matty, that's what I programmed. Not very well, but it does work.

Nate:

The sieve is super fast because it's doing what the computer does best:
count... (no division, just simple memory copying and very simple condition loops...)

Using elipitical curves is the fastest method for large numbers, but not small numbers.

Sieve of Eratosthenes is the fastest for most primes, but of course you have to have the memory to build the array.

For my program (badly written) it only needs:
Total Loops: 47,086,551

Your program needs:
Total Loops: 504,508,950

The only real fast way, is to create a table (bit database) of precalculated primes and just do a lookup. That is really the fastest way if your program is going to need to access primes very frequently. (Factoring for example.)

If your interested start here:
http://primes.utm.edu/prove/

And if you want to make your head hurt go here:
http://polymathprojects.org/category/finding-primes/

And of course you could also start with the obvious:
http://en.wikipedia.org/wiki/Integer_factorization

:D

Last edited 2010

Matty

(Posted 2010) [#22]

Sorry Shortwind, I didn't read all the posts (my bad...)

Jesse

(Posted 2010) [#23]

I doubt this will make much difference to anyone here as everyone seem to have pretty fast computers but by removing the second if statement in the for loop and removing one unnecessary calculation in the last demo, I was able to reduce the time on my wife's ancient imac 1.2ghz from 13 secs to 7.5 secs. I tried it in my netbook and it was prety much unchanged.

SuperStrict

Const stopme:Int=100000000
Local flags:Int[stopme+1]
Local i:Int
Local theprimes:Int[1000000]

For i=0 To stopme
	flags[i]=-1
Next

Local z:Int=0
Local endprime:Int=15485863
Local test:Int=2
Local x:Int=0
Local count:Int=0

z=MilliSecs()

For Local test:Int=2 To endprime

	If flags[test]<>-1 Then
		Continue
	EndIf
	
	x=test
	theprimes[count]=test
	count=count+1
	flags[test]=test

	Repeat
		x:+test
		flags[x]=test
	Until x>=endprime

Next

z=MilliSecs()-z

Print "Done in "+z/1000.0+" seconds."
Print "Total primes generated: "+count


For i=25 To 50
	Print theprimes[i]
Next

Last edited 2010

Shortwind

(Posted 2010) [#24]

LOL, thanks Jesse. I didn't see that until you pointed it out. :D

It's better now, but I still code for crap on the fly. Ha! Ha!

:D

zzz	(Posted 2010) [#25]

Thats nice :p Feels a little weird that 10x more iterations only doubles the time taken. Also, theres a sieve of atkin (sp?) around, which, if i dont remember wrong uses some clever way to calculate prime candidates instead of just working through the entire 2..n range like we are doing here.

zzz	(Posted 2010) [#26]

Wrote a simple sieve of atkin based on the wikipedia pseudocode;

SuperStrict
Framework brl.standardio
Import brl.math

Const LIMIT:Int=15500000
Local isprime:Int[LIMIT]

'timer
Local time:Int=MilliSecs()

Rem
########################################################
sieve of atkin starts here :P
End Rem
For Local i:Int=5 To isprime.length-1
	isprime[i]=False
Next
For Local x:Int=1 To Sqr(LIMIT)
	For Local y:Int=1 To Sqr(LIMIT)
	
		Local n:Int
		
		n=(4*x*x)+(y*y)
		If (n<=LIMIT) And ((n Mod 12 = 1) Or (n Mod 12 = 5)) Then
			isprime[n]=1-isprime[n]
		EndIf
		
		n=3*x*x+(y*y)
		If (n<=LIMIT) And (n Mod 12 = 7) Then
			isprime[n]=1-isprime[n]
		EndIf
		
		n=3*x*x-(y*y)
		If (x>y) And (n<=LIMIT) And (n Mod 12 = 11) Then
			isprime[n]=1-isprime[n]
		EndIf
	
	Next
Next
Rem
sieve out false primes
End Rem
For Local i:Int=5 To Sqr(LIMIT)
	If (isprime[i]) Then
		Local k:Int=i*i
		Repeat
			isprime[k]=False
			k:+(i*i)
		Until (k>LIMIT)
	EndIf
Next
Rem
done, put primes in a separate array
End Rem
Local primes:Int[LIMIT]
Local prime_index:Int=2
primes[0]=2
primes[1]=3
For Local i:Int=5 To isprime.length-1
	If (isprime[i]) Then
		primes[prime_index]=i
		prime_index:+1
	EndIf
Next
Rem
########################################################
DONE, write result
End Rem
time=MilliSecs()-time
Print time/1000.0+" seconds.."+prime_index+" primes.."
Rem
validate result against the original function
End Rem
Print "validating.."
Local ref_primes:Int[]=primes_01(prime_index)
For Local i:Int=0 To prime_index-1
	If (ref_primes[i]<>primes[i]) Then
		Print "mismatch.."
		Print ref_primes[i]+":"+primes[i]+":"+i
		Exit
	EndIf
Next
Print "Done!"
End
Rem
########################################################
original code by nate the great for reference
End Rem
Function primes_01:Int[](primes_count:Int)
	Local PR:Int=3
	Local primelist:Int[primes_count]
	Local prcnt:Int=1
	primelist[0]=2
	Local n:Int=0
	Local flag:Byte=False
	Repeat
		n=0
		flag=False
		While primelist[n]*primelist[n]<=PR
			If (PR/(primelist[n]*1.0))=Int(PR/primelist[n]) Then
				flag=True
				Exit
			EndIf
			n:+1
		Wend
		If Not flag
			primelist[prcnt]=pr
			prcnt:+1
		EndIf
		PR:+2
	Until prcnt=primes_count
	Return primelist
End Function

The original function runs in about 2,9 seconds, the "other sieve" :o takes about 1,4 seconds on average and this one does the work in just under 0,8 for me :D

Supposedly this one can be optimized dramatically also.

Shortwind

(Posted 2010) [#27]

I completely forgot about this one. This sieve only needs 16,911,909 loops to find 1,000,862 primes. Nice! And doesn't require the all the extra memory that sieve of Eratosthenes. Cool.

On one final note it should be known that all primes greater than 6 are of the form:

6k + or - 1

:D

Last edited 2010

Nate the Great

(Posted 2010) [#28]

0.4 seconds!

I cant believe how much faster its gotten. This is kind of a mix of the 4 ideas presented above by zzz, jessie and shortwind and me. I don't know what you would call it since I just combined ideas but its a sieve I guess.

Global num:Int = 15485864
Global sqrnum:Int = Sqr(num)+1
Global PrimeFlags:Byte[num]
Global Primes:Int[num]
Local tim:Int
Local mill:Int = MilliSecs()
For Local x:Int = 2 Until sqrnum
	If primeflags[x] = 0 Then
		tim = x
		Repeat
			tim:+x
			primeflags[tim] = True
		Until x+tim => num
	EndIf
Next 
Local count:Int
For Local y:Int = 2 Until num
	If primeflags[y] = False Then
		primes[count] = y
		count:+1
		
	EndIf
Next
mill = MilliSecs()-mill
Print mill/1000.0 +" seconds"
Rem
For Local i:Int = 0 Until count
	Print primes[i]
Next
EndRem
Print count + " primes found!"

Last edited 2010

Last edited 2010

Last edited 2010

zzz	(Posted 2010) [#29]

0.52 seconds here, Nate is on top again :p
EDIT: just figured out how i can make some massive improvements on the atkin sieve ;D

Last edited 2010

ima747

(Posted 2010) [#30]

8 core 2010 2.4ghz Mac pro (under OS X, under vista in bootcamp first post variant was actually a tiny bit slower, tough oddly the same install through parallels was the same as the mac... very wierd)

First post variant (current form with the first optimization update) 2.48600006 solidly

nate's release 2 posts from here up 0.186000004 to 0.189999998

running stuff in the background but that shouldn't matter...

would REALLY love to see a multi-threaded variant... 8 cores, with hyper-threading wants to know... :0)

Last edited 2010

Last edited 2010

ima747

(Posted 2010) [#31]

nevermind this post, was just OS randomness

Last edited 2010

zzz	(Posted 2010) [#32]

@im747, couldnt see any difference in the generated assembly so i guess its just a 'placebo' effect :)

zzz	(Posted 2010) [#33]

Got the aktin sieve down to 0.58 seconds with some basic optimizations, Nates latest version is still slightly faster :p

SuperStrict
Framework brl.standardio
Import brl.math

Const LIMIT:Int=15485864
Local isprime:Int[LIMIT]

'timer
Local time:Int=MilliSecs()

Rem
########################################################
sieve of atkin starts here :P
End Rem
For Local i:Int=5 To isprime.length-1
	isprime[i]=False
Next
Rem
'the original prime candidate generation code
For Local x:Int=1 To Sqr(LIMIT)
	Local xx:Int=x*x
	For Local y:Int=1 To Sqr(LIMIT)
	
		Local n:Int
		Local yy:Int=y*y
		
		n=4*xx+yy
		If (n<=LIMIT) Then
			If ((n Mod 12 = 1) Or (n Mod 12 = 5)) Then
				isprime[n]=1-isprime[n]
			EndIf
		EndIf
		
		n=3*xx+yy
		If (n<=LIMIT) And (n Mod 12 = 7) Then
			isprime[n]=1-isprime[n]
		EndIf
		
		n=3*xx-yy
		If (x>y) And (n<=LIMIT) And (n Mod 12 = 11) Then
			isprime[n]=1-isprime[n]
		EndIf
	
	Next
Next
End Rem
Rem
############################################################
optimized prime candidate code

# optimization 1, reduced x/y loop overhead

End Rem
Rem
4xx+yy
End Rem
For Local x:Int=1 To Sqr(LIMIT)
	Local xx:Int=x*x
	'only make y go up to values that keep 4xx+yy under LIMIT
	Local y_limit:Int=Sqr(LIMIT-(4*xx))
	For Local y:Int=1 To y_limit
		Local yy:Int=y*y
		Local n:Int=4*xx+yy
		'test if the n<=LIMIT check is still needed?
		'If (n<=LIMIT) And ((n Mod 12=1) Or (n Mod 12=5)) Then
		If (n Mod 12=1) Or (n Mod 12=5) Then
			isprime[n]=1-isprime[n]
		EndIf
	Next
Next
Rem
3xx+yy
End Rem
For Local x:Int=1 To Sqr(LIMIT)
	Local xx:Int=x*x
	'only make y go up to values that keep 3xx+yy under LIMIT
	Local y_limit:Int=Sqr(LIMIT-(3*xx))
	For Local y:Int=1 To y_limit
		Local yy:Int=y*y
		Local n:Int=3*xx+yy
		'If (n<=LIMIT) And (n Mod 12=7) Then
		If (n Mod 12=7) Then
			isprime[n]=1-isprime[n]
		EndIf
	Next
Next
Rem
3xx-yy
End Rem
For Local x:Int=1 To Sqr(LIMIT)
	Local xx:Int=x*x
	'special case here for -yy.. y should only go up to x-1 instead
	For Local y:Int=1 To x-1
		Local yy:Int=y*y
		Local n:Int=3*xx-yy
		If (x>y) And (n<=LIMIT) And (n Mod 12=11) Then
			isprime[n]=1-isprime[n]
		EndIf
	Next
Next
Rem
End Rem
Print (MilliSecs()-time)/1000.0
Rem
sieve out false primes
End Rem
For Local i:Int=5 To Sqr(LIMIT)
	If (isprime[i]) Then
		Local k:Int=i*i
		Repeat
			isprime[k]=False
			k:+(i*i)
		Until (k>LIMIT)
	EndIf
Next
Rem
End Rem
Print (MilliSecs()-time)/1000.0
Rem
done, put primes in a separate array
End Rem
Local primes:Int[LIMIT]
Local prime_index:Int=2
primes[0]=2
primes[1]=3
For Local i:Int=5 To isprime.length-1 Step 2
	If (isprime[i]) Then
		primes[prime_index]=i
		prime_index:+1
	EndIf
Next
Rem
########################################################
DONE, write result
End Rem
time=MilliSecs()-time
Print time/1000.0+" seconds.."+prime_index+" primes.."
Rem
validate result against the original function
End Rem
Print "validating.."
Local ref_primes:Int[]=primes_01(prime_index)
For Local i:Int=0 To prime_index-1
	If (ref_primes[i]<>primes[i]) Then
		Print "mismatch.."
		Print ref_primes[i]+":"+primes[i]+":"+i
		Exit
	EndIf
Next
Print "Done!"
End
Rem
########################################################
original code by nate the great for reference
End Rem
Function primes_01:Int[](primes_count:Int)
	Local PR:Int=3
	Local primelist:Int[primes_count]
	Local prcnt:Int=1
	primelist[0]=2
	Local n:Int=0
	Local flag:Byte=False
	Repeat
		n=0
		flag=False
		While primelist[n]*primelist[n]<=PR
			If (PR/(primelist[n]*1.0))=Int(PR/primelist[n]) Then
				flag=True
				Exit
			EndIf
			n:+1
		Wend
		If Not flag
			primelist[prcnt]=pr
			prcnt:+1
		EndIf
		PR:+2
	Until prcnt=primes_count
	Return primelist
End Function

Last edited 2010

zzz	(Posted 2010) [#34]

Slight improvement on nates latest code, keep finishing at around 0,49 seconds for me. Will probably run slower on older machines.

SuperStrict
Framework brl.standardio
Import brl.math

Global num:Int = 15485864
Global sqrnum:Int = Sqr(num)+1
Global PrimeFlags:Byte[num]
Global Primes:Int[num]
Local tim:Int
Local mill:Int = MilliSecs()
For Local x:Int = 2 Until sqrnum
	If primeflags[x] = 0 Then
		tim = x
		Repeat
			tim:+x
			primeflags[tim] = True
		Until x+tim => num
	EndIf
Next 
Local count:Int=0
primes[count]=2
count:+1
For Local y:Int = 3 Until num Step 2
	primes[count]=y
	count:+(primeflags[y]=False)
Next
mill = MilliSecs()-mill
Print mill/1000.0 +" seconds"
Rem
For Local i:Int = 0 Until count
	Print primes[i]
Next
EndRem
Print count + " primes found!"

Nate the Great

(Posted 2010) [#35]

first 10,000,000 primes in 6.0 seconds

SuperStrict
Framework brl.standardio
Import brl.math

Global num:Int = 179424674
Global sqrnum:Int = Sqr(num)+1
Global PrimeFlags:Byte[num]
Global Primes:Int[num]
Local tim:Int
Local mill:Int = MilliSecs()
For Local x:Int = 2 Until sqrnum
	If primeflags[x] = 0 Then
		tim = x
		Repeat
			tim:+x
			primeflags[tim] = True
		Until x+tim => num
	EndIf
Next 
Local count:Int=0
primes[count]=2
count:+1
For Local y:Int = 3 Until num Step 2
	primes[count]=y
	count:+(primeflags[y]=False)
Next
mill = MilliSecs()-mill
Print mill/1000.0 +" seconds"
Rem
For Local i:Int = 0 Until count
	Print primes[i]
Next
EndRem
Print count + " primes found!"

thanks for the optimization zzz

ima747

(Posted 2010) [#36]

Sweet

2.87800002 seconds
10000000 primes found!

zzz	(Posted 2010) [#37]

Nice :P Finished in 7.13600016 seconds for me

zzz	(Posted 2010) [#38]

SuperStrict
Framework brl.standardio
Import brl.math

Global num:Int = 179424674
'Global num:Int=15485864
Global sqrnum:Int = Sqr(num)+1
Global PrimeFlags:Byte[num]
Global Primes:Int[num]
Local tim:Int
Local mill:Int = MilliSecs()
'special case for 2,3
For Local x:Int=2 To 3
	If (primeflags[x]=0) Then
		Local tim:Int=x
		Repeat
			tim:+x
			primeflags[tim]=True
		Until x+tim>=num
	EndIf
Next
'Print (MilliSecs()-mill)/1000.0
'mark of the rest
Rem
For Local i:Int = 6 Until sqrnum-1 Step 6
	Local x:Int
	x=i-1
	If primeflags[x] = 0 Then
		tim = x
		Repeat
			tim:+x
			primeflags[tim] = True
		Until x+tim => num
	EndIf
	x=i+1
	If primeflags[x]=0 Then
		tim=x
		Repeat
			tim:+x
			primeflags[tim]=True
		Until x+tim>=num
	EndIf
Next 
End Rem
Local x:Int=5
Repeat
	If primeflags[x]=0 Then
		Local tim:Int=x
		Repeat
			tim:+x
			primeflags[tim]=True
		Until x+tim>=num
	EndIf
	x:+2
	If primeflags[x]=0 Then
		Local tim:Int=x
		Repeat
			tim:+x
			primeflags[tim]=True
		Until x+tim=>num
	EndIf
	x:+4
Until x>=sqrnum-1
'..
Local count:Int=0
primes[count]=2
count:+1
For Local y:Int = 3 Until num Step 2
	primes[count]=y
	count:+(primeflags[y]=False)
Next
mill = MilliSecs()-mill
Print mill/1000.0 +" seconds"
Rem
For Local i:Int = 0 Until count
	Print primes[i]
Next
EndRem
Print count + " primes found!"

Modified Nates a bit, finishes in around 6,85 for me now.

zzz	(Posted 2010) [#39]

SuperStrict

Framework brl.standardio
Import brl.math

'Local max_value:Int=179424674
Local max_value:Int=15485864

Local flags:Int[max_value/32+1]

Const mod32_mask:Int=%11111
Const div32_mask:Int=%11111111111111111111111111100000
Const even_mask:Int=%01010101010101010101010101010101

Local mask:Int[32]
For Local i:Int=0 To 31
	mask[i]=2^i
Next

'start
Local timer:Int=MilliSecs()

For Local i:Int=0 To flags.length-1
	flags[i]=even_mask
Next
flags[0]:~mask[2]

Local sqrnum:Int=Sqr(max_value)

For Local x:Int=3 To sqrnum Step 2
	Local index_byte:Int=(x & div32_mask) Shr 5
	Local index_bit:Int=1 Shl (x & mod32_mask)
	If ((flags[index_byte] & index_bit)=0) Then
		Local tmp:Int=x
		Repeat
			tmp:+x
			index_byte=(tmp & div32_mask) Shr 5
			index_bit=1 Shl (tmp & mod32_mask)
			flags[index_byte]:|index_bit
		Until (tmp+x>=max_value)
	EndIf
Next
Print (MilliSecs()-timer)/1000.0+" ..counting primes.."
Local count:Int=1
For Local k:Int=3 To max_value-1 Step 2
	Local index_byte:Int=(k & div32_mask) Shr 5
	Local index_bit:Int=1 Shl (k & mod32_mask)
	Rem
	If ((flags[index_byte] & index_bit)=0) Then
		count:+1
	EndIf
	End Rem
	count:+((flags[index_byte] & index_bit)=0)
Next

'finish
timer=MilliSecs()-timer
Print (timer/1000.0)+" seconds.. "+count+" primes found."

End

A real mess, i tend to not delete any code i write when "scribbling" like this, instead putting comments all over the place :P

Anyways, nates latest sieve inspired me to write this one, the speed increase is pretty much only from reduced memory overhead. Churns out a million primes in ~0,145 seconds and 10 million in ~4,9 seconds on my core2duo t6400 (2ghz).

EDIT: cleaned up the code
EDIT2: small change suggested by Czar
EDIT3: changed the prime counting code a little, now doing ~0,11s/1m primes for me.

Last edited 2010

Last edited 2010

Last edited 2010

Nate the Great

(Posted 2010) [#40]

Hey zzz, what are the masks for? I dont see how the masks are used to optimize it?

Czar Flavius

(Posted 2010) [#41]

Will making the masks constant improve speed more?

zzz	(Posted 2010) [#42]

I used the div/mod mask to get rid of division. cpus calculate the remainder as a byproduct when doing division so the div/idiv opcodes gets used for either operation, and theyre really slow.. And with a fixed known value its pretty straightforward to use bitwise operators instead.
EDIT: that doesnt really have much impact on the code, at least not in a positive way. the main speed increase comes from reducing the memory footprint by using bits for flags instead of bytes. The extra fiddling is outweighted by the reduced memory overhead by a lot :) (and thats also a big reason to why it isnt much faster on the 10m prime test)

The third 'mask' is just a quickie way to mark of all even numbers, 16 numbers at a time :o

@Czar: that works indeed, got it down to ~0,12seconds/1m primes for me :)

sample code:

SuperStrict

Framework brl.standardio

Const LIMIT:Int=500000000

Const MOD32:Int=%11111
Const MASK32:Int=%11111111111111111111111111100000

Local timer:Int=MilliSecs()
Local k:Int
Local count:Int

For Local i:Int=1 To LIMIT
	k=i Mod 32
	count:+(k=0)
Next

timer=MilliSecs()-timer
Print "standard mod: "+timer/1000.0+" sec, count: "+count

count=0
timer=MilliSecs()

For Local i:Int=1 To LIMIT
	k=i & mod32
	count:+(k=0)
Next

timer=MilliSecs()-timer
Print "bitwise and mod: "+timer/1000.0+" sec, count: "+count

End

Last edited 2010

Last edited 2010

Last edited 2010

Last edited 2010

Last edited 2010

zzz	(Posted 2010) [#43]

Its now doing 1m primes in under 0,1 seconds for me :P, getting steady 0,097/0,099 second results. 10m primes in ~4,5 seconds.

SuperStrict

Framework brl.standardio
Import brl.math

'Local max_value:Int=179424674
Local max_value:Int=15485864

Local flags:Int[max_value/32+1]

Const mod32_mask:Int=%11111
Const div32_mask:Int=%11111111111111111111111111100000
Const even_mask:Int=%01010101010101010101010101010101

Local mask:Int[32]
For Local i:Int=0 To 31
	mask[i]=2^i
Next

'start
Local timer:Int=MilliSecs()

For Local i:Int=0 To flags.length-1
	flags[i]=even_mask
Next
flags[0]:~mask[2]

Local sqrnum:Int=Sqr(max_value)

For Local x:Int=3 To sqrnum Step 2
	Local index_byte:Int=(x & div32_mask) Shr 5
	Local index_bit:Int=1 Shl (x & mod32_mask)
	If ((flags[index_byte] & index_bit)=0) Then
		Local tmp:Int=x
		Repeat
			tmp:+x
			index_byte=(tmp & div32_mask) Shr 5
			index_bit=1 Shl (tmp & mod32_mask)
			flags[index_byte]:|index_bit
		Until (tmp+x>=max_value)
	EndIf
Next
Print (MilliSecs()-timer)/1000.0+" ..counting primes.."

Local count:Int=0
Rem
For Local k:Int=3 To max_value-1 Step 2
	Local index_byte:Int=(k & div32_mask) Shr 5
	Local index_bit:Int=1 Shl (k & mod32_mask)
	'If ((flags[index_byte] & index_bit)=0) Then
	'	count:+1
	'EndIf
	count:+((flags[index_byte] & index_bit)=0)
Next
End Rem

'count primes
For Local ii:Int=0 To max_value/32-1
	
	Local b:Int=flags[ii]
	
	'unrolled loop..
	
	count:+(b & (1 Shl 1))=0
	count:+(b & (1 Shl 3))=0
	count:+(b & (1 Shl 5))=0
	count:+(b & (1 Shl 7))=0
	count:+(b & (1 Shl 9))=0
	count:+(b & (1 Shl 11))=0
	count:+(b & (1 Shl 13))=0
	count:+(b & (1 Shl 15))=0
	count:+(b & (1 Shl 17))=0
	count:+(b & (1 Shl 19))=0
	count:+(b & (1 Shl 21))=0
	count:+(b & (1 Shl 23))=0
	count:+(b & (1 Shl 25))=0
	count:+(b & (1 Shl 27))=0
	count:+(b & (1 Shl 29))=0
	count:+(b & (1 Shl 31))=0
	
Next
' check the last few candidates
Local m:Int=max_value/32
Local diff:Int=max_value-m*32

Local b:Int=flags[m]
For Local i:Int=1 To diff Step 2
	count:+(b & mask[i])=0
	'Print m*32+i
Next

'finish
timer=MilliSecs()-timer
Print (timer/1000.0)+" seconds.. "+count+" primes found."

Local primes:Int[count]
Local prime_index:Int=0


For Local ii:Int=0 To max_value/32
	For Local b:Int=0 To 31
		If ((flags[ii] & mask[b])=0) Then
			If (prime_index<count) Then
				primes[prime_index]=ii*32+b
				If (primes[prime_index]>1) Then
					prime_index:+1
				EndIf
			EndIf
		EndIf
	Next
Next

Print "done gathering primes.."
Print "getting ref primes.."

Local ref_primes:Int[]=primes_01(count)

For Local i:Int=0 To ref_primes.length-1
	If (ref_primes[i]<>primes[i]) Then
		Print ref_primes[i]+":"+primes[i]
		Exit
	EndIf
Next

Print "Done!"

End






Function primes_01:Int[](primes_count:Int)
	Local PR:Int=3
	Local primelist:Int[primes_count]
	Local prcnt:Int=1
	primelist[0]=2
	Local n:Int=0
	Local flag:Byte=False
	Repeat
		n=0
		flag=False
		While primelist[n]*primelist[n]<=PR
			If (PR/(primelist[n]*1.0))=Int(PR/primelist[n]) Then
				flag=True
				Exit
			EndIf
			n:+1
		Wend
		If Not flag
			primelist[prcnt]=pr
			prcnt:+1
		EndIf
		PR:+2
	Until prcnt=primes_count
	Return primelist
End Function

Last edited 2010

Shortwind

(Posted 2010) [#44]

This is very exciting watching this thread over the last few days.

zzz, Nate, super fun!

Still analyzing your last bit of code zzz. A pretty original line of thought there!

:D

Nate the Great

(Posted 2010) [#45]

wow zzz, 3.9 seconds, 10,000,000 primes found. I believe we have reached a sort of road block for efficiency of counting primes but zzz may prove me wrong. (it may be faster in asm hint hint haha) Anyway, I hadn't touched bmax for a while so this helped me really get into it and pick it back up. thanks!

zzz	(Posted 2010) [#46]

@Shortwind, i realized it was the way to go when i messed around with nates latest creation. changing the flags array to ints instead of bytes made it way slower, so figured it would work the other way around too.

@Nate, try doing 1m primes instead, it will be way faster since the 10m prime test re-introduces the memory issue i tried to get rid of.. :)

I dont think ill mess around much more with this. Assembly wont really help since its mostly a memory issue now, the algo touches memory all over the place over and over again so it is most likely racking up a huge number of chache/page mispredictions. A single L2 misprediction can cost over 200 cycles accordingly to google, and for reference i think an ordinary integer addition is a one cycle ordeal.

One could probably rewrite the code to work around this in some way, but i cant really come up with any ideas atm.

zzz	(Posted 2010) [#47]

Expanded a bit on the preliminary culling, the code gives a pretty solid average of 0,075 seconds / 1m primes :D

SuperStrict

Framework brl.standardio
Import brl.math
Import brl.retro

'Local max_value:Int=179424674
Local max_value:Int=15485864

Local flags:Int[max_value/32+1]

Const mod32_mask:Int=%11111
Const div32_mask:Int=%11111111111111111111111111100000
Const even_mask:Int=%01010101010101010101010101010101

Local mask:Int[32]
For Local i:Int=0 To 31
	mask[i]=2^i
Next

'start
Local timer:Int=MilliSecs()

'/////////////////////////////////////////
'build a premask array
Local premask_size:Int=10
Repeat
	If (premask_size Mod 2 = 0) And (premask_size Mod 3 = 0) And (premask_size Mod 5 = 0) And (premask_size Mod 7 = 0) Then
		If (premask_size Mod 32 = 0) Then
			Exit
		EndIf
	EndIf
	premask_size:+10
Forever
'build array
Local premask_array:Int[premask_size/32]
For Local index:Int=0 To premask_array.length-1
	Local b:Int=0
	For Local i:Int=0 To 31
		If ((index*32+i) Mod 2 = 0) Then
			b :| Int(2^i)
		EndIf
		If ((index*32+i) Mod 3 = 0) Then
			b :| Int(2^i)
		EndIf
		If ((index*32+i) Mod 5 = 0) Then
			b :| Int(2^i)
		EndIf
		If ((index*32+i) Mod 7 = 0) Then
			b :| Int(2^i)
		EndIf
	Next
	premask_array[index]=b
Next
'debug
Rem
For Local index:Int=0 To premask_array.length-1
	Print Bin(premask_array[index])
Next
End Rem
'apply mask
Rem
For Local i:Int=0 To flags.length-1
	flags[i]=even_mask
Next
flags[0]:~mask[2]
End Rem

'apply premask

For Local index:Int=0 To flags.length-1
	flags[index]=premask_array[index Mod premask_array.length]
Next

flags[0]:~mask[2]
flags[0]:~mask[3]
flags[0]:~mask[5]
flags[0]:~mask[7]

'/////////////////////////////////////////

Local sqrnum:Int=Sqr(max_value)

For Local x:Int=9 To sqrnum Step 2
	Local index_byte:Int=(x & div32_mask) Shr 5
	Local index_bit:Int=1 Shl (x & mod32_mask)
	If ((flags[index_byte] & index_bit)=0) Then
		Local tmp:Int=x
		Repeat
			tmp:+x
			index_byte=(tmp & div32_mask) Shr 5
			index_bit=1 Shl (tmp & mod32_mask)
			flags[index_byte]:|index_bit
		Until (tmp+x>=max_value)
	EndIf
Next
Print (MilliSecs()-timer)/1000.0+" ..counting primes.."

Local count:Int=0
Rem
For Local k:Int=3 To max_value-1 Step 2
	Local index_byte:Int=(k & div32_mask) Shr 5
	Local index_bit:Int=1 Shl (k & mod32_mask)
	'If ((flags[index_byte] & index_bit)=0) Then
	'	count:+1
	'EndIf
	count:+((flags[index_byte] & index_bit)=0)
Next
End Rem

'count primes
For Local ii:Int=0 To max_value/32-1
	
	Local b:Int=flags[ii]
	
	'unrolled loop..
	
	count:+(b & (1 Shl 1))=0
	count:+(b & (1 Shl 3))=0
	count:+(b & (1 Shl 5))=0
	count:+(b & (1 Shl 7))=0
	count:+(b & (1 Shl 9))=0
	count:+(b & (1 Shl 11))=0
	count:+(b & (1 Shl 13))=0
	count:+(b & (1 Shl 15))=0
	count:+(b & (1 Shl 17))=0
	count:+(b & (1 Shl 19))=0
	count:+(b & (1 Shl 21))=0
	count:+(b & (1 Shl 23))=0
	count:+(b & (1 Shl 25))=0
	count:+(b & (1 Shl 27))=0
	count:+(b & (1 Shl 29))=0
	count:+(b & (1 Shl 31))=0
	
Next
' check the last few candidates
Local m:Int=max_value/32
Local diff:Int=max_value-m*32

Local b:Int=flags[m]
For Local i:Int=1 To diff Step 2
	count:+(b & mask[i])=0
	'Print m*32+i
Next

'finish
timer=MilliSecs()-timer
Print (timer/1000.0)+" seconds.. "+count+" primes found."

Local primes:Int[count]
Local prime_index:Int=0


For Local ii:Int=0 To max_value/32
	For Local b:Int=0 To 31
		If ((flags[ii] & mask[b])=0) Then
			If (prime_index<count) Then
				primes[prime_index]=ii*32+b
				If (primes[prime_index]>1) Then
					prime_index:+1
				EndIf
			EndIf
		EndIf
	Next
Next

Print "done gathering primes.."
Print "getting ref primes.."

Local ref_primes:Int[]=primes_01(count)

For Local i:Int=0 To ref_primes.length-1
	If (ref_primes[i]<>primes[i]) Then
		Print ref_primes[i]+":"+primes[i]
		Exit
	EndIf
Next

Print "Done!"

End






Function primes_01:Int[](primes_count:Int)
	Local PR:Int=3
	Local primelist:Int[primes_count]
	Local prcnt:Int=1
	primelist[0]=2
	Local n:Int=0
	Local flag:Byte=False
	Repeat
		n=0
		flag=False
		While primelist[n]*primelist[n]<=PR
			If (PR/(primelist[n]*1.0))=Int(PR/primelist[n]) Then
				flag=True
				Exit
			EndIf
			n:+1
		Wend
		If Not flag
			primelist[prcnt]=pr
			prcnt:+1
		EndIf
		PR:+2
	Until prcnt=primes_count
	Return primelist
End Function

Last edited 2010

xcessive

(Posted 2010) [#48]

With the above code: 0.009200... seconds.

i7 extreme OC'd :D

Imphenzia

(Posted 2010) [#49]

0.128000006 seconds
1000000 primes found!

i7 950 @ original 3.06GHz (non-debug)

Grey Alien

(Posted 2010) [#50]

This thread is fantastic, well done all! Fascinating.

Grisu

(Posted 2010) [#51]

0.0540000014
1000000 primes found

i7 920 @ 2.98 GHZ (non-debug)

zzz	(Posted 2010) [#52]

Had a bunch of school-related stuff to finish today, so figured the sensible thing to do was to write yet another sieve to post here. Same as the last one pretty much. Continuing on the memory footprint philosophy i threw out even numbers completely.

EDIT: on my computer it does 30m primes in just under 10 seconds, where 10m takes just over 2 seconds, so it scales fairly well this time. Set limit to 573259392 if you want to try :)

EDIT2: updated the prime number counting code so now it should run considerably faster on the 1m prime test too. For me its slightly more than twice as fast as the last one i posted.


	SuperStrict
	
	Framework BRL.StandardIO
	Import BRL.Math
	Import BRL.Blitz
	
	'Const LIMIT:Int=1000
	'Const LIMIT:Int=982451654 '50m primes
	'Const LIMIT:Int=573259392 '30m primes
	'Const LIMIT:Int=179424674 '10m primes
	Const LIMIT:Int=15485864 '1m primes
	'Const LIMIT:Int=10000
	Const WORD_SIZE:Int=32
	
	Global wheel_mask:Int[]
	Global wheel_values:Int[]=[2,3,5,7,11,13]
	Global wheel_size:Int
	
	
	Global flags:Int[(LIMIT/2)/32+1]
	
	Const mod32:Int=%11111
	Const div32:Int=%11111111111111111111111111100000
	
	' ///////////////////////////////////////////////////////////////////////////
	
	InitWheel()
	
	WriteStdout "wheel bitmask initialized. "
	WriteStdout "size: "+(wheel_mask.length*WORD_SIZE)
	WriteStdout ", highest value: "+wheel_values[wheel_values.length-1]
	WriteStdout "~n"
	
	Delay 1000
	
	Print "Sieving values.."
	Local timer:Int=MilliSecs()
	
	For Local i:Int=0 To flags.length-1
		flags[i]=wheel_mask[i Mod wheel_mask.length]
	Next
	
	Local sqr_limit:Int=Sqr(LIMIT)/2+1
	For Local i:Int=1 To sqr_limit 'Sqr(LIMIT)/2+1
		Local index:Int=(i & div32) Shr 5
		Local bit:Int=1 Shl (i & mod32)
		If (flags[index] & bit=0) Then
			Local v:Int=1+(i*2)
			Local temp:Int=i+v
			Repeat
				Local index:Int=(temp & div32) Shr 5
				Local bit:Int=1 Shl (temp & mod32)
				flags[index]:|bit
				temp:+v			
			Until (temp>=(LIMIT/2))
		EndIf
	Next
	
	WriteStdout "Done! Took "+((MilliSecs()-timer)/1000.0)
	WriteStdout " seconds.. Counting primes.."
	WriteStdout "~n"
	
	' ----------------------------------------------------
		
	Local primes_count:Int=wheel_values.length
	For Local i:Int=0 To (LIMIT/2)/32-1
		Local b:Int=flags[i]
		Local count:Int=0
		'8
		count:+(b & (1 Shl 0) = 0)
		count:+(b & (1 Shl 1) = 0)
		count:+(b & (1 Shl 2) = 0)
		count:+(b & (1 Shl 3) = 0)
		count:+(b & (1 Shl 4) = 0)
		count:+(b & (1 Shl 5) = 0)
		count:+(b & (1 Shl 6) = 0)
		count:+(b & (1 Shl 7) = 0)
		'8
		count:+(b & (1 Shl 8) = 0)
		count:+(b & (1 Shl 9) = 0)
		count:+(b & (1 Shl 10) = 0)
		count:+(b & (1 Shl 11) = 0)
		count:+(b & (1 Shl 12) = 0)
		count:+(b & (1 Shl 13) = 0)
		count:+(b & (1 Shl 14) = 0)
		count:+(b & (1 Shl 15) = 0)
		'8
		count:+(b & (1 Shl 16) = 0)
		count:+(b & (1 Shl 17) = 0)
		count:+(b & (1 Shl 18) = 0)
		count:+(b & (1 Shl 19) = 0)
		count:+(b & (1 Shl 20) = 0)
		count:+(b & (1 Shl 21) = 0)
		count:+(b & (1 Shl 22) = 0)
		count:+(b & (1 Shl 23) = 0)
		'8
		count:+(b & (1 Shl 24) = 0)
		count:+(b & (1 Shl 25) = 0)
		count:+(b & (1 Shl 26) = 0)
		count:+(b & (1 Shl 27) = 0)
		count:+(b & (1 Shl 28) = 0)
		count:+(b & (1 Shl 29) = 0)
		count:+(b & (1 Shl 30) = 0)
		count:+(b & (1 Shl 31) = 0)
		' done
		primes_count:+count
	Next
	
	Local m:Int=(LIMIT/2)/32
	Local diff:Int=LIMIT-m*64
	diff:/2
	For Local i:Int=0 To diff-2
		primes_count:+(flags[m] & Int(2^i) = 0)
	Next		
	
	timer=MilliSecs()-timer
	
	Local primes:Int[primes_count]
	primes_count=0
	
	For Local i:Int=0 To wheel_values.length-1
		primes[primes_count]=wheel_values[i]
		primes_count:+1
	Next
	
	For Local i:Int=1 To LIMIT/2-1
		Local index:Int=(i & div32) Shr 5
		Local bit:Int=1 Shl (i & mod32)
		If (flags[index] & bit = 0) Then
			Local v:Int=1+(i*2)
			primes[primes_count]=v
			primes_count:+1
		EndIf	
	Next
	
	WriteStdout "Done! Took "+(timer/1000.0)+" seconds to find "
	WriteStdout "and count "+primes_count+" primes."
	
	Print ""
	Print "Validating result.."
	
	Local refprimes:Int[]=ReferencePrimes(primes_count)
	For Local i:Int=0 To primes_count-1
		'Print refprimes[i]
		If (refprimes[i]<>primes[i]) Then
			Print "!!"
			Print refprimes[i]+":"+primes[i]
			Exit
		EndIf
	Next
	
	Print "Done!"
	
	End
	
	' ///////////////////////////////////////////////////////////////////////////
	
	Function InitWheel()
		Local gcd:Int=wheel_values[0]
		For Local i:Int=1 To wheel_values.length-1
			gcd:*wheel_values[i]
		Next
		Local max_value:Int=gcd*32
		wheel_mask=New Int[gcd]
		For Local i:Int=0 To max_value-1
			Local value:Int=1+(i*2)
			For Local j:Int=0 To wheel_values.length-1
				If (value Mod wheel_values[j] = 0) Then
					Local index:Int=(i & div32) Shr 5
					Local bit:Int=1 Shl (i & mod32)
					wheel_mask[index]:|bit
				EndIf
			Next
		Next
		Return 
	End Function
	
	Function ReferencePrimes:Int[](primes_count:Int)
		Local PR:Int=3
		Local primelist:Int[primes_count]
		Local prcnt:Int=1
		primelist[0]=2
		Local n:Int=0
		Local flag:Byte=False
		Repeat
			n=0
			flag=False
			While primelist[n]*primelist[n]<=PR
				If (PR/(primelist[n]*1.0))=Int(PR/primelist[n]) Then
					flag=True
					Exit
				EndIf
				n:+1
			Wend
			If Not flag
				primelist[prcnt]=pr
				prcnt:+1
			EndIf
			PR:+2
		Until prcnt=primes_count
		Return primelist
	End Function

Last edited 2010

Last edited 2010

zzz	(Posted 2010) [#53]

Am I the only one who still have some interest in this? Nate, Jesse, Shortwind? :) You guys probably have better things to do (im studying ~2days a week atm so i dont), just wondering if I should continue to post updates etc.

I could try and clarify whats going on in the code a bit if anyone wants me to. I havent done that much more than obfuscating the code from the original sieve so its not hard to understand the logic in any way.

Last edited 2010

Shortwind

(Posted 2010) [#54]

Sorry, zzz, haven't lost interest, just been busy.

I believe the code we have in this thread is very usefull, and a good learning experience. What I'd suggest now is that you modify your code to quickly, with reasonable memory requirements, find the X'th prime number, the next prime, and the previous prime. (If you feel so inclined.)

In other words, if I tell your program that I want the 1,578th prime number, how fast can your program return that number?

If I give your program the number 95,456, how fast can it tell me the next or previous prime number?

And finally, (you'll have to start a new thread), what is the fastest method in BlitzMax for "factoring" a number into it's prime elements?

:D

zzz	(Posted 2010) [#55]

That functionality would come as a bonus if i could get a segmented sieve to work properly (and fast enough). Youd still have a bit of overhead but it would be able to quite quickly calculate enough primes to be able to answer those questions :)

Ive been thinking about upping the performance of the sieve by packing values.. The last piece of code excludes even numbers completely, so:

a) I get a smaller memory footprint since the algorithm needs one bit per two values in the 1-limit range when culling prime candidates.

b) The culling code can be made much more efficient. Initially i hade the problem of prime multiples not existing in the flag array. For example the first multiple of 5 is 10, which is even so it doesnt exist as far as the sieve is concerned. This one was pretty easy to figure out and had the nice side-effect of making the culling loop twice as fast since i could skip every second multiple completely.

Ive been tinkering a bit with excluding multiples of 2,3,5 and 7 in the same way, but this makes the 1-limit range non-linear as opposed to just excluding odd numbers.. It would give the flag array a value density of around 440% (ie, 4,4values per bit) tho if i figured it out, which would give a very nice performance boost :D

After i nail those two ill get a factorization thread going :P

Last edited 2010

zzz	(Posted 2010) [#56]

Proof of concept for a segmentet sieve. Seems to work properly so time to trim the code a bit i guess. Also opens up for multithreading which is always nice


	SuperStrict

	Framework BRL.StandardIO
	Import BRL.LinkedList
	
	Const MEMORY_LIMIT:Int=512 'kb
	Const SEGMENT_SIZE:Int=MEMORY_LIMIT*1000
	
	Global flags:Int[SEGMENT_SIZE]
	Global cull_values:TList=New TList
	Global current_segment:Int=0
	
	Const PRIME_INCREMENT:Int=64000 '64kb a time
	Global primes:Int[PRIME_INCREMENT]
	Global primes_count:Int=0
		
	'#####################################################################
		
	Const TARGET:Int=2000000
	Local time:Int=MilliSecs()
	
	Repeat
		Sieve()
	Until (primes_count>=TARGET)
	
	time=MilliSecs()-time
	Print "The "+TARGET+"th prime is "+primes[TARGET-1]+". ("+(time/1000.0)+" s, "+current_segment+" segments)"
	
	End
	
	'#####################################################################
		
	Type TCullValue
		
		Field _base:Int
		Field _current:Int
		
		Function Create:TCullValue(p_base:Int, p_current:Int)
			Local this:TCullValue=New TCullValue
			this._base=p_base
			this._current=p_current
			Return this
		End Function
	
	End Type
	
	Function Sieve:Int()
		For Local i:Int=0 To flags.length-1
			flags[i]=0
		Next
		If (current_segment=0) Then
			''''''''''''''''''''''''''''''''''''''''''''''''''
			'first segment, do the ordinary sieve stuff here
			'sieve
			Local value:Int=2
			For Local i:Int=2 To SEGMENT_SIZE-1
				If (flags[i]=0) Then
					Local temp:Int=value+value
					If (temp<SEGMENT_SIZE) Then
						Repeat
							flags[temp]=1
							temp:+value
						Until (temp>=SEGMENT_SIZE)
					EndIf
					cull_values.AddLast(TCullValue.Create(value, temp))
				EndIf
				value:+1	
			Next
			'gather primes
			For Local i:Int=2 To SEGMENT_SIZE-1
				If (flags[i]=0) Then
					If (primes_count>=primes.length) Then
						primes=primes[..primes.length+PRIME_INCREMENT]
					EndIf
					primes[primes_count]=i
					primes_count:+1
				EndIf
			Next
		Else
			''''''''''''''''''''''''''''''''''''''''''''''''''
			'arbitrary segment, cull with the values stored
			'from previous segments and then get the new
			'primes.
			Local value:Int=current_segment*SEGMENT_SIZE
			Local limit:Int=value+SEGMENT_SIZE
			'cull with stored values
			For Local o:TCullValue=EachIn cull_values
				While (o._current<limit)
					If (o._current>=value) Then
						flags[o._current-value]=1
					EndIf
					o._current:+o._base
				Wend
			Next
			'since no multiples can be found in the segment, count
			'new primes directly
			For Local i:Int=0 To SEGMENT_SIZE-1
				If (flags[i]=0) Then
					'add to cull values
					Local prime:Int=value+i
					cull_values.AddLast(TCullValue.Create(prime, prime+prime)) 
					'store prime
					If (primes_count>=primes.length) Then
						primes=primes[..primes.length+PRIME_INCREMENT]
					EndIf
					primes[primes_count]=prime
					primes_count:+1
				EndIf
			Next
		EndIf
		current_segment:+1
	End Function

	End

zzz	(Posted 2010) [#57]

Messed around a bit today with this and managed to get a segmented sieve going. The code needs some major cleaning up which i dont feel like doing now so heres some pretty numbers to look at instead:

edit: changed my mind, heres the code instead. Bit of warning that while it churns out primes pretty fast (100000000th prime took 12sec for me) it gobbles up lots of memory.. The TARGET constant is what prime the code will find so mess around with that.

edit2: small optimization. its currently doing 10m primes in just over 0,7 seconds for me :P

	SuperStrict

	Framework BRL.StandardIO
	Import BRL.LinkedList
	Import BRL.Math
	
	'GCD of 2,3,5,7 is 210, so segment size should be a multiple of this
	'expand this to 2,3,5,7,11,13 at least (GCD:30030, which is manageable)
	
	Const MEMORY_LIMIT:Int=500 'kb
	'Const SEGMENT_SIZE:Int=(MEMORY_LIMIT*210*4)
	Const SEGMENT_SIZE:Int=(2*3*5*7*11*13)
	
	Global flags:Byte[SEGMENT_SIZE]
	Global cull_values:TList=New TList
	Global current_segment:Int=0
	Global wheel:Byte[SEGMENT_SIZE]
	
	Global primes:TList=New TList
	Global primes_count:Int=0
	
	Global timer0:Int=0
	Global timer1:Int=0
	Global timer2:Int=0
	
	Global mem_usage:Int=0 'bytes
	mem_usage:+(SEGMENT_SIZE*2)
		
	'init wheel
	For Local i:Int=0 To SEGMENT_SIZE-1
		Local v:Int=1+(2*i)
		If (v Mod 2 = 0) Or (v Mod 3 = 0) Or (v Mod 5 = 0) Or (v Mod 7 = 0) Or (v Mod 11 = 0) Or (v Mod 13 = 0) Then
			wheel[i]=1
		EndIf
	Next
		
	'#####################################################################
		
	Delay 500
		
	Const TARGET:Int=10000000
	Local time:Int=MilliSecs()
	
	Repeat
		Sieve()
	Until (primes_count>=TARGET)

	Local targetprime:Int=0
	Local count:Int=6
	
	For Local o:TCullValues=EachIn primes
		For Local i:Int=0 To o._count-1
			count:+1
			If (count=TARGET) Then
				targetprime=o._primes[i,0]
				Exit
			EndIf
		Next
	Next
	
	time=MilliSecs()-time
	Print "The "+TARGET+"th prime is "+targetprime+". ("+(time/1000.0)+" s, "+current_segment+" segments)"
	Print (timer0/1000.0)
	Print (timer1/1000.0)
	Print (timer2/1000.0)
	Print ""
	Print "Total memory usage: "+(mem_usage/1000)+" kb"
	Print ""

	
	End
	
	'#####################################################################
		
	Type TCullValue
		
		Field _base:Int
		Field _current:Int
		
		Function Create:TCullValue(p_base:Int, p_current:Int)
			Local this:TCullValue=New TCullValue
			this._base=p_base
			this._current=p_current
			Return this
		End Function
	
	End Type
	
	Type TCullValues
	
		Global prime_ratio:Int=1 'expected prime ratio 
	
		Field _base:Int
		Field _limit:Int
		Field _count:Int=0
		Field _primes:Int[,]
		
		Function Create:TCullValues(p_base:Int, p_limit:Int)
			Local this:TCullValues=New TCullValues
			this._base=p_base
			this._limit=p_limit
			this._count=0
			this._primes=New Int[(p_limit-p_base)/TCullValues.prime_ratio,2]
			mem_usage:+this._primes.length*8
			Return this
		End Function
		
	End Type
	
	Function Sieve:Int()
	
		Rem
		For Local i:Int=0 To flags.length-1
			flags[i]=wheel[i]
			'flags[i]=0
		Next
		End Rem
		
		MemCopy(Varptr flags[0], Varptr wheel[0], SEGMENT_SIZE)
		
		Local cull:TCullValues=TCullValues.Create(current_segment*SEGMENT_SIZE, current_segment*SEGMENT_SIZE+SEGMENT_SIZE)
		cull._count=0
		If (current_segment=0) Then
			Local t0:Int=MilliSecs()
			''''''''''''''''''''''''''''''''''''''''''''''''''
			'first segment, do the ordinary sieve stuff here
			primes_count:+4 '2,3,5,7
			'Local tt:Int=MilliSecs()
			For Local i:Int=1 To SEGMENT_SIZE-1
				Local value:Int=1+(i*2)
				If (flags[i]=0) Then
					Local temp:Int=i+value 'temp is an index, not value!
					If (temp<SEGMENT_SIZE) Then
						Repeat
							flags[temp]=1
							temp:+value
						Until (temp>=SEGMENT_SIZE)
					EndIf
					'cull_values.AddLast(TCullValue.Create(value, temp))
					cull._primes[cull._count,0]=value
					cull._primes[cull._count,1]=temp
					cull._count:+1
					primes_count:+1
				EndIf
			Next
			'tt=MilliSecs()-tt
			'Print "sieve: "+(tt/1000.0)
			t0=MilliSecs()-t0
			timer0:+t0
		Else
			''''''''''''''''''''''''''''''''''''''''''''''''''
			'arbitrary segment, cull with the values stored
			'from previous segments and then get the new
			'primes.
			'Local tt:Int=MilliSecs()
			Local t:Int
			t=MilliSecs()
			Local value:Int=current_segment*SEGMENT_SIZE
			Local limit:Int=value+SEGMENT_SIZE

				
			'cull with stored values..		
			For Local o:TCullValues=EachIn primes
				Local done:Int=False
				
				For Local i:Int=0 To o._count-1
				
					Local base:Int=o._primes[i,0]
					Local cur:Int=o._primes[i,1]
					
					If (base*base>=limit*2) Then
						done=True
						Exit
					EndIf
					
					While (cur<value)
						cur:+base
					Wend
					
					While (cur<limit)
						flags[cur-value]=1
						cur:+base
					Wend
					
					o._primes[i,0]=base
					o._primes[i,1]=cur
					
				Next
				
				If (done=True) Then
					Exit
				EndIf
				
			Next
			
			
			
			t=MilliSecs()-t
			timer1:+t
			t=MilliSecs()
			
			'since no multiples can be found in the segment new primes can
			'be fetched right away
			For Local i:Int=0 To SEGMENT_SIZE-1
				If (flags[i]=0) Then
					'add to cull values
					Local prime:Int=1+2*(value+i)
					'Print prime
					'cull_values.AddLast(TCullValue.Create(prime, value+i)) 
					cull._primes[cull._count,0]=prime
					cull._primes[cull._count,1]=value+i
					cull._count:+1
					'Print cull._count
					primes_count:+1
				EndIf
			Next
			t=MilliSecs()-t
			timer2:+t
			'tt=MilliSecs()-tt
			'Print "sieve: "+(tt/1000.0)+" range: "+limit+"  primes: "+primes_count
		EndIf
		' update prime ratio and add current cull values to list
		TCullValues.prime_ratio=((cull._limit-cull._base)/cull._count)-1
		primes.AddLast(cull)
		current_segment:+1
	End Function

	End

Last edited 2010

Last edited 2010

Last edited 2010

zzz	(Posted 2011) [#58]

A bit of thread necromancy here.

I started refactoring the above code into something useful but lost interest soon afterwards.. Remembered the topic and figured id post it here anyways if someone is interested.

Feel free to try it and post some execution times :)


	Rem
	this version uses ~520mb of memory for 100m primes, while the
	bit array version uses ~140mb. The speed difference is quite big
	tho, 8s compared to ~15s, almost twice as fast..
	End Rem
	
	SuperStrict
	
	Framework brl.standardio
	Import brl.linkedlist
	Import brl.blitz
	Import brl.random
	Import brl.math

	Rem
	End Rem

	Print ""

	Local t:Int=MilliSecs()

	For Local i:Int=1 To 20
	
		Local index:Int=i*5000000
		Local prime:Int=TSoE.GetPrime(index)
		't=MilliSecs()-t
		Print "The "+index+"th prime is "+prime+". ( total running time: "+((MilliSecs()-t)/1000.0)+" seconds )"
		
	Next
	
	Input()
	
	End

	Rem
	######################################################################
	End Rem

	Type TSoE
	
		Rem
		End Rem
		
		Const SEGMENTSIZE:Int=2*3*5*7*11*13
		
		Global Init:Int=False
		
		Global arrSkipValues:Int[]=[2,3,5,7,11,13]
		Global arrWheel:Byte[SEGMENTSIZE]
	
		Global listCullValues:TList
		Global Primes:Int
		Global BiggestPrime:Int
		Global CurrentSegment:Int
		
		Rem
		End Rem
		
		Function Initialize:Int()
			If (TSoe.Init=False) Then
				'init wheel
				For Local i:Int=0 To TSoE.SEGMENTSIZE-1
					Local v:Int=1+(i Shl 1)
					'WriteStdout v
					For Local j:Int=0 To TSoe.arrSkipValues.length-1
						If (v Mod TSoE.arrSkipValues[j] = 0) Then
							'WriteStdout "!"
							TSoE.arrWheel[i]=1
							Continue
						EndIf
					Next
					'WriteStdout "~n"
				Next
				'init misc
				TSoE.listCullValues=New TList
				TSoE.Primes=0
				TSoE.BiggestPrime=0
				TSoE.CurrentSegment=0
				'done
				TSoE.Init=True
				Return True
			EndIf
			'no init done if object was already initialized
			Return False
		End Function
		
		Rem
		End Rem
		
		Function Sieve:TSoECull(segment:Int)
			'create cull object and flags
			Local cull:TSoECull=TSoECull.Create(segment*TSoE.SEGMENTSIZE, TSoE.SEGMENTSIZE)
			Local flags:Byte[TSoE.SEGMENTSIZE]
			'copy wheel to flags
			MemCopy(Varptr flags[0], Varptr TSoE.arrWheel[0], TSoE.SEGMENTSIZE)
			'sieve
			If (segment=0) Then
				'----------------------------------------------------
				'first segment
				TSoE.Primes:+TSoE.arrSkipValues.length
				For Local i:Int=1 To TSoE.SEGMENTSIZE-1
					Local value:Int=1+(i Shl 1)
					If (flags[i]=0) Then
						Local temp:Int=i+value
						If (temp<TSoE.SEGMENTSIZE) Then
							Repeat
								flags[temp]=1
								temp:+value
							Until (temp>=TSoE.SEGMENTSIZE)
						EndIf
						cull._flags[cull._primes]=value
						cull._primes:+1
					EndIf
				Next
			Else
				'----------------------------------------------------
				'arbitrary segment
				
				Local value:Int=segment*TSoE.SEGMENTSIZE
				Local limit:Int=value+TSoE.SEGMENTSIZE
				Local limit2:Int=limit Shl 1
				'limit2=Sqr(limit2)
				
				Local o:TSoECull
				For o=EachIn TSoE.listCullValues
					Exit
				Next
				
				'For Local o:TSoECull=EachIn TSoE.listCullValues
				
					'Local done:Int=False
					For Local i:Int=0 To o._primes-1
						Local base:Int=o._flags[i]
						If (base*base>=limit2) Then
						'If (base>limit2) Then
							Exit
						EndIf
						Local cur:Int=base Shr 1
						Local t:Int=value-cur
						cur:+(t-(t Mod base))
						cur:-value
						If (cur<0) Then
							cur:+base
						EndIf
						While (cur<TSoE.SEGMENTSIZE)
							flags[cur]=1
							cur:+base
						Wend
					Next
			
					'If (done=True) Then
					'	Exit
					'EndIf
				'Next
				
				'considering the limiting values above the code is a bit oddly
				'written, but nothing to gain from rewriting it really..
				
				For Local i:Int=0 To TSoE.SEGMENTSIZE-1
					If (flags[i]=0) Then
						Local prime:Int=1+((value+i) Shl 1)
						cull._flags[cull._primes]=prime
						cull._primes:+1
					EndIf
				Next
				
			EndIf
			'done
			Return cull
		End Function
		
		Rem
		End Rem
		
		Function GetPrime:Int(index:Int, threads:Int=1)
			TSoE.Initialize()
			While (TSoE.Primes<index)
				Local cull:TSoECull=TSoE.Sieve(TSoE.CurrentSegment)
				TSoE.Primes:+cull._primes
				TSoE.listCullValues.AddLast(cull)
				TSoE.CurrentSegment:+1
				TSoECull.prime_ratio=(cull._range/cull._primes)-1
			Wend
			'prime can now be found in cull values.. find and
			'return it.
			Local prime:Int=0
			Local count:Int=TSoE.arrSkipValues.length
			For Local o:TSoECull=EachIn TSoE.listCullValues
				If (count+o._primes<index) Then
					count:+o._primes
				Else
					For Local i:Int=0 To o._primes-1
						count:+1
						If (count=index) Then
							Return o._flags[i]
						EndIf
					Next
				EndIf
			Next
			'shouldnt get here..
			Return -1
		End Function
	
	End Type

	Rem
	######################################################################
	End Rem

	Type TSoECull

		Global prime_ratio:Int=1
		
		Field _base:Int
		Field _range:Int
		Field _primes:Int
		Field _flags:Int[]
		
		Function Create:TSoECull(base:Int, range:Int)
			Local this:TSoECull=New TSoECull
			this._base=base
			this._range=range
			this._primes=0
			this._flags=New Int[this._range/TSoECull.prime_ratio+1]
			Return this
		End Function

	End Type

Peter

(Posted 2011) [#59]

The 100000000th prime is 2038074743. ( total running time: 8.94799995 seconds )

from my company notebook (Intel Core 2 Duo, 2GHz, 2 GB RAM)

The 100000000th prime is 2038074743. ( total running time: 7.98600006 seconds )

from my Intel Core 2 6400 2,13 GHz, 2 GB RAM (growing old now)

-- Peter

Last edited 2011

_JIM	(Posted 2011) [#60]

The 100000000th prime is 2038074743. ( total running time: 5.26900005 seconds )

work pc: i3 3.06GHz

ima747

(Posted 2011) [#61]

The 100000000th prime is 2038074743. ( total running time: 6.32000017 seconds )

Mac Pro (mid 2010) 8 core (2x4 2.4ghz Xenon) with a few things in the background (shouldn't make much of a difference, extra cores and such)

Wonder if memory/bus speed would start having an impact with that high a load generated in such a short time...

virtlands

(Posted 2012) [#62]

Hi Nate the Grate!

It works just great on my laptop:

1.48599994 seconds
1000000 primes found!

My computer type is:
Intel Core 2 CPU T5200 @ 1.6 Ghz speed
(Centrino Duo) HP Pavilion DV6245us
4GB memory, 32-bit, Windows 7

PowerPC603

(Posted 2012) [#63]

The 100000000th prime is 2038074743. ( total running time: 3.68799996 seconds )

64-bit Quad Core i5-3450 @3.1GHz
12GB RAM
64-bit Windows 7 Home Premium SP 1
AMD Radeon HD7770 1GB

Using the edited code in first post:

0.194999993 seconds
1000000 primes found!

Last edited 2012