| #generate a list of small primes represented as a u64 array that fit into an OS page. | |
| #calculate their inverse mod u64 so we can test them without division | |
| import subprocess | |
| import os | |
| import math | |
| from typing import List, Tuple | |
| def get_page_size() -> int: | |
| """Get the system's page size""" | |
| return os.sysconf("SC_PAGE_SIZE") | |
| def is_prime_using_factor(n: int) -> bool: | |
| """Use GNU factor to check if a number is prime""" | |
| result = subprocess.run(['factor', str(n)], | |
| capture_output=True, | |
| text=True) | |
| # Parse output - format is "n: factors" | |
| factors = result.stdout.strip().split(': ')[1].split() | |
| # If there's only one factor and it equals n, it's prime | |
| return len(factors) == 1 and int(factors[0]) == n | |
| def extended_gcd(a: int, b: int) -> Tuple[int, int, int]: | |
| """Extended Euclidean Algorithm""" | |
| if a == 0: | |
| return b, 0, 1 | |
| gcd, x1, y1 = extended_gcd(b % a, a) | |
| x = y1 - (b // a) * x1 | |
| y = x1 | |
| return gcd, x, y | |
| def mod_inverse(a: int, m: int) -> int: | |
| """Calculate modular multiplicative inverse""" | |
| gcd, x, _ = extended_gcd(a, m) | |
| if gcd != 1: | |
| raise ValueError(f"Modular inverse does not exist for {a} mod {m}") | |
| return x % m | |
| def generate_primes_and_inverses(max_size: int) -> List[Tuple[int, int]]: | |
| """Generate primes and their mod inverses up to max_size bytes""" | |
| UINT64_MAX = (1 << 64) - 1 | |
| primes_and_inverses = [] | |
| total_size = 0 | |
| n = 3 # Start with 3 as we want odd primes | |
| while total_size < max_size: | |
| if is_prime_using_factor(n): | |
| try: | |
| inv = mod_inverse(n, UINT64_MAX + 1) | |
| # Each tuple is 16 bytes (8 bytes for prime + 8 bytes for inverse) | |
| size_increase = 16 | |
| if total_size + size_increase <= max_size: | |
| primes_and_inverses.append((n, inv)) | |
| total_size += size_increase | |
| else: | |
| break | |
| except ValueError: | |
| pass # Skip if no modular inverse exists | |
| n += 2 # Check only odd numbers | |
| return primes_and_inverses | |
| def main(): | |
| page_size = get_page_size() | |
| print(f"System page size: {page_size} bytes") | |
| # Generate primes and inverses to fill one page | |
| primes_and_inverses = generate_primes_and_inverses(page_size) | |
| # Print results in a format suitable for use in systems programming | |
| print("\nPrime number table with modular inverses (u64):") | |
| print("const PRIME_INVERSES: [(u64, u64)] = [") | |
| for prime, inv in primes_and_inverses: | |
| print(f" ({prime}, {inv}),") | |
| print("];") | |
| print(f"\nGenerated {len(primes_and_inverses)} prime-inverse pairs") | |
| print(f"Total size: {len(primes_and_inverses) * 16} bytes") | |
| # Verify correctness | |
| print("\nVerifying correctness...") | |
| UINT64_MAX = (1 << 64) - 1 | |
| for prime, inv in primes_and_inverses: | |
| # Check if (prime * inverse) mod 2^64 ≡ 1 | |
| if (prime * inv) % (UINT64_MAX + 1) != 1: | |
| print(f"ERROR: Invalid inverse for prime {prime}") | |
| break | |
| else: | |
| print("All modular inverses verified correct!") | |
| if __name__ == "__main__": | |
| main() |
Twitter thread https://x.com/SMT_Solvers/status/1871600292762181681
What is the best size table?
Seems you would want have a "batch" or "buffer" mode -b , which stores the sieve in memory and then re-uses it.
https://en.wikipedia.org/wiki/Cycle_detection
How do Floyd,Brent, array of (Floyd, Brent) perform emperically? Seems you coudld do some SIMD speedup as exhaustive testing up to u32 is feasable.
| from typing import List, Optional, Dict | |
| from dataclasses import dataclass | |
| from math import gcd | |
| @dataclass | |
| class Factors: | |
| """Track the prime factorization of a number""" | |
| factors: Dict[int, int] = None # factor -> exponent | |
| def __init__(self): | |
| self.factors = {} | |
| def insert(self, p: int, e: int = 1): | |
| """Insert a prime factor p with exponent e""" | |
| if p in self.factors: | |
| self.factors[p] += e | |
| else: | |
| self.factors[p] = e | |
| def miller_rabin(n: int, a: int) -> bool: | |
| """ | |
| Miller-Rabin primality test for a single base a | |
| """ | |
| if n == 2: | |
| return True | |
| if n % 2 == 0 or n <= 1: | |
| return False | |
| # Write n as 2^s * d | |
| s = 0 | |
| d = n - 1 | |
| while d % 2 == 0: | |
| s += 1 | |
| d //= 2 | |
| # Compute a^d % n using square-and-multiply | |
| x = pow(a, d, n) | |
| if x == 1 or x == n - 1: | |
| return True | |
| # Keep squaring x while we haven't tried all possibilities | |
| for _ in range(s - 1): | |
| x = (x * x) % n | |
| if x == n - 1: | |
| return True | |
| if x == 1: | |
| return False | |
| return False | |
| def is_prime(n: int) -> bool: | |
| """ | |
| Determine if n is prime using Miller-Rabin test with first few prime bases | |
| """ | |
| if n <= 1: | |
| return False | |
| if n <= 3: | |
| return True | |
| if n % 2 == 0: | |
| return False | |
| # First few prime bases give strong guarantees for numbers in practical ranges | |
| witnesses = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37] | |
| for a in witnesses: | |
| if n == a: | |
| return True | |
| if not miller_rabin(n, a): | |
| return False | |
| return True | |
| def pollard_rho(n: int, factors: Factors, a: int = 1) -> None: | |
| """ | |
| Pollard's rho algorithm with Montgomery arithmetic optimization | |
| Args: | |
| n: Number to factor | |
| factors: Factors object to store the factorization | |
| a: Parameter for the polynomial x^2 + a | |
| """ | |
| if n <= 1: | |
| return | |
| if is_prime(n): | |
| factors.insert(n) | |
| return | |
| # Find p using Pollard's rho algorithm | |
| def g(x: int) -> int: | |
| """The iteration function (polynomial) used for Pollard rho""" | |
| return (x * x + a) % n | |
| # Floyd's cycle-finding algorithm | |
| x = 2 | |
| y = 2 | |
| d = 1 | |
| while d == 1: | |
| x = g(x) | |
| y = g(g(y)) | |
| d = gcd(abs(x - y), n) | |
| if d == n: | |
| # If we found n itself as a factor, try again with different parameter | |
| pollard_rho(n, factors, a + 1) | |
| return | |
| # Factor each part recursively | |
| if is_prime(d): | |
| factors.insert(d) | |
| else: | |
| pollard_rho(d, factors, a + 1) | |
| n //= d | |
| if is_prime(n): | |
| factors.insert(n) | |
| else: | |
| pollard_rho(n, factors, a) | |
| def factor(n: int) -> Dict[int, int]: | |
| """ | |
| Factor a number into its prime factors | |
| Returns dict mapping prime factors to their exponents | |
| """ | |
| factors = Factors() | |
| # Handle small factors first | |
| for p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]: | |
| while n % p == 0: | |
| factors.insert(p) | |
| n //= p | |
| if n > 1: | |
| if is_prime(n): | |
| factors.insert(n) | |
| else: | |
| pollard_rho(n, factors) | |
| return factors.factors | |
| def main(): | |
| # Test the implementation with some examples | |
| test_numbers = [ | |
| 2**67 - 1, # Mersenne number | |
| (2**31 - 1) * (2**19 - 1), # Product of two Mersenne primes | |
| 341550071728321, # A large semiprime | |
| 2**32 + 1, # Fermat number | |
| ] | |
| for n in test_numbers: | |
| print(f"\nFactoring {n}:") | |
| factors = factor(n) | |
| # Print in factor^exponent format | |
| factorization = " * ".join(f"{p}^{e}" if e > 1 else str(p) | |
| for p, e in sorted(factors.items())) | |
| print(f"= {factorization}") | |
| if __name__ == "__main__": | |
| main() |
| # https://eprint.iacr.org/2021/933 | |
| import sympy as sp | |
| import numpy as np | |
| from math import log, e | |
| def create_rn_f_matrix(n, N): | |
| """Create the Rn,f matrix using sympy for prime generation""" | |
| # Get first n primes using sympy | |
| primes = list(sp.primerange(2, 1000))[:n] | |
| # Calculate N0 | |
| N0 = N ** (1/(n+1)) | |
| # Create matrix | |
| matrix = np.zeros((n+1, n+1)) | |
| # Set diagonal elements for first n rows | |
| for i in range(n): | |
| matrix[i][i] = i + 1 | |
| # Set last row | |
| for i in range(n): | |
| matrix[n][i] = N0 * log(primes[i]) | |
| matrix[n][n] = N0 * log(N) | |
| return matrix, primes | |
| def find_smooth_relations(N, primes, bound): | |
| """Find relations using sympy's factorint""" | |
| relations = [] | |
| for a in range(2, bound): | |
| val = (a * a) % N | |
| factors = sp.factorint(val) | |
| # Check if smooth over our prime base | |
| if all(p in primes for p in factors.keys()): | |
| relations.append((a, val, factors)) | |
| if len(relations) >= len(primes) + 1: | |
| break | |
| return relations | |
| def solve_for_factors(relations, N): | |
| """Try to find factors using the relations""" | |
| for i, (a, val, _) in enumerate(relations): | |
| # Try simple cases first | |
| gcd = sp.gcd(a + 1, N) | |
| if 1 < gcd < N: | |
| return gcd, N // gcd | |
| # If simple cases don't work, try combining relations | |
| n_rels = len(relations) | |
| for mask in range(1, 1 << n_rels): | |
| x = 1 | |
| y = 1 | |
| for i in range(n_rels): | |
| if mask & (1 << i): | |
| a, _, _ = relations[i] | |
| x = (x * a) % N | |
| y = (y * ((a * a) % N)) % N | |
| # Check if we found factors | |
| y_root = sp.sqrt_mod(y, N) | |
| if y_root is not None: | |
| gcd1 = sp.gcd(x + y_root, N) | |
| if 1 < gcd1 < N: | |
| return gcd1, N // gcd1 | |
| gcd2 = sp.gcd(x - y_root, N) | |
| if 1 < gcd2 < N: | |
| return gcd2, N // gcd2 | |
| return None | |
| def lattice_factor(N, n=47): | |
| """ | |
| Main factoring function using lattice methods with sympy | |
| N: number to factor | |
| n: dimension parameter (47 for N ≈ 2^400, 95 for N ≈ 2^800) | |
| """ | |
| if N < 2: | |
| return None | |
| if N % 2 == 0: | |
| return (2, N//2) | |
| # Quick check for prime | |
| if sp.isprime(N): | |
| return None | |
| # Create matrix and get prime base | |
| matrix, primes = create_rn_f_matrix(n, N) | |
| # Find relations | |
| bound = int(sp.sqrt(N)) | |
| relations = find_smooth_relations(N, primes, bound) | |
| if len(relations) < n + 1: | |
| return None | |
| # Try to find factors using relations | |
| return solve_for_factors(relations, N) | |
| # Example usage: | |
| def test_factorization(): | |
| # Test with a smaller number first | |
| N = 15347 | |
| print(f"Attempting to factor {N}") | |
| factors = lattice_factor(N, n=10) | |
| if factors: | |
| print(f"Found factors: {factors}") | |
| # Verify | |
| a, b = factors | |
| print(f"Verification: {a} * {b} = {a * b}") | |
| print(f"Original N: {N}") | |
| else: | |
| print("Failed to factor") | |
| if __name__ == "__main__": | |
| test_factorization() |