Primality test

A primality test decides whether a number is prime: greater than 1, with no divisors other than 1 and itself. Trial division, the test in this article, needs at most about 10^6 tries for one number up to 10^{12}. To find every prime in a range, use the sieve of Eratosthenes instead.

Explanation

Is 13 prime? Try every number from 2 to 12 and look at the remainder when 13 is divided by it. None of the 11 remainders is 0, so no number between 1 and 13 divides 13, and 13 is prime. Is 15 prime? Try 2: the remainder is 1. Try 3: the remainder is 0, so 15 = 3 · 5. That already proves 15 is not prime, and there is no need to try 4 or anything after it.

Fig. 1. Trying divisors of 13 and 15. Red divides, grey does not, faded numbers are never tried.

A number greater than 1 that is not prime is composite.

Problems often give a number up to about 2 \cdot 10^9. Take 1 999 999 973, the largest prime below 2 \cdot 10^9. No number from 2 to 1 999 999 972 divides it, so all of them are tried: almost 2 \cdot 10^9 numbers. By the rule of thumb from time complexity, about 10^8 simple steps per second, that is about 20 seconds.

Most composite numbers have a small divisor: every second number is divisible by 2, every third by 3. A natural shortcut is to try only 2 to 10, which is 9 tries for any number. Is 143 prime? It is odd, and none of 2, 3, …, 10 divides it, so the shortcut answers "prime". But 143 = 11 · 13.

If 143 is a product of two numbers greater than 1, the smaller one is at most 11. Two numbers that are both 12 or more multiply to at least 12 · 12 = 144, which is already more than 143. So trying 2 to 11 is enough.

Call 11 the limit for 143: the largest number whose square is at most 143, which is the square root of 143 rounded down. The same holds for every number: if it is a · b with 1 < a ≤ b, then a · a is at most a · b, the number itself, so a is at most the limit. If nothing from 2 to the limit divides the number, no such pair exists and the number is prime.

Fig. 2. The smaller divisor of 143 is left of the dashed line. The shortcut stops at 10 and misses it.

For a prime near 2 \cdot 10^9 the limit is 44 721, because 44 721 · 44 721 = 1 999 967 841 and 44 722 · 44 722 = 2 000 057 284. That is 44 720 tries instead of 2 \cdot 10^9.

This method is called trial division.

Rule out 0 and 1: 139 is at least 2, so go on.

Try the next number, 2: 2 · 2 = 4 ≤ 139, and 139 % 2 = 1, not 0. The yellow cell is the number tried.

Try the next number, 3: 3 · 3 = 9 ≤ 139, and 139 % 3 = 1, not 0.

Try the next number, 4: 4 · 4 = 16 ≤ 139, and 139 % 4 = 3, not 0.

Try the next number, 5: 5 · 5 = 25 ≤ 139, and 139 % 5 = 4, not 0.

Try the next number, 6: 6 · 6 = 36 ≤ 139, and 139 % 6 = 1, not 0.

Try the next number, 7: 7 · 7 = 49 ≤ 139, and 139 % 7 = 6, not 0.

Try the next number, 8: 8 · 8 = 64 ≤ 139, and 139 % 8 = 3, not 0.

Try the next number, 9: 9 · 9 = 81 ≤ 139, and 139 % 9 = 4, not 0.

Try the next number, 10: 10 · 10 = 100 ≤ 139, and 139 % 10 = 9, not 0.

Try the next number, 11: 11 · 11 = 121 ≤ 139, and 139 % 11 = 7, not 0.

The next number is 12: 12 · 12 = 144 > 139, so the loop stops without trying it.

No divisor found among 2 to 11: 139 is prime, and is_prime returns true.

1/13

Fig. 3. Step through with the buttons or the arrow keys. Grey numbers leave a remainder.

Implementation

The number tested is n, and the number tried is i. "Its square is at most n" is the loop condition i * i <= n, and "divides" is n % i == 0. The test goes into a function bool is_prime(long long n) that returns from inside the loop as soon as it knows the answer.

1. Rule out 0 and 1: if n < 2, return false.

2. Try the next number: i = 2, 3, … while i * i <= n.

3. Stop at the first divisor: if n % i == 0, return false.

4. No divisor found: return true.

The input is one integer n, 0 \le n \le 10^{12}. The program prints YES if n is prime and NO otherwise. For 143 the input is:

143

#include <iostream>
using namespace std;

bool is_prime(long long n) {
    if (n < 2) return false;                 // rule out 0 and 1
    for (long long i = 2; i * i <= n; i++)   // try the next number
        if (n % i == 0) return false;        // stop at the first divisor
    return true;                             // no divisor found
}

int main() {
    long long n;
    cin >> n;
    cout << (is_prime(n) ? "YES" : "NO") << '\n';
}

The condition i * i <= n uses only integers, so it is exact, while sqrt(n) works with floating point. The <= matters: for the square of a prime, such as 25 = 5 · 5 or 121 = 11 · 11, the limit itself is the only divisor to find, and i * i < n never tries it.

The loop makes at most about \sqrt{n} tries, so the test takes O(\sqrt{n}) time, about 10^6 tries for a prime near 10^{12}, and $O(1)$ memory.

If the function returns i instead of false, and n instead of true, it gives the smallest divisor of n greater than 1, which is n itself when n is prime.

Trying 2 and then only odd i halves the number of tries, but the time stays O(\sqrt{n}). If you rule out even numbers with if (n % 2 == 0) return false;, put if (n == 2) return true; before it, or 2 comes out composite.

Numbers up to 10^18

For a prime near 10^{18} trial division needs 10^9 tries, about 10 seconds by the rule of thumb, and some problems ask about thousands of such numbers. The Miller–Rabin test answers with 12 modular powers instead, each about 60 squarings with binary exponentiation. The function below gives the exact answer for every n below 2^{64}; copy it as it is. It uses __int128 for the products, which GCC and Clang support. The Miller–Rabin test explains how it works.

typedef unsigned long long ull;

ull power(ull a, ull e, ull m) {             // a^e mod m, binary exponentiation
    ull r = 1;
    a %= m;
    while (e > 0) {
        if (e & 1) r = (__int128)r * a % m;
        a = (__int128)a * a % m;
        e >>= 1;
    }
    return r;
}

bool is_prime(ull n) {
    if (n < 2) return false;
    for (ull p : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37})
        if (n % p == 0) return n == p;
    ull d = n - 1;
    int s = 0;
    while (d % 2 == 0) { d /= 2; s++; }
    for (ull a : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
        ull x = power(a, d, n);
        if (x == 1 || x == n - 1) continue;
        bool composite = true;
        for (int r = 1; r < s; r++) {
            x = (__int128)x * x % n;
            if (x == n - 1) { composite = false; break; }
        }
        if (composite) return false;
    }
    return true;
}

Use it like the first is_prime: read n as unsigned long long and call is_prime(n). You don't need to follow the lines inside to use it.

Similar algorithms

The sieve of Eratosthenes marks every prime up to about 10^7 at once in O(n \log \log n) time, so it is the tool when many numbers from one range are asked. Trial division on every number up to 10^7 makes about 1.7 \cdot 10^9 tries, about 17 seconds.

Prime factorization runs the same loop up to the limit, but divides out each divisor it finds and goes on; the primality test stops at the first one.

Listing the divisors of a number also runs the same loop, collecting each divisor and its partner instead of stopping.

Practice