Pythonism

joy in simplicity – platonic world of code and math

How to find prime numbers with Python

with 27 comments

One of the first things a beginning programmer should try is code for testing the primality of a given number, or getting a list of primes. A prime is a number only divisible by itself and one. (Note that divisibility here means divisible without a remainder…)

Searching for big primes can get very subtle. This piece will only scratch the surface of primes, but it’s still a good challenge which sharpens our thinking. The code we will use will evolve through the article, so please read sequentially from the start.

There are two naive methods to find primes. First the Divisor Test which is useful for checking the primality of a given number, or small collection of numbers, second the Prime Sieve which is faster to generate lots of primes.

The Divisor Test

If we have a number, n, then the first immediate way to test if it is prime is to divide every number less than n into n and if none divide then answer Yes. If any divide then answer No. This is a candidate for using a loop. What about:


startnumber=7
startnumber*=1.0
prime = True
for divisor in range(2,startnumber):
    if startnumber/divisor==int(startnumber/divisor):
        prime=False
if prime==False:
    print startnumber," is not prime"
else:
    print startnumber, "is prime"

When Python divides two numbers of type integer it returns an integer which is the result of the division rounded down. To get the answer in decimal we have to set the type of at least one of the numbers to a float type by multiplying by 1.0.

output:

>>>
7.0 is prime
>>>

So that works for 7 but there is a lot of refining still to do. We can tighten up the loop as follows: Look at each number the loop tries as the divisor: 2,3,4,5,6… It is apparent that there’s not much point in dividing 7 by 6 since we know that 6 is too big to divide 7. In fact n-1 will only ever divide n when n=2. so what do we want the loop to count up to ? The square root of n. If there are no divisors below this value there aren’t any at all. This will make the code quicker, especially when we plug in some big numbers later.


startnumber=7
startnumber*=1.0
prime = True
for divisor in range(2,startnumber**0.5+1):
    if startnumber/divisor==int(startnumber/divisor):
        prime=False
if prime==False:
    print startnumber," is not prime"
else:
    print startnumber, "is prime"

The +1 in the range statement is necessary because python would stop counting at int(n**0.5)-1 otherwise.

Time for more improvements. Let’s turn the code into a function so we can combine it with other code more easily.


def isprime(startnumber):
    startnumber*=1.0
    prime = True
    for divisor in range(2,startnumber**0.5+1):
        if startnumber/divisor==int(startnumber/divisor):
            prime=False
    return prime

…this can be shortened too. What’s this business of creating the variable, “prime” ? It isn’t necessary. Like removing excess weight from a racing car we can chop out the padding and have:


def isprime(startnumber):
    startnumber*=1.0
    for divisor in range(2,int(startnumber**0.5)+1):
        if startnumber/divisor==int(startnumber/divisor):
            return False
    return True

This is getting nicely streamlined. As soon as a divisor is found we break out of the loop, if no divisors are ever found we know n is prime so the function returns True once the loop has terminated. Now we can call this function from a loop and generate a list of primes by adding


for a in range(2,100):
    if isprime(a):
        print a

which gives:

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

OK, fine. This will work in reasonable time for quite big numbers, which could be fun. To find all the primes under 10,000 took about 20 seconds on my pc. The Prime Sieve is faster for lots of primes so read on.

Looking again at what that loop is doing we see that it is using even divisors to test. Well if n isn’t divisible by 2 it cant be divisible by any even number so let’s stop it testing those.

what about:


def isprime(startnumber):
    startnumber*=1.0
    if startnumber/2==int(startnumber/2) and startnumber!=2:
        return False
    for divisor in range(3,int(startnumber**0.5)+1,2):
        if startnumber/divisor==int(startnumber/divisor):
            return False
    return True

We had to add a new line testing solely for 2, but its outside the loop so won’t slow us down much. Also this function might return False for n == 2, so we’ve added another condition that stops this. Good, but there’s more.

The division operation inside the loop could be trimmed down. At the moment, within the loop, there are two division operations and an int(n) operation. All divisions in the function could be changed to single modulus operations, using Python’s % operator.


def isprime(startnumber):
    startnumber*=1.0
    if startnumber%2==0 and startnumber!=2:
        return False
    for divisor in range(3,int(startnumber**0.5)+1,2):
        if startnumber%divisor==0:
            return False
    return True

Whatever code is inside the loop will be repeated lots of times. Thus trimming and speeding up this region of the code is essential to shorten runtime. It is fun to evolve code like this, don’t you think? This is more how a real programmer might work, happy the person who can write down perfect code immediately, for me it tends to evolve, as I have tried to show here in these simple examples.

Our final evolutionary step is as follows. We still have redundant calculations going on in our code, for instance although we only check odd numbers > 2 now we still check for divisibility by a lot of numbers that have prime factors in common, such as checking for divisibility by 5 which is good but later checking for 15 which is a wasted task. We start by observing that numbers of the form 6n+1 or 6n-1 are usually prime. In fact this formula generates all the primes > 3 as well as a few composites. Fewer composites than above though so its still an improvement. We check for divisibility separately by 2 then 3, then by 6n+1 and 6n-1 for all n until 6n+1>= root n


def isprime(n):
    n*=1.0
    if n%2==0 and n!=2 or n%3==0 and n!=3:
        return False
    for b in range(1,int((n**0.5+1)/6.0+1)):
        if n%(6*b-1)==0:
            return False
        if n %(6*b+1)==0:
           return False
    return True

Note that the best way of testing speed is to use the time module in Python. The shortest code isn’t always the fastest either, this last code is quicker than our earlier attempts but slightly longer. If the program is only occupying a few bytes of memory (space complexity) then speed (time complexity) is clearly more important…

The Prime Sieve

If you want lots of primes then there is a faster way, though. Start with a collection of all the natural numbers and remove all the even ones. Then remove all multiples of 3, then 5 and so on through the earlier primes. After you have taken out all possible multiples what remains must be the primes, since they are not multiples. The code for this works quicker when you need a collection of the first n primes, for large n. Here’s a naive yet short Sieve of Eratosthenes in python. This one removes the composites from the list in situ until there are none left and only primes remain.

A flaw is that many numbers are checked more than once.

a=list(range(2,1000))
for x in a:
    n=x
    while x*n <= max(a):
        if x*n in a:
            a.remove(x*n)
        n+=1
print(a)

How can we eliminate the redundancy ? The next one i tried uses a slightly different approach, using two lists, one of composites which is populated in the first loop and one of the primes which are every integer not in the composites, populated in the second loop.

comp=[]
for a in range(2,1000):
    if a not in comp:
        n=a
        while n*a <=1000:
            comp.append(n*a)
            n+=1
primes=[]
for c in range(2,1000):
    if c not in comp:
        primes.append(c)
print(primes)

but this next one is even better since it seems to avoid any repetition. For each starting factor a temporary list of composites is generated, these are then removed from the primes list, and then the logic moves to the next prime factor in the list and repeats the same process until it’s done.

a=list(range(2,1000))
for x in a:
    comp=[]
    for y in a:
        comp.append(x*y)
    for z in comp:
        if z in a:
            a.remove(z)
print(a)

aha ! but there was redundancy in that one too. The first unchecked composite is n*n, but we were counting n*2,n*3, n*5 etc up to n*n, this can be remedied. Also we need to avoid checking the numbers when they are too large, so here is our final version.

a=list(range(2,1000))
for x in a:
    comp=[]
    for y in range(a.index(x),len(a)):
        q=x*a[y]
        if q<=max(a):
            comp.append(q)
        else:
            break
    for z in comp:
        a.remove(z)

A pretty picture

Finally, here is a plot of the first thousand primes with a vertical line for each prime and blank lines for composites. It gives you a loose idea of the density of primes. Looks a bit like a giant barcode, as Arka says.

The large primes used for encryption have hundreds of digits. With what we have done here it would probably take thousands of years of computer time to verify the primality of such numbers, for those you use clever math, see for example this wikipedia page. Brute force is usually the easiest to code, but as our skills develop we will try and use insight to work in harmony with the computer. Use the wetware inside your skull for the bits it is best for, and leave the silicon to its specialty.

About these ads

Written by Trip Technician

May 4, 2008 at 6:36 pm

Posted in Math, Primes

Tagged with ,

27 Responses

Subscribe to comments with RSS.

  1. [...] finding more than one prime. It’s possibly the simplest to code but you should also look at exploring primes since the function there is faster as I timed it. For really large numbers even these tests can [...]

    • This helped me with a program that counts perfect numbers, but it stopped after the eighth one, probably because it was getting too large.

      Dylan

      August 3, 2012 at 2:44 pm

  2. muy bueno me sirvio de ayuda

    #126.- Es_primo.py…

    num = 7
    primo = Trae

    for divisor in range(2, num):
    if ((num % divisor) == 0):
    primo = False
    if primo:
    print “El número” , num, “es primo”
    else:
    print “El número”, num, “No es primo”

    #128.- Primo.py…

    primo = Trae
    divisor = 2
    num = int(raw_input(“Dame un número: ”))
    while divisor < num and primo:
    if ((num % divisor) == 0):
    primo = False
    divisor += 1

    if primo:
    print “El número”, num, “es primo”
    else:
    print “El número”, num, “No es primo”

    yasix

    January 4, 2009 at 3:15 am

  3. bien hecho … seguir adelante y vamos a hacer un espacio de piolot que todavía!

    pythonisms

    January 5, 2009 at 8:38 pm

  4. [...] the prime number page [...]

  5. Sorry to disent, but the sieve of Erathostenes is much, much more efficient than the code you propose… To calculate all primes tor 10,000 it may be there is no significant difference, but try it again for 100,000, or 1,000,0000 and you are in for a big shock.

    Jaime

    March 4, 2009 at 6:28 pm

  6. can you write the program for a particular prime number to find example i want 1000 prime is 7919

    mallikarjun

    June 16, 2010 at 3:18 am

  7. how can i find the 1000th prime number , plz reply urgently

    poornima

    June 16, 2010 at 4:51 am

  8. that it is possible to buy exactly 50, 51, 52, 53, 54, and 55 McNuggets, by finding solutions to the Diophantine equation. You can solve this in your head, using paper and pencil, or writing a program. However you chose to solve this problem, list the combinations of 6, 9 and 20 packs of McNuggets you need to buy in order to get each of the exact amounts. Given that it is possible to buy sets of 50, 51, 52, 53, 54 or 55 McNuggets by combinations of 6, 9 and 20 packs, show that it is possible to buy 56, 57,…, 65 McNuggets. In other words, show how, given solutions for 50-55, one can derive solutions for 56-65.

    fhiroz

    June 22, 2010 at 6:22 am

  9. Numpy can help implement some of these sieves… take a look at my blog post:
    http://rebrained.com/?p=458

    nolfonzo

    September 3, 2010 at 6:28 pm

  10. take a look at this site:
    http://stackoverflow.com/questions/2068372/fastest-way-to-list-all-primes-below-n-in-python/3035188#3035188

    and associated timings:
    http://ideone.com/oDg2Y

    The sieve implementations are very fast.

    I posted an implementation on my blog rebrained.com. It uses numpy.

    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)
    

    nolfonzo

    September 4, 2010 at 3:08 pm

  11. yes

    bhgathsingh

    October 5, 2010 at 3:14 pm

  12. i want to know how get prime number with python

    mallesh@786

    January 29, 2011 at 4:12 pm

  13. That plot for the first 5000 prime numbers looks like some bar code ! ….. LOL !

    Arkapravo

    July 13, 2011 at 4:09 pm

  14. [...] finally upgraded my isPrime function, as my old one was inefficient. Here’s a good blog on Pythonism that walks through how to find prime [...]

  15. could you explain * and **. I am not sure of their functions. I am guessing ** is ‘to the power of’? but what is the single star as in startnumber*=1.0 ?

    lwclasses

    August 31, 2011 at 9:53 pm

    • ** is indeed ‘to the power of’. a*=1.0 means the same as a=a*1.0 ie just multiply a by 1.0.

      this is done to convert a to a decimal otherwise python will do integer division which doesn’t give a remainder. Try it in the interpreter:

      >>> a=5
      >>> a/2
      2
      >>> a*=1.0
      >>> a/2
      2.5

      I am guessing by your question that you are new to Python. I recommend obtaining it on your own machine and playing !! fun awaits.

      Trip Technician

      September 1, 2011 at 9:42 am

  16. Great, great run-through! I am new to Python and am crap at math, so this is really helpful.

    Anthony Arroyo

    February 26, 2012 at 6:37 pm

  17. 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

    gooeystuffbsdoc

    August 20, 2012 at 8:20 pm

    • what you are talking about in this problem is known as the “prime sieve” or “sieve of Eratosthenes”. It is a method of finding primes by starting with all numbers… crossing out multiples of 2… then 3 then 4 and going through untill what you have left cannot be divisible by 2… 3… 4 etc and thus must be a list of primes. If you look at the info on sieves below on this page it might help.

      You’re doing well but keep plugging at it. I can’t write the solution for you because you have to learn, but good luck. Think of making a list and eliminating any even number – you’ve done this easily in your comment. It is not a great leap then to see a way to cut out every multiple of 3… something like

      if number % 3 == 0:
      list.remove(number)

      etc. I wouldn’t worry about the other code because using lambda and numpy is more advanced and not worth getting into for now. go for it and post your results… I’ll respond unless I have gone to sleep !

      Trip Technician

      August 20, 2012 at 8:29 pm

  18. Hi Trip,

    I just need do understand how to get the divisors in the reverse order.

    Regards,
    babsdoc

    babsdoc

    August 20, 2012 at 8:32 pm

    • if you have a list like [0,1,2,3,4,5] then to iterate through the list you might do:

      for x in my_list:
      do something

      if you want to go through in the opposite direction then do the same thing but on the reversed list. You get the reversed list by doing

      for x in my_list[::-1]:
      do something

      Trip Technician

      August 20, 2012 at 8:39 pm

  19. Python style:

    def is_prime(n):
    if n<2:
    return False
    p=n
    n=range(n)
    n.remove(0)
    n.remove(1)
    for i in n:
    if p%i==0:
    return False
    return True

    qiubox

    October 29, 2012 at 3:38 am

  20. Hey nice work, I used your code in my example of solving some math problems:
    http://deadendmath.com/all-primes-ending-in-1/

    obitus

    February 11, 2013 at 12:41 am

  21. […] online I came across this one person's website that had exactly what I am looking for. If you look here there is a nice program that tells us if a number is prime, which is the first key in our quest […]

  22. Cool article. Very detailed. I wrote a similar article a while back about finding prime numbers in python.

  23. Help with a program that determines if a number is a prime number then finds how many factors it has

    JujuBuoy

    November 5, 2013 at 7:59 am


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 45 other followers

%d bloggers like this: