A fast approach to prime number sieving (non-threaded array)

Question

I am posting this question for input on comparisons between other sieving methods and opinions on what I have discovered. I believe this might be the fastest (singular, non-threaded) array implementation of prime sieving to date, correct me if I am wrong.

My rudimentary benchmark tests stand as follows:

1. My optimized code: 4 seconds N = 1,000,000,000

2. My original code: 7 seconds N = 1,000,000,000

3. Standard Sieve-of-Eratosthenes: 16 seconds N = 1,000,000,000

I have seen PrimeSieve and as I understand it utilizes better means to compute faster, but as it seems - it must dump the data somewhere and retrieve it to utilize it again, because it cycles the computations through the cache only.

My original thread with same pseudo-code is here, however, my original Java code did not match the pseudo-code I developed (but essentially performed similar functions at a lesser speed and easier complexity to code).

Pseudo-code as follows:

1. Initialize list with initial values \$(2, 3)\$.
2. Use the formula \$6k ± 1\$ where \$k = 1\$ and must increment by \$1\$. Store values to list.
3. Compute step two up to a product \$n\$, where \$n > 7\$ if \$k = 1\$.
4. Starting with \$x = 5\$, remove all products of \$x \times y\$ (up to \$n\$), where \$y\$ starts at \$x\$ and
   increments to the next number in the list as \$x\$ remains the same.
5. Repeat step four by assigning \$x\$ to the next number in the list. Stop when \$x \times y > n\$.
6. Final step: We must remove all factors of the square of every prime number \$(25, 49)\$

//-----------------------------------------------------------------
//  Class:   SixesSieve.java
//  Author:  Alex Lieberman
//  Purpose: Finds all prime numbers up to a specified
//           maxNumber (exclusive) and marks them as true in the boolean array.
//-----------------------------------------------------------------

public class SixesSieve {    
  public static void main (String[] args) { 

   int maxNumber = 1000000000;
   boolean[] isPrime = new boolean [maxNumber];
   isPrime[2] = true;
   isPrime[3] = true;

   for (int i = 5, j = 7; i < maxNumber; i+=6, j+=6) {
     isPrime[i] = true;
     if (j < maxNumber)       
      isPrime[j] = true;
    }

   for (int a = 5; a*a <= maxNumber; a+=2) {
     if (isPrime[a]) {

       for (int b = a, z = a*a; b*a <= maxNumber; ) {
        if(b*a == z && z < maxNumber/2) {
         while (true){
           if (z <= maxNumber) {
            isPrime[z] = false; 
            z = z+(a*a*2);
           }
           else
            {b+=2; break;}
           }
         }   
        else {
         isPrime[b*a] = false;
         b+=2;
         while (!isPrime[b])
         {b+=2;}    
        }  
       }
      }   
    }

    for (int i = 0; i < maxNumber; i++) {
      if (isPrime[i])
       System.out.print(i + " ");
    }
  }
 }

Answer 1

Just some minor comments to the code style. Let's start with your fileheader...in these times and days it should be unnecessary. First it doesn't make sense to have the file/classname in the comment, it's only something that will get out of date because refactoring can't reach it. Second, the other parts should be in the JavaDoc of the class.

---

Variable names. I know it is hard to find good names for such...mathematical stuff, but try to find them. 'a', 'b', 'z' are meaningless. It is hard to read the code if variables have meaningless names.

---

The code formatting is partly a bloody mess. If you use an IDE, make use of these formatting features. It is so much easier to read code if it is well formatted.

First, don't omit braces, even if it is just a one line 'if' or 'for' block. Second, be consistent about your brace style, and be consequent.

---

System.out.print(i + " ");

I still don't understand why the Java compiler does allow this. This is not a "cast", this is a hack in the best of cases.

System.out.print(Integer.toString(i) + " ");

Answer 2

Some comments on the code:

---

Don't write a big main, this makes the code unusable for anything but the current demo. Whenever main gets more than a few lines, extract a method. Btw., Java loves short methods and optimizes them way better than long ones.

---

Move constants out of methods. Actually, int maxNumber = 1000; should be an argument. Then a comment like

// maxNumber must be an even number

is obviously not enough. You need a Javadoc. And then don't be afraid of

if (maxNumber % 2 != 0) {
    throw new IllegalArgumentException("maxNumber must be an even");
}

as trying to process wrong input is waste of time at best.

---

Try to separate the output from the computation, even for small programs. This is sometimes more work, but you get a reusable result ad it helps you write cleaner code.

---

If you're interested in optimization, be assured that it's hard. And also fun. In Java it's especially complicated and usually you need some benchmarking framework (they're free).

---

The code is multiplying possibly big numbers, but I can't see any guarantee that they don't overflow. When I increased the limit to 1.3e9, problems came. That's not that bad as it's close enough to the maximum possible input (so adding a check would do).

The code is badly missing a test. The idea is good (though I still don't understand it fully), but it may produce wrong results. Testing each prime below Integer.MAX_VALUE via BigInteger.isProbablePrime should be doable.

---

Some code pieces should be improved, e.g.

secondPrime += 2;
while (!isPrime[secondPrime]) {
    secondPrime += 2;
}

is a do-while loop. This piece

while (true) {
    if (square <= maxNumber) {
        isPrime[square] = false;
        square += (prime * prime * 2);
    }
    else {
        secondPrime += 2;
        break;
    }
}

is nothing but

while (square <= maxNumber) {
    isPrime[square] = false;
    square += (prime * prime * 2);
}
secondPrime += 2;

---

After a few similar changes and some debugging, I finally understood what it's all about.

It's pretty fast and clever. The speed comes from testing only

• products p*q with p being a prime, isPrime[q] being true¹, and p<q
• multiples of p**2 for a prime p

Finding this out from the code alone is a like the twelve labours of Hercules. No offense, this happens quite often with optimized code. That's what CR is for.

---

¹ Because of q being a either prime or not yet detected composite number.

Answer 3

The revisions you've already made (variable names, and especially consistent indentation of {, }, and the things between them, helps a lot.

It's still better to organize the code in several functions instead of one main() function. Among other things, this helps reduce the required indentation of some of the code. But mainly it makes functionality a lot easier to see. At a minimum, the code that does the sieving and the code that prints the results should each be in independent functions.

Consider this code fragment:

                    secondPrime += 2;
                    while (!isPrime[secondPrime]) {
                        secondPrime += 2;
                    }

This appears to serve the function of changing secondPrime to the next larger prime. You could put something like this in the body of a function which you would then call like this:

                    secondPrime = getPrimeAfter(secondPrime);

In fact, maybe you could use this to write

            for (int secondPrime = prime, square = prime * prime; 
                 secondPrime * prime <= maxNumber;
                 secondPrime = getPrimeAfter(secondPrime)) {

which would make the for statement a bit better organized. (I am always uncomfortable when the "advance to the next iteration value" code is not between the second ; and the ) of the for loop control structure. The only question about this is that in the other main branch of logic within that same for loop, you wrote

                            secondPrime += 2;

... that is, in this case you don't increase secondPrime to the next prime, but rather just to the next odd number. The fact that you chose secondPrime as a name of this variable, however, suggests that you didn't actually mean to stop at the next odd number when that number isn't prime, so maybe you do really want to increase secondPrime in the same way on every iteration of the for (int secondPrime loop. Defining and using a function such as getPrimeAfter(int) will help avoid this error (if it is an error) and make it clearer what is actually supposed to happen.

It looks like the variable named square is actually multipleOfSquare; that is, it takes on the values prime*prime, 3*prime*prime, 5*prime*prime, and so forth.

Consider this line:

                if (secondPrime * prime == square && square < maxNumber / 2) {

If secondPrime is, in fact, a prime, then this condition can be true only on the first iteration, when secondPrime == prime. That would suggest that it might be better to execute the code of this if branch just once before the loop and only do the else code during the loop. (On the other hand, since the code as it stands at this moment does not appear to guarantee that secondPrime is prime, it's hard to determine exactly what would happen if you reorganized that loop as I just suggested.)

These three comments have me a bit confused:

                // if the product of 'prime' and 'secondPrime' equals
                // 'square'

                    // searches for all multiples of square that may exist in
                    // array

                // if 'prime' and 'secondPrime' are not equal to square
                // then search for all multiples of 'prime' by 'secondPrime'

The first comment is pretty useless, since it just says exactly what the first part of the if condition below says. In the second comment, we know that secondPrime * prime == square, so actually this is searching for all multiples of secondPrime*prime. But that's just what the third comment says. So do the if and else branches of if (secondPrime * prime == square (the branch with the second comment and the branch with the third comment) actually do anything different from each other? Do you really need two branches? If so, that's what you should explain in the comments.

You might be able to save some computation here:

    for (int prime = 5; prime * prime <= maxNumber; prime += 2) {
        // the if-statement below skips the testing of all non-prime numbers
        if (isPrime[prime]) {

If you compute the square root of maxNumber just once at the very start of the sieve, you could write prime <= maxSquareRoot instead of prime * prime <= maxNumber, and save a multiplication on every iteration. Another question is whether you could use prime = getPrimeAfter(prime) instead of prime += 2. Before doing that, you should show that getPrimeAfter(prime) will always terminate before its internal iterator runs off the end of the isPrime array. But come to think of it, you really should also establish that fact in the place where you write

                    while (!isPrime[secondPrime]) {

deep within the bowels of your sieve. On the other hand, it might slow your code down a tiny bit if you change that line to

                    while (secondPrime <= maxNumber && !isPrime[secondPrime]) {

but it would guarantee no running off the end of the array, and it would take far less development time to achieve that guarantee. Alternatively, you could allocate an additional element in the isPrime array that would act as a "guard byte": that is, you allocate space for isPrime[maxNumber + 1] and set isPrime[maxNumber + 1] = true, and rely on your other logic to ensure that you will never set isPrime[maxNumber + 1] = false (and also never actually use maxNumber + 1 as a prime number).

Answer 4

Regarding "I believe this might be the fastest (singular, non-threaded) array implementation of prime sieving to date", I disagree, at least for C implementations.

First, never use printing for benchmarking unless that is specifically what you want to test. For good sieves, the I/O will take more time than the sieving. A custom written buffered syswrite can be far faster than printf or cout, and I imagine there is something analogous in Java. I can beat primesieve by over 2x if we include printing, merely because that isn't something primesieve has optimized at all (it's irrelevant for reasonable amounts of output so there doesn't seem like there is a need).

Counting can be also be done faster than you're doing it. Clearly using advanced prime counting methods are orders of magnitude faster than any sieve, but I mean that you can use methods to count bytes or words which can be faster, and most sieve examples trying to show fast performance do this.

What other sieves did you test? How did you test them to reach this conclusion?

It's possible I messed something up in the translation to C, but here are my results for counting to 10¹¹, using the June 26th code, all single threaded:

Monolithic:

423s   Lieberman
416s   Trivial odd SoE
234s   Mathisen 1998
147s   MPU 0.36 monolithic SoE
 96s   Pari 2.6.2 monolithic

Segmented:

208s   Walisch byte segment example
104s   Oliveira e Silva v1 example
 92s   Walisch bit segment example
 81s   tverniquet
 79s   Flammenkamp
 78s   Mario Roy algorithm 3
 55s   Oliveira e Silva v2 example
 43s   Math::Prime::Util
 34s   primegen 0.97 (tuned cache size)
 23s   primesieve 5.0
 20s   yafu

Note that this is just at one size, and starting from 2. Doing ranges will be vastly more efficient with the segmented sieves, and doing counts to larger values starts to choke the monolithic sieves. The segmented sieves of course use less memory. At 10¹² yafu and primesieve are at 5 minutes; primegen, MPU, and OeS's v2 are all at 10 minutes; while the others are at 16+ minutes or worse.

Answer 5

Updated the code with as many suggestions as I felt optimized.

I tested it a lot of different ways and some ways it was slower - even with more methods.

In addition - I modified the initializing loop to skip all products of five.

import java.util.*;
import java.lang.*;
import java.math.*;

public class SixSieveFinal {

  static int maxNumber = 1000000000;
  static boolean[] isPrime = new boolean[maxNumber];
  static int limit = (int) Math.sqrt(maxNumber);
  static long timeStart, timeEnd, timeDifference;
  static double timeElapsed;
  static int totalPrimes;

    public static void main(String[] args) {

      isPrime[2] = true;
      isPrime[3] = true;

      System.out.println("Starting Computation.");
      timeStart = System.currentTimeMillis();

      getPotentialPrimes();
      eliminateProducts();

      timeEnd = System.currentTimeMillis();
      timeDifference = timeEnd - timeStart;
      timeElapsed = timeDifference / 1000.000000;

      System.out.println("Finished.");
      System.out.println("Time elapsed: " + timeElapsed);

      totalPrimes = countPrimes(totalPrimes);
      System.out.println("Total primes: " + totalPrimes);

      queryPrint();

    }

    static void getPotentialPrimes() {

      for (int countBySix = 6, countToSkip = 2; countBySix < maxNumber;) {

        if (countToSkip == 5) {

            isPrime[(countBySix) - 1] = true;
            countToSkip++;

        } else if (countToSkip == 7) {

            isPrime[(countBySix) + 1] = true;
            countToSkip = 3;

        } else {

            isPrime[(countBySix) - 1] = true;
            isPrime[(countBySix) + 1] = true;
            countToSkip++;

        }

        countBySix += 6;

      }

    }

    static void eliminateProducts() {

      int productOfSquare;
      int possiblePrime;

      for (int primeOne = 7; primeOne * primeOne <= maxNumber; primeOne += 2) {

        if (isPrime[primeOne]) {

          if (primeOne <= limit) {

            productOfSquare = primeOne * primeOne;

            while (productOfSquare < maxNumber && productOfSquare > 0) {

              isPrime[productOfSquare] = false;
              productOfSquare += (primeOne * primeOne * 2);

            }

          }

          possiblePrime = primeOne;

          do {

          possiblePrime += 2;

          } while (!isPrime[possiblePrime]);

          while (possiblePrime * primeOne <= maxNumber) {

            isPrime[possiblePrime * primeOne] = false;

            do {

            possiblePrime += 2;

            } while (!isPrime[possiblePrime]);
          }

        }
      }

    }

    static void queryPrint() {

      System.out.println("Would you like to print the primes (Y/N)?");

      Scanner inputQuery = new Scanner(System.in);
      char printResponse;

      do {

        printResponse = inputQuery.next().charAt(0);
        if (printResponse == 'Y') {

          printPrimes();

        } else {

          System.out.println("Exiting program.");
          return;

        }

      } while (printResponse != 'N' || printResponse != 'N');

    }

    static void printPrimes() {

      for (int primeIndex = 0; primeIndex < maxNumber; primeIndex++) {

        if (isPrime[primeIndex]) {

          System.out.print(primeIndex + " ");

        }

      }

      return;

    }

    static int countPrimes(int countPrimes) {

      for (int primeIndex = 0; primeIndex < maxNumber; primeIndex++) {

        if (isPrime[primeIndex]) {

          countPrimes++;

        }

    }

    return countPrimes;

  }

}