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):
    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

  1. 7 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
        possible_primes.insert(0,2)
        return [x for x in possible_primes if x > 0]
    
    

    By Peter on Sep 2, 2010

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

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

    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

    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 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))

    Regards,
    babsdoc

    By babsdoc on Aug 20, 2012

Post a Comment