Fastest way to find all primes under 4 billion

Question

I'm trying to print every prime number under 2**32. Right now I'm using a bool vector to build a sieve and then print out the primes after making the sieve. It takes 4 minutes just to print out the primes upto 1 billion. Is there a faster way to do this?? Here is my code

#include <iostream>
#include <cstdlib>
#include <vector>
#include <math.h>

using namespace std;

int main(int argc, char **argv){
  long long limit = atoll(argv[1]);
  //cin >> limit;
  long long sqrtlimit = sqrt(limit);

  vector<bool> sieve(limit+1, false);

  for(long long n = 4; n <= limit; n += 2)
    sieve[n] = true;

  for(long long n=3; n <= sqrtlimit; n = n+2){
    if(!sieve[n]){
      for(long long m = n*n; m<=limit; m=m+(2*n))
        sieve[m] = true;
    }
  }

  long long last;
  for(long long i=limit; i >= 0; i--){
    if(sieve[i] == false){
      last = i;
      break;
    }
  }
  cout << last << endl;

  for(long long i=2;i<=limit;i++)
  {
    if(!sieve[i])
      if(i != last)
        cout<<i<<",";
      else
        cout<<i;
  }
  cout<<endl;

Answer 1

I discuss the problem of generating large numbers of primes at my blog, where I find the sum of the first billion primes is 11138479445180240497. I describe four different methods:

1. Brute force, testing each number starting from 2 using trial division.

2. Generate candidates using a 2,3,5,7-wheel, then test primality with strong pseudoprime tests to bases 2, 7, and 61; this method works only up to 2^32, which was insufficient for me to sum the first billion primes but will be sufficient for you.

3. An algorithm due to Melissa O'Neill that uses sieve embedded in a priority queue, which is quite slow.

4. A segmented sieve of Eratosthenes, which is very fast but requires space to store both the sieving primes and the sieve itself.

Answer 2

This will probably speed it up a bit:

#include <algorithm>
#include <iostream>
#include <iterator>
#include <vector>

int main()
{
    std::vector<unsigned long long> numbers;
    unsigned long long maximum = 4294967296;
    for (unsigned long long i = 2; i <= maximum; ++i)
    {
        if (numbers.empty())
        {
            numbers.push_back(i);
            continue;
        }

        if (std::none_of(numbers.begin(), numbers.end(), [&](unsigned long long p)
        {
            return i % p == 0;
        }))
        {
            numbers.push_back(i);
        }

    }

    std::cout << "Primes:  " << std::endl;
    std::copy(numbers.begin(), numbers.end(), std::ostream_iterator<int>(std::cout, " "));

    return 0;
}

It is kind of the inverse of the Sieve of Eratosthenes (instead of starting with every number under the limit and eliminating multiples, it starts at 2 and ignores multiples up to the limit).

Answer 3

The fastest way would probably be to take a pregenerated list.

http://www.bigprimes.net/ has the first 1.4 billion primes available for download, which should include every prime under 30 billion or so.

I suppose loading the binary might take too long when it is a few gigabytes in size.

Answer 4

Have you benchmarked what is taking the most time? Is it the sieve itself, or the writing of the output?

One quick way to speed up the sieve is to stop worrying about all the even numbers. There's only one even number that's a prime and you can hard-code it. That cuts the size of your array in half which will help immensely if you're bumping up against the limits of physical memory.

vector<bool> sieve((limit+1)/2, false);
...
  for(long long m = n*n/2; m<=limit/2; m=m+n)
    sieve[m] = true;

As for the output itself, cout is notoriously inefficient. It might be more efficient to call itoa or some equivalent yourself, then use cout.write to output it. You might even go old school and use fwrite with stdout.

Answer 5

I wrote a fast way in C to output primes up to 4 billion in under three minutes using a single thread on my Ryzen 9 3900x. If you output it to a file, it ends up being 2.298GB and I think it uses about 20GB of memory to complete.

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

#define ARRSIZE 4000000000
#define MAXCALC ARRSIZE/2

int main() {
    long int z;
    long int *arr = (long int*) malloc((ARRSIZE) * sizeof(long int));
    for (long int x=3;x <= MAXCALC; x++) {
        if (x % 10 == 3 || x % 10 == 7 || x % 10 == 9) {
            for (long int y=3; y < MAXCALC; y++){
                    z = x * y;
                    if (z < ARRSIZE)
                        arr[z] = 1;
                    else
                        break;
            }
        }
    }
    printf("2 3 5 ");
    for (long int x=7; x < ARRSIZE; x++) {
        if (x % 2 != 0 && x % 10 != 5)
            if (arr[x] != 1)
                printf("%ld ", x);
    }
    printf("\n");

    free(arr);
    return EXIT_SUCCESS;
}

Answer 6

I wrote code i python that outputs all primes less than a billion in 8.7 seconds. But I was unsure if you ment the first 4 billion primes or allan primes less than that. Anyways, here is my solution:

import numpy as np
from math import isqrt

def all_primes_less_than(n):
    is_prime = np.full(n,True)
    is_prime[:2] = False
    is_prime[4::2] = False
    for p in range(3,isqrt(n),2):
        if is_prime[p]: is_prime[p*p::p] = False
    return is_prime.nonzero()[0]