Primes: trial division, sieves and factorisation

A prime is a whole number greater than 1 whose only divisors are 1 and itself: 2, 3, 5, 7, 11 and so on. Three tools cover almost every prime question in DSA. To test one number n, trial division tries divisors up to √n, which costs O(√n), because a divisor larger than √n always comes paired with one smaller than √n. To list every prime up to n, the sieve of Eratosthenes crosses out the multiples of each prime, starting at its square, in O(n log log n) time and O(n) memory. To factorise many numbers up to n, a smallest-prime-factor sieve stores each number’s smallest prime factor, and any number then splits into primes with one lookup per factor, O(log n) each.

Trial division up to the square root

If n = d × e with d ≤ e, then d × d ≤ n, so d ≤ √n: every composite number has a divisor no larger than its square root. Trial division can therefore stop at √n, and after 2 it needs to try only odd divisors, about √n / 2 divisions for a prime:

import math


def is_prime(n):
    """Trial division by 2 and by odd numbers up to the integer square root. Returns (answer, divisions made)."""
    if n < 2:
        return False, 0
    if n % 2 == 0:
        return n == 2, 1
    divisions = 1
    for d in range(3, math.isqrt(n) + 1, 2):
        divisions += 1
        if n % d == 0:
            return False, divisions
    return True, divisions


for n in [97, 91, 1_000_003, 999_999_937, 10**12 + 39, 10**12 + 41]:
    prime, divisions = is_prime(n)
    print(f"{n:>17,}  {'prime' if prime else 'composite':<9}  {divisions:>7,} divisions")

print()
print("bound for 10^16 - 1:      int(math.sqrt(n)) =", int(math.sqrt(10**16 - 1)), " math.isqrt(n) =", math.isqrt(10**16 - 1))
print("bound for (2^53 + 1)^2:   int(math.sqrt(n)) =", int(math.sqrt((2**53 + 1) ** 2)), " math.isqrt(n) =", math.isqrt((2**53 + 1) ** 2))

               97  prime            5 divisions
               91  composite        4 divisions
        1,000,003  prime          500 divisions
      999,999,937  prime       15,811 divisions
1,000,000,000,039  prime      500,000 divisions
1,000,000,000,041  composite        2 divisions

bound for 10^16 - 1:      int(math.sqrt(n)) = 100000000  math.isqrt(n) = 99999999
bound for (2^53 + 1)^2:   int(math.sqrt(n)) = 9007199254740992  math.isqrt(n) = 9007199254740993

A prime near 10¹² takes 500,000 divisions, while most composites stop early: 10¹² + 41 is divisible by 3, so the second division settles it. The bound itself has to be exact. int(math.sqrt(n)) turns n into a float first, which gave 100,000,000 for 10¹⁶ − 1, one too large, and 2⁵³ for (2⁵³ + 1)², one too small: for the square of a prime that large, the loop would stop one short of the only divisor it has to try. math.isqrt returns the exact integer square root of any whole number.

The sieve of Eratosthenes

To list the primes up to n, mark every number from 2 to n as possibly prime. The smallest number still marked, p, is prime, since no smaller prime divides it. Cross out its multiples from p² on: each smaller multiple k × p with k < p has a prime factor smaller than p and is already crossed out. Move to the next marked number, and stop once p² is larger than n; every number still marked is prime.

This program sieves up to 10⁷ with a bytearray, one byte a number, and counts how many times a number is crossed out:

import math


def sieve(n):
    """is_prime[i] is 1 exactly when i is prime, for 0 <= i <= n (the sieve of Eratosthenes)."""
    is_prime = bytearray([1]) * (n + 1)
    is_prime[0:2] = b"\x00\x00"
    for p in range(2, math.isqrt(n) + 1):
        if is_prime[p]:
            # p is prime: cross out p*p, p*p + p, ..., up to n
            is_prime[p * p :: p] = bytes(len(range(p * p, n + 1, p)))
    return is_prime


is_prime = sieve(100)
print("primes up to 100:", [i for i in range(101) if is_prime[i]])

is_prime = sieve(10**7)
print(f"primes up to 10^7: {sum(is_prime):,}")
print()
print("         n   crossings   crossings / n   ln ln n")
for n in (10**4, 10**5, 10**6, 10**7):
    crossings = sum((n - p * p) // p + 1 for p in range(2, math.isqrt(n) + 1) if is_prime[p])
    print(f"{n:>10,}  {crossings:>10,}  {crossings / n:>14.3f}  {math.log(math.log(n)):>8.3f}")

primes up to 100: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
primes up to 10^7: 664,579

         n   crossings   crossings / n   ln ln n
    10,000      16,981           1.698     2.220
   100,000     193,078           1.931     2.443
 1,000,000   2,122,048           2.122     2.626
10,000,000  22,850,051           2.285     2.780

It finds the 664,579 primes up to 10⁷. While n grows a thousandfold, the crossings per number grow only from 1.70 to 2.29, in step with ln ln n, which is where the O(n log log n) comes from. The memory is one byte a number, 10 MB for 10⁷.

Smallest prime factors: factorising many numbers

A sieve can record more than yes or no. Keep for every number x its smallest prime factor, spf[x], set by the first prime that crosses x out. Then x splits into primes by repeated division: spf[x] is one factor, and x / spf[x] holds the rest. A number up to n has at most log₂ n prime factors, so each factorisation takes O(log n) lookups:

import math


def smallest_prime_factors(n):
    """spf[i] is the smallest prime that divides i, for 2 <= i <= n."""
    spf = list(range(n + 1))
    for p in range(2, math.isqrt(n) + 1):
        if spf[p] == p:  # nothing smaller divides p, so p is prime
            for m in range(p * p, n + 1, p):
                if spf[m] == m:
                    spf[m] = p
    return spf


def factorise(x, spf):
    """The prime factors of x, smallest first, with one table lookup per factor."""
    factors = []
    while x > 1:
        factors.append(spf[x])
        x //= spf[x]
    return factors


def trial_factorise(x):
    """The same factors by trial division, and the number of divisions it made."""
    factors, divisions, d = [], 0, 2
    while d * d <= x:
        divisions += 1
        if x % d == 0:
            factors.append(d)
            x //= d
        else:
            d += 1
    if x > 1:
        factors.append(x)
    return factors, divisions


N = 1_000_000
spf = smallest_prime_factors(N)
for x in [360, 999_983, 524_288, 735_134, 997_002]:
    print(f"{x:>7,} = {' × '.join(map(str, factorise(x, spf)))}")

lookups = divisions = 0
for x in range(2, 100_001):
    by_table = factorise(x, spf)
    by_trial, made = trial_factorise(x)
    assert by_table == by_trial
    lookups += len(by_table)
    divisions += made
print(f"factorising every number from 2 to 100,000: {lookups:,} table lookups, {divisions:,} trial divisions")

    360 = 2 × 2 × 2 × 3 × 3 × 5
999,983 = 999983
524,288 = 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2 × 2
735,134 = 2 × 31 × 71 × 167
997,002 = 2 × 3 × 3 × 3 × 37 × 499
factorising every number from 2 to 100,000: 343,614 table lookups, 6,956,686 trial divisions

The most lookups below 10⁶ go to 524,288 = 2¹⁹, nineteen of them. Factorising every number from 2 to 100,000 took 343,614 lookups with the table and 6,956,686 trial divisions, about twenty times as many.

Choosing a method

Trial division needs no memory and suits a handful of questions. The sieve answers a primality question in one lookup after one pass, and the factor table answers factorisation questions the same way. Three further tools reach where these stop. A segmented sieve finds the primes in a far window, such as 10¹⁴ to 10¹⁴ + 10⁶, with memory for the window only. The Miller–Rabin test decides whether one large number is prime with a few modular powers. A linear sieve crosses out every composite exactly once.

Exercise and self-check

Count the primes in up to 100,000 ranges of numbers up to 10⁷, sieving only once.

def count_primes(queries):
    """For each (l, r), the number of primes p with l <= p <= r."""
    raise NotImplementedError("write count_primes")

import random

from ranges import count_primes


def slow_is_prime(n):
    """Trial division: right, and too slow for large ranges."""
    if n < 2:
        return False
    d = 2
    while d * d <= n:
        if n % d == 0:
            return False
        d += 1
    return True


def test_single_numbers():
    """one-number ranges: 1 is not prime, 2 and 97 are, 91 is not"""
    assert count_primes([(1, 1), (2, 2), (97, 97), (91, 91)]) == [0, 1, 1, 0]


def test_tiny_inputs():
    """a call whose largest r is 1 or 2 sieves only that far"""
    assert count_primes([(1, 1)]) == [0]
    assert count_primes([(2, 2)]) == [1]
    assert count_primes([(1, 2), (1, 1)]) == [1, 0]


def test_small_ranges():
    """ranges that start at 1 or 2, and a range without primes"""
    assert count_primes([(1, 10), (2, 10), (1, 100), (90, 100), (114, 126)]) == [4, 4, 25, 1, 0]


def test_order_is_kept():
    """the answers come back in the order of the queries"""
    assert count_primes([(1, 100), (1, 10), (11, 20)]) == [25, 4, 4]


def test_random_against_trial_division():
    """300 random ranges below 3,000, compared with trial division"""
    rng = random.Random(9)
    queries = []
    for _ in range(300):
        l = rng.randint(1, 3000)
        queries.append((l, rng.randint(l, 3000)))
    expected = [sum(1 for x in range(l, r + 1) if slow_is_prime(x)) for l, r in queries]
    assert count_primes(queries) == expected


def test_up_to_a_million():
    """78,498 primes up to 10^6, and the ones in its last thousand numbers"""
    assert count_primes([(1, 10**6), (999_001, 10**6)]) == [78_498, 65]


def test_ten_thousand_ranges():
    """10,000 ranges below 10^6: their answers add up to a known total"""
    rng = random.Random(4)
    queries = []
    for _ in range(10_000):
        l = rng.randint(1, 10**6)
        queries.append((l, min(10**6, l + rng.randint(0, 5_000))))
    answers = count_primes(queries)
    assert len(answers) == 10_000
    assert sum(answers) == 1_977_492

/** For each [l, r], the number of primes p with l <= p <= r. */
export function countPrimes(queries) {
  throw new Error('write countPrimes');
}

import { test, assert } from 'mysmartcopilot:test';
import { countPrimes } from './ranges.mjs';

/** Trial division: right, and too slow for large ranges. */
function slowIsPrime(n) {
  if (n < 2) return false;
  for (let d = 2; d * d <= n; d++) if (n % d === 0) return false;
  return true;
}

/** A small seeded generator (Park and Miller), so every run builds the same ranges. */
function generator(seed) {
  const random = () => (seed = (seed * 48271) % 2147483647) / 2147483647;
  return (k) => Math.floor(random() * k);
}

test('one-number ranges: 1 is not prime, 2 and 97 are, 91 is not', () => {
  assert.deepEqual(countPrimes([[1, 1], [2, 2], [97, 97], [91, 91]]), [0, 1, 1, 0]);
});

test('a call whose largest r is 1 or 2 sieves only that far', () => {
  assert.deepEqual(countPrimes([[1, 1]]), [0]);
  assert.deepEqual(countPrimes([[2, 2]]), [1]);
  assert.deepEqual(countPrimes([[1, 2], [1, 1]]), [1, 0]);
});

test('ranges that start at 1 or 2, and a range without primes', () => {
  assert.deepEqual(countPrimes([[1, 10], [2, 10], [1, 100], [90, 100], [114, 126]]), [4, 4, 25, 1, 0]);
});

test('the answers come back in the order of the queries', () => {
  assert.deepEqual(countPrimes([[1, 100], [1, 10], [11, 20]]), [25, 4, 4]);
});

test('300 random ranges below 3,000, compared with trial division', () => {
  const pick = generator(9);
  const queries = [];
  for (let i = 0; i < 300; i++) {
    const l = 1 + pick(3000);
    queries.push([l, l + pick(3001 - l)]);
  }
  const expected = queries.map(([l, r]) => {
    let c = 0;
    for (let x = l; x <= r; x++) if (slowIsPrime(x)) c++;
    return c;
  });
  assert.deepEqual(countPrimes(queries), expected);
});

test("78,498 primes up to 10^6, and the ones in its last thousand numbers", () => {
  assert.deepEqual(countPrimes([[1, 1000000], [999001, 1000000]]), [78498, 65]);
});

test('10,000 ranges below 10^6: their answers add up to a known total', () => {
  const pick = generator(4);
  const queries = [];
  for (let i = 0; i < 10000; i++) {
    const l = 1 + pick(1000000);
    queries.push([l, Math.min(1000000, l + pick(5001))]);
  }
  const answers = countPrimes(queries);
  assert.equal(answers.length, 10000);
  assert.equal(answers.reduce((s, x) => s + x, 0), 1963500);
});

Why the sieve is right and what it costs, proved and counted Pro: opens with a Pro, Premium or Ultimate pass

Primes in a far window, and numbers too large to sieve Pro: opens with a Pro, Premium or Ultimate pass

The factor table in five languages, the 32-bit trap and a practice set Pro: opens with a Pro, Premium or Ultimate pass

Interview questions

Key takeaways

• Trial division stops at the integer square root: a composite always has a divisor no larger than √n. Compute the bound with math.isqrt, never through a float.
• The sieve of Eratosthenes crosses out the multiples of each prime from p² and stops at √n: O(n log log n) time and one byte (or one bit) per number.
• A smallest-prime-factor table factorises any number up to n with one lookup per prime factor, at most log₂ n.
• Choose by the workload: a few numbers, trial division; many primality questions, a sieve; many factorisations, a factor table.