133. 循环单位非因子(Repunit Nonfactors)

完全由数字 1 组成的数称为循环单位（repunit）。定义 \(R(k)\) 为长度为 \(k\) 的循环单位；例如，\(R(6)=111111\)。

考虑形如 \(R(10^n)\) 的循环单位。

虽然 \(R(10)\)、\(R(100)\) 或 \(R(1000)\) 都不能被 \(17\) 整除，但 \(R(10000)\) 能被 \(17\) 整除。而对于 \(19\)，无论 \(n\) 取何值，\(R(10^n)\) 都不能被 \(19\) 整除。事实上，值得注意的是，\(11\)、\(17\)、\(41\) 和 \(73\) 是一百以内仅有的四个能成为 \(R(10^n)\) 因子的质数。

求所有小于十万且永远不能成为 \(R(10^n)\) 因子的质数之和。

分析：本题的核心是判断某个质数 \(p\) 能否整除 \(R(10^n)\)。由 \(R(k)=\frac{10^k-1}{9}\) 可得，对于 \(p\neq 3\)，\(p\mid R(k)\) 当且仅当 \(10^k\equiv 1\pmod p\)。引入 \(A(p)=\operatorname{ord}_p(10)\) 为 \(10\) 模 \(p\) 的乘法阶，问题转化为：\(A(p)\mid 10^n=2^n\cdot5^n\) 是否有解。这意味着 \(A(p)\) 的质因子必须全部属于 \(\{2,5\}\)。利用 \(A(p)\mid(p-1)\)，只需提取 \(p-1\) 中因子 \(2\) 和 \(5\) 的幂次构成 \(M\)，再检查 \(10^M\equiv 1\pmod p\) 即可。对 \(p=2,3,5\) 需要单独处理。

一、数学背景

循环单位 \(R(k)\) 定义为 \(R(k)=\frac{10^k-1}{9}\)。对于质数 \(p\)，欲知是否存在 \(n\) 使得 \(p\mid R(10^n)\)，分两种情况讨论。

\(p\neq 3\) 的情形：此时 \(\gcd(9,p)=1\)，\(9\) 在模 \(p\) 下可逆，因此

设 \(A(p)=\operatorname{ord}_p(10)\) 为 \(10\) 模 \(p\) 的乘法阶，即满足 \(10^d\equiv1\pmod p\) 的最小正整数 \(d\)。由费马小定理，\(A(p)\mid(p-1)\)。于是

因此 \(A(p)\mid 10^n\) 成立当且仅当 \(A(p)\) 的所有质因子均来自 \(\{2,5\}\)。

\(p=3\) 的情形：此时 \(9\) 与 \(3\) 不互质，等价的 \(10^k\equiv1\pmod3\) 不成立。直接计算：\(10\equiv1\pmod3\)，\(R(k)=\sum_{i=0}^{k-1}10^i\equiv\sum_{i=0}^{k-1}1=k\pmod3\)。当 \(k=10^n\) 时，\(k\equiv1\pmod3\)，故 \(R(10^n)\equiv1\not\equiv0\pmod3\)。因此 \(p=3\) 永远不整除 \(R(10^n)\)。

\(p=2,5\) 的情形：\(\gcd(10,p)\neq1\)，\(10^k\not\equiv1\pmod p\) 恒成立。由 \(R(k)\) 的奇偶性和末位数字亦可知 \(p=2,5\) 永远不整除任何 repunit。

判定方法：对 \(p\neq3\)，设 \(v_2,v_5\) 分别为 \(p-1\) 中因子 \(2\) 和 \(5\) 的幂次，令 \(M=2^{v_2}\cdot5^{v_5}\)。则

若 \(10^M\equiv1\pmod p\)，则 \(p\) 可整除某个 \(R(10^n)\)；否则 \(p\) 永远不能。

二、算法设计

对于小于 \(10^5\) 的所有质数 \(p\)：

• 若 \(p=3\)，直接判定为 nonfactor（永远不整除）。
• 否则提取 \(p-1\) 中因子 \(2\) 和 \(5\) 的幂次，构造 \(M\)。
• 使用快速模幂 pow(10, M, p) 计算 \(10^M\bmod p\)。
• 若结果等于 \(1\)，则 \(p\) 是 potential factor，跳过；否则累加 \(p\)。

筛法与判定均可在极短时间内完成。

三、复杂度分析

• 埃拉托色尼筛法生成 \(10^5\) 以内的质数，时间复杂度 \(O(N\log\log N)\)，空间 \(O(N)\)。
• 对每个质数 \(p\)，提取因子和一次模幂均为 \(O(\log p)\)。
• 质数个数 \(\pi(10^5)\approx9592\)，总时间远小于 \(60\) 秒（实测约 \(70\) 毫秒）。

四、代码实现与说明

"""
Project Euler Problem 133: Repunit Nonfactors

Find the sum of all primes below 100,000 that will never be a factor of R(10^n),
where R(k) = (10^k - 1) / 9 is the repunit of length k.
"""

from math import isqrt


def sieve_primes(limit: int) -> list[int]:
    """Return all primes below `limit` using the Sieve of Eratosthenes."""
    is_prime = [True] * limit
    is_prime[0] = is_prime[1] = False
    for i in range(2, isqrt(limit) + 1):
        if is_prime[i]:
            step = i
            start = i * i
            is_prime[start:limit:step] = [False] * ((limit - start - 1) // step + 1)
    return [i for i, p in enumerate(is_prime) if p]


def can_be_factor(p: int) -> bool:
    """
    Return True if prime p can divide R(10^n) for some integer n >= 0.

    For p == 3, the formula (10^k-1)/9 is not equivalent to 10^k \equiv 1 (mod 3)
    because 9 is not invertible mod 3. R(10^n) \equiv 1 (mod 3) always.
    """
    if p == 3:
        return False

    t = p - 1
    m = 1
    while t % 2 == 0:
        t //= 2
        m *= 2
    while t % 5 == 0:
        t //= 5
        m *= 5

    return pow(10, m, p) == 1


def solve(limit: int = 100000) -> int:
    """Return the sum of all primes below `limit` that never divide R(10^n)."""
    primes = sieve_primes(limit)
    total = 0
    for p in primes:
        if not can_be_factor(p):
            total += p
    return total


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

代码说明（按执行顺序）：

导入与常量

• 从 math 导入 isqrt（整数平方根），用于筛法上限判定。

sieve_primes

• 初始化布尔数组 is_prime，下标 \(0\) 和 \(1\) 置为 False，其余为 True。
• 外层循环 \(i\) 从 \(2\) 到 \(\lfloor\sqrt{N}\rfloor\)：若 \(i\) 为质数，则用切片赋值将其所有倍数标记为 False。
• 切片步进 step = i，起始 start = i * i（更小的倍数已被更小的质因子筛去）。
• 切片长度由 (limit - start - 1) // step + 1 计算，确保覆盖到 limit-1。
• 返回 True 对应的下标列表。

can_be_factor

• 参数 \(p\) 是质数，返回该质数能否整除某个 \(R(10^n)\)。
• 若 \(p=3\)：直接返回 False。因为 \(R(10^n)\equiv1\pmod3\) 恒成立，\(3\) 永不整除。
• 从 \(t = p-1\) 开始，反复除以 \(2\) 直至不能整除，同时将 \(m\) 乘以 \(2\)。
• 同理处理因子 \(5\)。
• 最终 \(m = 2^{v_2(p-1)}\cdot5^{v_5(p-1)}\)。
• 调用 pow(10, m, p) 计算 \(10^m\bmod p\)，若等于 \(1\) 则 \(p\) 可整除某个 \(R(10^n)\)。

solve

• 获取小于 limit 的全部质数。
• 累加所有 not can_be_factor(p) 的 \(p\)，即所有 nonfactor。
• 返回累加和。

main 区块

• 以 \(10^5\) 为上限调用 solve 并打印结果。