Primes generation

The classic approach at generating primes less than a certain value \( N \) is the Sieve of Eratosthenes. The algorithm simply proceeds by dropping all the numbers generated by repeatedly adding a number (thus multiples of the first number).

#include <iostream>
#include <cmath>
#include <vector>
using namespace std;

void printPrimes(unsigned int N) {
  // Generate all prime numbers less than or equal to N

  vector<bool> isPrime(N + 1, true);
  isPrime[0] = false;
  isPrime[1] = false; // By definition

  for (unsigned int i = 4; i<isPrime.size(); i += 2) {
    isPrime[i] = false;
  }
  unsigned int sqr = static_cast<unsigned int>(sqrt(N));
  for (unsigned int i = 3; i <= sqr; i += 2) {
    unsigned int val = i + i;
    while (val < isPrime.size()) {
      isPrime[val] = false;
      val += i;
    }
  }
  for (unsigned int i = 0; i <= N; ++i) {
    if (isPrime[i] == true)
      cout << i << " ";
  }
}

int main() {
  const int N = 100;
  printPrimes(N);
  return 0;
}

The code above features two optimization:

• Even numbers are computed beforehand while odd numbers proceed skipping the even ones
• Only numbers below \( \sqrt{N} \) are considered since any other greater number would be useless to compute by recursively adding it multiple times

A quick computational analysis shows that the running time is

and since by harmonic number

we can conclude that a loose upper-bound is \( O(N\log{N}) \) (a better analysis would yield \( O(N\log{\log{N}}) \)).

Interval prime generation using \( 6k\pm1 \)

If we’re prompted with the task of generating primes in a given interval \( [x;y] \) with \( x \) and \( y \) potentially large integers but close to each other (i.e. a small interval on a high-range,) a faster approach could be exploiting the \( 6k\pm1 \) property (a variation on Euclid’s proof of infinite primes). The argument shows that any prime \( P \) is either of the form \( 6k+1 \) or \( 6k-1 \).

A prime number greater than 3 cannot have the following forms:

• \( 6k \) since it would be divisible by 2 and 3
• \( 6k\pm2 \) since it would be divisible by 2
• \( 6k\pm3 \) since it would be divisible by 3
• \( 6k\pm4 \) since it would be divisible by 2

so it is either of the form \( 6k\pm1 \) or \( 6k\pm5 \). The latter form is equivalent to the first one since increasing the index \( k \) yields \( 6(k+1)-5 = 6k+1 \). The same holds for \( 6k-1 \).

Therefore if we need to generate primes into a given interval we can simply reach a \( 6k\pm1 \) number in that interval and proceed with the primality testing. Testing for a \( 6k\pm1 \) number to be a prime is a necessary check since all primes are of the form \( 6k\pm1 \) but not all \( 6k\pm1 \) numbers are prime (e.g. \( 6k-1 \) with \( k=6 \)).

bool isPrime(unsigned int num) {
  if (num == 0 || num == 1)
    return false;
  if (num == 2 || num == 3)
    return true;
  unsigned int sq = sqrt(num);
  if (num % 2 == 0)
    return false;
  if (sq % 2 == 0)
    --sq;
  for (unsigned int i = sq; i >= 2; i -= 2) {
    if (num % i == 0)
      return false;
  }
  return true;
}

void printPrimesInInterval(unsigned int A, unsigned int B) {
  if (2 >= A && 2 <= B)
    cout << 2 << endl;
  if (3 >= A && 3 <= B)
    cout << 3 << endl;

  for (int i = 1;; ++i) {
    unsigned int value = 6 * i;
    if (value - 1 >= A && value - 1 <= B) {
      if (isPrime(value - 1))
        cout << value - 1 << endl;
    } else if (value - 1 > B)
      break;
    if (value + 1 >= A && value + 1 <= B) {
      if (isPrime(value + 1))
        cout << value + 1 << endl;
    }
  }
}

The code features many optimization spots and has been kept overly verbose for readability purposes. Complexity is \( O((y-x)\sqrt{y}) \) being the primality test the most expensive operation.