When using trial division to test a number n for primality, the highest potential divisor that needs to be checked is floor(sqrt(n)). The first number that must not be checked is n itself, because that would result in a false verdict of non-primality (seeing that every number is divisible by itself). But that is precisely what happens for the code posted by user5486241 in the case n == 2:
for (int i = 2; i < std::sqrt(n)+1; ++i)
if (n % i == 0)
return false;
Recomputing the square root during each iteration carries a heavy runtime cost, though, unless the compiler is able to quietly correct this programming mistake. This requires the compiler to prove that the expression std::sqrt(n) is loop-invariant, which in turn requires it to prove or know that std::sqrt() is idempotent (free of observable side-effects).
The point about only scanning up to the square root is valid and important, though. Scanning up to n instead of sqrt(n) leads to an exponential increase in runtime, which is why the OP's original code is not included in the benchmarks. These benchmarks show the performance of user5486241's code (shown as 'v0') vs. the following cleaned-up version 'v1':
typedef std::uint32_t num_t;
bool v1::is_prime (num_t n)
{
if (n < 2)
return false;
for (num_t i = 2, sqrt_n = num_t(std::sqrt(n)); i <= sqrt_n; ++i)
if (n % i == 0)
return false;
return true;
}
The num_t typedef has no deeper purpose here, I introduced it only to facilitate experimentation without turning everything into a template (and making it thus unreadable).
The benchmarks show the time in milliseconds for scanning a maximum-width range according to the tasks SPOJ PRIME1 (n - m <= 100,000, n <= 1,000,000,000, timeout 6 s) and SPOJ PRINT (n - m <= 1,000,000, n <= INT32_MAX, timeout 1.223 s) respectively.
The first column shows the value for n, which is the upper end of the range and in the first row it is also equal to the constant n - m. The variants v2 through v6 are explained later.
*** BC++ 5.5.1 (free) 32 ***
PRIME1 10^5 v0 v1 v2 v3 v4 v5 v6 #primes
---------------------------------------------------------------------
100000 75 11 5 4 4 3 1 9592
1000000 253 33 16 11 11 7 1 7224
10000000 667 84 40 27 27 13 1 6134
100000000 1813 228 107 73 71 30 1 5454
1000000000 4964 624 291 196 191 69 1 4832
PRINT 10^6 v0 v1 v2 v3 v4 v5 v6 #primes
---------------------------------------------------------------------
1000000 1834 236 114 80 79 55 6 78498
10000000 6519 798 386 267 264 134 7 62090
100000000 17948 2237 1083 719 709 291 6 54332
1000000000 49234 6093 2863 1955 1890 693 6 47957
2147483647 69528 8749 4078 2741 2624 939 7 46603
With bcc64 the overhead of v0 is less pronounced than with bcc32, because the 64-bit compiler can inline the square root computation with three XMM instructions instead of having to call a library function. The timings for variants v1 through v6 are noticeably worse, though (showing only the last line):
*** BC++ 7.10 (RX/Seattle) 64 *** CLANG 3.3.1 (35465.f352ad3.17344af)
2147483647 20007 9512 4781 3373 3195 1710 22 46603
Now, gcc and VC++ as examples for highly optimising compilers that can quietly correct programming mistakes. Their timings - regardless of bitness - are exactly the same, within the limits of jitter.
*** g++ 4.8.1 32 *** on MAGNUM (3.0 GHz Haswell)
PRIME1 10^5 v0 v1 v2 v3 v4 v5 v6 #primes
---------------------------------------------------------------------
100000 10 10 5 4 3 3 0 9592
1000000 31 31 16 11 9 7 0 7224
10000000 82 82 41 27 24 14 0 6134
100000000 216 219 110 73 64 31 1 5454
1000000000 580 580 293 195 170 71 1 4832
PRINT 10^6 v0 v1 v2 v3 v4 v5 v6 #primes
---------------------------------------------------------------------
1000000 226 222 113 77 68 52 5 78498
10000000 781 780 391 263 228 133 5 62090
100000000 2104 2128 1065 706 608 294 5 54332
1000000000 5761 5824 2935 1942 1685 706 6 47957
2147483647 8221 8162 4108 2743 2379 952 6 46603
v0 ... user5486241's code
v1 ... sqrt computation hoisted out of the loop
v2 ... skipping multiples of 2
v3 ... skipping multiples of 2 and 3
v4 ... skipping multiples of 2, 3 and 5
v5 ... trial-dividing only by primes (SoE for the factors)
v6 ... simple windowed Sieve of Eratosthenes
Comparing the timings for v0 and v1 shows clearly that these two compilers managed to hoist the square root computation out of the loop all by themselves.
All these timings were taken on a 3 GHz Haswell laptop. The gcc32 timings are for a P4 build - no advanced instructions - so that I could run the same executable on a 2.8 GHz Northwood (P4), which is the kind of iron that SPOJ ran on when the limits for PRIME1 and PRINT were set:
*** g++ 4.8.1 32 *** on HELIOS (2.8 GHz Northwood)
PRIME1 10^5 v0 v1 v2 v3 v4 v5 v6 #primes
---------------------------------------------------------------------
1000000000 2378 2324 1165 778 625 266 2 4832
PRINT 10^6 v0 v1 v2 v3 v4 v5 v6 #primes
---------------------------------------------------------------------
2147483647 32806 32679 16344 10933 8800 3549 107 46603
It can be clearly seen that division instructions are about four times as fast on current machines compared to back when the tests limits were set, since v0 through v5 are essentially limited by division performance. The sieve code (v6) profits even more from the many improvements in current processors, across the board from memory system performance to superscalar out-of-order execution.
The performance of simple, division-based primality tests is interesting because they are a perfectly adequate solution for many applications where the superior bulk performance of sieves is not necessary. Euler problem #35 (circular primes) comes to mind, where only 5460 candidates need to be checked for covering the range up to 1,000,000. The advantage of the division-based code - besides its simplicity - is that it does not need auxiliary structures like factor sieves.
One simple improvement is to skip the number 2 and its multiples, which immediately doubles performance:
bool v2::is_prime (num_t n)
{
if (n <= 2)
return n == 2;
if ((n & 1) == 0)
return false;
num_t sqrt_n = num_t(std::sqrt(n));
for (num_t i = 3; i <= sqrt_n; i += 2)
if (n % i == 0)
return false;
return true;
}
The performance is actually slightly more than doubled because multiples of 2 get rejected faster now: before the call to sqrt() and before the loop is set up, and without requiring a division. This code already shows the complications necessary for wheel-based operations:
o special handling of the wheel primes before operations commence
o special handling for the multiples of the wheel primes
o special rules for striding through the search space, so that wheel prime multiples are skipped
Skipping the remaining multiples of 3 cuts the workload by another third. I'm showing it here because it makes the approach behind v4 easier to understand:
bool v3::is_prime (num_t n)
{
if (n < 4)
return n == 2 || n == 3;
if ((n & 1) == 0 || n % 3 == 0)
return false;
num_t sqrt_n = num_t(std::sqrt(n));
for (num_t i = 5, d = 2; i <= sqrt_n; i += d, d ^= 6)
if (n % i == 0)
return false;
return true;
}
The step sequence for avoiding multiples of 3 as well as multiples of 2 is 2, 4, 2, 4... (when starting at 5) which I produce here with a simple xor cheat.
The step sequence gets longer when 5 is added to the wheel; it has 8 steps which are 4, 2, 4, 2, 4, 6, 2, 6 (when starting at 7). In other words, the mod 30 wheel has 8 prime-bearing spokes and numbers above are the distances between non-dud spokes when the wheel is currently positioned at 7. I couldn't find an elegant way of creating the sequence with a few arithmetic operations, so I simply stuffed it into the nibbles of a 32-bit word that gets rotated as the algorithm chugs along. For compatibility I did the rotation with a workaround instead of a non-standard compiler intrinsic or library function. This hurts performance a bit (increasing the instruction count from one to three, with attendant increase in register pressure), which is why the performance increase is slightly less than the expected 20%.
I'm also showing a different way of dealing with the wheel primes: instead of testing for equality with each of them, I use a constant as a bitmap. This allows me to deal with 32 numbers in one go. The performance of this branch is insignificant because it is taken very rarely, but the tighter code means less pressure on the instruction cache.
A similar approach could be used for discarding all multiples of 2, 3 and 5 after a single division by 30 (if (EMPTY_SPOKES & (1 << (n % 30))) return false) but this does not result in an improvement because it rejects multiples of 2 slightly more slowly than the code shown in v4, which is why performance-wise it is pretty much a wash. Hence the simpler, more pedestrian code.
template<unsigned SHIFT> // replace with _lrotl() or similar
inline
num_t rotate_right (num_t n)
{
return num_t( (n >> SHIFT) | (n << (sizeof(num_t) * CHAR_BIT - SHIFT)) );
}
enum
{
PRIMES_7_TO_29 = (1<<7) | (1<<11) | (1<<13) | (1<<17) | (1<<19) | (1<<23) | (1<<29),
PRIMES_2_TO_31 = (1<<2) | (1<< 3) | (1<< 5) | PRIMES_7_TO_29 | (1<<31)
};
bool v4::is_prime (num_t n)
{
if (n < 32)
return (PRIMES_2_TO_31 >> n) & 1; // normalised to 0/1; good compilers recognise this
if ((n & 1) == 0 || n % 3 == 0 || n % 5 == 0)
return false;
num_t sqrt_n = num_t(std::sqrt(n));
for (num_t i = 7, w = 0x62642424; i <= sqrt_n; i += (w & 15), w = rotate_right<4>(w))
if (n % i == 0)
return false;
return true;
}
v5 needs a vector of potential prime factors up to the square root of the upper limit. For a limit of INT32_MAX the maximum factor not exceeding the square root is 46340, and even the simplest, most-straightforward code can sieve this in less than a millisecond:
// intended to go no higher than sqrt(2^31) == 46340, hence no need for segmentation or wheeling
// (despite exceeding the L1 cache not slower than vector<bool> for gcc and lots faster for bcc)
void sieve_small_primes_up_to (unsigned limit, std::vector<num_t> &result)
{
result.resize(0);
unsigned max_factor = unsigned(std::sqrt(limit)); // = 215
std::vector<char> is_composite(limit + 1); // +1 because we're indexing with numbers
for (unsigned i = 2; i <= max_factor; ++i) // * 214
if (!is_composite[i]) // * 47
for (unsigned j = i * i; j <= limit; j += i)
is_composite[j] = true;
for (unsigned i = 2; i <= limit; ++i) // * 46339
if (!is_composite[i]) // * 4792
result.push_back(num_t(i));
}
Obviously, this initialisation must be arranged to happen before the primality test function is called. In order to make v5 directly comparable to v4 I decided to handle the primes 2 to 5 in exactly the same way as v4 does. This gives v5 the same advantage as v4 has, which is fast rejection of the multiples of 2, 3 and 5 (73.3% of all numbers). Hence the wheel primes need to be removed from the vector:
sieve_small_primes_up_to(num_t(std::sqrt(INT32_MAX)), small_factors_gt_5);
// erase the wheel primes
small_factors_gt_5.erase(small_factors_gt_5.begin(), small_factors_gt_5.begin() + 3);
Normally this would all be encapsulated in a class that houses the function and its auxiliary data, guarding its invariants. Here I'm showing only the raw primality test code:
bool v5::is_prime (num_t n)
{
if (n < 32)
return (PRIMES_2_TO_31 >> n) & 1; // normalised to 0/1; good compilers recognise this
if ((n & 1) == 0 || n % 3 == 0 || n % 5 == 0)
return false;
num_t sqrt_n = num_t(std::sqrt(n));
for (auto prime: small_factors_gt_5)
{
if (prime > sqrt_n)
return true;
if (n % prime == 0)
return false;
}
assert("too few small_factors_gt_5[] (not initialised properly)");
return false;
}
Older compilers and the classic (non-clang) versions of the Emborlandero compiler need some #ifdefery for replacing the range-based for loop with a normal iterator-based for loop, but I've edited that out here in order to improve readability.
And finally, some code that can beat even PRINT without breaking a sweat, with several orders of magnitude reserve regarding time and memory limits. Most of the bits and pieces are already in place; all that remains is replacing iterated calls to is_prime() with a small Eratosthenean sieve. The only tricky bit here is computing exactly where in the sieve window the sequence of hops for a given prime starts. Here I formulated this in a way that anchors the computation at the slot before the start of the sieve window (hence the apparently unmotivated -1); this has the advantage that it avoids additional branching or a double mod.
Since the code is orders of magnitude faster than the SPOJ requirements, it can simply sieve its factor primes on every call.
void v6::get_primes (num_t m, num_t n, std::vector<num_t> &result)
{
std::vector<num_t> small_factors;
sieve_small_primes_up_to(unsigned(std::sqrt(n)), small_factors);
num_t window_size = n - m + 1;
std::vector<char> is_composite(window_size);
for (num_t prime: small_factors)
{
num_t first_multiple = prime * prime; // absolute, may lie outside the window
if (first_multiple > n)
break;
num_t stride = prime; // no need to get cute with wheely things and whatnot
num_t i;
if (first_multiple >= m)
i = first_multiple - m;
else
i = stride - (m - 1 - first_multiple) % stride - 1;
for ( ; i < window_size; i += stride)
is_composite[i] = true;
}
result.resize(0);
// 0 and 1 aren't composites, but they aren't primes either -> skip them if they're there
for (num_t i = num_t((m < 1) + (m < 2)); i < window_size; ++i)
if (!is_composite[i])
result.push_back(m + i);
}
Nis prime if no smaller prime divides into it exactly (ie you don't need to try all the numbers less thanNjust the primes less thanN). If you keep track of all the primes you have found that makes it easier. You can seed your prime list (you don't need to calculate all of them). primes.utm.edu/lists/small/10000.txt If you want to do it fast just get a list of pre-generated primes and print them: primes.utm.edu/lists/small/millions \$\endgroup\$ – Martin York Feb 5 '16 at 1:36