Segmented Sieve of Eratosthenes in Java

As promised, here's a simple implementation of the segmented sieve. The segmented sieve is a very straight forward algorithm. The objective is to find all prime numbers between L and R. Typically L and R are very large (1e16-1e18), while R-L is much smaller (1e8-1e10). The odd prime numbers that have multiples in this interval are all smaller than sqrt(R) can be found by using a basic sieve. The sieve works in exactly the same way as before, except for each prime p we start sifting from its first multiple contained in the given interval. To find it, we need to find the first number n that satisfies the equation n mod p = 0. Now be q = -L mod p, then n = L + q. In fact, (L + p) mod p = (L - L) mod p = 0. L + q is the first number to be discarded; the rest of the multiples are found by repeatedly adding p to this number.
For example, if L = 100 and R = 200, all the odd primes smaller than sqrt(200) are 3, 5, 7, 11 and 13. -100 mod 3 = 2 and 102 is the first multiple of 3 contained in the interval. Likewise, -100 mod 5 = 0, -100 mod 7 = 5 and so on.
This implementation of the segmented sieve makes use of the basic sieve from my previous post. Note that for very big numbers (e.g., R=2^63, sqrt(R)>2^31) the basic sieve might not work as expected (int should be replaced by long).

 public class SegmentedPrimeSieve {  
      private final byte sieve[];  
      private final long start;  
      /**  
       * Creates a segmented sieve in the interval defined by the values of start  
       * and end.  
       *   
       * @param start  
       * @param end  
       */  
      public SegmentedPrimeSieve(long start, long end) {  
           // if the starting value is not odd, choose the next one  
           start = start % 2 == 0 ? ++start : start;  
           // length of the byte array  
           int length = (int) ((end - start) / 16 + 1);  
           sieve = new byte[length];  
           // finally, let's compute the extended range  
           end = start + length * 16 - 2;  
           // find all the primes up to sqrt(end)  
           int maxPrime = (int) Math.floor(Math.sqrt(end));  
           PrimeSieve baseSieve = new PrimeSieve((int) maxPrime);  
           System.out.println("Sieving numbers between " + start + " and " + end);  
           // maximum value of k to sift multiples of primes in the form 2*k+1  
           int maxK = maxPrime / 2;  
           long intervalHalfSize = 8 * length;  
           // let's assume primes in the form 2*k+1 starting from k=1  
           for (int k = 1; k <= maxK; k++) {  
                // if the number is marked as a prime in the basic sieve start  
                // sifting all of its multiples in the given interval  
                if (baseSieve.get(k)) {  
                     final int p = 2 * k + 1;  
                     // This is the initial offset to start sifting from (-start%p)  
                     int offset = (int) ((p - (start % p)) % p);  
                     // if the offset is odd, start+offset is even, skip it because  
                     // we don't have even numbers in the sieve. divide by two for  
                     // the same reason. Note that this step is crucial!  
                     if (offset % 2 == 1)  
                          offset += p;  
                     offset /= 2;  
                     for (; offset < intervalHalfSize; offset += p) {  
                          sieve[offset >> 3] |= (1 << (offset & 7));  
                     }  
                }  
           }  
           this.start = start;  
      }  
      public boolean isPrime(long n) {  
           if (n < start)
                throw new RuntimeException("The number " + n
                          + " is too small for the values in the sieve.");
           if (n == 2)  
                return true;  
           if (n == 1 || n % 2 == 0)  
                return false;  
           int dn = (int) (n - start);  
           int i = dn / 16;  
           if (i >= sieve.length)  
                throw new RuntimeException("The number " + n  
                          + " exceeds the values in the sieve.");  
           return ((sieve[i] >> ((dn / 2) & 7)) & 1) == 0;  
      }  
 }  

Eratosthenes Prime Sieve: a Java Implementation

One of the classical problems in the Sphere Online Judge consists in generating prime numbers. Such a problem can be solved using a prime sieve, a technique that allows finding new prime numbers by discarding all the multiples of known prime numbers. For example, we know that 2 is a prime, so we would automatically discard all even numbers, then we would go ahead and discard the multiples of 3, 5, 7 and so on. Eratosthenes came up with this algorithm roughly 2250 years ago (Wikipedia provides a good description of the sieve).
I wrote down a very straightforward implementation of Eratosthenes sieve in java. Since the algorithm can be used to find a considerable amount of prime numbers, usually it's implemented so that the amount of memory used is as small as possible. A single bit of information is used to store whether or not a number is prime. Thus each byte contains information about 8 numbers. Discarding all even numbers beforehand cuts in half the number of bits needed. Other techniques can be used to reduce the amount of memory used, but I decided to stop here for simplicity. Here's my source code:

 public class PrimeSieve {  
      private final byte sieve[];  
      /**  
       * Creates a sieve of integers up to n.  
       *   
       * @param n  
       */  
      public PrimeSieve(int n) {  
           // using one bit per number, skipping even numbers  
           int sieveSize = n / 16 + 1;  
           // round up to the next multiple of 16  
           n = sieveSize * 16;  
           System.out.println("Sieving numbers up to " + n);  
           // initialize the array of bytes. Each bit corresponds to an odd integer  
           // between 1 and n, starting with the rightmost bit of the first byte.  
           // If the bit is 0, the number is prime. Initially, all numbers are  
           // assumed to be prime and some will be sifted out.  
           sieve = new byte[sieveSize];  
           // 1 is composite  
           sieve[0] = 0x01;  
           // this is the maximum starting number to search for primes in the form  
           // 2*k+1  
           int maxK = (int) Math.floor(Math.sqrt(n / 2));  
           int nHalf = n / 2;  
           // loop on numbers of the form 2*k+1  
           for (int k = 1; k <= maxK; k++) {  
                // if 2*k+1 is marked as prime, sift all its multiples  
                if (get(k)) {  
                     // start from (2*k+1)^2: must divide this by two since the array  
                     // doesn't contain multiples of two. Thus the starting number is  
                     // (2*k+1)^2/2 = 2*k*(k+1). Note that this is odd.  
                     // the increment is 2*k+1 (since the sieve contains only odd  
                     // numbers, using this increment automatically skips to the next  
                     // odd multiple).  
                     final int increment = 2 * k + 1;  
                     for (int composite = 2 * k * (k + 1); composite < nHalf; composite += increment)  
                          // the index in the array is obtained by discarding the  
                          // rightmost 3 bits (or divide by 8); likewise, the position  
                          // in the byte is obtained by right shifting one as many  
                          // time as the number represented by the same 3 bits.  
                          // Note that the function get() is implemented similarly.  
                          sieve[composite >> 3] |= (1 << (composite & 7));  
                }  
           }  
      }  
      /**  
       * Checks if the number 2*n+1 is marked as prime.  
       *   
       * @param n  
       *      An integer number.  
       * @return True if 2*n+1 is marked as prime, false otherwise.  
       */  
      boolean get(int n) {  
           return ((sieve[n >> 3] >> (n & 7)) & 1) == 0;  
      }  
      /**  
       *   
       * @param n  
       *      An integer.  
       * @return True if n is prime.  
       */  
      public boolean isPrime(int n) {  
           if (n == 2)  
                return true;  
           if (n == 1 || n % 2 == 0)  
                return false;  
           int i = n / 16;  
           if (i >= sieve.length)  
                throw new RuntimeException("The number " + n  
                          + " exceeds the values in the sieve.");  
           return ((sieve[i] >> ((n / 2) & 7)) & 1) == 0;  
      }  
 }  

Sometimes it's useful to find prime numbers in a small interval instead of all of them up to a certain value. In cases like this, it's inconvenient to use the full sieve, since its speed and memory requirements always depend on the upper primes bound. A segmented sieve allows sifting the composite numbers in a given interval, ignoring most of the previous numbers. In this case speed and memory requirements will depend (mostly) on the length of the interval. I'll give an implementation of the segmented sieve in my next blog post.