Prime number calculation fun

Question

We're having a bit of fun here at work. It all started with one of the guys setting up a Hackintosh and we were wondering whether it was faster than a Windows Box of (nearly) same specs that we have. So we decided to write a little test for it. Just a simple Prime number calculator. It's written in Java and tells us the time it takes to calculate the first n Prime numbers.

Optimised version below - now takes ~6.6secs

public class Primes {

    public static void main(String[] args) {
     int topPrime = 150000;
     int current = 2;
     int count = 0;
     int lastPrime = 2;

     long start = System.currentTimeMillis();

     while (count < topPrime) {

      boolean prime = true;

      int top = (int)Math.sqrt(current) + 1;

      for (int i = 2; i < top; i++) {
       if (current % i == 0) {
        prime = false;
        break;
       }
      }

      if (prime) {
       count++;
       lastPrime = current;
      }
      if (current == 2) {
       current++;
      } else {
       current = current + 2;
      }
     }

     System.out.println("Last prime = " + lastPrime);
     System.out.println("Total time = " + (double)(System.currentTimeMillis() - start) / 1000);
    } 
}

We've pretty much lost the plot of the whole Hackintosh vs PC thing and are just having some fun with optimising it. First attempt with no optimisations (the above code has a couple) ran around 52.6min to find the first 150000 prime numbers. This optimisation is running around 47.2mins.

If you want to have a go and post your results, then stick em up.

Specs for the PC I'm running it on are Pentium D 2.8GHz, 2GB RAM, running Ubuntu 8.04.

Best Optimisation so far has been the square root of current, first mentioned by Jason Z.

Answer 1

cough http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes cough

Answer 2

Keeping in mind that there are better ways to look for primes...

I think that you can take this loop:

for (int i = 2; i < top; i++)

and make it so that your counter variable i goes from 3 and only tries to do the mod on odd numbers, since all primes other than 2 are never divisible by any even numbers.

Answer 3

Well I see a couple of quick optimizations that can be done. First you don't have to try each number up to half of the current number.

Instead you only have try up to the square root of the current number.

And the other optimization was what BP said with a twist: Instead of

int count = 0;
...
for (int i = 2; i < top; i++)
...
if (current == 2)
  current++;
else
  current += 2;

use

int count = 1;
...
for (int i = 3; i < top; i += 2)
...
current += 2;

This should speed things up quite a lot.

EDIT:
More optimization courtesy of Joe Pineda:
Remove the variable "top".

int count = 1;
...
for (int i = 3; i*i <= current; i += 2)
...
current += 2;

If this optimization indeed increases speed is up to java.
Calculating the square root takes a lot of time compared to multiplying two numbers. However since we move the multiplication into the for loop this is done every single loop. So this COULD slow things down depending on how fast the square root algorithm in java is.

Answer 4

That's a bit worse than my sieve did on a 8 Mhz 8088 in turbo pascal in 1986 or so. But that was after optimisations :)

Answer 5

Since you're searching for them in ascending order, you could keep a list of the primes you've already found and only check for divisibility against them, since all non-prime numbers can be reduced to a list of lesser prime factors. Combine that with the previous tip about not checking for factors over the square root of the current number, and you'll have yourself a pretty darn efficient implementation.

Answer 6

In C#:

class Program
{
    static void Main(string[] args)
    {
        int count = 0;
        int max = 150000;
        int i = 2;

        DateTime start = DateTime.Now;
        while (count < max)
        {
            if (IsPrime(i))
            {
                count++;
            }

            i++;

        }
        DateTime end = DateTime.Now;

        Console.WriteLine("Total time taken: " + (end - start).TotalSeconds.ToString() + " seconds");
        Console.ReadLine();
    }

    static bool IsPrime(int n)
    {
        if (n < 4)
            return true;
        if (n % 2 == 0)
            return false;

        int s = (int)Math.Sqrt(n);
        for (int i = 2; i <= s; i++)
            if (n % i == 0)
                return false;

        return true;
    }
}

Output:

Total time taken: 2.087 seconds

Answer 7

Posting this as a logical explanation as I don't know java.

All prime numbers above 3 are of the form 6n − 1 or 6n + 1, because all other numbers are divisible by 2 or 3.

Also, if there is no divisor up until the square root of a number, it's prime.

so, if we're looking for x, where 6n < sqrt(x):

if x % 6n±1 == 0 : not prime.

until you get to the the square root, add 1 to n. so let's see if 147 is prime:

sqrt(147) ==  12.1244
147%5 !=0, 147%7 == 0
> not a prime

sqrt(2027) == 45.02
2027 % 5 !=0, 2027 % 7 !=0 (this is 6±1, obviously)
2027 % 12±1 != 0
2027 % 18±1 != 0 
2027 % 24±1 != 0
2027 % 30±1 != 0
2027 % 36±1 != 0
2027 % 42±1 != 0
 -- over the square root, number is prime! --

Then, all you need to do is iterate over a list of numbers.

You can't find all primes only with this as you need to exclude non-primes under a certain threshold, which is not too hard to find.

If you can manage to fit this logic inside a prime validation routine, you should get any number you want real fast (should take sqrt(n)/6 iterations for an integer n).

Hope this helps.

Answer 8

C#

Enhancement to Aistina's code:

This makes use of the fact that all primes greater than 3 are of the form 6n + 1 or 6n - 1.

This was about a 4-5% speed increase over incrementing by 1 for every pass through the loop.

class Program
{       
    static void Main(string[] args)
    {
        DateTime start = DateTime.Now;

        int count = 2; //once 2 and 3

        int i = 5;
        while (count < 150000)
        {
            if (IsPrime(i))
            {
                count++;
            }

            i += 2;

            if (IsPrime(i))
            {
                count++;
            }

            i += 4;
        }

        DateTime end = DateTime.Now;

        Console.WriteLine("Total time taken: " + (end - start).TotalSeconds.ToString() + " seconds");
        Console.ReadLine();
    }

    static bool IsPrime(int n)
    {
        //if (n < 4)
        //return true;
        //if (n % 2 == 0)
        //return false;

        int s = (int)Math.Sqrt(n);
        for (int i = 2; i <= s; i++)
            if (n % i == 0)
                return false;

        return true;
    }
}

Answer 9

Does the re-declaration of the variable prime

        while (count < topPrime) {

            boolean prime = true;

within the loop make it inefficient? (I assume it doesn't matter, since I would think Java would optimize this)

boolean prime;        
while (count < topPrime) {

            prime = true;

Answer 10

My take at optimization, avoiding too cryptic tricks. I use the trick given by I-GIVE-TERRIBLE-ADVICE, which I knew and forgot... :-)

public class Primes
{
  // Original code
  public static void first()
  {
    int topPrime = 150003;
    int current = 2;
    int count = 0;
    int lastPrime = 2;

    long start = System.currentTimeMillis();

    while (count < topPrime) {

      boolean prime = true;

      int top = (int)Math.sqrt(current) + 1;

      for (int i = 2; i < top; i++) {
        if (current % i == 0) {
          prime = false;
          break;
        }
      }

      if (prime) {
        count++;
        lastPrime = current;
//      System.out.print(lastPrime + " "); // Checking algo is correct...
      }
      if (current == 2) {
        current++;
      } else {
        current = current + 2;
      }
    }

    System.out.println("\n-- First");
    System.out.println("Last prime = " + lastPrime);
    System.out.println("Total time = " + (double)(System.currentTimeMillis() - start) / 1000);
  }

  // My attempt
  public static void second()
  {
    final int wantedPrimeNb = 150000;
    int count = 0;

    int currentNumber = 1;
    int increment = 4;
    int lastPrime = 0;

    long start = System.currentTimeMillis();

NEXT_TESTING_NUMBER:
    while (count < wantedPrimeNb)
    {
      currentNumber += increment;
      increment = 6 - increment;
      if (currentNumber % 2 == 0) // Even number
        continue;
      if (currentNumber % 3 == 0) // Multiple of three
        continue;

      int top = (int) Math.sqrt(currentNumber) + 1;
      int testingNumber = 5;
      int testIncrement = 2;
      do
      {
        if (currentNumber % testingNumber == 0)
        {
          continue NEXT_TESTING_NUMBER;
        }
        testingNumber += testIncrement;
        testIncrement = 6 - testIncrement;
      } while (testingNumber < top);
      // If we got there, we have a prime
      count++;
      lastPrime = currentNumber;
//      System.out.print(lastPrime + " "); // Checking algo is correct...
    }

    System.out.println("\n-- Second");
    System.out.println("Last prime = " + lastPrime);
    System.out.println("Total time = " + (double) (System.currentTimeMillis() - start) / 1000);
  }

  public static void main(String[] args)
  {
    first();
    second();
  }
}

Yes, I used a labeled continue, first time I try them in Java...
I know I skip computation of the first few primes, but they are well known, no point to recompute them. :-) I can hard-code their output if needed! Beside, it doesn't give a decisive edge anyway.

Results:

-- First
Last prime = 2015201
Total time = 4.281

-- Second
Last prime = 2015201
Total time = 0.953

Not bad. Might be improved a bit, I suppose, but too much optimization can kill good code.

Answer 11

You should be able to make the inner loop twice as fast by only evaluating the odd numbers. Not sure if this is valid Java, I'm used to C++, but I'm sure it can be adapted.

            if (current != 2 && current % 2 == 0)
                    prime = false;
            else {
                    for (int i = 3; i < top; i+=2) {
                            if (current % i == 0) {
                                    prime = false;
                                    break;
                            }
                    }
            }

Answer 12

I decided to try this in F#, my first decent attempt at it. Using the Sieve of Eratosthenes on my 2.2Ghz Core 2 Duo it runs through 2..150,000 in about 200 milliseconds. Each time it calls it self it's eliminated the current multiples from the list, so it just gets faster as it goes along. This is one of my first tries in F# so any constructive comments would be appreciated.

let max = 150000
let numbers = [2..max]
let rec getPrimes sieve max =
    match sieve with
    | [] -> sieve
    | _ when sqrt(float(max)) < float sieve.[0] -> sieve
    | _ -> let prime = sieve.[0]
           let filtered = List.filter(fun x -> x % prime <> 0) sieve // Removes the prime as well so the recursion works correctly.
           let result = getPrimes filtered max
           prime::result        // The filter removes the prime so add it back to the primes result.

let timer = System.Diagnostics.Stopwatch()
timer.Start()
let r = getPrimes numbers max
timer.Stop()
printfn "Primes: %A" r
printfn "Elapsed: %d.%d" timer.Elapsed.Seconds timer.Elapsed.Milliseconds

Answer 13

I bet Miller-Rabin would be faster. If you test enough contiguous numbers it becomes deterministic, but I wouldn't even bother. Once a randomized algorithm reaches the point that its failure rate is equal to the likelihood that a CPU hiccup will cause a wrong result, it just doesn't matter any more.

Answer 14

It takes us under a second (2.4GHz) to generate the first 150000 prime numbers in Python using Sieve of Eratosthenes:

#!/usr/bin/env python
def iprimes_upto(limit):
    """Generate all prime numbers less then limit.

    http://stackoverflow.com/questions/188425/project-euler-problem#193605
    """
    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

def sup_prime(n):
    """Return an integer upper bound for p(n).

       p(n) < n (log n + log log n - 1 + 1.8 log log n / log n)

       where p(n) is the n-th prime. 
       http://primes.utm.edu/howmany.shtml#2
    """
    from math import ceil, log
    assert n >= 13
    pn = n * (log(n) + log(log(n)) - 1 + 1.8 * log(log(n)) / log(n))
    return int(ceil(pn))

if __name__ == '__main__':
    import sys
    max_number_of_primes = int(sys.argv[1]) if len(sys.argv) == 2 else 150000
    primes = list(iprimes_upto(sup_prime(max_number_of_primes)))
    print("Generated %d primes" % len(primes))
    n = 100
    print("The first %d primes are %s" % (n, primes[:n]))

Example:

$ python primes.py

Generated 153465 primes
The first 100 primes are [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, 101, 103, 107, 109, 113, 
127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197,
199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379,
383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
467, 479, 487, 491, 499, 503, 509, 521, 523, 541]

Answer 15

Here's my solution... its fairly fast... it calculates the primes between 1 and 10,000,000 in 3 seconds on my machine (core i7 @ 2.93Ghz) on Vista64.

My solution is in C, but I am not a professional C programmer. Feel free to criticize the algorithm and the code itself :)

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<time.h>

//5MB... allocate a lot of memory at once each time we need it
#define ARRAYMULT 5242880 


//list of calculated primes
__int64* primes; 
//number of primes calculated
__int64 primeCount;
//the current size of the array
__int64 arraySize;

//Prints all of the calculated primes
void PrintPrimes()
{
    __int64 i;
    for(i=0; i<primeCount; i++)
    {
     printf("%d ", primes[i]);
    }

}

//Calculates all prime numbers to max
void CalcPrime(__int64 max)
{
    register __int64 i;
    double square;
    primes = (__int64*)malloc(sizeof(__int64) * ARRAYMULT);
    primeCount = 0;
    arraySize = ARRAYMULT;

    //we provide the first prime because its even, and it would be convenient to start
    //at an odd number so we can skip evens.
    primes[0] = 2;
    primeCount++;



    for(i=3; i<max; i+=2)
    {
     int j;
     square = sqrt((double)i);

     //only test the current candidate against other primes.
     for(j=0; j<primeCount; j++)
     {
      //prime divides evenly into candidate, so we have a non-prime
      if(i%primes[j]==0)
       break;
      else
      {
       //if we've reached the point where the next prime is > than the square of the
       //candidate, the candidate is a prime... so we can add it to the list
       if(primes[j] > square)
       {
        //our array has run out of room, so we need to expand it
        if(primeCount >= arraySize)
        {
         int k;
         __int64* newArray = (__int64*)malloc(sizeof(__int64) * (ARRAYMULT + arraySize));

         for(k=0; k<primeCount; k++)
         {
          newArray[k] = primes[k];
         }

         arraySize += ARRAYMULT;
         free(primes);
         primes = newArray;
        }
        //add the prime to the list
        primes[primeCount] = i;
        primeCount++;
        break;

       }
      }

     }

    }


}
int main()
{
    int max;
    time_t t1,t2;
    double elapsedTime;

    printf("Enter the max number to calculate primes for:\n");
    scanf_s("%d",&max);
    t1 = time(0);
    CalcPrime(max);
    t2 = time(0);
    elapsedTime = difftime(t2, t1);
    printf("%d Primes found.\n", primeCount);
    printf("%f seconds elapsed.\n\n",elapsedTime);
    //PrintPrimes();
    scanf("%d");
    return 1;
}

Answer 16

@ Mark Ransom - not sure if this is java code

They will moan possibly but I wished to rewrite using paradigm I have learned to trust in Java and they said to have some fun, please make sure they understand that spec says nothing that effects ordering on the returned result set, also you would cast result set dot values() to a list type given my one-off in Notepad before taking a short errand

=============== begin untested code ===============

package demo;

import java.util.List;
import java.util.HashSet;

class Primality
{
  int current = 0;
  int minValue;
  private static final HashSet<Integer> resultSet = new HashSet<Integer>();
  final int increment = 2;
  // An obvious optimization is to use some already known work as an internal
  // constant table of some kind, reducing approaches to boundary conditions.
  int[] alreadyKown = 
  {
     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, 101, 103, 107, 109, 113, 
     127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197,
     199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
     283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379,
     383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
     467, 479, 487, 491, 499, 503, 509, 521, 523, 541
  };
  // Trivial constructor.

  public Primality(int minValue)
   {
      this.minValue = minValue;
  }
  List calcPrimes( int startValue )
  {
      // eliminate several hundred already known primes 
      // by hardcoding the first few dozen - implemented 
      // from prior work by J.F. Sebastian
      if( startValue > this.minValue )
      {
          // Duh.
          current = Math.abs( start );
           do
           {
               boolean prime = true;
               int index = current;
               do
               {
                  if(current % index == 0)
                  {
                      // here, current cannot be prime so break.
                      prime = false;
                      break;
                   }
               while( --index > 0x00000000 );

               // Unreachable if not prime
               // Here for clarity

               if ( prime )
               {     
                   resultSet dot add ( or put or whatever it is )
                           new Integer ( current ) ;
               }
           }
           while( ( current - increment ) > this.minValue );
           // Sanity check
           if resultSet dot size is greater that zero
           {
              for ( int anInt : alreadyKown ) { resultSet.add( new Integer ( anInt ) );}
             return resultSet;
           }
           else throw an exception ....
      }

=============== end untested code ===============

Using Hash Sets allows searching results as B-Trees, thus results could be stacked up until the machine begins to fail then that starting point could be used for another block of testing == the end of one run used as a Constructor value for another run, persisting to disk work already accomplished and allowing incremental feed-forward designs. Burnt out right now, loop logic needs analysis.

patch (plus add sqrt) :

  if(current % 5 == 0 )
  if(current % 7 == 0 )
  if( ( ( ( current % 12 ) +1 ) == 0) || ( ( ( current % 12 ) -1 ) == 0) ){break;}
  if( ( ( ( current % 18 ) +1 ) == 0) || ( ( ( current % 18 ) -1 ) == 0) ){break;}
  if( ( ( ( current % 24 ) +1 ) == 0) || ( ( ( current % 24 ) -1 ) == 0) ){break;}
  if( ( ( ( current % 36 ) +1 ) == 0) || ( ( ( current % 36 ) -1 ) == 0) ){break;}
  if( ( ( ( current % 24 ) +1 ) == 0) || ( ( ( current % 42 ) -1 ) == 0) ){break;}


// and - new work this morning:


package demo;

/**
 *
 * Buncha stuff deleted for posting .... duh.
 *
 * @author  Author
 * @version 0.2.1
 *
 * Note strings are base36
 */
public final class Alice extends java.util.HashSet<java.lang.String>
{
    // prints 14551 so it's 14 ½ seconds to get 40,000 likely primes
    // using Java built-in on amd sempron 1.8 ghz / 1600 mhz front side bus 256 k L-2
    public static void main(java.lang.String[] args)
    {
        try
        {
            final long start=System.currentTimeMillis();
            // VM exhibits spurious 16-bit pointer behaviour somewhere after 40,000
            final java.lang.Integer upperBound=new java.lang.Integer(40000);
            int index = upperBound.intValue();

            final java.util.HashSet<java.lang.String>hashSet
            = new java.util.HashSet<java.lang.String>(upperBound.intValue());//
            // Arbitraily chosen value, based on no idea where to start.
            java.math.BigInteger probablePrime
            = new java.math.BigInteger(16,java.security.SecureRandom.getInstance("SHA1PRNG"));
            do
            {
                java.math.BigInteger nextProbablePrime = probablePrime.nextProbablePrime();
                if(hashSet.add(new java.lang.String(nextProbablePrime.toString(Character.MAX_RADIX))))
                {
                    probablePrime = nextProbablePrime;
                    if( ( index % 100 ) == 0x00000000 )
                    {
                        // System.out.println(nextProbablePrime.toString(Character.MAX_RADIX));//
                        continue;
                    }
                    else
                    {
                        continue;
                    }
                }
                else
                {
                    throw new StackOverflowError(new String("hashSet.add(string) failed on iteration: "+
                    Integer.toString(upperBound.intValue() - index)));
                }
            }
            while(--index > 0x00000000);
            System.err.println(Long.toString( System.currentTimeMillis() - start));
        }
        catch(java.security.NoSuchAlgorithmException nsae)
        {
            // Never happen
            return;
        }
        catch(java.lang.StackOverflowError soe)
        {
            // Might happen
            System.out.println(soe.getMessage());//
            return;
        }
    }
}// end class Alice

Answer 17

Here is my take on it. The program is writtern in C and takes 50 milliseconds on my laptop(Core 2 Duo, 1 GB Ram). I am keeping all the calculated primes in an array and trying divisibility only till sqrt of number. Of course, this doesnt work when we need very large number of primes(tried with 100000000) as array grows too big and gives seg fault.

/*Calculate the primes till TOTALPRIMES*/
#include <stdio.h>
#define TOTALPRIMES 15000

main(){
int primes[TOTALPRIMES];
int count;
int i, j, cpr;
char isPrime;

primes[0] = 2;
count = 1;

for(i = 3; count < TOTALPRIMES; i+= 2){
    isPrime = 1;

    //check divisiblity only with previous primes
    for(j = 0; j < count; j++){
        cpr = primes[j];
        if(i % cpr == 0){
            isPrime = 0;
            break;
        }
        if(cpr*cpr > i){
            break;
        }
    }
    if(isPrime == 1){
        //printf("Prime: %d\n", i);
        primes[count] = i;
        count++;
    }


}

printf("Last prime = %d\n", primes[TOTALPRIMES - 1]);
}

$ time ./a.out 
Last prime = 163841
real    0m0.045s
user    0m0.040s
sys 0m0.004s

Answer 18

I found this code somewhere on my machine when I started reading this blog entry about prime numbers. The code is in C# and the algorithm I used came from my head although it is probably somewhere on Wikipedia. ;) Anyway, it can fetch the first 150000 prime numbers in about 300ms. I discovered that the sum of the n first odd numbers is equal to n^2. Again, there is probably a proof of this somewhere on wikipedia. So knowing this, I can write an algorithm that wil never have to calculate a square root but I have to calculate incrementally to find the primes. So if you want the Nth prime, this algo will have to find the (N-1) preceding primes before! So there it is. Enjoy!

//
// Finds the n first prime numbers.
//
//count: Number of prime numbers to find.
//listPrimes: A reference to a list that will contain all n first prime if getLast is set to false.
//getLast: If true, the list will only contain the nth prime number.
//
static ulong GetPrimes(ulong count, ref IList listPrimes, bool getLast)
{
    if (count == 0)
        return 0;
    if (count == 1)
    {
        if (listPrimes != null)
        {
            if (!getLast || (count == 1))
                listPrimes.Add(2);
        }

        return count;
    }

    ulong currentSquare = 1;
    ulong nextSquare = 9;
    ulong nextSquareIndex = 3;
    ulong primesCount = 1;

    List dividers = new List();

    //Only check for odd numbers starting with 3.
    for (ulong curNumber = 3; (curNumber  (nextSquareIndex % div) == 0) == false)
                dividers.Add(nextSquareIndex);

            //Move to next square number
            currentSquare = nextSquare;

            //Skip the even dividers so take the next odd square number.
            nextSquare += (4 * (nextSquareIndex + 1));
            nextSquareIndex += 2;

            //We may continue as a square number is never a prime number for obvious reasons :).
            continue;
        }

        //Check if there is at least one divider for the current number.
        //If so, this is not a prime number.
        if (dividers.Exists(div => (curNumber % div) == 0) == false)
        {
            if (listPrimes != null)
            {
                //Unless we requested only the last prime, add it to the list of found prime numbers.
                if (!getLast || (primesCount + 1 == count))
                    listPrimes.Add(curNumber);
            }
            primesCount++;
        }
    }

    return primesCount;
}