|
|
楼主 |
发表于 2026-8-19 11:18
|
显示全部楼层
什么是素数先要厘清：
#!/usr/bin/env python3
"""
╔══════════════════════════════════════════════════════════════════╗
║ PrimeGenerator v3.0 (Ultimate) ║
║ 核心哲学：乘法定位合数 → 加法标记合数 → 质数自然浮现 ║
╚══════════════════════════════════════════════════════════════════╝
"""
from __future__ import annotations
import math
import random
import time
from typing import Iterator, List, Optional, Tuple
from dataclasses import dataclass
from functools import lru_cache
# ============================================================
# 全局小素数缓存（启动即加载，加速所有试除/预筛逻辑）
# ============================================================
def _load_small_primes(limit: int = 10000) -> List[int]:
"""快速生成小素数表，用于后续预筛"""
sieve = bytearray(b"\x00") * (limit + 1)
primes = []
for i in range(2, limit + 1):
if not sieve:
primes.append(i)
for m in range(i * i, limit + 1, i):
sieve[m] = 1
return primes
SMALL_PRIMES = _load_small_primes(10000)
SMALL_PRIME_SET = frozenset(SMALL_PRIMES)
# ============================================================
# 错误类型
# ============================================================
class PrimeError(Exception):
"""素数模块基础异常"""
pass
class NotFoundError(PrimeError):
"""在指定范围内未找到素数"""
pass
# ============================================================
# 性能统计
# ============================================================
@dataclass
class PrimeStats:
"""性能统计"""
count: int = 0
elapsed_ms: float = 0.0
method: str = ""
@property
def rate(self) -> float:
return self.count / self.elapsed_ms * 1000 if self.elapsed_ms > 0 else 0.0
# ============================================================
# 主类
# ============================================================
class PrimeGenerator:
"""
工业级素数生成器
功能覆盖：
✅ 基础埃氏筛（bytearray 内存优化）
✅ 奇数-only 筛（内存减半）
✅ 线性筛（O(n) 时间，每个合数只被标记一次）
✅ 分段筛（超大区间，缓存友好）
✅ 迭代器 / 生成器接口
✅ 米勒-拉宾确定性素性测试
✅ 大素数生成（密码学级）
✅ 第 N 个素数 / 素数计数 π(x)
✅ 质因数分解（Pollard-Rho）
✅ 孪生素数 / 梅森素数检测
"""
# ----------------------------------------------------------
# ① 经典埃氏筛（简洁版）
# ----------------------------------------------------------
@staticmethod
def sieve(n: int) -> List[int]:
"""
经典埃拉托斯特尼筛法
返回 [2, n] 内所有素数
时间 O(n log log n)，空间 O(n)
"""
if n < 2:
return []
if n == 2:
return [2]
# bytearray: 0=素数候选, 1=合数
composite = bytearray(n + 1)
composite[0] = composite[1] = 1
primes: List[int] = []
sqrt_n = math.isqrt(n)
for i in range(2, sqrt_n + 1):
if not composite:
primes.append(i)
# 乘法定位：i*i 起跳，加法步进标记
composite[i*i : n+1 : i] = b"\x01" * ((n - i*i) // i + 1)
# 收集 sqrt(n) 以上的素数
for i in range(max(sqrt_n + 1, 3), n + 1):
if not composite:
primes.append(i)
return primes
# ----------------------------------------------------------
# ② 奇数筛（内存减半）
# ----------------------------------------------------------
@staticmethod
def sieve_odd(n: int) -> List[int]:
"""
仅筛奇数：索引 k 对应数值 2k+1
内存占用约为经典版的 1/2
"""
if n < 2:
return []
if n == 2:
return [2]
half = (n + 1) // 2 # 奇数个数
composite = bytearray(half) # composite[k] → 2k+1 是否为合数
primes = [2]
for k in range(1, math.isqrt(n) // 2 + 1):
if not composite[k]:
p = 2 * k + 1
primes.append(p)
# 起始：p*p → 索引 = (p*p-1)//2
start = (p * p - 1) // 2
# 步长 p 个"奇数位" = 数值步长 2p
composite[start : half : p] = b"\x01" * ((half - start - 1) // p + 1)
# 剩余未标记的奇数
for k in range(math.isqrt(n) // 2 + 1, half):
if not composite[k]:
primes.append(2 * k + 1)
return primes
# ----------------------------------------------------------
# ③ 线性筛（Euler 筛，O(n) 最优）
# ----------------------------------------------------------
@staticmethod
def sieve_linear(n: int) -> List[int]:
"""
Euler 线性筛
每个合数仅被其最小质因子标记一次
时间复杂度严格 O(n)
"""
if n < 2:
return []
if n == 2:
return [2]
composite = bytearray(n + 1)
primes: List[int] = []
for i in range(2, n + 1):
if not composite:
primes.append(i)
for p in primes:
if i * p > n:
break
composite[i * p] = 1
if i % p == 0: # 保证每个合数只被最小质因子标记
break
return primes
# ----------------------------------------------------------
# ④ 生成器 / 迭代器（惰性逐个产出）
# ----------------------------------------------------------
@staticmethod
def iter_primes(n: int) -> Iterator[int]:
"""惰性生成 [2, n] 内素数，几乎不占额外内存"""
if n < 2:
return
yield 2
if n < 3:
return
half = (n + 1) // 2
composite = bytearray(half)
for k in range(1, math.isqrt(n) // 2 + 1):
if not composite[k]:
yield 2 * k + 1
start = (2 * k + 1) ** 2
step = 2 * (2 * k + 1)
for m in range(start, n + 1, step):
idx = m // 2
if idx < half:
composite[idx] = 1
for k in range(math.isqrt(n) // 2 + 1, half):
if not composite[k]:
yield 2 * k + 1
# ----------------------------------------------------------
# ⑤ 分段筛（Segmented Sieve）
# ----------------------------------------------------------
@staticmethod
def segmented_sieve(low: int, high: int) -> List[int]:
"""
分段筛：生成 [low, high] 内所有素数
空间复杂度 O(√high + 区间宽度)，适合超大范围
"""
if high < 2:
return []
low = max(low, 2)
if low > high:
return []
sqrt_h = math.isqrt(high) + 1
base_primes = PrimeGenerator.sieve(sqrt_h)
result: List[int] = []
segment_size = max(32768, sqrt_h) # 缓存友好块大小
while low <= high:
upper = min(low + segment_size - 1, high)
size = upper - low + 1
block = bytearray(size) # 0=素数候选
for p in base_primes:
# 找到 >= low 的最小 p 倍数
start = ((low + p - 1) // p) * p
for m in range(start, upper + 1, p):
block[m - low] = 1
for i in range(size):
if not block:
val = low + i
if val >= 2:
result.append(val)
low = upper + 1
return result
# ----------------------------------------------------------
# ⑥ 米勒-拉宾确定性素性测试
# ----------------------------------------------------------
@staticmethod
def is_prime(n: int) -> bool:
"""
确定性米勒-拉宾素性测试
- n < 2^64 ：固定基底，100% 正确
- n ≥ 2^64 ：40 轮随机，错误概率 < 4^(-40)
"""
if n < 2:
return False
if n in SMALL_PRIME_SET:
return True
if n % 2 == 0 or n % 3 == 0 or n % 5 == 0:
return False
# 快速试除小素数（排除大量合数）
for p in SMALL_PRIMES:
if p * p > n:
break
if n % p == 0:
return False
# 分解 n-1 = d * 2^s
d = n - 1
s = 0
while d & 1 == 0:
d >>= 1
s += 1
# ---- 确定性基底 ----
if n < 2_152_302_898_747:
bases = [2, 3, 5, 7, 11]
elif n < 3_474_749_660_383:
bases = [2, 3, 5, 7, 11, 13]
elif n < 341_550_071_728_321:
bases = [2, 3, 5, 7, 11, 13, 17]
elif n < 2**64:
bases = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
else:
# 超大数：随机 40 轮
for _ in range(40):
a = random.randrange(2, n - 1)
if not PrimeGenerator._miller_witness(a, d, s, n):
return False
return True
for a in bases:
if a % n == 0:
continue
if PrimeGenerator._miller_witness(a, d, s, n):
return False
return True
@staticmethod
def _miller_witness(a: int, d: int, s: int, n: int) -> bool:
"""米勒-拉宾单次 witness 检测，返回 True 表示 n 一定是合数"""
x = pow(a, d, n)
if x == 1 or x == n - 1:
return False
for _ in range(s - 1):
x = (x * x) % n
if x == n - 1:
return False
return True
# ----------------------------------------------------------
# ⑦ 大素数生成（密码学级）
# ----------------------------------------------------------
@staticmethod
def generate_prime(bits: int, safe: bool = False, timeout: float = 30.0) -> int:
"""
生成 bits 位素数
safe=True → 生成安全素数 p = 2q+1（q 也是素数）
timeout → 超时（秒），防止死循环
"""
if bits < 2:
raise ValueError("bits 必须 ≥ 2")
start = time.monotonic()
while True:
if time.monotonic() - start > timeout:
raise TimeoutError(f"{timeout}s 内未生成 {bits} 位素数")
# 最高位 + 最低位（奇数）置 1
n = random.getrandbits(bits)
n |= (1 << (bits - 1)) | 1
if PrimeGenerator.is_prime(n):
if not safe:
return n
# 安全素数：q = (n-1)/2 也必须是素数
q = (n - 1) // 2
if PrimeGenerator.is_prime(q):
return n
@staticmethod
def generate_prime_pair(bits: int) -> Tuple[int, int]:
"""生成 RSA 用的素数对 (p, q)"""
p = PrimeGenerator.generate_prime(bits)
q = PrimeGenerator.generate_prime(bits)
while q == p:
q = PrimeGenerator.generate_prime(bits)
return p, q
# ----------------------------------------------------------
# ⑧ 第 N 个素数（带近似定位）
# ----------------------------------------------------------
@staticmethod
def nth_prime(n: int) -> int:
"""
返回第 n 个素数（从 1 开始）
利用素数定理近似上界，再筛出
"""
if n < 1:
raise ValueError("n 必须 ≥ 1")
if n <= len(SMALL_PRIMES):
return SMALL_PRIMES[n - 1]
# 素数定理近似上界：n ln n + n ln ln n
upper = int(n * (math.log(n) + math.log(math.log(n))) + 50)
upper = max(upper, 100)
while True:
primes = PrimeGenerator.sieve(upper)
if len(primes) >= n:
return primes[n - 1]
upper *= 2 # 不够就翻倍
# ----------------------------------------------------------
# ⑨ 质因数分解（Pollard's Rho）
# ----------------------------------------------------------
@staticmethod
def factorize(n: int) -> List[int]:
"""
Pollard's Rho + 米勒拉宾
返回素因子列表（含重复）
"""
if n < 2:
return []
if PrimeGenerator.is_prime(n):
return [n]
factors: List[int] = []
def _rho(x: int) -> int:
"""Pollard Rho 找因子"""
if x % 2 == 0:
return 2
if PrimeGenerator.is_prime(x):
return x
while True:
c = random.randrange(1, x)
f = lambda a: (pow(a, 2, x) + c) % x
a, b, d = 2, 2, 1
while d == 1:
a = f(a)
b = f(f(b))
d = math.gcd(abs(a - b), x)
if d != x:
return d
stack = [n]
while stack:
v = stack.pop()
if v == 1:
continue
if PrimeGenerator.is_prime(v):
factors.append(v)
continue
d = _rho(v)
stack.append(d)
stack.append(v // d)
return sorted(factors)
@staticmethod
def factorize_unique(n: int) -> List[Tuple[int, int]]:
"""返回 [(素因子, 指数), ...] 格式"""
from collections import Counter
f = PrimeGenerator.factorize(n)
return sorted(Counter(f).items())
# ----------------------------------------------------------
# ⑩ 特殊素数检测
# ----------------------------------------------------------
@staticmethod
def is_twin(p: int) -> bool:
"""判断 p 是否属于孪生素数对（p 和 p+2 或 p-2 同为素数）"""
return PrimeGenerator.is_prime(p) and (
PrimeGenerator.is_prime(p + 2) or PrimeGenerator.is_prime(p - 2)
)
@staticmethod
def is_mersenne(p: int) -> bool:
"""判断 p 是否为梅森素数（2^p - 1 是素数）"""
if not PrimeGenerator.is_prime(p):
return False
m = (1 << p) - 1 # 2^p - 1
return PrimeGenerator.is_prime(m)
@staticmethod
def pi(x: int) -> int:
"""素数计数函数 π(x)：返回 ≤ x 的素数个数"""
if x < 2:
return 0
return len(PrimeGenerator.sieve(x))
# ----------------------------------------------------------
# ⑪ 性能基准测试
# ----------------------------------------------------------
@staticmethod
def benchmark(n: int = 10_000_000) -> str:
"""对比三种筛法性能"""
import time as t
methods = [
("经典埃氏筛", PrimeGenerator.sieve),
("奇数筛", PrimeGenerator.sieve_odd),
("线性筛", PrimeGenerator.sieve_linear),
]
lines = [f"{'='*55}", f" 素数筛性能基准 (n={n:,})", f"{'='*55}"]
for name, func in methods:
t0 = t.perf_counter()
primes = func(n)
elapsed = (t.perf_counter() - t0) * 1000
lines.append(f" {name:12s} → {len(primes):>7,} 个素数 | {elapsed:>8.1f} ms")
lines.append(f"{'='*55}")
return "\n".join(lines)
# ============================================================
# 便捷函数（模块级，方便直接调用）
# ============================================================
def primes_upto(n: int) -> List[int]:
return PrimeGenerator.sieve(n)
def is_prime(n: int) -> bool:
return PrimeGenerator.is_prime(n)
def nth_prime(n: int) -> int:
return PrimeGenerator.nth_prime(n)
def factorize(n: int) -> List[int]:
return PrimeGenerator.factorize(n)
# ============================================================
# 演示 & 自测
# ============================================================
if __name__ == "__main__":
pg = PrimeGenerator()
print("╔══════════════════════════════════════╗")
print("║ PrimeGenerator v3.0 Demo ║")
print("╚══════════════════════════════════════╝\n")
# 1. 基础筛
print("▶ 100 以内素数:")
print(f" {pg.sieve(100)}\n")
# 2. 迭代器
print("▶ 50 以内素数（迭代器）:")
print(f" {list(pg.iter_primes(50))}\n")
# 3. 分段筛
print("▶ [1000, 1050] 区间素数:")
print(f" {pg.segmented_sieve(1000, 1050)}\n")
# 4. 第 N 个素数
print("▶ 第 100 / 1000 / 10000 个素数:")
for n in [100, 1000, 10000]:
p = pg.nth_prime(n)
print(f" π({n}) = {p}")
print()
# 5. 素性测试
print("▶ 米勒-拉宾素性测试:")
test_nums = [982451653, 9999999967, 1234567891, 999999999989]
for v in test_nums:
print(f" is_prime({v:,}) = {pg.is_prime(v)}")
print()
# 6. 大素数生成
print("▶ 生成 128 位素数:")
big = pg.generate_prime(128)
print(f" {big} (bit_length={big.bit_length()})")
print(f" 素性验证: {pg.is_prime(big)}\n")
# 7. 安全素数
print("▶ 生成 64 位安全素数 (p=2q+1):")
safe_p = pg.generate_prime(64, safe=True)
q = (safe_p - 1) // 2
print(f" p = {safe_p} (prime: {pg.is_prime(safe_p)})")
print(f" q = {q} (prime: {pg.is_prime(q)})\n")
# 8. 质因数分解
print("▶ Pollard-Rho 质因数分解:")
demo_nums = [123456789012345, 987654321098765, 2**10 * 3**5 * 7**3]
for v in demo_nums:
fac = pg.factorize_unique(v)
fac_str = " × ".join(f"{p}^{e}" if e > 1 else str(p) for p, e in fac)
print(f" {v:,} = {fac_str}")
print()
# 9. 孪生素数
print("▶ 孪生素数对 (前 10 对):")
twins = []
for p in pg.iter_primes(1000):
if pg.is_twin(p) and p + 2 <= 1000 and pg.is_prime(p + 2):
twins.append((p, p + 2))
if len(twins) >= 10:
break
for a, b in twins:
print(f" ({a}, {b})")
print()
# 10. 梅森素数
print("▶ 梅森素数 (p 使 2^p-1 为素数, p≤31):")
for p in [2, 3, 5, 7, 13, 17, 19, 31]:
if pg.is_mersenne(p):
m = (1 << p) - 1
print(f" M_{p} = {m:,}")
print()
# 11. 性能基准
print(pg.benchmark(5_000_000))
print()
# 12. 素数定理验证
print("▶ 素数定理验证 π(x) ≈ x/ln(x):")
for x in [1000, 10000, 100000, 1000000]:
actual = pg.pi(x)
approx = x / math.log(x)
err = abs(actual - approx) / actual * 100
print(f" π({x:,}) = {actual:,} | x/ln(x) ≈ {approx:,.1f} | 误差 {err:.1f}%")
print()
print("✅ 全部演示完成")
|
|
IP卡
狗仔卡
发表于 2026-8-19 11:09
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
显身卡