146. 质数模式研究(Investigating a Prime Pattern)

使得 \(n^2+1\)、\(n^2+3\)、\(n^2+7\)、\(n^2+9\)、\(n^2+13\) 和 \(n^2+27\) 为连续质数的最小正整数 \(n\) 是 \(10\)。一百万以下所有这样的整数 \(n\) 之和为 \(1242490\)。

求一亿五千万以下所有这样的整数 \(n\) 之和。

分析：\(n\) 必须为偶数（否则 \(n^2+1\) 为 \(>2\) 的偶数），且由模 \(3,5,7,\dots\) 等小质数约束，大量 \(n\) 被筛除。用模乘积 \(P=9699690\)（\(2\cdot3\cdot5\cdot7\cdot11\cdot13\cdot17\cdot19\)）预计算有效剩余类，仅约 24024 个。对每个候选 \(n\)，用 Miller-Rabin 检验 6 个目标数与 2 个中间数。

一、数学背景

要使 \(n^2+1, n^2+3, n^2+7, n^2+9, n^2+13, n^2+27\) 为连续质数，需要： 1. 这 6 个数都是质数。 2. 它们之间的所有奇数（\(n^2+5,11,15,17,19,21,23,25\)）都是合数。

\(n\) 的模约束：

• \(n\) 必须为偶数。若 \(n\) 为奇数，\(n^2\) 为奇数，\(n^2+1\) 为偶数 \(>2\)，必为合数。
• \(n\equiv\pm1\pmod3\)。若 \(n\equiv0\pmod3\)，则 \(n^2+3, n^2+9, n^2+27\) 均被 \(3\) 整除。
• \(n\equiv0\pmod5\)。若 \(n\not\equiv0\pmod5\)，则 \(n^2\equiv1\) 或 \(4\pmod5\)，分别导致 \(n^2+9\equiv0\) 或 \(n^2+1\equiv0\pmod5\)。
• 综合得 \(n\equiv10,20\pmod{30}\)。

更进一步，对质数 \(p=7,11,13,17,19,\dots\) 筛选 \(n\bmod p\) 的合法取值，使 \(n^2+k\) 对目标偏移量 \(k\in\{1,3,7,9,13,27\}\) 均不被 \(p\) 整除。

中间数的自动合数性：

在 \(n\equiv\pm1\pmod3\) 且 \(n\equiv0\pmod5\) 的条件下： - \(n^2+5, n^2+11, n^2+17, n^2+23\) 均被 \(3\) 整除（\(n^2\equiv1\pmod3\)）。 - \(n^2+15, n^2+25\) 均被 \(5\) 整除（\(n^2\equiv0\pmod5\)）。 - 仅 \(n^2+19\) 和 \(n^2+21\) 不自动为合数，需显式检测。

二、算法设计

筛法生成有效剩余类

选取小质数集合 \(S=\{2,3,5,7,11,13,17,19\}\)，乘积 \(P=9699690\)。用 bytearray 构建长度为 \(P\) 的布尔筛，初始全为真。对每个 \(p\in S\)： - 枚举 \(r=0,\dots,p-1\)，若 \(r^2+k\equiv0\pmod p\) 对任意目标 \(k\) 成立，则 \(r\) 为禁止剩余类。 - 在筛中将所有满足 \(n\equiv r\pmod p\) 的位置标记为假。

最终筛中为真的位置即为有效剩余类。

候选迭代与质数检验

对每个有效剩余类 \(r\)，以 \(P\) 为步长枚举 \(n=r, r+P, r+2P, \dots < 1.5\times10^8\)。对每个候选 \(n\)： - 用 Miller-Rabin 检测 6 个目标数是否均为质数。若任一不是，跳过。 - 用 Miller-Rabin 检测 \(n^2+19\) 和 \(n^2+21\) 是否为合数。若任一是质数，跳过。 - 通过全部检验则累加 \(n\)。

Miller-Rabin 使用确定性 7 基组 \(\{2,325,9375,28178,450775,9780504,1795265022\}\)，覆盖全部 64 位整数。

三、复杂度分析

• 筛法耗时 \(O(P\cdot|S|)\approx 8\times10^7\) 字节操作，约 0.5 秒。
• 有效剩余类 24024 个，每个产生约 15.5 个候选，总计约 37 万次检验。
• 每次检验最多 8 次 Miller-Rabin（6 目标 + 2 中间数），每个 MR 为 \(O(\log^3 n)\) 模幂。
• 总耗时约 7-8 秒，远低于 60 秒限制。

四、代码实现与说明

"""
Project Euler Problem 146: Investigating a Prime Pattern

Find the sum of all integers n < 150,000,000 such that:
n^2+1, n^2+3, n^2+7, n^2+9, n^2+13, n^2+27 are consecutive primes.
"""

from math import isqrt

TARGET_OFFSETS = (1, 3, 7, 9, 13, 27)
EXTRA_CHECK_OFFSETS = (19, 21)

FILTER_PRIMES = (2, 3, 5, 7, 11, 13, 17, 19)

MR_BASES = (2, 325, 9375, 28178, 450775, 9780504, 1795265022)


def is_prime_mr(n: int) -> bool:
    """Miller-Rabin deterministic primality test for 64-bit integers."""
    if n < 2:
        return False
    if n in (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37):
        return True
    if n % 2 == 0:
        return False

    d = n - 1
    s = 0
    while d % 2 == 0:
        d //= 2
        s += 1

    for a in MR_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 = (x * x) % n
            if x == n - 1:
                break
        else:
            return False
    return True


def solve() -> int:
    N = 150_000_000

    P = 1
    for p in FILTER_PRIMES:
        P *= p

    sieve = bytearray(b'\x01') * P
    for p in FILTER_PRIMES:
        forbidden = bytearray(p)
        for r in range(p):
            r2 = (r * r) % p
            for off in TARGET_OFFSETS:
                if (r2 + off) % p == 0:
                    forbidden[r] = 1
                    break
        for f in range(p):
            if forbidden[f]:
                sieve[f:P:p] = b'\x00' * len(sieve[f:P:p])

    total = 0
    for r in range(P):
        if not sieve[r]:
            continue
        n = r
        while n < N:
            n2 = n * n

            all_prime = True
            for off in TARGET_OFFSETS:
                if not is_prime_mr(n2 + off):
                    all_prime = False
                    break
            if not all_prime:
                n += P
                continue

            bad = False
            for off in EXTRA_CHECK_OFFSETS:
                if is_prime_mr(n2 + off):
                    bad = True
                    break
            if bad:
                n += P
                continue

            total += n
            n += P

    return total


if __name__ == "__main__":
    result = solve()
    print(result)

代码说明：

常量区

• TARGET_OFFSETS：需要为质数的 6 个偏移量。
• EXTRA_CHECK_OFFSETS：需验证为合数的 2 个中间偏移量（其余 6 个中间数由模 3 或 5 自动确定为合数）。
• FILTER_PRIMES：用于筛选的小质数集合，乘积 \(P=9699690\)。
• MR_BASES：Miller-Rabin 确定性基组，可覆盖全部 \(<2^{64}\) 的整数。

is_prime_mr

• 对 64 位整数执行确定性 Miller-Rabin 质数判定。
• \(n<2\) 返回假；小质数直接返回真；偶数返回假。
• 分解 \(n-1 = d \cdot 2^s\)。
• 对每个基 \(a\)：若 \(a\) 是 \(n\) 的倍数则跳过；计算 \(x = a^d \bmod n\)；若 \(x=1\) 或 \(x=n-1\) 则本轮通过；否则重复平方 \(s-1\) 次，期望某次得到 \(n-1\)。若始终未出现，返回假。
• 全部基通过后返回真。

solve

• 计算 \(P = \prod_{p\in S} p\)。
• 构建 bytearray 筛子，初始化所有位置为 \(1\)（有效）。

筛法

• 对外层 \(p \in S\)：创建 forbidden 列表记录 \(p\) 的禁止剩余类。
• 枚举 \(r=0,\dots,p-1\)，计算 \(r^2 \bmod p\)，检查与 6 个目标偏移量之和是否被 \(p\) 整除。
• 对每个禁止剩余类 \(f\)，使用切片赋值 sieve[f:P:p] = b'\x00' * count 高效标记全部 \(n\equiv f\pmod p\) 的位置为无效。

候选迭代

• 遍历 \(r=0,\dots,P-1\)，跳过筛子标记为无效的 \(r\)。
• 内层 while n < N：以 \(P\) 为步长枚举所有 \(n\equiv r\pmod P\)。
• 计算 \(n^2\)，对 6 个目标偏移量依次执行 MR 判定。若某数非质数，all_prime 置为假并跳出。
• 若 6 个目标数均为质数，再对 \(n^2+19\) 和 \(n^2+21\) 执行 MR 判定。若任一是质数，bad 置为真。
• 全部通过后累加 \(n\) 到 total。
• 返回 total。

主程序

• 调用 solve() 并打印结果。