""" 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()