Is there a simple algorithm that can determine if X is prime?

Question

I have been trying to work my way through Project Euler, and have noticed a handful of problems ask for you to determine a prime number as part of it.

1. I know I can just divide x by 2, 3, 4, 5, ..., square root of X and if I get to the square root, I can (safely) assume that the number is prime. Unfortunately this solution seems quite klunky.

2. I've looked into better algorithms on how to determine if a number is prime, but get confused fast.

Is there a simple algorithm that can determine if X is prime, and not confuse a mere mortal programmer?

Thanks much!

Answer 1

The first algorithm is quite good and used a lot on Project Euler. If you know the maximum number that you want you can also research Eratosthenes's sieve.

If you maintain the list of primes you can also refine the first algo to divide only with primes until the square root of the number.

With these two algoritms (dividing and the sieve) you should be able to solve the problems.

Edit: fixed name as noted in comments

Answer 2

To generate all prime numbers less than a limit Sieve of Eratosthenes (the page contains variants in 20 programming languages) is the oldest and the simplest solution.

In Python:

def iprimes_upto(limit):
    is_prime = [True] * limit
    for n in range(2, limit):
        if is_prime[n]:
           yield n
           for i in range(n*n, limit, n): # start at ``n`` squared
               is_prime[i] = False

Example:

>>> list(iprimes_upto(15))
[2, 3, 5, 7, 11, 13]

Answer 3

I see that Fermat's primality test has already been suggested, but I've been working through Structure and Interpretation of Computer Programs, and they also give the Miller-Rabin test (see Section 1.2.6, problem 1.28) as another alternative. I've been using it with success for the Euler problems.

Answer 4

Keeping in mind the following facts (from MathsChallenge.net):

• All primes except 2 are odd.
• All primes greater than 3 can be written in the form 6k - 1 or 6k + 1.
• You don't need to check past the square root of n

Here's the C++ function I use for relatively small n:

bool isPrime(unsigned long n)
{
    if (n == 1) return false; // 1 is not prime
    if (n < 4) return true; // 2 and 3 are both prime
    if ((n % 2) == 0) return false; // exclude even numbers
    if (n < 9) return true; //we have already excluded 4, 6, and 8.
    if ((n % 3) == 0) return false; // exclude remaining multiples of 3

    unsigned long r = floor( sqrt(n) );
    unsigned long f = 5;
    while (f <= r)
    {
        if ((n % f) == 0)  return false;
        if ((n % (f + 2)) == 0) return false;
        f = f + 6;
    }
    return true; // (in all other cases)
}

You could probably think of more optimizations of your own.

Answer 5

Here's a simple optimization of your method that isn't quite the Sieve of Eratosthenes but is very easy to implement: first try dividing X by 2 and 3, then loop over j=1..sqrt(X)/6, trying to divide by 6*j-1 and 6*j+1. This automatically skips over all numbers divisible by 2 or 3, gaining you a pretty nice constant factor acceleration.

Answer 6

I'd recommend Fermat's primality test. It is a probabilistic test, but it is correct surprisingly often. And it is incredibly fast when compared with the sieve.

Answer 7

For reasonably small numbers, x%n for up to sqrt(x) is awfully fast and easy to code.

Simple improvements:

test 2 and odd numbers only.

test 2, 3, and multiples of 6 + or - 1 (all primes other than 2 or 3 are multiples of 6 +/- 1, so you're essentially just skipping all even numbers and all multiples of 3

test only prime numbers (requires calculating or storing all primes up to sqrt(x))

You can use the sieve method to quickly generate a list of all primes up to some arbitrary limit, but it tends to be memory intensive. You can use the multiples of 6 trick to reduce memory usage down to 1/3 of a bit per number.

I wrote a simple prime class (C#) that uses two bitfields for multiples of 6+1 and multiples of 6-1, then does a simple lookup... and if the number i'm testing is outside the bounds of the sieve, then it falls back on testing by 2, 3, and multiples of 6 +/- 1. I found that generating a large sieve actually takes more time than calculating primes on the fly for most of the euler problems i've solved so far. KISS principle strikes again!

I wrote a prime class that uses a sieve to pre-calculate smaller primes, then relies on testing by 2, 3, and multiples of six +/- 1 for ones outside the range of the sieve.

Answer 8

For Project Euler, having a list of primes is really essential. I would suggest maintaining a list that you use for each problem.

I think what you're looking for is the Sieve of Eratosthenes.

Answer 9

Your right the simples is the slowest. You can optimize it somewhat.

Look into using modulus instead of square roots. Keep track of your primes. you only need to divide 7 by 2, 3, and 5 since 6 is a multiple of 2 and 3, and 4 is a multiple of 2.

Rslite mentioned the eranthenos sieve. It is fairly straight forward. I have it in several languages it home. Add a comment if you want me to post that code later.

---

Here is my C++ one. It has plenty of room to improve, but it is fast compared to the dynamic languages versions.

// Author: James J. Carman
// Project: Sieve of Eratosthenes
// Description: I take an array of 2 ... max values. Instead of removeing the non prime numbers,
// I mark them as 0, and ignoring them.
// More info: http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
#include <iostream>

int main(void) {
        // using unsigned short.
        // maximum value is around 65000
        const unsigned short max = 50000;
        unsigned short x[max];
        for(unsigned short i = 0; i < max; i++)
                x[i] = i + 2;

        for(unsigned short outer = 0; outer < max; outer++) {
                if( x[outer] == 0)
                        continue;
                unsigned short item = x[outer];
                for(unsigned short multiplier = 2; (multiplier * item) < x[max - 1]; multiplier++) {
                        unsigned int searchvalue = item * multiplier;
                        unsigned int maxValue = max + 1;
                        for( unsigned short maxIndex = max - 1; maxIndex > 0; maxIndex--) {
                                if(x[maxIndex] != 0) {
                                        maxValue = x[maxIndex];
                                        break;
                                }
                        }
                        for(unsigned short searchindex = multiplier; searchindex < max; searchindex++) {
                                if( searchvalue > maxValue )
                                        break;
                                if( x[searchindex] == searchvalue ) {
                                        x[searchindex] = 0;
                                        break;
                                }
                        }
                }
        }
        for(unsigned short printindex = 0; printindex < max; printindex++) {
                if(x[printindex] != 0)
                        std::cout << x[printindex] << "\t";
        }
        return 0;
}

I will throw up the Perl and python code I have as well as soon as I find it. They are similar in style, just less lines.

Answer 10

Here is a simple primality test in D (Digital Mars):

/** 
 * to compile:
 * $ dmd -run prime_trial.d
 * to optimize:
 * $ dmd -O -inline -release prime_trial.d 
 */
module prime_trial;

import std.conv : to;  
import std.stdio : w = writeln;

/// Adapted from: http://www.devx.com/vb2themax/Tip/19051 
bool 
isprime(Integer)(in Integer number) 
{
  /* manually test 1, 2, 3 and multiples of 2 and 3 */
  if (number == 2 || number == 3)
    return true;
  else if (number < 2 || number % 2 == 0 || number % 3 == 0)
    return false;

  /* we can now avoid to consider multiples 
   * of 2 and 3. This can be done really simply 
   * by starting at 5 and incrementing by 2 and 4 
   * alternatively, that is: 
   *    5, 7, 11, 13, 17, 19, 23, 25, 29, 31, 35, 37, ...    
   * we don't need to go higher than the square root of the number */
  for (Integer divisor = 5, increment = 2; divisor*divisor <= number; 
       divisor += increment, increment = 6 - increment) 
    if (number % divisor == 0)
      return false;

  return true;  // if we get here, the number is prime
}

/// print all prime numbers less then a given limit
void main(char[][] args) 
{
  const limit = (args.length == 2) ? to!(uint)(args[1]) : 100;
  for (uint i = 0; i < limit; ++i) 
    if (isprime(i))
      w(i);
}

Answer 11

I am working thru the Project Euler problems as well and in fact just finished #3 (by id) which is the search for the highest prime factor of a composite number (the number in the ? is 600851475143).

I looked at all of the info on primes (the sieve techniques already mentioned here), and on integer factorization on wikipedia and came up with a brute force trial division algorithm that I decided would do.

So as I am doing the euler problems to learn ruby I was looking into coding my algorithm and stumbled across the mathn library which has a Prime class and an Integer class with a prime_division method. how cool is that. i was able to get the correct answer to the problem with this ruby snippet:

require "mathn.rb"
puts 600851475143.prime_division.last.first

this snippet outputs the correct answer to the console. of course i ended up doing a ton of reading and learning before i stumbled upon this little beauty, i just thought i would share it with everyone...

Answer 12

Is not optimized but it's a very simple function.

    function isprime(number){

    if (number == 1)
        return false;

    var times = 0;

    for (var i = 1; i <= number; i++){
        if(number % i == 0){
            times ++;
        }
    }
        if (times > 2){
            return false;
        }

    return true;
    }

Answer 13

Maybe this implementation in Java can be helpful:

public class SieveOfEratosthenes {

    /**
     * Calling this method with argument 7 will return: true true false false true false true false
     * which must be interpreted as : 0 is NOT prime, 1 is NOT prime, 2 IS prime, 3 IS prime, 4 is NOT prime
     * 5 is prime, 6 is NOT prime, 7 is prime.
     * Caller may either revert the array for easier reading, count the number of primes or extract the prime values
     * by looping.
     * @param upTo Find prime numbers up to this value. Must be a positive integer.
     * @return a boolean array where index represents the integer value and value at index returns
     * if the number is NOT prime or not.
     */
    public static boolean[] isIndexNotPrime(int upTo) {
        if (upTo < 2) {
            return new boolean[0];
        }

        // 0-index array, upper limit must be upTo + 1
        final boolean[] isIndexNotPrime = new boolean[upTo + 1];

        isIndexNotPrime[0] = true; // 0 is not a prime number.
        isIndexNotPrime[1] = true; // 1 is not a prime number.

        // Find all non primes starting from 2 by finding 2 * 2, 2 * 3, 2 * 4 until 2 * multiplier > isIndexNotPrime.len
        // Find next by 3 * 3 (since 2 * 3 was found before), 3 * 4, 3 * 5 until 3 * multiplier > isIndexNotPrime.len
        // Move to 4, since isIndexNotPrime[4] is already True (not prime) no need to loop..
        // Move to 5, 5 * 5, (2 * 5 and 3 * 5 was already set to True..) until 5 * multiplier > isIndexNotPrime.len
        // Repeat process until i * i > isIndexNotPrime.len.
        // Assume we are looking up to 100. Break once you reach 11 since 11 * 11 == 121 and we are not interested in
        // primes above 121..
        for (int i = 2; i < isIndexNotPrime.length; i++) {
            if (i * i >= isIndexNotPrime.length) {
                break;
            }
            if (isIndexNotPrime[i]) {
                continue;
            }
            int multiplier = i;
            while (i * multiplier < isIndexNotPrime.length) {
                isIndexNotPrime[i * multiplier] = true;
                multiplier++;
            }
        }

        return isIndexNotPrime;
    }

    public static void main(String[] args) {
        final boolean[] indexNotPrime = SieveOfEratosthenes.isIndexNotPrime(7);
        assert !indexNotPrime[2]; // Not (not prime)
        assert !indexNotPrime[3]; // Not (not prime)
        assert indexNotPrime[4]; // (not prime)
        assert !indexNotPrime[5]; // Not (not prime)
        assert indexNotPrime[6]; // (not prime)
        assert !indexNotPrime[7]; // Not (not prime)
    }
}

Answer 14

I like this python code.

def primes(limit) :
    limit += 1
    x = range(limit)
    for i in xrange(2,limit) :
        if x[i] ==  i:
            x[i] = 1
            for j in xrange(i*i, limit, i) :
                x[j] = i
    return [j for j in xrange(2, limit) if x[j] == 1]

A variant of this can be used to generate the factors of a number.

def factors(limit) :
    limit += 1
    x = range(limit)
    for i in xrange(2,limit) :
        if x[i] == i:
            x[i] = 1
            for j in xrange(i*i, limit, i) :
                x[j] = i
    result = []
    y = limit-1
    while x[y] != 1 :
        divisor = x[y]
        result.append(divisor)
        y /= divisor
    result.append(y)
    return result

Of course, if I were factoring a batch of numbers, I would not recalculate the cache; I'd do it once and do lookups in it.

Answer 15

The AKS prime testing algorithm:

Input: Integer n > 1  


if (n is has the form ab with b > 1) then output COMPOSITE  

r := 2  
while (r < n) {  
    if (gcd(n,r) is not 1) then output COMPOSITE  
    if (r is prime greater than 2) then {  
        let q be the largest factor of r-1  
        if (q > 4sqrt(r)log n) and (n(r-1)/q is not 1 (mod r)) then break  
    }  
    r := r+1  
}  

for a = 1 to 2sqrt(r)log n {  
    if ( (x-a)n is not (xn-a) (mod xr-1,n) ) then output COMPOSITE  
}  

output PRIME;   

Answer 16

another way in python is:

import math

def main():
    count = 1
    while True:
        isprime = True

        for x in range(2, int(math.sqrt(count) + 1)):
            if count % x == 0: 
                isprime = False
                break

        if isprime:
            print count


        count += 2


if __name__ == '__main__':
    main()