Lucy_Hedgehog prime counting in Python
Computes pi(x) by running the sieve of Eratosthenes on the counts S(v) for the values v = floor(x/k) only. It takes O(x^(3/4) / log x) time and O(x^(1/2)) memory.
- File
lucy_hedgehog_prime_counting.py, 54 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/lucy-hedgehog-prime-counting/. The authors are listed in prov/KWQ8J3.jsonl.
from math import isqrt
def lucy_hedgehog_prime_counting(x: int) -> int:
"""Return pi(x), the number of primes p <= x.
Precondition: x is an integer (x < 2 gives 0).
Cost: O(x^(3/4) / log x) word operations and O(x^(1/2)) words of memory.
S(v) is stored in small[v] for v <= r = isqrt(x) and in large[k] for v = x // k, 1 <= k <= r.
"""
if x < 2: # @line L1
return 0
r = isqrt(x) # @line L2
small = [0] + [v - 1 for v in range(1, r + 1)] # @line L3
large = [0] + [x // k - 1 for k in range(1, r + 1)]
for p in range(2, r + 1): # @line L4
if small[p] == small[p - 1]: # @line L5
continue # p is composite
c = small[p - 1] # @line L6
p2 = p * p
for k in range(1, min(r, x // p2) + 1): # @line L7
kp = k * p
u = large[kp] if kp <= r else small[x // kp]
large[k] -= u - c # @line L8
for v in range(r, p2 - 1, -1):
small[v] -= small[v // p] - c
return large[1] # @line L9
def lucy_prime_sum(x: int) -> int:
"""Return the sum of the primes p <= x, by the same recurrence with weights p.
Precondition: x is an integer. Cost: as lucy_hedgehog_prime_counting.
"""
if x < 2:
return 0
r = isqrt(x)
small = [0] + [v * (v + 1) // 2 - 1 for v in range(1, r + 1)]
large = [0] + [(x // k) * (x // k + 1) // 2 - 1 for k in range(1, r + 1)]
for p in range(2, r + 1):
if small[p] == small[p - 1]:
continue
c = small[p - 1]
p2 = p * p
for k in range(1, min(r, x // p2) + 1):
kp = k * p
u = large[kp] if kp <= r else small[x // kp]
large[k] -= p * (u - c)
for v in range(r, p2 - 1, -1):
small[v] -= p * (small[v // p] - c)
return large[1]
test_lucy_hedgehog_prime_counting.py 79 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/lucy-hedgehog-prime-counting/. The authors are listed in prov/KWQ8J3.jsonl.
import random
from lucy_hedgehog_prime_counting import lucy_hedgehog_prime_counting, lucy_prime_sum
def is_prime_by_trial_division(m: int) -> bool:
if m < 2:
return False
d = 2
while d * d <= m:
if m % d == 0:
return False
d += 1
return True
LIMIT = 400000
PI = [0] * (LIMIT + 1)
SUM = [0] * (LIMIT + 1)
def _fill() -> None:
flags = bytearray([1]) * (LIMIT + 1)
flags[0] = flags[1] = 0
p = 2
while p * p <= LIMIT:
if flags[p]:
flags[p * p :: p] = bytes(len(range(p * p, LIMIT + 1, p)))
p += 1
c = s = 0
for i in range(LIMIT + 1):
if flags[i]:
c += 1
s += i
PI[i] = c
SUM[i] = s
_fill()
def test_table_against_trial_division() -> None:
c = 0
for m in range(5001):
c += is_prime_by_trial_division(m)
assert PI[m] == c
def test_every_x_up_to_2000() -> None:
for x in range(-3, 2001):
assert lucy_hedgehog_prime_counting(x) == PI[max(x, 0)]
assert lucy_prime_sum(x) == SUM[max(x, 0)]
def test_squares_of_primes_and_neighbours() -> None:
for q in (2, 3, 5, 7, 11, 13, 31, 97, 101, 211, 307, 541, 631):
for x in (q * q - 1, q * q, q * q + 1, q**4 - 1, q**4, q**4 + 1):
if 0 <= x <= LIMIT:
assert lucy_hedgehog_prime_counting(x) == PI[x]
assert lucy_prime_sum(x) == SUM[x]
def test_random() -> None:
rng = random.Random(2013)
for _ in range(80):
x = rng.randint(2001, LIMIT)
assert lucy_hedgehog_prime_counting(x) == PI[x]
assert lucy_prime_sum(x) == SUM[x]
def test_known_values() -> None:
assert lucy_hedgehog_prime_counting(100) == 25
assert lucy_hedgehog_prime_counting(10**6) == 78498
assert lucy_hedgehog_prime_counting(10**8) == 5761455
assert lucy_hedgehog_prime_counting(10**9) == 50847534
assert lucy_prime_sum(2 * 10**6) == 142913828922
What the next level needs
- the tests did not run in this build