Linear sieve in Python
List the primes up to n and the least prime factor lp(m) of every m up to n in Θ(n) time. Each composite m is written once, as p times i with p = lp(m) and p at most lp(i).
- File
linear_sieve.py, 38 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/linear-sieve/. The authors are listed in prov/0C3CNA.jsonl.
def linear_sieve(n: int) -> tuple[list[int], list[int]]:
"""Return (lp, primes): lp[m] is the least prime factor of m for 2 <= m <= n, lp[0] = lp[1] = 0,
and primes is the increasing list of primes up to n.
Precondition: n >= 0.
Cost: Theta(n) word operations; n + 1 words for lp and pi(n) words for the primes.
"""
lp = [0] * (n + 1) # @line L1
primes: list[int] = []
for i in range(2, n + 1): # @line L2
if lp[i] == 0: # @line L3
lp[i] = i # @line L4
primes.append(i)
for p in primes: # @line L5
if p > lp[i] or p * i > n: # @line L6
break
lp[p * i] = p # @line L7
return lp, primes # @line L8
def factorize(m: int, lp: list[int]) -> list[tuple[int, int]]:
"""Return the prime factorization of m as (prime, exponent) pairs in increasing order.
Precondition: 1 <= m < len(lp), and lp comes from linear_sieve.
Cost: O(log m) lookups and divisions.
"""
factors: list[tuple[int, int]] = []
while m > 1:
p, e = lp[m], 0
while m % p == 0:
m //= p
e += 1
factors.append((p, e))
return factors
test_linear_sieve.py 86 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/linear-sieve/. The authors are listed in prov/0C3CNA.jsonl.
import random
from linear_sieve import factorize, linear_sieve
def least_prime_factor_by_trial_division(m: int) -> int:
d = 2
while d * d <= m:
if m % d == 0:
return d
d += 1
return m
def factor_by_trial_division(m: int) -> list[tuple[int, int]]:
out, d = [], 2
while d * d <= m:
e = 0
while m % d == 0:
m //= d
e += 1
if e:
out.append((d, e))
d += 1
if m > 1:
out.append((m, 1))
return out
LIMIT = 40000
ORACLE_LP = [0, 0] + [least_prime_factor_by_trial_division(m) for m in range(2, LIMIT + 1)]
def check(n: int) -> None:
lp, primes = linear_sieve(n)
assert lp == ORACLE_LP[: n + 1]
assert primes == [m for m in range(2, n + 1) if ORACLE_LP[m] == m]
def test_example() -> None:
lp, primes = linear_sieve(20)
assert primes == [2, 3, 5, 7, 11, 13, 17, 19]
assert lp[2:] == [2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 17, 2, 19, 2]
def test_edges() -> None:
assert linear_sieve(0) == ([0], [])
assert linear_sieve(1) == ([0, 0], [])
assert linear_sieve(2) == ([0, 0, 2], [2])
for n in (3, 4, 8, 9, 10, 24, 25, 26, 48, 49, 120, 121, 168, 169, 289):
check(n)
def test_every_n_up_to_2000() -> None:
for n in range(2001):
check(n)
def test_random_larger_n() -> None:
rng = random.Random(1978)
for _ in range(10):
check(rng.randint(2001, LIMIT))
def test_each_composite_written_once() -> None:
n = 5000
lp, primes = linear_sieve(n)
writes = 0
for i in range(2, n + 1):
for p in primes:
if p > lp[i] or p * i > n:
break
writes += 1
assert writes == n - 1 - len(primes)
def test_factorize() -> None:
lp, _ = linear_sieve(LIMIT)
assert factorize(360, lp) == [(2, 3), (3, 2), (5, 1)]
assert factorize(1, lp) == []
rng = random.Random(87)
for m in [2, 3, 4, 39601, 39989, LIMIT] + [rng.randint(2, LIMIT) for _ in range(500)]:
assert factorize(m, lp) == factor_by_trial_division(m)
What the next level needs
- the tests did not run in this build