From 0ec6d144d6535a31caae75af0260351345d653cf Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 10 Jul 2026 18:11:16 +0000 Subject: [PATCH 1/2] Replace pure trial division with Miller-Rabin + Pollard rho - is_prime: deterministic Miller-Rabin (exact for all n < 3.3e24) instead of odd trial division - factorize: trial division strips factors below 10^4, then Pollard rho (Floyd cycle detection) handles the remaining cofactor; 18-digit semiprimes factor in milliseconds - tests: is_prime vs a sieve below 10^4, strong-pseudoprime and Mersenne cases, Fermat number F6, large semiprimes and prime powers, and the semiprime collapse identity at 10^18 scale - README: update the factorization note in Notes Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01898JobGjDfsgm93YW6jH6s --- README.md | 5 +-- lambda_ratio_explorer.py | 74 +++++++++++++++++++++++++++++++++------- tests/test_core.py | 38 +++++++++++++++++++++ 3 files changed, 102 insertions(+), 15 deletions(-) diff --git a/README.md b/README.md index f983053..f8c06e8 100644 --- a/README.md +++ b/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 diff --git a/lambda_ratio_explorer.py b/lambda_ratio_explorer.py index 531e298..c2d54df 100644 --- a/lambda_ratio_explorer.py +++ b/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,32 @@ def is_carmichael(n: int) -> bool: return all((n - 1) % (p - 1) == 0 for p in factors) +# Deterministic Miller-Rabin witness set: correct for all n < 3.3 * 10^24. +_MILLER_RABIN_WITNESSES = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37) + + 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 diff --git a/tests/test_core.py b/tests/test_core.py index 73bc098..9fe930e 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -47,6 +47,44 @@ 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) + + +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): From 4473b526643c458c70f8f08202ecf99541db4f7f Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 10 Jul 2026 18:37:44 +0000 Subject: [PATCH 2/2] Add base 41 to Miller-Rabin witnesses for the stated 3.3e24 bound The 12-base set (primes <= 37) is only deterministic below psi_12 = 318665857834031151167461; the advertised 3.3e24 bound requires the 13-base set including 41 (Sorenson & Webster 2015). Adds psi_12 as a regression test. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01898JobGjDfsgm93YW6jH6s --- lambda_ratio_explorer.py | 6 ++++-- tests/test_core.py | 3 +++ 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/lambda_ratio_explorer.py b/lambda_ratio_explorer.py index c2d54df..a828c61 100644 --- a/lambda_ratio_explorer.py +++ b/lambda_ratio_explorer.py @@ -347,8 +347,10 @@ def is_carmichael(n: int) -> bool: return all((n - 1) % (p - 1) == 0 for p in factors) -# Deterministic Miller-Rabin witness set: correct for all n < 3.3 * 10^24. -_MILLER_RABIN_WITNESSES = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37) +# 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: diff --git a/tests/test_core.py b/tests/test_core.py index 9fe930e..50b0d20 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -65,6 +65,9 @@ def test_large_primes_and_composites(self): 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: