## Prime numbers and Numpy – Python

August 27, 2010 – 4:27 pmI’ve been looking at the element-wise operations in Numpy arrays. The code below generates 100K primes in about 1.3 secs on my laptop.

import numpy import math def prime(upto=100): return filter(lambda num: (num % numpy.arange(2,1+int(math.sqrt(num)))).all(), range(2,upto+1))

>>> prime() [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97] >>>

Similiar code that uses a for loop runs about 5 times longer for 100K primes, about 7 secs:

import numpy import math def prime2(upto=100): return filter(lambda num: numpy.array([num % factor for factor in range(2,1+int(math.sqrt(num)))]).all(), range(2,upto+1))

By not including even numbers, the original code can be optimized to run in half the time again, about .7 secs for 100K primes.

import numpy import math def prime3(upto=100): return [2]+filter(lambda num: (num % numpy.arange(3,1+int(math.sqrt(num)),2)).all(), range(3,upto+1,2))

The above code is not the fastest implementation, but illustrative of using the element-wise operations in Numpy arrays rather than looping through a list.

We can write faster code by breaking after finding the first modulo equal zero condition. The code below runs in about .4 seconds for 100K primes:

import numpy import math def prime4(upto=100): primes=[2] for num in range(3,upto+1,2): isprime=True for factor in range(3,1+int(math.sqrt(num)),2): if not num % factor: isprime=False; break if isprime: primes.append(num) return primes

Ok Scruffy Pete has laid down the challenge in the comments section so I need to get serious about optimization. The code below uses a few Numpy tricks for a zippy implementation and generates 1 million primes in about .03 secs, about 6X faster than Scruffy Pete’s.

import math import numpy def prime5(upto): primes=numpy.arange(2,upto+1) isprime=numpy.ones(upto-1,dtype=bool) for factor in primes[:int(math.sqrt(upto))]: if isprime[factor-2]: isprime[factor*2-2::factor]=0 return primes[isprime]

Scruffy Pete… Ball’s in your court matey!

UPDATE 9/3 – did a half sieve style implementation of the above (see comments) which runs a little faster – under .02 secs for 1 million primes:

import math import numpy def prime6(upto): primes=numpy.arange(3,upto+1,2) isprime=numpy.ones((upto-1)/2,dtype=bool) for factor in primes[:int(math.sqrt(upto))]: if isprime[(factor-2)/2]: isprime[(factor*3-2)/2::factor]=0 return numpy.insert(primes[isprime],0,2)

I posted this on the stack overflow site (see comments section). They run timings which you can see on:

http://ideone.com/oDg2Y

My impl is there listed under nolfonzo_prime6

## 12 Responses to “Prime numbers and Numpy – Python”

The Sieve Of Eratosthenes is faster. One possible implementation:By

Peteron Sep 2, 2010how about this one…. primesfrom2to(n)

http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188

By

bobon Sep 3, 2010Very interesting – thanks Bob.

I notice that the fastest pure python solution on that post runs at about .15 secs on my machine.

My Numpy solution is about 3X better than that.

The half/third-sieve Numpy solutions primesfrom3to(n) and primesfrom2to(n) from your link are almost twice as fast as mine, which is probably expected.

At some point I’ll try to turn mine into a half sieve and get proper timings.

By

nolfonzoon Sep 3, 2010ok – made the change for the 1/2 sieve (see update). The code now runs comparable to the 1/2 sieve impl in bob’s link. Though dare I say my impl is a little more elegant?

I thought about having a shot at the 1/3 sieve but after a few minutes thought it may be time to stop the madness.

By

nolfonzoon Sep 3, 2010I posted my 1/2 sieve implementation on

http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188

They run automated tests against the submitted prime scripts, and you can see that mine (nolfonzo_prime6) wasn’t the fastest but held its own against the other numpy implementations clocking in at 12ms for 1 million primes.

You can look at the scripts and results here: http://ideone.com/oDg2Y

Some of the faster implementations dispensed with the initial prime array I had, and backed into it from the sieve, or indexes, directly at the end. I guess that’s the way you’re supposed to do a sieve but it does get a little more complicated – more memory efficient though I’m guessing.

Interesting stuff.

By

nolfonzoon Sep 4, 2010Pete’s is the fatest my Mac. But I love your first numpy pythonic solution without loops.

By

Tom Luoon Nov 2, 2010Hi,

I’m not a programmer, but I’m interested in learning python. I’m trying to write a program in python. Can you please help me out?

I need to generate primes by using 2 lists,

the first list will be have only one value [2].

The second list will have odd numbers lets say from 3 to 20 [3,5,7,9,11,13,15,17,19]

The code should will take the first number 3 and check if 3 % 2 ==0, if it isnt, then 3 should be added to the first list [2,3] the code will then go to the next number 5 and check divisibility in reverse, if 5 % 3 ==0 if not then check if 5 % 2==0 if not then add 5 to the list [2,3,5] and so on.

May seem funny, but I’m finding it hard to program.

Also, kindly explain what is happening in the following line of your code:

return filter(lambda num: (num % numpy.arange(2,1+int(math.sqrt(num)))).all(), range(2,upto+1))

Regards,

babsdoc

By

babsdocon Aug 20, 2012This content is very interesting but it took me a long time to find it in google.

I found it on 13 spot, you should focus on quality backlinks building, it will help you

to increase traffic. And i know how to help you, just search in google – k2 seo tips

By

Heidion Jul 4, 2014I read a lot of interesting articles here. Probably you spend

a lot of time writing, i know how to save you a lot of time, there is an online tool that creates high quality, google

friendly posts in seconds, just search in google – laranitas free content source

By

Octavioon Sep 11, 2014Apologies if this comment is useless but I tried my own routine for calculating primes under a certain number. It might not be the best but it’s far from the slowest and is easy to understand. It works by doing sieve on a boolean vector and converts all True values to indices. It’s odd because it avoids all division in the code.

def ajs_primes3a(upto):

mat = np.ones((upto), dtype=bool)

mat[0] = False # 0 is not a prime

mat[1] = False # 1 is (probably) not a prime

mat[4::2] = False # neither is anything divisible by 2 except 2

for idx in range(3, int(sqrt(upto)), 2): # go through odd numbers from 3 up

mat[idx*2::idx] = False # sieve them out

return np.where(mat == True) # produce a vector where elements = True which gives primes!

I also wondered whether using a sparse matrix would be quicker.

By

Alan J Salmonion Feb 12, 2015