Replace pure trial division with Miller–Rabin + Pollard rho

README.md

@@ -206,8 +206,9 @@ identity at runtime, so they double as tests for the theorem in
 
 ## Notes
 
-- Factorization is trial division. Practical up to `q` of order 10^5.
-  Replace with Pollard rho if you want to push further.
+- Factorization strips small factors by trial division, then switches to
+  Miller-Rabin primality testing plus Pollard rho. Comfortable with
+  18-digit semiprimes; primality is deterministic below `3.3 * 10^24`.
 - The legacy ratio `lambda(q) / log(n)` is retained in `Row` and the CLI
   for backward compatibility, but `log(n)` is just a constant scalar
   and does not enter the structural story. The interesting metrics are

lambda_ratio_explorer.py

@@ -25,6 +25,7 @@
 import argparse
 import csv
 import math
+import random
 from dataclasses import dataclass
 from functools import reduce
 from math import gcd, lcm
@@ -57,21 +58,54 @@ def carmichael_lambda_bruteforce(n: int) -> int:
     raise RuntimeError(f"bruteforce search failed for n={n}")
 
 
+# Factors below this bound are stripped by trial division; anything left
+# is handled by Miller-Rabin + Pollard rho.
+_TRIAL_DIVISION_LIMIT = 10_000
+
+
+def _pollard_rho(n: int) -> int:
+    """Return a nontrivial factor of odd composite n (Floyd cycle detection)."""
+    while True:
+        c = random.randrange(1, n)
+        x = y = random.randrange(2, n)
+        d = 1
+        while d == 1:
+            x = (x * x + c) % n
+            y = (y * y + c) % n
+            y = (y * y + c) % n
+            d = gcd(abs(x - y), n)
+        if d != n:
+            return d
+
+
+def _factor_hard(n: int) -> list[int]:
+    """Prime factors (with multiplicity) of n, which has no factor below
+    _TRIAL_DIVISION_LIMIT."""
+    if n == 1:
+        return []
+    if is_prime(n):
+        return [n]
+    d = _pollard_rho(n)
+    return _factor_hard(d) + _factor_hard(n // d)
+
+
 def factorize(n: int) -> dict[int, int]:
-    """Trial-division factorization, good enough for toy scans up to moderate q."""
+    """Prime factorization: trial division for small factors, then
+    Pollard rho for the remaining cofactor. Practical well beyond the
+    old pure-trial-division limit (e.g. 18-digit semiprimes)."""
     if n < 1:
         raise ValueError("n must be positive")
 
     factors: dict[int, int] = {}
     d = 2
-    while d * d <= n:
+    while d * d <= n and d <= _TRIAL_DIVISION_LIMIT:
         while n % d == 0:
             factors[d] = factors.get(d, 0) + 1
             n //= d
         d += 1 if d == 2 else 2  # 2, then odd candidates only
-    if n > 1:
-        factors[n] = factors.get(n, 0) + 1
-    return factors
+    for p in _factor_hard(n):
+        factors[p] = factors.get(p, 0) + 1
+    return dict(sorted(factors.items()))
 
 
 def carmichael_prime_power(p: int, k: int) -> int:
@@ -313,18 +347,34 @@ def is_carmichael(n: int) -> bool:
     return all((n - 1) % (p - 1) == 0 for p in factors)
 
 
+# Deterministic Miller-Rabin witness set: the first 13 primes are exact
+# for all n < 3,317,044,064,679,887,385,961,981 ~ 3.3 * 10^24
+# (Sorenson & Webster 2015).
+_MILLER_RABIN_WITNESSES = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41)
+
+
 def is_prime(n: int) -> bool:
+    """Deterministic Miller-Rabin, exact for all n < 3.3 * 10^24."""
     if n < 2:
         return False
-    if n == 2:
-        return True
-    if n % 2 == 0:
-        return False
-    d = 3
-    while d * d <= n:
-        if n % d == 0:
+    for p in _MILLER_RABIN_WITNESSES:
+        if n % p == 0:
+            return n == p
+    d = n - 1
+    r = 0
+    while d % 2 == 0:
+        d //= 2
+        r += 1
+    for a in _MILLER_RABIN_WITNESSES:
+        x = pow(a, d, n)
+        if x == 1 or x == n - 1:
+            continue
+        for _ in range(r - 1):
+            x = x * x % n
+            if x == n - 1:
+                break
+        else:
             return False
-        d += 2
     return True
 
 

tests/test_core.py

@@ -47,6 +47,47 @@ def test_reconstructs_n(self, n):
         assert product == n
 
 
+class TestIsPrime:
+    def test_matches_sieve_below_10000(self):
+        limit = 10000
+        sieve = [True] * (limit + 1)
+        sieve[0] = sieve[1] = False
+        for i in range(2, int(limit ** 0.5) + 1):
+            if sieve[i]:
+                for j in range(i * i, limit + 1, i):
+                    sieve[j] = False
+        for n in range(limit + 1):
+            assert is_prime(n) == sieve[n]
+
+    def test_large_primes_and_composites(self):
+        assert is_prime(2 ** 61 - 1)  # Mersenne prime
+        assert is_prime(1000000007)
+        assert not is_prime(1000000007 * 1000000009)
+        # strong pseudoprime to base 2, composite
+        assert not is_prime(3215031751)
+        # psi_12: smallest strong pseudoprime to all 12 prime bases <= 37;
+        # base 41 in the witness set must catch it
+        assert not is_prime(318665857834031151167461)
+
+
+class TestFactorizeLarge:
+    def test_fermat_number_f6(self):
+        # F6 = 2^64 + 1 = 274177 * 67280421310721, both prime
+        assert factorize(2 ** 64 + 1) == {274177: 1, 67280421310721: 1}
+
+    def test_18_digit_semiprime(self):
+        p, q = 1000000007, 1000000009
+        assert factorize(p * q) == {p: 1, q: 1}
+
+    def test_large_prime_power(self):
+        p = 1000003
+        assert factorize(p ** 3) == {p: 3}
+
+    def test_semiprime_collapse_identity_large(self):
+        p, q = 1000000007, 1000000009
+        assert collapse_index(p * q) == gcd(p - 1, q - 1)
+
+
 class TestCarmichaelLambda:
     @pytest.mark.parametrize("n", range(1, 151))
     def test_fast_matches_bruteforce(self, n):