Today, a fellow student claimed that it would take much time to check the first 1,000,000 numbers for primes. I claimed that it would be a matter of seconds to do so for the first 1,000,000,000 numbers.
So, lets prove my claim.
Trivial approach
#include <iostream> // cin, cout
#include <cmath> // sqrt
#include <vector>
using namespace std;
vector<unsigned long> primeList;
bool isPrime(unsigned long n) {
vector<unsigned long>::iterator myIt;
unsigned long root = (unsigned long) sqrt(n);
for(myIt=primeList.begin(); myIt != primeList.end(); myIt++){
if(n%(*myIt) == 0) {
return false;
} else if ((*myIt) > root) {
return true;
}
}
return true;
}
int main(int argc, char* argv[]) {
primeList.push_back(2);
cout << 2 << endl;
unsigned long long max = (unsigned long long) atoi(argv[1]);
for (unsigned long i=3; i<max; i+=2) {
if(isPrime(i)) {
primeList.push_back(i);
cout << i << endl;
}
}
return 0;
}
Execute it for 100,000,000:
time ./generate-prime-list.out 100000000 > primes.txt
real 0m57.274s
user 0m41.855s
sys 0m15.229s
So 41 seconds for all primes not bigger than 100,000,000.
Now, lets test it for 1,000,000,000:
time ./generate-prime-list.out 1000000000 > primes.txt
real 18m18.205s
user 15m56.904s
sys 2m12.500s
16 minutes ... not exactly "seconds". But please keep in mind that I also wrote the result to a txt-file. This txt file is 502.0 MB big. It takes quite a lot of time to write such an amount of data from memory to disk.
Sieve of Eratosthenes
A first try
#include <iostream> // cin, cout
#include <vector>
using namespace std;
void sieveOfEratosthenes(unsigned long n) {
vector<bool> primesEratosthenes (n+1, true);
cout << 2 << endl;
for (unsigned long i=3; i<n; i+=2) {
if (primesEratosthenes[i] == true) {
cout << i << endl;
for (unsigned long j=2; j*i<=n; j++) {
primesEratosthenes[j*i] = false;
}
}
}
}
int main(int argc, char* argv[]) {
unsigned long long n = (unsigned long long) atoi(argv[1]);
sieveOfEratosthenes(n);
return 0;
}
This one takes 1 minute 5 seconds:
make;time ./generate-prime-list.out 1000000000 > testPrimes.txt
g++ -std=c++0x -Wall -pedantic -O3 -D NDEBUG generate-prime-list.cpp -o generate-prime-list.out
real 3m20.436s
user 1m4.908s
sys 2m11.748s
I'm getting closer to "seconds". ☺
ofstream, endl, \n and buffers
Writing 502.0 MB takes some time. It's not getting better when I pipe this through bash. So lets write it directly:
#include <fstream> // ofstream
#include <vector>
using namespace std;
void sieveOfEratosthenes(unsigned long n) {
ofstream myfile;
myfile.open ("huge-prime-list.txt");
vector<bool> primesEratosthenes (n+1, true);
myfile << 2 << endl;
for (unsigned long i=3; i<n; i+=2) {
if (primesEratosthenes[i] == true) {
myfile << i << endl;
for (unsigned long j=2; j*i<=n; j++) {
primesEratosthenes[j*i] = false;
}
}
}
myfile.close();
}
int main(int argc, char* argv[]) {
unsigned long long n = (unsigned long long) atoi(argv[1]);
sieveOfEratosthenes(n);
return 0;
}
time ./generate-prime-list.out 1000000000
real 3m21.249s
user 1m4.016s
sys 2m12.332s
Ok, no real change :-/
Another idea is to replace endl by \n (see explanation)
That was a good try. Now it needs only 46 seconds:
time ./generate-prime-list.out 1000000000
real 0m49.539s
user 0m46.619s
sys 0m0.920s
Another reason why this might be slow could be too many system calls. So I could buffer some and write them in blocks. I could also just write blocks of unsigned long numbers instead of a string representation. This might lead to a much smaller file size which in consequence is faster to write:
#include <stdio.h> // fopen
#include <iostream> // atoi
#include <vector>
using namespace std;
void sieveOfEratosthenes(unsigned long n) {
FILE* pFile;
pFile = fopen("huge-prime-list.txt", "wb");
vector<bool> primesEratosthenes (n+1, true);
unsigned long tmp = 2;
fwrite(&tmp, sizeof(unsigned long),1, pFile);
for (unsigned long i=3; i<n; i+=2) {
if (primesEratosthenes[i] == true) {
fwrite(&i, sizeof(unsigned long),1, pFile);
for (unsigned long j=2; j*i<=n; j++) {
primesEratosthenes[j*i] = false;
}
}
}
fclose(pFile);
}
int main(int argc, char* argv[]) {
unsigned long long n = (unsigned long long) atoi(argv[1]);
sieveOfEratosthenes(n);
return 0;
}
And execute it:
time ./generate-prime-list.out 1000000000
real 0m39.700s
user 0m38.546s
sys 0m0.640s
38 seconds for all primes from 2 to 1,000,000,000. The file size is now only 203.4 MB.
By the way, simply setting this with setbuf(pFile, NULL); to unbuffered resulted in 50 seconds execution time.
You can get the list with this snippet:
#include <iostream>
#include <fstream>
using namespace std;
int main(int argc, char* argv[]) {
if (argc != 2) {
cout << "You have to specify a file name" << endl;
} else {
FILE* pFile;
pFile = fopen(argv[1], "rb");
long long x;
size_t read;
while (!feof(pFile)) {
read = fread(&x, sizeof(long long), 1, pFile);
(void) read;
if (feof(pFile)){
break; // otherwise it duplicates the last entry
}
cout << x << endl;
}
fclose(pFile);
}
}
Improve sieving
The following change was suggested by Niklas B.. Thanks!
Take a look at the inner for loop. This one does the sieving, so it gets executed very often. In this loop, you have to calculate ji for checking the condition of the loop and again for setting it to false. You can get rid of one of those operations. Additionally, you don't have to start sieving at 2p, but you can start at p*p as you already sieved out all multiples of the first, second, ..., current-1-th prime.
#include <stdio.h> // fopen
#include <iostream> // atoi
#include <vector>
using namespace std;
void sieveOfEratosthenes(long long n) {
FILE* pFile;
pFile = fopen("huge-prime-list.bin", "wb");
vector<bool> primesEratosthenes (n+1, true);
long long tmp = 2;
fwrite(&tmp, sizeof(long long),1, pFile);
for (long long i=3; i<n; i+=2) {
if (primesEratosthenes[i]) {
fwrite(&i, sizeof(long long), 1, pFile);
for (long long j=i*i; j<=n; j+=i) {
primesEratosthenes[j] = false;
}
}
}
fclose(pFile);
}
int main(int argc, char* argv[]) {
if (argc != 2) {
cout << "You have to specify n" << endl;
} else {
long long n = (long long) atoi(argv[1]);
sieveOfEratosthenes(n);
}
return 0;
}
Execute:
time ./generate-prime-list.out 1000000000
real 0m24.222s
user 0m23.485s
sys 0m0.436s
primesieve
Primesieve is a free software program and C++ library that generates prime numbers and prime k-tuplets (twin primes, prime triplets, ...) $< 2^{64}$ using a highly optimized implementation of the sieve of Eratosthenes.
According to the GUI, it finds all 50,847,534 primes below 1,000,000,000 in 0.16 seconds. But write them to a file...
Sieve of Atkin
Arthur Oliver Lonsdale Atkin (July 31, 1925 – December 28, 2008) was a British mathematician who invented this sieve. I've implemented it according to the pseudocode provided in Wikipedia. A very long explanation of Atkins sieve is on StackOverflow
My first implementation
#include <stdio.h> // fopen
#include <iostream> // atoi
#include <vector>
#include <cmath> // sqrt, ceil,
using namespace std;
void sieveOfAtkin(long long limit) {
FILE* pFile;
pFile = fopen("huge-prime-list.bin", "wb");
long long root = ceil(sqrt(limit));
// initialize the sieve
vector<bool> is_prime(limit, false);
// put in candidate primes:
// integers which have an odd number of
// representations by certain quadratic forms
for (long long x = 1; x <= root; x++) {
long long xSquare = x*x;
for (long long y = 1; y <= root; y++) {
long long n = (4*xSquare)+(y*y);
if (n <= limit && (n % 12 == 1 || n % 12 == 5)) {
is_prime[n] = !is_prime[n];
}
n = (3*xSquare)+(y*y);
if (n <= limit && n % 12 == 7) {
is_prime[n] = !is_prime[n];
}
n = (3*xSquare)-(y*y);
if (x > y && n <= limit && n % 12 == 11) {
is_prime[n] = !is_prime[n];
}
}
}
// eliminate composites by sieving
for (long long n = 5; n <= root; n++) {
if (is_prime[n]) {
long long add = n*n;
for (long long k = add; k < limit; k += add) {
// n is prime, omit multiples of its square; this is
// sufficient because composites which managed to get
// on the list cannot be square-free
is_prime[k] = false;
}
}
}
// Output primes
long long primTmp = 2;
fwrite(&primTmp, sizeof(long long), 1, pFile);
primTmp = 3;
fwrite(&primTmp, sizeof(long long), 1, pFile);
for (long long n = 5; n < limit; n++) {
if (is_prime[n]) {
fwrite(&n, sizeof(long long), 1, pFile);
}
}
fclose(pFile);
}
int main(int argc, char* argv[]) {
if (argc != 2) {
cout << "You have to specify n" << endl;
} else {
long long n = (long long) atoi(argv[1]);
sieveOfAtkin(n);
}
return 0;
}
Atkins sieve has a much worse performance than Sieve of Eratosthenes:
time ./generate-prime-list.out 1000000000
real 1m6.001s
user 1m4.724s
sys 0m0.604s
Primegen
Primegen is an implementation by Daniel J. Bernstein. It's a little bit cluttered with 80 files and 3854 LOC in total.
When I have some time, I'll update this article and create a new version of my sieve with ideas from Primegen.
Primegen is fast:
time ./primes > primes.txt
real 0m11.677s
user 0m10.649s
sys 0m0.708s
About 11 seconds for writing all primes between 2 and 1,000,000,000 to a txt file.
See also
You might want to try the following Project Euler problem sets:
- Project Euler
- Problem 3: What is the largest prime factor of the number 600851475143 ?
- Problem 7: What is the 10 001st prime number?
- Problem 60
- SPOJ
- Segmented sieve of Eratosthenes
- Fun with prime numbers: This looks almost like my article. It might be interesting, because the author thought about finding primes above $10^9$ which I didn't.
Finally
The last script that took 23 seconds for calculating and writing all primes in 2 to 1,000,000,000 seems to be the best one. Do you know anything else that could get improved? Please provide a working example (e.g. as a gist)