Printing all the prime numbers between two bounds

Question

I keep getting a time limit exceeded for this solution to SPOJ Prime1. I need to print all primes within a given interval (with intervals as large as 10⁵, and numbers as large as 10⁹).

I've tried everything I could think of to get this code to be faster but nothing works. What can I do?

#include <iostream>

bool isprime(int n) {
    if (n < 2) return false;
    for (int i = 2; i < n; ++i) {
        if (n % i == 0) return false;
    }
    return true;
}

int main() {
    std::ios_base::sync_with_stdio(false);
    std::cin.tie(0);
    // Code to to print all the prime numbers between two bounds. n is the
    // number of test cases.
    int n, j, k;
    std::cin >> n;

    for (int i = 0; i < n; ++i) {
        std::cin >> j >> k;
        while (j-- != k) {
            if (isprime(j)) std::cout << j << '\n';
        }
    }
}

Answer 1

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);
}

Answer 2

When checking if a number is prime, you only need to check up to square root of the number. This is because their are pairs of factors, and one is below the square root, and the other is above the square root. To improve your code use

for (int i = 2; i < std::sqrt(n)+1;++i) {
    if (n % i == 0) return false;
}