Sieve of Eratosthenes

[Gauri Shankar Iyer]

My first attempt at writing the algo for finding primes in a range using the Sieve of Eratosthenes.
(Wikipedia: http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes).

I'm wondering what are the bad practices you can see in this small piece of code, and any optimizations you see fit? Also, looking to reduce the number of lines of code in any way possible. Some criticism will be very helpful. Thanks.

Some observations you might find interesting:

- For the range (1, 500000), it takes about 5 seconds with the code below, on my AMD dual core machine.
- On line 33, changing the Math.pow(variable, 2) to (variable * variable) reduces the running time to 2 seconds!
- The reduction of lines 59, 60, 61 into line 62 has no effect whatsoever. Removing explicit allocations doesn't necessarily mean implicit VM allocations will be removed too.

import java.util.ArrayList;
import java.util.Iterator;
import java.util.List;


public class Sieve {

	public static void main(String[] args) {
		
		/* small test */
		int[] primes = sieve(1, 100);
		for (int i : primes) {
			System.out.println(i);
		}
	}
	
	private static int[] sieve(int start, int end) {
		
		/* some basic range checks */
		if(start < 1 || end < start || start == end) {
			throw new ArithmeticException("Messed up input");
		}
		
		/* generate ints within range */
		int[] primes = new int[end-start+1];
		for (int j = 0; j < end - start + 1; j++) {
			primes[j] = start + j;
		}
		
		/* init running prime to start, and increment until
		 * running prime squared is greater than the end
		 */
		for (int runningPrime = (start == 1 ? 2: start); end > Math.pow(runningPrime, 2); runningPrime++) {
			for (int i = runningPrime; i < primes.length; i++) {
				if(-1 != primes[i]) {
					if(primes[i] % runningPrime == 0) {
						primes[i] = -1;
					}
				}
			}
		}
		
		/* handling that twat */
		if(primes[0] == 1) {
			primes[0] = -1;
		}
		
		/* list primes */
		List list = new ArrayList();
		for (int i = 0; i < primes.length; i++) {
			if(-1 != primes[i])
				list.add(primes[i]);
		}
		
		/* create the return int array */
		int[] primesFound = new int[list.size()];
		int k = 0;
		for (Iterator iterator = list.iterator(); iterator.hasNext();) {
			//Integer i = (Integer) iterator.next();
			//primesFound[k] = i.intValue();
			//k++;
			primesFound[k++] = ((Integer) iterator.next()).intValue();
		}
		
		return primesFound;
	}
}

[Hunter McMillen]

any sort of built-in operators (*,/,+,-) are always going to be less expensive to use than recursive functions like Math.pow, because functions have to be called, execute their own code and possibly call other functions, then they have to return. And for a smaller program you can never tell the difference. But when you are trying to do 500,000 calculations, that overhead adds up quickly.


Hunter

[Campbell Ritchie]

Always use x * x rather than Math.pow(x, 2). It is more precise, as well as faster.

[Lester Burnham]

"end < start || start == end" is the same as "end <= start"

[Campbell Ritchie]

It is a long time since I wrote a sieve application, but yours looks much more complicated than I remember mine.

[Matthew Brown]

Interesting comment on line 43 .

[Gauri Shankar Iyer]

Matthew Brown wrote:Interesting comment on line 43 .

Ha ha, isn't that a remarkable trait of most programmers?

[Campbell Ritchie]

I wrote mine with a boolean[] array. Set all the values to true (except 0 and 1, which are implicitly false) and change all the places found with the sieve to false. Start looping at 2. It is easier to have an array 1 size larger than your range, so you can use the index to represent the number.

[Paul Beckett]

see problem 10 "Calculate the sum of all the primes below two million" on projecteuler.net. If you can solve the problem then a solution pdf is provided giving hints (including pseudo code) for how to efficiently do a sieve (solves in much less than a second on a modest PC). Your code should be more than sufficient to solve the problem and is way more efficient than my first attempt where I didn't even use a sieve. The solution pdf will explain things much more eloquently than I could here.

Note: Campbell is correct to use a boolean[] but this can be optimized even further thus reducing the size of the array.
(i) no even number other than 2 can be a prime
(ii) can you reduce the max number you need to check for?

[Campbell Ritchie]

I tried it with a boolean[] array yesterday, taking about 3 seconds to sieve from 2 to 100000000 (1e8), but it took about 1 minute to display all those numbers on screen. For 2 to 1000000000 (1e9) I suffered an OutOfMemoryError.

[Campbell Ritchie]

Paul Beckett wrote: . . . Campbell is correct . . .

Can I have that in writing please

can you reduce the max number you need to check for?

I think the squaring allows for that. But surely it is simpler to use the square root?

[Paul Beckett]

I think the squaring allows for that. But surely it is simpler to use the square root?

ok, I'd missed that bit
Using the square root would be more efficient in this case. So work out the square root of the end number (crosslimit) before the loop and then the termination expression becomes "end<crosslimit".

Using a boolean array you can change your increment expression to increment the loop by the number you are currently removing multiples of (runningPrime in the original code).

My algorithm takes <1.5 seconds to work out all the primes 1e8 but also gets an OutOfMemoryError with 1e9. Obviously the time taken is relative to the performance of the machine.

My version of the sieve:

 		int limit = maxNumber;
		int sieveBound = limit / 2 // Paul Beckett: corrected from (limit - 1) / 2;
		boolean[] sieve = new boolean[sieveBound];
		int crossLimit = (int) Math.floor(Math.sqrt(limit));

		for (int i = 1; i < crossLimit; i++) {
			if (!sieve[i]) {
				for (int j = 2 * i * (i + 1); j < sieveBound; j = j
						+ (2 * i + 1)) {
					sieve[j] = true;
				}
			}
		}

[Campbell Ritchie]

Paul Beckett wrote: . . . the time taken is relative to the performance of the machine. . . .

I think it runs in quadratic time.

Are you using false == prime? If so, what about 0, which is an exact multiple of every number, but is not a "natural" number? And what about 1, which doesn't divide by any other numbers, but is not usually considered prime?
I seemed to get a different result from yours when I printed the numbers out.

~/java$ java Eratosthenes
Time for sieving 0.001s Please push "enter" to continue.
Primes found:
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

Total of all those primes = 1060

Duration by Paul's method = 0.002s
Please push "enter" to continue.
Primes found by Paul's method:
4 7 10 12 13 16 17 19 22 24 25 27 28 31 32 34 37 38 40 42 43 45 46 47

Total of all those primes = 659

import static java.lang.System.*;
import java.util.Scanner;

public class Eratosthenes
{
    public static void main(String[] args)
    {
        int i = args.length == 0 ? 100 : Integer.parseInt(args[0]);
        Eratosthenes er = new Eratosthenes(i);
        er.printPrimes();
        er.new Sieve(i).printPrimes();
    }
    
    private final boolean[] primes;
    
    public Eratosthenes(int largest)
    {
        long duration = currentTimeMillis();
        primes = new boolean[largest + 1];
        for (int i = 2; i < primes.length; i++)
        {
            primes[i] = true;
        }
        sieve(2);
        for (int i = 3; i <= Math.sqrt(primes.length); i += 2)
        {
            if (primes[i])
            {
                sieve(i);
            }
        }
        out.print("Time for sieving ");
        {
            long i = currentTimeMillis();
            if (i > duration)
            {
                duration = i - duration;
                out.printf("%d.%03ds ", duration / 1000, duration % 1000);
            }
            else
            {
                out.print("< 1ms. ");
            }
        }
        out.print("Please push \"enter\" to continue.");
        new Scanner(in).nextLine();
    }
    
    private void sieve(int i)
    {
        for (int j = i * 2; j < primes.length; j += i)
        {
            primes[j] = false;
        }
    }
    
    public void printPrimes()
    {
        long total = 0;
        out.println("Primes found:");
        for (int i = 0; i < primes.length; i++)
        {
            if (primes[i])
            {
                out.printf("%8d", i);
                total += i;
            }
        }
        out.printf("%n%nTotal of all those primes = %d%n%n", total);
    }
    
    /**
     * Paul Beckett's technique, I believe copied correctly from this forum
     */
    private class Sieve
    {
        private boolean[] sieve;
        
        public Sieve(int maxNumber)
        {
            long duration = currentTimeMillis();
     		int limit = maxNumber;
            int sieveBound = (limit - 1) / 2;
            sieve = new boolean[sieveBound];
            int crossLimit = (int) Math.floor(Math.sqrt(limit));
            for (int i = 1; i < crossLimit; i++) {
            if (!sieve[i]) {
                for (int j = 2 * i * (i + 1); j < sieveBound; j = j
                        + (2 * i + 1)) {
                    sieve[j] = true;
                }
            }
        }
            duration -= currentTimeMillis();
            out.printf("Duration by Paul's method = %d.%03ds%n",
                    -duration / 1000, -duration % 1000);
            out.print("Please push \"enter\" to continue.");
            new Scanner(in).nextLine();
        }
    
        public void printPrimes()
        {
            long total = 0;
            out.println("Primes found by Paul's method:");
            for (int i = 0; i < sieve.length; i++)
            {
                if (sieve[i])
                {
                    out.printf("%8d", i);
                    total += i;
                }
            }
            out.printf("%n%nTotal of all those primes = %d%n%n", total);// As suggested by Euler website.
        }
    }
}

[Kurt Van Etten]

Campbell Ritchie wrote:I tried it with a boolean[] array yesterday, taking about 3 seconds to sieve from 2 to 100000000 (1e8), but it took about 1 minute to display all those numbers on screen. For 2 to 1000000000 (1e9) I suffered an OutOfMemoryError.

Well, that raises some interesting questions... How much memory does the Java VM use to store a boolean value? I'm guessing one byte on most platforms. (Does Java have anything similar to the C sizeof operator? It would be nice to be able to estimate available resources before attempting a memory allocation.) In that case, you could increase the range of your sieve by using bitwise operators to store 8 0/1 values per byte, with a bit of a performance trade-off.

Of course that would only increase your range by a factor of 8. To get a truly big sieve, you'd need to write the values out to disk. With some judicious buffering and separate threads to handle the I/O, the performance should still be reasonable.

[Paul Beckett]

Campbell Ritchie wrote:Are you using false == prime? If so, what about 0, which is an exact multiple of every number, but is not a "natural" number? And what about 1, which doesn't divide by any other numbers, but is not usually considered prime?
I seemed to get a different result from yours when I printed the numbers out.

I didn't include it before but my summing is done like so:

		long sum = 2;
		for (int i = 1; i < sieveBound; i++) {
			if (!sieve[i]) {	
				sum = sum + (2 * i + 1);
			}
		}
		System.out.println("sum=" + sum);    

So a prime is false, don't think it matters which way round you do it. With prime = false then I don't need to explicitly set every value in the array to true before sieving.
The array doesn't contain even numbers so have to handle 2 as a special case. Other than that the value in the array at index i represents the number 2i+1, so index 3 represents the number 7.
Where maxNumber=12:
sieve: [false, false, false, false, true]
index: [ 0, 1, 2, 3, 4]
numbers: [ 1, 3, 5, 7, 9]
Obviously 1 is not a prime number so notice that when iterating the sieve always start at index 1.

Also my sieveBound should have been limit/2 instead of "(limit-1)/2".

[Campbell Ritchie]

So what did I copy wrongly that I got your code to produce that series of numbers?

[Paul Beckett]

you didn't copy anything incorrectly. The implementation you have used of my version of the sieve is correct I think. It is the printPrimes method that is incorrect. You are printing out the index of any "true" in the array. Whereas in my implementation of the sieve, "false" is a prime. And then instead of the index of the "false" elements you need to use 2i+1 to find the real number and start iterating at index 1.

Obviously you have identified a big flaw in the method I've used - its difficult to understand what is going on!!! However this is a java implementation of the method from Project Euler problem 10 solution pdf (section 10.2.2). I prefer to make code readable as a first priority and then only reduce readability to optimise if absolutely necessary but for project euler execution speed seems to be the priority over all others.

[Campbell Ritchie]

Please post working code, which we can actually execute and understand the output of. You ought to have posted the printing method as well, so we can see a set of prime numbers. If I had seen an output of all non-prime numbers, I would quickly have worked out that I was printing the wrong half of the array.
I don't think Project Euler do want speed above all else. I think they want us to think up a correct solution.

[Paul Beckett]

Campbell Ritchie wrote:Please post working code, which we can actually execute and understand the output of. You ought to have posted the printing method as well, so we can see a set of prime numbers

The sieve I posted in my second post is a fully functioning sieve minus the method declarations.
My version of the sieve was for solving the project euler problem 10 and doesn't need a printing method for that purpose. In my third post I did provide the complete summing functionality which shows how to convert an array element to a prime number. Also I provided a "diagram" of the sieve, the array indexes, and the corresponding number the array index represented.

[Carey Brown]

I'll just add my 2cents. Using a BitSet should help with memory issues and appears to be plenty fast. On my machine (your milage may vary) I ran 500,000 numbers in 0.02 seconds.

import java.util.BitSet;


public class Primes
{
	public static void main( String[] args )
	{
		long			tmp = Long.parseLong( args[0] );

		if( tmp > Integer.MAX_VALUE )
		{
			System.err.println( "TOO BIG - maximum must be less than or equal to " + Integer.MAX_VALUE );
			System.exit( 1 );
		}

		int				max = (int)tmp;
		int				sqrt = (int)Math.sqrt( max );
		BitSet			bitSet = new BitSet( max + 1 );
		int				i, j;
		long			start, stop;

		start = System.currentTimeMillis();

		//+++++ the meat and potatoes +++++

		i = 2;
		for( j=i*2 ; j <= max ; j+=i )
			bitSet.set( j );

		for( i=3 ; i <= sqrt ; i+=2 )
			if( ! bitSet.get(i) )
				for( j=i*2 ; j <= max ; j+=i )
					bitSet.set( j );

		//+++++++++++++++++++++++++++++++++

		stop = System.currentTimeMillis();

		System.out.format( "N=%,d sec=%.2f\n", max, ((stop - start)/1000.0) );

		for( i=2 ; i < 200 ; i++ )
			if( ! bitSet.get(i) )
				System.out.format( "%d ", i );
	}
}

[Mike Simmons]

Carey Brown wrote:I'll just add my 2cents. Using a BitSet should help with memory issues and appears to be plenty fast.

Yeah - boolean[] arrays seem to often be implemented using 8 bits for each boolean. The boolean array is faster, but not by enough to justify using eight times as much memory (IMO).

Carey Brown wrote: On my machine (your milage may vary) I ran 500,000 numbers in 0.02 seconds

On my machine, it ran in 0.02 sec as well. And the same code can run 500,000,000 numbers in 18-19 sec. (Going to 1,000,000,000 would require changing the standard JVM memory settings, so I stopped there.) I got it down to 11-12 sec by replacing the second for loop with this:

    for (i = 3; i <= sqrt; i += 2)
      if (!bitSet.get(i)) {
        int inc = i * 2;
        for (j = i * 3; j <= max; j += inc) 
          bitSet.set(j);
      }

As long as you've already handled all the multiples of 2, we might as well only deal with the odd multiples of i at this point.

A next step might be to eliminate all multiples of 2 or 3. My guess is that might reduce the time by about 1/6 (since the only new skips are for multiples of 3 which aren't already multiples of 2, like 9, 15, 21, 27, 33... etc). Making the code somewhat more complex to boot. Not sure it's worth it. But I thought I'd toss it out there; maybe someone can use the idea to good effect.

[Kurt Van Etten]

Carey Brown wrote:I'll just add my 2cents. Using a BitSet should help with memory issues and appears to be plenty fast. On my machine (your milage may vary) I ran 500,000 numbers in 0.02 seconds.

I was curious about how big of a time hit this would make, so I switched my previous version to use a BitSet instead of an array (keeping everything else the same), and it took about 2.5 times longer to run. So I suppose you have to weigh whether the time or memory is more important to you. Now I wonder how much of the time hit is due to the overhead for the function calls. I ought to try implementing the bit operations directly in an array and see how much of an improvement that makes.

[Kurt Van Etten]

Mike Simmons wrote:A next step might be to eliminate all multiples of 2 or 3. My guess is that might reduce the time by about 1/6 (since the only new skips are for multiples of 3 which aren't already multiples of 2, like 9, 15, 21, 27, 33... etc). Making the code somewhat more complex to boot. Not sure it's worth it. But I thought I'd toss it out there; maybe someone can use the idea to good effect.

I tried this out. Instead of adding 2 times the current prime each time you mark an element in the sieve, you alternate adding 4 times and 2 times the current prime. Unfortunately the time was essentially unchanged--even though I needed to mark 1/3 fewer elements in the sieve for each prime greater than 3, the time savings was apparently offset by the slight increase in complexity.

[Mike Simmons]

Yeah, I was getting very similar results, with the same skip-2/skip-4 idea. Seems like the added calculations to figure out which elements to hit offset the savings of hitting fewer elements. Annoying. I was looking for cleaner ways to refactor it, but not seeing much improvement. Oh well.

[Paul Beckett]

The BitSet is more efficient in terms of memory usage - no argument there. However both BitSet and boolean[] will hit the limit soon enough of being indexed by an int.

Arbitrarily allocating Xmx768m, my boolean[] can handle 1e9 in 15s with the BitSet in 36s. The point being 1e9 is approximately half of Integer.MAX_VALUE. So why did I stop at 1e9 and not go on to Integer.MAX_VALUE? My implementation doesn't contain even numbers in the sieve so requires an i*2 calculation to work out the real number from the index (see earlier posts) so 1e9 is approximately the max value I can handle without recoding. So with recoding then memory requirement would be fulfilled by Xms1536m. I'm not advocating always chucking more memory at a problem but in this instance IMHO the inherent Integer.MAX_VALUE limit is more of a concern than the memory efficiency.