Prime numbers and Numpy – Python

August 27, 2010 – 4:27 pm

I’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):
    for num in range(3,upto+1,2):
        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):
    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):
    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:

My impl is there listed under nolfonzo_prime6

  1. 14 Responses to “Prime numbers and Numpy – Python”

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

    def FastPrimeSieve(max):
        possible_primes =  range(3,max+1, 2)
        curr_index = -1
        max_index = len(possible_primes)
        for latest_prime in possible_primes:
            curr_index +=1
            if not latest_prime : continue
            for index_variable_not_named_j in xrange((curr_index+latest_prime),max_index, latest_prime): possible_primes[index_variable_not_named_j]=0
        return [x for x in possible_primes if x > 0]

    By Peter on Sep 2, 2010

  3. how about this one…. primesfrom2to(n)

    By bob on Sep 3, 2010

  4. Very 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 nolfonzo on Sep 3, 2010

  5. ok – 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 nolfonzo on Sep 3, 2010

  6. I posted my 1/2 sieve implementation on

    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:

    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 nolfonzo on Sep 4, 2010

  7. Pete’s is the fatest my Mac. But I love your first numpy pythonic solution without loops.

    By Tom Luo on Nov 2, 2010

  8. Hi,
    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))


    By babsdoc on Aug 20, 2012

  9. This 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 Heidi on Jul 4, 2014

  10. I 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 Octavio on Sep 11, 2014

  11. Apologies 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 Salmoni on Feb 12, 2015

  12. What’s up to every body, it’s my first go to see of
    this web site; this website includes remarkable and in fact
    fine material in support of readers.

    By hoodies for women on Aug 18, 2015

  13. Every weekend i used to go to see this website, as i want enjoyment, as this this website conations truly nice funny information too.

    By Pearline on Aug 20, 2015

  1. 2 Trackback(s)

  2. Dec 26, 2014: Fixed Fastest way to list all primes below N #dev #it #asnwer | Good Answer
  3. Jan 8, 2015: Fixed: Fastest way to list all primes below N #development #solution #it | IT Info

Post a Comment