Segmented Sieve of Eratosthenes?

It’s easy enough to make a simple sieve:

for (int i=2; i<=N; i++){
    if (sieve[i]==0){
        cout << i << " is prime" << endl;
        for (int j = i; j<=N; j+=i){
            sieve[j]=1;
        }
    }
    cout << i << " has " << sieve[i] << " distinct prime factors/n";
}

But what about when N is very large and I can’t hold that kind of array in memory? I’ve looked up segmented sieve approaches and they seem to involve finding primes up until sqrt(N) but I don’t understand how it works. What if N is very large (say 10^18)?

Segmented sieve of eratosthenes: sieving composite numbers by their indexes?

I’m trying to code a prime number finder that prints primes between two given values. I’ve got no problem with coding traditional sieve, but when it comes segmented, my python knowledge is coming up s

Segmented Sieve of Atkin, possible?

I am aware of the fact that the Sieve of Eratosthenes can be implemented so that it finds primes continuosly without an upper bound (the segmented sieve). My question is, could the Sieve of Atkin/Bern

about Sieve of Eratosthenes

I’ve got some trouble with the sieve of eratosthenes. So i’ve got the mathmatics calculations for the sieve from a book called Schaum’s outlines but i think the book has programmed to code wrong…

Sieve of Eratosthenes implementation

I am trying to implement algorithm for Sieve of Eratosthenes but I don’t know why this program crashes for larger programs. Initially I was using vector but now I am implementing this using dynamic me

Sieve of Eratosthenes

I read up on the sieve of Eratosthenes while solving a question on Project Euler. I’m sure you guys know which question im talking about… So here’s the thing. My code manages to show all the primes

Sieve of Eratosthenes Scheme

I’ve been searching the web for an implementation of the Sieve of Eratosthenes in scheme and although I came up with a lot of content, none of them seemed to have made it like I need it to be done. Th

Sieve of Eratosthenes malfunction

I am trying to create a Sieve of Eratosthenes in Java but the code I wrote seems to be flawed. I was wondering if anyone else can spot my mistake since I cannot. The output I get is simply [2], which

Sieve of Eratosthenes Halting

I’m trying to implement the Sieve of Eratosthenes via the Wikipedia page, for some reason this code halts and doesn’t finish. I’m a beginner at C, so please explain If I misused anything. I’m not sure

Help understanding eratosthenes sieve implementation

I found this LINQ implementation of the eratosthenes sieve on this website. I understand the basic concept of the sieve, but there’s one detail I don’t get. What is the purpose of the first Enumerable

Spoj PRIME1 using sieve of eratosthenes (in C)

I’m trying to solve the problem PRIME1 using segmented sieve of Eratosthenes. My program works correctly with the normal sieve that is up to NEW_MAX. But there is some problem with cases n > NEW_MA

Answers

There’s a version of the Sieve based on priority queues that yields as many primes as you request, rather than all of them up to an upper bound. It’s discussed in the classic paper “The Genuine Sieve of Eratosthenes” and googling for “sieve of eratosthenes priority queue” turns up quite a few implementations in various programming languages.

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.

The size of n doesn’t matter, except that a larger n will take longer to sieve than a smaller n; the size that matters is the size of the segment, which should be as large as convenient (say, the size of the primary memory cache on the machine).

You can see a simple implementation of a segmented sieve here. Note that a segmented sieve will be very much faster than O’Neill’s priority-queue sieve mentioned in another answer; if you’re interested, there’s an implementation here.

EDIT: I wrote this for a different purpose, but I’ll show it here because it might be useful:

Though the Sieve of Eratosthenes is very fast, it requires O(n) space. That can be reduced to O(sqrt(n)) for the sieving primes plus O(1) for the bitarray by performing the sieving in successive segments. 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 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, the smallest multiple of each sieving prime is the multiple that ended the sieving in the prior segment, and so the sieving continues until finished.

Consider the example of sieve from 100 to 200 in segments of 20. The five sieving primes are 3, 5, 7, 11 and 13. In the first segment from 100 to 120, the bitarray has ten 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 primesRange takes arguments lo, hi and delta; lo and hi must be even, with lo < hi, and lo must be greater than sqrt(hi). The segment size is twice delta. Ps is a linked list containing the sieving primes less than sqrt(hi), with 2 removed since even numbers are ignored. Qs is a linked list containing the offest 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 + 2i + 1.

function primesRange(lo, hi, delta)
    function qInit(p)
        return (-1/2 * (lo + p + 1)) % p
    function qReset(p, q)
        return (q - delta) % p
    sieve := makeArray(0..delta-1)
    ps := tail(primes(sqrt(hi)))
    qs := map(qInit, ps)
    while lo < hi
        for i from 0 to delta-1
            sieve[i] := True
        for p,q in ps,qs
            for i from q to delta step p
                sieve[i] := False
        qs := map(qReset, ps, qs)
        for i,t from 0,lo+1 to delta-1,hi step 1,2
            if sieve[i]
                output t
        lo := lo + 2 * delta

When called as primes(100, 200, 10), the sieving primes ps are [3, 5, 7, 11, 13]; 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.

You can see this program in action at http://ideone.com/iHYr1f. And in addition to the links shown above, if you are interested in programming with prime numbers I modestly recommend this essay at my blog.

You ask how you could find the number of distinct prime factors for each number n up to 1018.

Segmented sieve is not applicable for this.

The sieve as shown in your question is marking each entry with the number of distinct prime factors. Thus you have to start the marking from prime itself and not its square, and thus will have to maintain the core primes reaching into N / 2, not the sqrt(N) as is usual with the Boolean sieves marking just the primality of each entry.

So instead of maintaining 10^18 8-bit ints (the multiplicities) at the same time, as is with your current scheme (i.e. 1,000,000 TB of memory), you’ll have to maintain up to 5*10^17 bits (the primes) at the same time, while sieving the segments one after another. Even wheel-packed by 30 numbers per byte, that’s still about 16,666 TB of memory.

That’s not a solution at all.

And choosing to maintain the pi(5*10^17) = 1.23*10^16 64-bit core primes themselves would take up about 100,000 TB of memory.

Your only option seems to be the various O(1)-space factorization algorithms, like e.g. deterministic Miller-Rabin with a specific set of witnesses, to actually find the prime factorization for each n.

It’s just that we are making segmented with the sieve we have. The basic idea is let’s say we have to find out prime numbers between 85 and 100. We have to apply the traditional sieve,but in the fashion as described below:

So we take the first prime number 2 , divide the starting number by 2(85/2) and taking round off to smaller number we get p=42,now multiply again by 2 we get p=84, from here onwards start adding 2 till the last number.So what we have done is that we have removed all the factors of 2(86,88,90,92,94,96,98,100) in the range.

We take the next prime number 3 , divide the starting number by 3(85/3) and taking round off to smaller number we get p=28,now multiply again by 3 we get p=84, from here onwards start adding 3 till the last number.So what we have done is that we have removed all the factors of 3(87,90,93,96,99) in the range.

Take the next prime number=5 and so on……………… Keep on doing the above steps.You can get the prime numbers (2,3,5,7,…) by using the traditional sieve upto sqrt(n).And then use it for segmented sieve.

Based on Swapnil Kumar answer I did the following algorithm in C. It was built with mingw32-make.exe.

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

int main()
{
    const int MAX_PRIME_NUMBERS = 5000000;//The number of prime numbers we are looking for
    long long *prime_numbers = malloc(sizeof(long long) * MAX_PRIME_NUMBERS);
    prime_numbers[0] = 2;
    prime_numbers[1] = 3;
    prime_numbers[2] = 5;
    prime_numbers[3] = 7;
    prime_numbers[4] = 11;
    prime_numbers[5] = 13;
    prime_numbers[6] = 17;
    prime_numbers[7] = 19;
    prime_numbers[8] = 23;
    prime_numbers[9] = 29;
    const int BUFFER_POSSIBLE_PRIMES = 29 * 29;//Because the greatest prime number we have is 29 in the 10th position so I started with a block of 841 numbers
    int qt_calculated_primes = 10;//10 because we initialized the array with the ten first primes
    int possible_primes[BUFFER_POSSIBLE_PRIMES];//Will store the booleans to check valid primes
    long long iteration = 0;//Used as multiplier to the range of the buffer possible_primes
    int i;//Simple counter for loops
    while(qt_calculated_primes < MAX_PRIME_NUMBERS)
    {
        for (i = 0; i < BUFFER_POSSIBLE_PRIMES; i++)
            possible_primes[i] = 1;//set the number as prime

        int biggest_possible_prime = sqrt((iteration + 1) * BUFFER_POSSIBLE_PRIMES);

        int k = 0;

        long long prime = prime_numbers[k];//First prime to be used in the check

        while (prime <= biggest_possible_prime)//We don't need to check primes bigger than the square root
        {
            for (i = 0; i < BUFFER_POSSIBLE_PRIMES; i++)
                if ((iteration * BUFFER_POSSIBLE_PRIMES + i) % prime == 0)
                    possible_primes[i] = 0;

            if (++k == qt_calculated_primes)
                break;

            prime = prime_numbers[k];
        }
        for (i = 0; i < BUFFER_POSSIBLE_PRIMES; i++)
            if (possible_primes[i])
            {
                if ((qt_calculated_primes < MAX_PRIME_NUMBERS) && ((iteration * BUFFER_POSSIBLE_PRIMES + i) != 1))
                {
                    prime_numbers[qt_calculated_primes] = iteration * BUFFER_POSSIBLE_PRIMES + i;
                    printf("%d/n", prime_numbers[qt_calculated_primes]);
                    qt_calculated_primes++;
                } else if (!(qt_calculated_primes < MAX_PRIME_NUMBERS))
                    break;
            }

        iteration++;
    }

    return 0;
}

It set a maximum of prime numbers to be found, then an array is initialized with known prime numbers like 2, 3, 5…29. So we make a buffer that will store the segments of possible primes, this buffer can’t be greater than the power of the greatest initial prime that in this case is 29.

I’m sure there are a plenty of optimizations that can be done to improve the performance like parallelize the segments analysis process and skip numbers that are multiple of 2, 3 and 5 but it serves as an example of low memory consumption.