233 lines
6.9 KiB
Python
233 lines
6.9 KiB
Python
"""
|
||
The first two consecutive numbers to have two distinct prime factors are:
|
||
|
||
14 = 2 * 7
|
||
15 = 3 * 5
|
||
|
||
The first three consecutive numbers to have three distinct prime factors are:
|
||
|
||
644 = 2^2 * 7 * 23
|
||
645 = 3 * 5 * 43
|
||
646 = 2 * 17 * 19
|
||
|
||
Find the first four consecutive integers to have four distinct prime factors each.
|
||
What is the first of these numbers?
|
||
"""
|
||
|
||
import time
|
||
from typing import Dict, List
|
||
|
||
|
||
def timer(func):
|
||
def wrapper(*args, **kwargs):
|
||
start_time = time.time()
|
||
result = func(*args, **kwargs)
|
||
end_time = time.time()
|
||
print(f"{func.__name__} took {end_time - start_time:.6f} seconds to run.")
|
||
return result
|
||
|
||
return wrapper
|
||
|
||
|
||
def is_probable_prime(n: int) -> bool:
|
||
"""
|
||
Miller-Rabin素性测试
|
||
对于n < 3,317,044,064,679,887,385,961,981 (3e24)，使用确定性基底集
|
||
"""
|
||
if n < 2:
|
||
return False
|
||
# 小素数快速检查
|
||
small_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
|
||
for p in small_primes:
|
||
if n % p == 0:
|
||
return n == p
|
||
|
||
# 写成 d * 2^s 的形式
|
||
d = n - 1
|
||
s = 0
|
||
while d % 2 == 0:
|
||
d //= 2
|
||
s += 1
|
||
|
||
# 确定性基底集（对64位整数足够）
|
||
# 根据研究，测试 [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
|
||
# 可以覆盖 < 2^64 的所有数
|
||
test_bases = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
|
||
|
||
for a in test_bases:
|
||
if a % n == 0:
|
||
continue
|
||
x = pow(a, d, n)
|
||
if x == 1 or x == n - 1:
|
||
continue
|
||
for _ in range(s - 1):
|
||
x = pow(x, 2, n)
|
||
if x == n - 1:
|
||
break
|
||
else:
|
||
return False
|
||
return True
|
||
|
||
|
||
def factorize(n: int) -> Dict[int, int]:
|
||
"""
|
||
优化版试除法因数分解
|
||
返回: {质因数: 指数, ...}
|
||
|
||
优化策略：
|
||
1. 快速排除小因子（2, 3）
|
||
2. 6k±1优化（跳过所有3的倍数）
|
||
3. 动态更新上界（随着因子被移除，sqrt(n)减小）
|
||
4. 对大余数进行素性预检
|
||
"""
|
||
if n < 2:
|
||
return {}
|
||
|
||
factors = {}
|
||
remaining = n
|
||
|
||
# 处理负数
|
||
if remaining < 0:
|
||
factors[-1] = 1
|
||
remaining = -remaining
|
||
|
||
# 步骤1: 快速处理小素数（2和3）
|
||
for p in [2, 3]:
|
||
if remaining % p == 0:
|
||
count = 0
|
||
while remaining % p == 0:
|
||
remaining //= p
|
||
count += 1
|
||
factors[p] = count
|
||
|
||
# 步骤2: 6k±1优化
|
||
# 所有大于3的素数都形如 6k±1
|
||
# 即依次检查 5, 7, 11, 13, 17, 19...
|
||
# 步长模式: +2 (到6k+1), +4 (到6k+5, 即下一个6(k+1)-1)
|
||
divisor = 5
|
||
step = 2 # 交替使用 2 和 4
|
||
|
||
# 动态计算上界：只需要检查到 sqrt(remaining)
|
||
# 随着remaining减小，上界也减小
|
||
while divisor * divisor <= remaining:
|
||
if remaining % divisor == 0:
|
||
count = 0
|
||
while remaining % divisor == 0:
|
||
remaining //= divisor
|
||
count += 1
|
||
factors[divisor] = count
|
||
# 更新上界（优化关键！）
|
||
|
||
divisor += step
|
||
step = 6 - step # 2 -> 4 -> 2 -> 4...
|
||
|
||
# 步骤3: 如果remaining > 1，说明remaining本身是质数
|
||
# 使用Miller-Rabin确认（对大数避免误判）
|
||
if remaining > 1:
|
||
# 对于小余数直接确认，大余数用素性测试
|
||
if remaining < 1_000_000 or is_probable_prime(remaining):
|
||
factors[remaining] = factors.get(remaining, 0) + 1
|
||
else:
|
||
# 极小概率：remaining是合数但试除法未找到因子
|
||
# 此时remaining必为两个大质数的乘积，且都 > sqrt(original_n)
|
||
# 对于这种情况，可回退到Pollard's Rho（可选）
|
||
factors[remaining] = 1
|
||
|
||
return factors
|
||
|
||
|
||
def factorize_list(n: int) -> List[int]:
|
||
"""返回展开形式的质因数列表，如 12 -> [2, 2, 3]"""
|
||
factors = factorize(n)
|
||
result = []
|
||
for p, exp in sorted(factors.items()):
|
||
result.extend([p] * exp)
|
||
return result
|
||
|
||
|
||
def distinct_prime_factors_sieve(limit: int) -> List[int]:
|
||
"""
|
||
返回一个列表pf，其中pf[i]是i的不同质因数个数，i从0到limit。
|
||
使用筛法计算，复杂度O(limit log log limit)。
|
||
"""
|
||
pf = [0] * (limit + 1)
|
||
for p in range(2, limit + 1):
|
||
if pf[p] == 0: # p是质数
|
||
for multiple in range(p, limit + 1, p):
|
||
pf[multiple] += 1
|
||
return pf
|
||
|
||
|
||
@timer
|
||
def main(limit: int = 4) -> None:
|
||
n = 1155
|
||
keep_ok = False
|
||
res = []
|
||
while True:
|
||
tl = set(factorize_list(n))
|
||
if len(tl) == limit:
|
||
res.append(n)
|
||
keep_ok = True
|
||
if len(res) == limit and keep_ok:
|
||
print(res)
|
||
break
|
||
else:
|
||
res = []
|
||
keep_ok = False
|
||
n += 1
|
||
|
||
|
||
@timer
|
||
def main_key(limit: int = 4) -> None:
|
||
if limit < 2:
|
||
raise ValueError("limit must be at least 2")
|
||
lease = [2, 3]
|
||
# 预计算一个足够大的范围，这里假设答案不会超过10^7，但为了效率，动态扩展
|
||
# 我们使用一个缓存，按需扩展
|
||
_pf_cache = [] # 缓存数组，索引i对应数字i的质因数个数
|
||
_pf_cache_limit = 0
|
||
|
||
def ensure_cache(upto: int) -> None:
|
||
nonlocal _pf_cache, _pf_cache_limit
|
||
if upto <= _pf_cache_limit:
|
||
return
|
||
# 扩展缓存，每次至少扩展到upto，或者按一定步长扩展
|
||
new_limit = max(upto, _pf_cache_limit * 2 if _pf_cache_limit else upto + 10000)
|
||
# 重新计算整个缓存，或者增量更新？为了简单，重新计算整个范围
|
||
# 但增量更新较复杂，我们重新计算到new_limit
|
||
_pf_cache = distinct_prime_factors_sieve(new_limit)
|
||
_pf_cache_limit = new_limit
|
||
|
||
for i in range(2, limit + 1):
|
||
if i in [2, 3]:
|
||
n = lease[i - 2]
|
||
else:
|
||
n = lease[i - 3] + lease[i - 4]
|
||
start = 10 ** (n - 1)
|
||
# 确保缓存至少覆盖 start + 一个估计的窗口大小，比如10000
|
||
ensure_cache(start + 10000)
|
||
keep_ok = False
|
||
res = []
|
||
current = start
|
||
while True:
|
||
# 如果当前数字超出缓存，扩展缓存
|
||
if current > _pf_cache_limit:
|
||
ensure_cache(current + 10000)
|
||
if _pf_cache[current] == i:
|
||
res.append(current)
|
||
keep_ok = True
|
||
if len(res) == i and keep_ok:
|
||
print(f"{i} - {res}")
|
||
if i > 3:
|
||
lease.append(len(str(max(res))))
|
||
break
|
||
else:
|
||
res = []
|
||
keep_ok = False
|
||
current += 1
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|
||
main_key()
|