173 lines
4.9 KiB
Python
173 lines
4.9 KiB
Python
"""
|
||
The prime 41, can be written as the sum of six consecutive primes:
|
||
41 = 2 + 3 + 5 + 7 + 11 + 13
|
||
This is the longest sum of consecutive primes that adds to a prime below one-hundred.
|
||
|
||
The longest sum of consecutive primes below one-thousand that adds to a prime, contains 21 terms,
|
||
and is equal to 953.
|
||
|
||
Which prime, below one-million, can be written as the sum of the most consecutive primes?
|
||
"""
|
||
|
||
import time
|
||
from math import isqrt, log, sqrt
|
||
|
||
from bitarray import bitarray
|
||
|
||
_prime_cache = set()
|
||
|
||
|
||
def timer(func):
|
||
def wrapper(*args, **kwargs):
|
||
start_time = time.perf_counter()
|
||
result = func(*args, **kwargs)
|
||
end_time = time.perf_counter()
|
||
elapsed_time = end_time - start_time
|
||
print(f"{func.__name__} time: {elapsed_time:.6f} seconds")
|
||
return result
|
||
|
||
return wrapper
|
||
|
||
|
||
def n_primes(n: int = 4000) -> list[int]:
|
||
if n < 0:
|
||
raise ValueError("n must be a positive integer")
|
||
if n == 0:
|
||
return []
|
||
if n == 1:
|
||
return [2]
|
||
|
||
# 使用Dusart方法估计上限
|
||
limit = int(n * (log(n) + log(log(n)) - 0.5))
|
||
|
||
# 初始化全1（假设都是素数），0和1置为0
|
||
is_prime = bitarray(limit + 1)
|
||
is_prime.setall(True)
|
||
is_prime[:2] = False # 0和1不是素数
|
||
|
||
# 只需筛到 sqrt(n)
|
||
imax = isqrt(limit) + 1
|
||
|
||
for i in range(2, imax):
|
||
if is_prime[i]:
|
||
# 步长i，从i*i开始（小于i*i的已被更小的素数筛过）
|
||
is_prime[i * i : limit + 1 : i] = False
|
||
|
||
# 提取结果
|
||
result = []
|
||
for i, val in enumerate(is_prime):
|
||
if val:
|
||
result.append(i)
|
||
if len(result) >= n:
|
||
break
|
||
|
||
return result
|
||
|
||
|
||
def is_prime(num: int) -> bool:
|
||
if num in _prime_cache:
|
||
return True
|
||
if num < 2:
|
||
return False
|
||
if num in (2, 3):
|
||
_prime_cache.add(num)
|
||
return True
|
||
if num % 2 == 0 or num % 3 == 0:
|
||
return False
|
||
|
||
# 检查6k±1形式的因子
|
||
i = 5
|
||
w = 2
|
||
while i * i <= num:
|
||
if num % i == 0:
|
||
return False
|
||
i += w
|
||
w = 6 - w # 在2和4之间切换，实现6k±1的检查
|
||
|
||
_prime_cache.add(num)
|
||
return True
|
||
|
||
|
||
def get_bound_dynanic(N, iterations=3):
|
||
# 初始猜测（基于简化公式）
|
||
k = sqrt(2 * N / log(N)) if N > 100 else 10
|
||
|
||
# 牛顿法迭代（3次足够）
|
||
for _ in range(iterations):
|
||
ln_k = log(k)
|
||
# f(k) = k^2 * ln_k - 2N, f'(k) = k * (2*ln_k + 1)
|
||
k = (k**2 * (1 + ln_k) + 2 * N) / (k * (2 * ln_k + 1))
|
||
|
||
return int(k) # 返回安全上界（如562）
|
||
|
||
|
||
def get_bound(limit: int = 10**6) -> int:
|
||
"""
|
||
估算素数个数上限。
|
||
"""
|
||
return int(isqrt(int(2 * limit / log(limit))) * 1.55)
|
||
|
||
|
||
@timer
|
||
def max_primes_sum_best(limit: int = 10**6) -> tuple[int, int] | None:
|
||
# === 1. 确定搜索范围 ===
|
||
bound = get_bound_dynanic(limit)
|
||
print(f"bound: {bound}")
|
||
primes = n_primes(bound)
|
||
|
||
# === 2. 构建前缀和数组 ===
|
||
# prefix[i] 表示前 i 个素数的和（primes[0] 到 primes[i-1]）
|
||
prefix = [0]
|
||
current_sum = 0
|
||
max_possible_length = 0
|
||
|
||
# 计算从2开始连续累加，不超过 limit 的最大素数个数
|
||
for p in primes:
|
||
current_sum += p
|
||
if current_sum >= limit:
|
||
break
|
||
prefix.append(current_sum)
|
||
max_possible_length += 1
|
||
|
||
# === 3. 搜索最长序列（倒序遍历 + 剪枝）===
|
||
best_length = 0
|
||
best_prime = 0
|
||
|
||
# 外层：右端点从大到小遍历（优先尝试更长序列）
|
||
for right in range(max_possible_length, 0, -1):
|
||
# 关键剪枝1：如果右端点 <= 当前最优长度，不可能找到更长序列
|
||
# 因为即使从左端点0开始，长度也只有 right，而 right <= best_length
|
||
if right <= best_length:
|
||
break
|
||
|
||
# 关键剪枝2：左端点只需考虑到 (right - best_length - 1)
|
||
# 因为我们需要的长度是 (right - left)，必须满足 > best_length
|
||
# 即 left < right - best_length
|
||
max_left = right - best_length
|
||
|
||
for left in range(max_left):
|
||
# 计算连续素数之和：primes[left] 到 primes[right-1]
|
||
consecutive_sum = prefix[right] - prefix[left]
|
||
|
||
if consecutive_sum >= limit:
|
||
continue # 跳过，但继续尝试更大的 left（sum 会变小）
|
||
|
||
if is_prime(consecutive_sum):
|
||
length = right - left
|
||
|
||
if length > best_length:
|
||
best_length = length
|
||
best_prime = consecutive_sum
|
||
|
||
return best_prime, best_length
|
||
|
||
|
||
def main():
|
||
limit = int(input("limit:"))
|
||
max_sum = max_primes_sum_best(limit)
|
||
print(f"max primt: {max_sum}")
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|