迭代式 · BigInt · Lucas

组合数的计算与溢出

$C(60, 30) = 118264581564861424$ ，约 $1.18 \times 10^{17}$ ；而定义式分子上的 $60!$ 有 82 位十进制数字，量级在 $10^{81}$ 。中间量比答案本身大出六十多个数量级，这段落差在双精度浮点里无处安放。本页量出四种求值实现各自失真的位置，并给出把 $n$ 放大到任意规模时仍然可算的模素数路线。

1 · 阶乘定义式的量级落差

JavaScript 的 Number 是 IEEE 754 双精度，有效位 53 位，可精确表示的最大整数 Number.MAX_SAFE_INTEGER 为 $9007199254740991$ ，约 $9.01 \times 10^{15}$ 。越过这个界，相邻可表示整数之间开始出现空隙，整数加乘不再封闭。

阶乘越界得很早。流行的说法是 $21!$ 起就不再精确，理由是 $21! \approx 5.11 \times 10^{19}$ 已超过 $2^{53}$ 。实测把这个界推后了两位： $21!$ 与 $22!$ 在 Number 里仍与精确值逐位相符，第一个存不下的阶乘是 $23!$ 。原因在于阶乘的二进制表示尾部带着大量零位： $21!$ 含 $2^{18}$ 这个因子，剔掉 18 个零位后有效部分只剩 48 位； $23!$ 的零位是 19 个，有效部分升到 56 位，这才装不下。

警示 · 阶乘式的失效方式不止「差几个数」。Number 能表示的最大值约 $1.8 \times 10^{308}$ ， $170!$ 尚在其内， $171!$ 溢出为 Infinity。此时 $C(171, 1)$ 按定义式算得 Infinity，而答案只是 $171$ ； $C(200, 100)$ 的分子与分母同为 Infinity，商是 NaN。

2 · 乘除交替的迭代式

避开大阶乘的办法是让乘与除交替进行， $C(n, k)$ 的展开式（见组合 C(n, k)）改写成 $k$ 个因子的连乘：

C(n, k) = \prod_{i=0}^{k-1} \frac{n-i}{i+1}

按这个次序求值，每一步的中间量都是整数，除法从不留余数。

定理 2.1 记 $R_j = \prod_{i=0}^{j-1} \frac{n-i}{i+1}$ ，则 $R_j = C(n, j)$ 对 $0 \le j \le k$ 成立。

证明 $R_j$ 展开即 $\frac{n(n-1)\cdots(n-j+1)}{j!}$ ，分子是 $j$ 个连续整数之积。 $j$ 个连续整数之积必被 $j!$ 整除，商恰为 $C(n, j)$ 。逐步看则有 $R_j = R_{j-1} \cdot \frac{n-j+1}{j}$ ，其中 $R_{j-1} = C(n, j-1)$ 已是整数，而 $C(n, j-1)(n-j+1) = C(n, j) \cdot j$ 保证乘出的分子被 $j$ 整除。∎

取 $k \le n/2$ 之后 $R_j$ 单调递增，中间量的峰值就是 $C(n, k)$ 本身。把 Number 换成 BigInt、把 / 换成整数除法，整条链精确无损，代价是 $k$ 次大整数乘除。

3 · 只用加法的递推

第三条路径绕开除法：Pascal 递推 $C(n, k) = C(n-1, k-1) + C(n-1, k)$ 只做加法（见 Pascal 三角与二项式定理）。滚动一行数组、每行从高位往低位覆写即可， $O(nk)$ 次加法、 $O(k)$ 空间。它不能跳算：求单个 $C(n, k)$ 也得把前 $n$ 行推一遍。

换来的是干净的误差行为。加法在 Number 上只要两个加数与和都不越过 $2^{53}$ 就是精确的，而递推路径上出现的全是 $C(i, j)$ 本身，没有比答案更大的量。

4 · 浮点实现的失真边界

IEEE 754 double · 53-bit · MAX_SAFE_INTEGER = 9007199254740991

图 4-1 · 阶乘式、迭代式、Pascal 递推与 BigInt 精确值的同台对照。可拖动

n

与

k

观察哪几行先与精确值分道，并切换素数看模

p

的余数。

三种 Number 实现的首处失真位置各不相同，扫描 $n \le 100$ 的全部 $(n, k)$ 与 BigInt 精确值逐对比对，结果如下：

实现	首处失真	该处精确值	是否越过 $2^{53}$
阶乘式 · `Number`	$C(55, 24)$	2488589544741300	否
迭代式 · `Number`	$C(56, 23)$	3167295784216200	否
Pascal 递推 · `Number`	$C(57, 25)$	9929472283517787	是
迭代式 · `BigInt`	无	不适用	不适用

Pascal 递推的失真点与「第一个越过 $2^{53}$ 的组合数」重合： $9929472283517787$ 正是所有 $C(n, k)$ 中首个超过 $9007199254740991$ 的值。扫描 $n \le 300$ 内全部未越界的 $(n, k)$ ，递推实现一处不差——它撑到了信息论允许的最后一格。乘除式则不然。

警示 · 本仓库 combinatorics/core/combinatorics.ts 的 comb() 走的是浮点迭代式加 Math.round。它与精确值的第一处分歧在 $C(56, 23)$ ：精确值 $3167295784216200$ ，comb() 给出 $3167295784216201$ ，差 $1$ ，相对误差 $3.2 \times 10^{-16}$ 。这个精确值只有 $3.17 \times 10^{15}$ ，仍在安全整数以内，失真来自中间量的舍入累积而非结果装不下。 $n \le 80$ 范围内共 23 个 $(n, k)$ 落在这类分歧上。本系列既有页面用到的 $n$ 都在 10 以内，远在这条界线之内；此处记的是适用范围。

5 · 模素数的组合数与 Lucas 定理

Lucas · C(n, k) ≡ ∏ C(nᵢ, kᵢ) (mod p)

许多场合只需要 $C(n, k) \bmod p$ 。 $p$ 为素数时，费马小定理给出 $a^{p-1} \equiv 1 \pmod p$ ，于是 $a^{p-2}$ 就是 $a$ 的逆元（modular inverse），除以 $k!$ 改写为乘上 $(k!)^{p-2} \bmod p$ ，全程在 $[0, p)$ 内运算。这条路要求 $n < p$ ：否则 $n!$ 含因子 $p$ ，模 $p$ 得零，逆元不存在。

$n$ 超过 $p$ 时接上 Lucas 定理。把 $n$ 与 $k$ 写成 $p$ 进制， $n = \sum_i n_i p^i$ 、 $k = \sum_i k_i p^i$ ，则

C(n, k) \equiv \prod_i C(n_i, k_i) \pmod p

每位的 $n_i$ 与 $k_i$ 都小于 $p$ ，逐位取组合数再相乘即可，位数只有 $\lfloor \log_p n \rfloor + 1$ 。某一位上 $k_i > n_i$ 时该位的组合数为零，整个乘积随之归零——这也给出了 $p$ 整除 $C(n, k)$ 的一个判据（Kummer 定理的特例）。

图 5-1 ·

n

与

k

的

p

进制分解与逐位组合数。可切换素数

p

观察位数如何伸缩，某一位上

k

的数字超过

n

的数字时乘积归零。

6 · 参考文献

Binomial coefficient. Wikipedia. 组合数的多种计算方法与整除性质。https://en.wikipedia.org/wiki/Binomial_coefficient
Lucas's theorem. Wikipedia. Lucas 定理的陈述与 $p$ 进制证明。https://en.wikipedia.org/wiki/Lucas%27s_theorem
Fermat's little theorem. Wikipedia. 费马小定理及其在模素数求逆上的用法。https://en.wikipedia.org/wiki/Fermat%27s_little_theorem
Double-precision floating-point format. Wikipedia. IEEE 754 双精度的 53 位有效位与安全整数范围。https://en.wikipedia.org/wiki/Double-precision_floating-point_format