154. 探索帕斯卡金字塔(Exploring Pascal's Pyramid)

用球形小球构造一个三角形金字塔,使得每个小球都恰好放在下一层的三个小球之上。

然后,我们计算从顶点通向每个位置的路径数量:

一条路径从顶点开始,并向下移动到当前位置正下方的三个球中的任意一个。

因此,到达某个位置的路径数量等于它正上方那些数字之和(根据位置不同,正上方最多有三个数字)。

其结果就是帕斯卡金字塔,并且第 \(n\) 层中的数字是三项式展开式 \((x + y + z)^n\) 的系数。

\((x + y + z)^{200000}\) 的展开式中,有多少个系数是 \(10^{12}\) 的倍数?

一、题意概述

\(n\) 层的每个位置对应一个三元组 \((a,b,c)\),满足:

$$a + b + c = n,\quad a,b,c \ge 0$$

对应的三项式系数为:

$$\binom{n}{a,b,c} = \frac{n!}{a!b!c!}$$

本题要求统计当 \(n=200000\) 时,有多少个有序三元组 \((a,b,c)\) 使得该系数能被 \(10^{12}\) 整除。由于 \(10^{12}=2^{12}5^{12}\),条件等价于:

$$v_2\left(\binom{n}{a,b,c}\right) \ge 12,\quad v_5\left(\binom{n}{a,b,c}\right) \ge 12$$

其中 \(v_p(m)\) 表示质数 \(p\)\(m\) 中的指数。

二、数学背景

对任意质数 \(p\),Legendre 公式给出:

$$v_p(m!) = \frac{m - s_p(m)}{p-1}$$

其中 \(s_p(m)\)\(m\)\(p\) 进制下的数位和。因此:

$$ \begin{aligned} v_p\left(\binom{n}{a,b,c}\right) &= v_p(n!) - v_p(a!) - v_p(b!) - v_p(c!) \\ &= \frac{s_p(a)+s_p(b)+s_p(c)-s_p(n)}{p-1} \end{aligned} $$

这也正是 Kummer 定理在多项式系数上的形式:该指数等于用 \(p\) 进制把 \(a,b,c\) 相加得到 \(n\) 时产生的进位总数。

对本题的两个质数:

直接枚举所有 \((a,b,c)\) 需要约 \(\frac{(n+1)(n+2)}{2}\) 次检查,数量级约为 \(2\times 10^{10}\),不可行。

三、算法分析

候选方案一:直接枚举加剪枝

可以只枚举 \(a\le b\le c\),再按三元组是否有相等元素乘以 \(1,3,6\) 的对称系数。这能把枚举量缩小约 \(6\) 倍,但仍然是数十亿级别,不适合 Python。

继续按五进制数位和剪枝,能减少许多无效点,但仍需要在大量区间内逐点扫描,性能不稳定。

候选方案二:进位 DP 加位集校正

本解法采用该方案。

先只考虑 \(5^{12}\)。因为五进制位数很少,可以用进位 DP 一次性统计所有有序三元组中满足 \(v_5\ge 12\) 的数量。状态为:

$$\mathrm{dp}[\text{carry}][\text{carry\_sum}]$$

逐位处理 \(n\) 的五进制数字,枚举当前位上 \(a,b,c\) 的三个数字,更新下一位进位与累计进位数。累计进位数只需要截断到 \(12\),因为更大的值在本题中没有区别。

接着扣除“已经满足 \(5^{12}\),但不满足 \(2^{12}\)”的三元组。二进制不足条件为:

$$\operatorname{popcount}(a)+\operatorname{popcount}(b)+\operatorname{popcount}(c) \le \operatorname{popcount}(n)+11$$

这是一个较强的限制。为了避免逐点枚举,构造按 \((\operatorname{popcount}(x), s_5(x))\) 分组的位集:

固定最小元素 \(a\),只考虑 \(a\le b\le c\),即:

$$a \le b \le \frac{n-a}{2}$$

对每组候选的 \(b\),用右移后的反向位集表示满足条件的 \(c=n-a-b\)。随后通过位与和 bit_count() 一次统计整段区间内的合法数量。最后按三元组元素相等情况把无序计数换成有序计数。

最终答案为:

$$ \#(v_5\ge 12) - \#(v_5\ge 12 \text{ 且 } v_2<12) $$

四、复杂度分析

五进制进位 DP 的状态很小,复杂度约为 \(O(\log_5 n \cdot 12 \cdot 5^3)\)

位集校正部分按 \(a\)\(0\) 遍历到 \(\lfloor n/3\rfloor\)。每次查询不是逐个枚举 \(b\),而是在 Python 大整数位集上进行批量位运算;设分组计划数为 \(G\)、机器字宽为 \(w\),可视为约 \(O(\frac{n}{3}\cdot G\cdot \frac{n}{w})\) 的底层字操作。实际 \(n=200000\) 时分组数量较小,运行时间低于 \(60\) 秒。

空间上主要保存若干按特征分组的位集,数量由可能的 popcount 和五进制数位和决定,空间复杂度约为 \(O(C\cdot n)\) 位,其中 \(C\) 是非空特征组数量。

五、代码实现与说明

完整代码位于 scripts/pe154.py。核心函数按代码顺序如下。

基础数位工具

def digits_base(number: int, base: int) -> list[int]:
    """Return the little-endian digits of number in the given base."""
    if number == 0:
        return [0]

    digits: list[int] = []
    while number:
        digits.append(number % base)
        number //= base
    return digits

digits_base 将整数转换成低位在前的进制数字列表,方便从最低位开始模拟加法进位。digit_sum_basebuild_digit_sum_table 分别计算单个整数和整段区间的进制数位和。

五进制进位 DP

def count_prime_carry_at_least(level: int, prime: int, threshold: int) -> int:
    """Count ordered triples whose base-prime carry sum is at least threshold."""
    if threshold <= 0:
        return (level + 1) * (level + 2) // 2

    transitions: list[list[list[tuple[int, int]]]] = []
    for digit in range(prime):
        by_carry: list[list[tuple[int, int]]] = []
        for carry_in in range(3):
            counts: dict[int, int] = {}
            for first in range(prime):
                for second in range(prime):
                    for third in range(prime):
                        total = first + second + third + carry_in
                        if total % prime == digit:
                            carry_out = total // prime
                            counts[carry_out] = counts.get(carry_out, 0) + 1
            by_carry.append(list(counts.items()))
        transitions.append(by_carry)

这一段预先枚举一位上的所有数字组合。对于给定的当前位目标数字 digit 和输入进位 carry_in,记录每种输出进位 carry_out 出现多少次。

    dp = [[0] * (threshold + 1) for _ in range(3)]
    dp[0][0] = 1

    for digit in digits_base(level, prime):
        next_dp = [[0] * (threshold + 1) for _ in range(3)]
        for carry_in, row in enumerate(dp):
            for carry_sum, ways in enumerate(row):
                if ways == 0:
                    continue
                for carry_out, transition_count in transitions[digit][carry_in]:
                    next_sum = min(threshold, carry_sum + carry_out)
                    next_dp[carry_out][next_sum] += ways * transition_count
        dp = next_dp

    return dp[0][threshold]

dp[carry][carry_sum] 表示处理到当前位时,输出进位为 carry、累计进位数为 carry_sum 的方案数。carry_sum 被截断到 threshold,因此所有已经达标的方案会汇总在最后一格。所有位处理完后,只接受最终进位为 0 的方案。

位集构造

def build_feature_bitsets(
    level: int,
) -> tuple[list[int], list[int], list[list[int]], list[list[int]]]:
    """Build exact and reversed bitsets grouped by binary and quinary digit sums."""
    popcounts = [value.bit_count() for value in range(level + 1)]
    quinary_sums = build_digit_sum_table(level, 5)

这一段预计算每个整数的二进制 popcount 与五进制数位和。之后的校正逻辑只使用这两个特征。

    for value in range(level + 1):
        popcount = popcounts[value]
        quinary_sum = quinary_sums[value]
        exact[popcount][quinary_sum] |= 1 << value
        reversed_exact[popcount][quinary_sum] |= 1 << (level - value)

exact 用第 value 位表示该值属于某个特征组;reversed_exact 把位置反向为 level - value。固定 a 后,reversed_exact >> a 的第 b 位就对应 \(c=level-a-b\)

校正计划与短缺计数

def make_correction_plan(
    popcount_budget: int,
    quinary_sum_needed: int,
    exact_nonempty: Iterable[tuple[int, int, int]],
    reversed_unions: list[list[int]],
) -> list[tuple[int, int]]:
    """Group bitsets needed for one correction query."""

校正计划把同一类 c 约束对应的所有 b 位集合并起来,减少主循环中的位运算次数。popcount_budget 是在固定 \(a\) 后,\(b+c\) 还能使用的 popcount 上限;quinary_sum_needed\(b+c\) 还需要达到的五进制数位和下限。

def count_power_of_two_shortfall(level: int, threshold: int) -> int:
    """Count ordered triples with enough factors of 5 but too few factors of 2."""

该函数统计需要从五进制 DP 结果中扣除的数量。它只枚举 \(a\le b\le c\) 的代表元:

    max_first = level // 3
    high = level // 2
    range_mask = (1 << (high + 1)) - 1

range_mask 表示当前合法的 \(b\) 区间。每次 a 增加后,低端和高端只会移动一位,因此代码用异或增量更新掩码,避免重复构造长区间位集。

                subtotal = 0
                for bits_c, bits_b in plan:
                    subtotal += (bits_b & (bits_c >> first) & range_mask).bit_count()
                correction += 6 * subtotal

这里完成核心统计:bits_b 给出候选 \(b\)bits_c >> first\(c=level-first-b\) 映射到 \(b\) 坐标,range_mask 限制 \(a\le b\le c\)。三者取交后用 bit_count() 得到本轮代表元数量,先按一般情况乘以 \(6\)

随后代码分别处理 \(a=b<c\)\(a<b=c\)\(a=b=c\) 三种边界。一般情况乘以 \(6\),两个相等时实际只有 \(3\) 个排列,三个全相等时实际只有 \(1\) 个排列,因此需要做小幅修正。

总求解入口

def solve(level: int = DEFAULT_LEVEL, threshold: int = DEFAULT_THRESHOLD) -> int:
    """Return the count of trinomial coefficients divisible by 10^threshold."""
    enough_fives = count_prime_carry_at_least(level, 5, threshold)
    two_shortfall = count_power_of_two_shortfall(level, threshold)
    return enough_fives - two_shortfall

solve 先统计所有满足 \(5^{12}\) 的系数,再扣除其中 \(2^{12}\) 不足的部分,即得到同时满足 \(2^{12}\)\(5^{12}\) 的系数数量。

文件末尾提供了 brute_forcerun_self_tests。它们只用于小规模对照:直接枚举所有三元组,并用数位和公式判断 \(v_2\)\(v_5\),从而验证优化算法的正确性。