sieving primes up to 10^12

Question

I'm working on something that requires me to generate all prime numbers up to 10^12.

Since I have never needed this many primes before, I normally just implement the algorithm on this webpage here

The problem here is, of course, that 10^12 is larger than the max value for an integer, and as a result I cannot make an array of that size.

I am unfamiliar with the methods one would use to generate so many primes efficiently, and was wondering if anyone could shed some light on the situation.

Answer 1

You will need to use a segmented sieve.

The basic idea of a segmented sieve is to choose the sieving primes less than the square root of n, choose a reasonably large segment size that nevertheless fits in memory, and then sieve each of the segments in turn, starting with the smallest. At the first segment, the smallest multiple of each sieving prime that is within the segment is calculated, then multiples of the sieving prime are marked as composite in the normal way; when all the sieving primes have been used, the remaining unmarked numbers in the segment are prime. Then, for the next segment, for each sieving prime you already know the first multiple in the current segment (it was the multiple that ended the sieving for that prime in the prior segment), so you sieve on each sieving prime, and so on until you are finished.

Consider the example of sieving from 100 to 200 in segments of 20; the 5 sieving primes are 3, 5, 7, 11 and 13. In the first segment from 100 to 120, the bitarray has 10 slots, with slot 0 corresponding to 101, slot k corresponding to 100 + 2k + 1, and slot 9 corresponding to 119. The smallest multiple of 3 in the segment is 105, corresponding to slot 2; slots 2+3=5 and 5+3=8 are also multiples of 3. The smallest multiple of 5 is 105 at slot 2, and slot 2+5=7 is also a multiple of 5. The smallest multiple of 7 is 105 at slot 2, and slot 2+7=9 is also a multiple of 7. And so on.

Function primes takes arguments lo, hi and delta; lo and hi must be even, with lo < hi, and lo must be greater than the square root of hi. The segment size is twice delta. Array ps of length m contains the sieving primes less than the square root of hi, with 2 removed since even numbers are ignored, calculated by the normal Sieve of Eratosthenes. Array qs contains the offset into the sieve bitarray of the smallest multiple in the current segment of the corresponding sieving prime. After each segment, lo advances by twice delta, so the number corresponding to an index i of the sieve bitarray is lo + 2 i + 1.

function primes(lo, hi, delta)
    sieve := makeArray(0..delta-1)
    ps := tail(primes(sqrt(hi)))
    m := length(ps)
    qs := makeArray(0..m-1)
    for i from 0 to m-1
        qs[i] := (-1/2 * (lo + ps[i] + 1)) % ps[i]
    while lo < hi
        for i from 0 to delta-1
            sieve[i] := True
        for i from 0 to m-1
            for j from qs[i] to delta step ps[i]
                sieve[j] := False
            qs[i] := (qs[i] - delta) % ps[i]
        for i from 0 to delta-1
            t := lo + 2*i + 1
            if sieve[i] and t < hi
                output t
        lo := lo + 2*delta

For the sample given above, this is called as primes(100, 200, 10). In the example given above, qs is initially [2,2,2,10,8], corresponding to smallest multiples 105, 105, 105, 121 and 117, and is reset for the second segment to [1,2,6,0,11], corresponding to smallest multiples 123, 125, 133, 121 and 143.

The value of delta is critical; you should make delta as large as possible so long at it fits in cache memory, for speed. Use your language's library for the bitarray, so that you only take a single bit for each sieve location. If you need a simple Sieve of Eratosthenes to calculate the sieving primes, this is my favorite:

function primes(n)
    sieve := makeArray(2..n, True)
    for p from 2 to n step 1
        if sieve(p)
            output p
            for i from p * p to n step p
                sieve[i] := False

Those functions are both in pseudocode; you'll have to translate to Java, with appropriate integer data types. Where the pseudocode says output you can print the prime, or collect the primes in an array, whatever you want to do with them.

I've done a lot of work with primes at my blog, including the essay Programming with Prime Numbers which includes a segmented sieve on the last page.