148. 探索帕斯卡三角(Exploring Pascal's Triangle)

我们可以轻松验证，帕斯卡三角形的前七行中没有任何一项能被 \(7\) 整除：

(帕斯卡三角前七行图示)

然而，如果我们检查前一百行，会发现 \(5050\) 项中仅有 \(2361\) 项不被 \(7\) 整除。

求帕斯卡三角形前 \(10^9\) 行中不被 \(7\) 整除的项数。

分析：由 Lucas 定理，第 \(n\) 行不被质数 \(p\) 整除的项数为 \(f(n)=\prod(d_i+1)\)，其中 \(d_i\) 是 \(n\) 的 \(p\) 进制各位。问题转化为对所有 \(n=0\) 到 \(N-1\) 求和 \(\sum f(n)\)，可通过 digit DP 递归实现。

一、数学背景

Lucas 定理：对于质数 \(p\)，将 \(n\) 和 \(k\) 按 \(p\) 进制展开为 \(n=\sum n_i p^i\)、\(k=\sum k_i p^i\)，则

因此 \(\binom{n}{k} \not\equiv 0 \pmod{p}\) 当且仅当对所有 \(i\) 有 \(k_i \le n_i\)。

对于固定的 \(n\)，每一位 \(n_i\) 有 \(n_i+1\) 种 \(k_i\) 的合法选择（\(0,1,\dots,n_i\)），故第 \(n\) 行不被 \(p\) 整除的项数：

其中 \(n_i\) 是 \(n\) 的 \(p\) 进制各位数码（\(0\le n_i < p\)）。

二、算法设计

设 \(S(N) = \sum_{n=0}^{N-1} f(n)\)，目标是求 \(S(10^9)\)，其中 \(p=7\)。

考虑所有 \(m\) 位 \(p\) 进制数（首位可以为 \(0\)），其 \(f(n)\) 之和为：

对 \(p=7\)，每位的贡献和为 \(\frac{7\cdot8}{2}=28\)，故 \(m\) 位全范围的 \(f(n)\) 和为 \(28^m\)。

对一般 \(N\)，从最高位逐位处理。设当前位值为 \(d\)（位于 \(p^k\) 位），前缀累积乘积为 \(\text{pref}\)：

• 该位选 \(0,1,\dots,d-1\) 时，低位任意取 \(p^{k}\) 种可能，贡献为 \(\text{pref}\cdot\frac{d(d+1)}{2}\cdot 28^k\)。
• 该位选 \(d\) 时，前缀乘积更新为 \(\text{pref}\cdot(d+1)\)，继续处理低位。

基：\(S(0)=0\)。

三、复杂度分析

\(10^9\) 的 7 进制最多 11 位，递归层数 \(O(\log_7 N)\approx 11\)，每层 \(O(1)\) 运算。总时间在微秒级。

四、代码实现与说明

"""
Project Euler Problem 148: Exploring Pascal's Triangle

Find the number of entries in the first 10^9 rows of Pascal's triangle
that are NOT divisible by 7.
"""


def to_base(n: int, base: int) -> list[int]:
    """Return the digits of n in the given base, from most to least significant."""
    if n == 0:
        return [0]
    digits = []
    while n > 0:
        digits.append(n % base)
        n //= base
    return digits[::-1]


def sum_not_divisible(N: int, base: int = 7) -> int:
    """
    Return the sum of f(n) for n = 0 to N-1, where f(n) = prod(d_i + 1)
    and d_i are the base-b digits of n.
    """
    if N <= 0:
        return 0

    digits = to_base(N, base)
    k = len(digits) - 1

    result = 0
    prefix_prod = 1

    for i, d in enumerate(digits):
        pos = k - i
        full_block_sum = (base * (base + 1) // 2) ** pos

        contrib = prefix_prod * (d * (d + 1) // 2) * full_block_sum
        result += contrib

        prefix_prod *= (d + 1)

    return result


def solve() -> int:
    """Return the answer for N = 10^9."""
    N = 10 ** 9
    return sum_not_divisible(N)


if __name__ == "__main__":
    assert sum_not_divisible(100) == 2361, "Sample verification failed!"
    result = solve()
    print(result)

代码说明：

to_base

• 将整数 \(n\) 转换为指定进制，返回从高位到低位的数码列表。
• 处理 \(n=0\) 的特殊情况，返回 [0]。
• 循环提取 \(n \bmod b\) 作为当前位，然后 \(n \gets \lfloor n/b \rfloor\)。
• 最后翻转得到从高到低的顺序。

sum_not_divisible

• 参数 \(N\)：计算 \(n = 0\) 到 \(N-1\) 的 \(f(n)\) 之和。
• 参数 \(\text{base}=7\)。
• 若 \(N\le 0\)，返回 \(0\)。

主循环

• 将 \(N\) 转为 7 进制数码列表，\(k\) 为最高位的幂次。
• result 累加贡献，prefix_prod 跟踪已固定高位对乘积的累积因子。
• 对每一位：
• pos = k - i 是该位对应的 \(7^{\text{pos}}\) 位置。
• full_block_sum = 28^pos 是低位 \(pos\) 位全组合的 \(f\) 之和。
• 该位取 \(0\) 到 \(d-1\) 的贡献：前缀乘积 × 该位平均贡献 \(\frac{d(d+1)}{2}\) × 低位全组合和。
• 累加到 result。
• 更新 prefix_prod *= (d + 1)，固定当前位为 \(d\)，继续处理下一位。

solve

• 调用 sum_not_divisible(10**9) 得到答案。

主程序

• 断言 \(S(100)=2361\) 验证正确性。
• 输出 \(S(10^9)\)。