Prime sieve implementation (using trial division) in C++

How can I optimize this?

#include <iostream>
using namespace std;
bool isprime(long nur);
void printAll();

long primes[10000000];
long primeCount = 0;

int main() {
long limit;
cout << "Enter the limit ";
cin >> limit;

for (long nur=2; nur<= limit; nur++)
{
if (isprime(nur)==true)
{
cout<<nur<<endl;
primes[primeCount] = nur;
primeCount++;

if(primeCount == 10000000)
{
cout << "Increase Buffer to store primes" << endl;
break;
}
}
}
cout << "Final " << primes[primeCount-1] << endl;
//    printAll();
return 0;
}

bool isprime(long nur){
for(long itr =0; itr < primeCount; ++itr)
{
if(nur % primes[itr]==0)
{
return false;
}
}
return true;
}

void printAll()
{
cout << "All Primes" << endl;
for(long itr =0; itr < primeCount; ++itr)
{
cout << "\t" << primes[itr] << endl;
}
}

• This is not (even similar to) the sieve of Eratosthenes. The first (and probably biggest) optimization would be to implement an actual sieve. For a sieve, you typically start with odd numbers from 3. You output 3 as being prime. Then you "cross out" all multiples of three. Then you find the next number that isn't crossed out, output it as a prime, and cross out its multiples. Continue until bored. – Jerry Coffin May 22 '16 at 23:01
• So the point of the sieve of Eratosthenes is to cancel out all multiples of primes and whatever is left over is prime. Is that right? – Neeraj Mula May 22 '16 at 23:05
• Yes, that's about it. – Jerry Coffin May 22 '16 at 23:09
• @NeerajMula: Yes. But note because you are crossing out non primes. Once you have done a prime you scan up to the next prime valid prime (not crossed out number). Then that is prime and you can cross out all its multiples. A good example: codereview.stackexchange.com/questions/80353/… good explanation algolist.net/Algorithms/Number_theoretic/Sieve_of_Eratosthenes – Martin York May 23 '16 at 2:36
• The simplest optimization is to ignore multiples of 2 and 3 automatically. Start at 5. Add 2 (check). Add 4 (check). Add 2 (check). Add 4 (check). etc... – Martin York May 23 '16 at 2:37

One remark up front, before I go more into detail: if you put the processing logic into a separate function instead of inlining everything in main() then it becomes much easier to work with - like doing experiments, benchmarking (perhaps comparing different versions of the code) and to run automated tests.

As it is you cannot even test your code without manually entering a limit and then scrutinising the printed primes, unless you write a separate program that pipes input into your program and checks its output.

Also, I must concur with Jerry Coffin: your code is to the Sieve of Eratosthenes as a slingshot is to a MAC-10. There is some similarity in purpose but fundamentally different operating principles and vastly different throughput.

The Achilles Heel of trial division is that it needs to perform many expensive, slow divisions.

On my notebook, the CPU manages to perform about 300,000 divisions per millisecond but in the same time it can sieve more than 2 million numbers in the manner of Eratosthenes. Also, in order to prove the primality of a candidate, trial division has to divide by all primes up to the square root of the candidate. That's a lot of divisions: up to 6,542 for a 32-bit number, and up to 203,280,221 for a 64-bit number.

The upshot is that the Sieve of Eratosthenes beats trial division by three to four orders of magnitude as far as bulk sieving is concerned. There are still some niches where trial division is a good choice, though:

• infrequent primality tests of smallish integers (not much bigger than 32 bits)
• primality tests for candidates scattered over unsievably wide ranges
• sifting tiny windows, as in next_prime(n) (gaps between primes up to 32 bit are no wider than 335)
• eliminating multiples of small primes before doing heavy-duty primality tests like Miller-Rabin

There are two basic approaches for improving the performance of trial division:

• make division faster by replacing it with multiplication
• eliminate multiples of small primes so that fewer divisions need be performed

Replace division with multiplication for a seven-fold speedup

To make a long story short, a divisibility test like

if (n % 23 == 0)
...


can be replaced with

if (n * 3921491879 <= 186737708)
...


On my laptop this clocks 2.2 million divisibility tests per millisecond (that's with the unpredictable if), more than seven times as fast as division.

The mighty gcc does this already for constant divisors/moduli, whereas most other compilers can use only a more general scheme where a modulo op involves two multiplications and several lesser arithmetic ops. The full story can be read in these two papers whose authors include Torbjörn Granlund and Peter L. Montgomery (yes, the Peter L. Montgomery):

Compilers use these tricks only for constant divisors/moduli but nothing keeps us from computing the data at runtime, for all primes up to the square root of some limit like 2^32 (6541 odd primes).

The multiplier is the multiplicative inverse (mod 2^32) of the prime, and the threshold being compared to is 2^32 divided by the prime. Extension to 64 bits is trivial. The modular inverse could be computed via the extended GCD algorithm, but there's a simpler and much faster way:

uint32_t ModularInverse (uint32_t n)
{
uint32_t x = 2 - n;

x *= 2 - x * n;
x *= 2 - x * n;
x *= 2 - x * n;
x *= 2 - x * n;

return x;
}


That's the Newton-Raphson method, a.k.a. Hensel lifting. The details are in

Eliminate multiples of small primes from consideration

This has already been written about extensively here on Code Review and on Stack Overflow. For example, in this answer elsewhere I discussed and benchmarked the skipping of small primes in connection with trial division. This approach is especially effective for applications like next_prime(n), since the primes skipped when enumerating candidates can also be excluded from the trial divisions itself. E.g. if enumeration of candidates skips multiples of 2 and 3 then the trial division can start with the prime 5.

Extension beyond a tiny handful of primes is difficult because wheels get prohibitively big very quickly:

[]

A wheel that skips all odd primes less than 255 - i.e. 3 through 251 - would have a circumference (modulus) of about 3.2 * 10^100, and it would have 6.4 * 10^99 prime-bearing spokes (or deltas analogous to the 4,2,4,2... of the mod 6 wheel). Here are the actual numbers in order to really drive the lesson home:

modulus = 32133165458954322436165317614053356655440093295804604057122379340449075183940351576262600371617210115
spokes = 6449338146298816039189204226752375763972945771408415527222676870616270485927336345600000000000000000


Such a wheel cannot be pre-computed, of course. However, if only a tiny window of numbers needs to be scanned (like the gap between two primes, which cannot be wider than 335 for 32-bit numbers) then one can easily compute the tiny section of the wheel that covers the range of interest.

An approach inspired by the deque sieve (a.k.a. 'incremental sieve' or 'sliding sieve') would need only a tiny array of 335+1 bytes that can easily be allocated on the stack, and an approach inspired by the windowed Sieve of Eratosthenes would be similarly effective.

The crucial point in both cases is to limit the primes used to some set of tiny primes (e.g. the byte-sized ones between 3 and 251) instead of going all the way to the square root. The latter would be prohibitive for numbers much bigger than 32 bits, and for ranges where it's not prohibitive one might as well use a windowed sieve directly (C#):

UInt32 first_prime_not_less_than (UInt32 n)
{
return primes_between(n, n + Math.Min(335, UInt32.MaxValue - n)).First();
}


The divisions needed for computing the starting offsets within the window for each prime are not lost, only borrowed from the trial division which can henceforth skip them. Needless to say, these divisions can also be replaced with multiplications.