Welcome, guest | Sign In | My Account | Store | Cart

This is a fast prime number list generator using sieve algorithm. This function return a list of prime numbers which <= argument.

Python, 54 lines
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54``` ```def primes(n): if n==2: return  elif n<2: return [] s=range(3,n+1,2) mroot = n ** 0.5 half=(n+1)/2-1 i=0 m=3 while m <= mroot: if s[i]: j=(m*m-3)/2 s[j]=0 while j
##### This program is almost 8 times faster than the other recipe(fast prime number list creator), and it's more than 2 time faster than Alex's following code:

sieve = range(max+1) def dosieve(sieve, max): for i in range(2, math.sqrt(max)+1): if not sieve[i]: continue for j in range(i*i, max+1, i): sieve[j] = None dosieve(sieve, max)

##### primes=map(str,filter(None, sieve[2:]))

one reason it's faster is that is it operate only on odd number. Another reason is n**0.5 is slightly faster than math.sqrt(n) :) Wensheng Wang (author) 18 years, 7 months ago

Oops! The second reason I gave is not correct. If the reviewer see this, please delete that second reason. Also please indent Alex Martelli's code. thanks. Kazuo Moriwaka 18 years, 7 months ago

another reason. Thanks for making progress. I think another reason why your code is faster is old (my) recipe was create lists each steps, but your code just rewrite numbers to 0. Pierre Quentel 18 years, 7 months ago

Improvement. Very nice recipe, thanks. A slight improvement : the function raises exceptions for some odd numbers (primes(13) for instance). It is safer if you change the value of "half" to len(s) Michael Spencer 18 years, 7 months ago

Nice speed, but the indexing is a little opaque... Here's the same algorithm with simpler indexing:

``````def primes2(n):
if nHere's the same algorithm with simpler indexing:

<pre>
def primes2(n):
if n
``````

</pre> Michael Spencer 18 years, 7 months ago

How do you add code in the comments. I've tried surrounding it in pre tags. I give up... Jimmy Hoeks 18 years, 5 months ago

Out by 1? You say the function returns primes that are less than the argument. Did you mean less than or equal to the argument?

In fact, as it stands, the function can return a prime higher than the argument (try it for 3000, for instance). Changing 's=range(3,n+2,2)' to 's=range(3,n+1,2)' would ensure you get all primes up to and including the argument. Paul Miller 18 years, 5 months ago

Seems to be a bug. For some reason, I get an IndexError for certain prime or prime-power values of n. For example, try n=25, n=251. Giuliano Castelli 18 years, 3 months ago

purpose for better implementation. Would it be possible to implement an analogous so efficient algorithm as primes(n), but with two inputs, primes(s,n), where s is the starting number from which the primes begin to be calculate? Instead of beginning always from 2. Thank you! Simon Gomizelj 17 years, 11 months ago

bug fix. as many of you coders have noticed that with this code there are some values that fail, such as 13, 25, 251, and 3000. there is a simple solution to this problem, the line reads

half = (n+1) / 2

half = (n+1) / 2 - 1

this variable seems to be the size of the list considering it should be 1/2 of the original list's size because we are working only with odd numbers. however, the coder forgot to take into account that he skipped the odd value of 1 and started straight at three. so half really was the size of the list + 1. most numbers, when looping, had increments that we greater than this + 1, but sometime it fell in it, hence giving an element doesnt exist error. Wensheng Wang (author) 17 years, 3 months ago

thank all. Thank all the commentors and your suggestions. I have made the correction: s=range(3,n+1,2) and half=(n+1)/2-1

Thanks. Hugh Brown 16 years, 1 month ago

More pythonesque. This works the same but has a cleaner loop:

``````for i in range(0, int(mroot/2)) :
m=2*i+3
if s[i]:
for j in range((m*m-3)/2, half, m) :
s[j]=0
`````` Hugh Brown 16 years, 1 month ago

More pythonesque. This works the same but has a cleaner loop:

``````for i in range(0, int(mroot/2)) :
m=2*i+3
if s[i]:
for j in range((m*m-3)/2, half, m) :
s[j]=0
`````` Dr. Drang 15 years, 10 months ago

Twice as fast? In my benchmarks (on an Intel iMac with n=1,000,000, n=10,000,000, and n=20,000,000), the recipe code is only about 10% faster than the "Alex" code, although its memory use is probably much smaller.

The code below runs about 50% faster than the "Alex" code and isn't too obscure, I think:

``````nroot = int(sqrt(n))
sieve = range(n+1)
sieve = 0

for i in xrange(2, nroot+1):
if sieve[i] != 0:
m = n/i - i
sieve[i*i: n+1:i] =  * (m+1)

sieve = [x for x in sieve if x !=0]
``````

The slice assignment seems to be the big time saver. The three terms that define the slice are the same as the three terms that define the range for j in the "Alex" code, which is why I think the code is still readable. m is a little tricky; it's the number of multiples of i that are greater than i**2 and less than or equal to n. The integer division used to calculate m will have to be changed for Python 3.

Since zero is false, I suppose the "!= 0" parts could be omitted from the conditionals, but I think it's easier to read with them in. larsspam 14 years, 6 months ago

I have made the correction: s=range(3,n+1,2) and half=(n+1)/2-1

Actually I still encountered the IndexError, even with these corrections. In my case, n was 100000.

half was 49999.5 and j was 49999.

I'm using Python 3. I guess the problem is that I need to use the new // integer division operator for computing both half and j, as Dr Drang mentioned.

Also, for others who may be using Python 3, you'll have to change the initial "s=range(3,n+1,2)" to "s=list(range(3,n+1,2))".

At least, that's what worked for me, a Python novice. Would be nice if someone more knowledgeable would say something definitive. Shfif py 14 years, 5 months ago

actually 209 is a prime number , but it isn't showing on your primes(3000) result list

Someone asked for a function that accepts a lower bound:

``````def primes(n,low=1):
'''return a list of prime numbers between LOW and N'''
low = low / 2 * 2 + 1   # convert LOW to odd number
lNum = range(low,n+1,2) # create a list of all the odd numbers between LOW and N
iRoot= n ** 0.5
iMid = len(lNum)        # number of numbers in the list
i = 0
m = 3                   # 2 is removed from the list because we didn't include the even numbers, so we start with removing
# the 3 multiples (9,15,21 ....)
while m < iRoot:
if lNum[i] != 0:      # if we haven't removed the number and its multiplies already
# I don't like this line, I think it might work in the original algorithm, but might make this one skip
# some numbers. If it does skip numbers, remove this 'if' statement
j = (m*m - low) / 2   # j is the _index_ of the first multiple of m to remove
while (j<iMid):

if (j >= 0):        # ignore indexes below zero (for obvious reasons)
lNum[j] = 0
j+=m

i += 1
m += 2

return [x for x in lNum if x != 0]
`````` Shfif py 14 years, 5 months ago

my mistake 209 / 11 is 19. Patrick Quinn 14 years, 3 months ago

For generating a large number of primes, the following improvement to Dr. Drang's code is helpful: def primes_up_to(n): nroot = int(sqrt(n)) sieve = [True] * (n+1) sieve = False sieve = False

``````    for i in xrange(2, nroot+1):
if sieve[i]:
m = n/i - i
sieve[i*i: n+1:i] = [False] * (m+1)

return [i for i in xrange(n+1) if sieve[i]]
``````

It's the same code, except it uses an array of booleans instead of ints for the sieve. For n = 50 million, it uses about 1/4 the memory and executes 25% faster (9 seconds instead of 12) on my machine. Patrick Quinn 14 years, 3 months ago

Sorry, I'm new here and haven't yet learned to use the preview button...

``````def primes_up_to(n):
nroot = int(sqrt(n))
sieve = [True] * (n+1)# range(n+1)
sieve = False
sieve = False

for i in xrange(2, nroot+1):
if sieve[i]:
m = n/i - i
sieve[i*i: n+1:i] = [False] * (m+1)

return [i for i in xrange(n+1) if sieve[i]]
`````` Storm 13 years, 5 months ago

Here is a clean starting point for an even faster sieve..

``````def primes(n):
s=range(0,n+1)
s=0
bottom=2
top=n//bottom
while (bottom*bottom<=n):
while (bottom<=top):
if s[top]:
s[top*bottom]=0
top-=1
bottom+=1
top=n//bottom
return [x for x in s if x]

print primes(1311)
``````

No sqrt isnt needed. duplicate work doesn't need to be done. And yes, this could be sped up even more by only doing odd numbers with 2 being a special case.. And faster yet with a 30/8 boolean-bit representation (2X3X5=30, 8 possible primes in any given section 30-60, 60-90...) offset(1,7 11, 13, 17, 19, 23, 29)

This base-code will never double-dip. so any given composite is removed once and only once. Giannis Fysakis 13 years, 3 months ago

Indeed it is fast :-) Junjie 13 years, 1 month ago

To Storm:

def myprime(n): t1=datetime.now() s=range(0,n+1) s=0 bottom=2 top=n//bottom while (bottombottom<=n): while (bottom<=top): if s[top]: s[topbottom]=0 top-=1 bottom+=1

###### Here is where modified
``````        while s[bottom]==0:
bottom+=1
``````
##### #
``````        top=n//bottom
print "cost : ",datetime.now()-t1
return [x for x in s if x]

>>> sum(prime.myprime(10000000))
Time cost :  0:00:08.536719
3203324994356L
>>> sum(prime.primes(10000000))
Time cost :  0:00:25.071403
``````

Program seems to run three times faster :) 3203324994356L rocks_smith 11 years, 6 months ago

with true false the timings are less for N<10000 but for higher values of N putting integers is faster davi mote 8 years, 2 months ago

I am new to programming, but wouldn't it be fastest to only check using previous primes up to half the value?

def primes(n): from datetime import datetime startTime = datetime.now() prime= list() x=5 prime.append(2) prime.append(3)

``````while x<n:
test=0
a=x/2
for i in prime:
if x%i==0:
test=1
break
if i > a:
break
if test==0:
prime.append(x)
x=x+1
print (prime)
print (datetime.now() - startTime)
`````` davi mote 8 years, 2 months ago

even faster

``````while x<n:
a=x**0.5
test=0
for i in prime:
if x%i==0:
test=1
break
if i>a:
break
if test==0:
prime.append(x)
x=x+1
print (prime)
print (datetime.now() - startTime)
`````` Created by Wensheng Wang on Mon, 7 Feb 2005 (PSF)

### Required Modules

• (none specified)