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;
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.