二項係数 mod Prime Power

$\binom{n}{m} \pmod{p^e}$ を高速に解く。($n, m, p^e < 10^{300}$, $pe < 10^6$, $p$: 素数

$\newcommand{\stirlingfirst}[2]{\genfrac{[}{]}{0pt}{}{#1}{#2}}$Andrew Granville さんの Binomial coefficients modulo prime powers (BinCoeff.pdf) とは違った計算方法。

前計算後の $\binom{n}{m} \pmod{p^e}$ の計算量(乗算・剰余算の回数)を、$O\left(e^2 \left(\frac{\log{n}}{\log{p}} + \min(\log{n}, p \log {p}) \right) \right)$ (?) から $O(e \log{n})$ にできる。

($\stirlingfirst{n}{m}$: 第 1 種 Stirling 数)

計算方法

$(n!)_p$ を $\gcd(i, p) = 1$ を満たす $i$ ($1 \le i \le n$) の積と定義するとき、$(n!)_p$ $\pmod{p^e}$ が求められれば十分。

まず $u, v$ ($0 \le v < p$) を用いて $n = up + v$ とかく。このとき、
$$
\begin{align*}
((up + v)!)_p
&= \left(\prod_{i=0}^{u-1} \prod_{j=1}^{p-1} (ip + j)\right) \cdot \prod_{j=1}^{v} (up+j) \\
&\equiv \left(\prod_{i=0}^{u-1} \left( \sum_{k=0}^{e-1} (ip)^k \stirlingfirst{p}{k+1} \right) \right) \cdot \left(\sum_{k=0}^{e-1} (up)^k \stirlingfirst{v+1}{k+1}\right) \pmod{p^e} \\
&= \stirlingfirst{p}{1}^u \left(\prod_{i=0}^{u-1} \left(1 + \sum_{k=1}^{e-1} \frac{\stirlingfirst{p}{k+1}}{\stirlingfirst{p}{1}} (ip)^k \right) \right) \cdot \left(\sum_{k=0}^{e-1} (up)^k \stirlingfirst{v+1}{k+1}\right) \\
&\equiv \stirlingfirst{p}{1}^u f_{p,e}(u) \left( \sum_{k=0}^{e-1} (up)^k \stirlingfirst{v+1}{k+1} \right) \pmod{p^e}
\end{align*}
$$
が成り立つ。ここで、$f_{p,e}(x)$ は適当な $2e-2$ 次の多項式である(後で証明)。よって、$f_{p,e}(0), \ldots, f_{p,e}(2e-2)$ を計算しておくことで $O(e \log{p})$ で $f_{p,e}(u)$ の値を補間できる。

したがって、適当な前計算を行っておけば、$(n!)_p$ は $O(e \log{p})$ で計算できる。

実装方法と計算量

1. $\stirlingfirst{n}{m}$ ($1 \le n \le p$, $1 \le m \le \min(p, e)$) の前計算: $O(p \cdot \min(p, e))$
2. $f_{p,e}(0), ..., f_{p,e}(2e-2)$ や多項式補間用の前計算: $O(e \cdot \min(p, e) + e \log{p})$
3. 必要な $(n_i!)_p$ の計算: $O(\log_p{n} \cdot (e \log{p})) = O(e \log{n})$

合計で $O(pe + e \log n)$ 回の乗算・剰余算で、$\binom{n}{m} \pmod{p^e}$ を計算できる。

多項式となる証明

任意に $e, p \in \mathbb{Z}$ ($e \ge 1$, $p \ge 2$) をとり、$a_k \in \mathbb{Q}$ で $a_k$ の分母は $p$ と互いに素とする。このとき、
$$
g(x, u) := \prod_{i=0}^{u-1} \left(1 + \sum_{k=1}^{e-1} a_k (xi)^k\right)
$$
とすると、適当な $b_k \in \mathbb{Q}$ を用いて、
$$
\begin{align*}
\log(g(x, u)) &= \sum_{i=0}^{u-1} \log{\left(1 + \sum_{k=1}^{e-1} a_k (xi)^k \right)} \\
&= \sum_{i=0}^{u-1} \sum_{k=1}^{\infty} b_k x^k i^k \\
&= \sum_{k=1}^{\infty} b_k x^k s_k(u), \\
g(x, u) &= \exp{\left(\sum_{k=1}^{\infty} b_k x^k s_k(u)\right)}\\
&= 1 + \sum_{k=1}^{\infty} t_{k}(u) \cdot x^k
\end{align*}
$$
とかける。ここで、$s_k(x)$ は $s_k(x) := \sum_{i=0}^{x-1} i^k$ なる $(k+1)$-次の多項式で、$t_k(x)$ は適当な $(2k)$-次の多項式となる。

最初の $g(x, u)$ の展開式 と 最後の $g(x, u)$ の $x^k$ の係数は一致することから、$t_k(u) \in \mathbb{Q}$ かつ $t_k(u)$ の分母は $p$ と互いに素であり、
$$
g(x, u)|_{x=p} \equiv 1 + \sum_{k=1}^{e-1} t_k(u) \cdot p^k \pmod{p^e}
$$
とかける。よって $g(x, u)|_{x=p}$ は $\bmod p^e$ で $u$ の $2e-2$ 次の多項式となる。

その他

Andrew さんの方法は多倍長整数なしで実装しようとすると手間がかかるのだけど、この方法だと多倍長整数を持たなくても実装しやすい。

  • $p$ が素数のとき、$\stirlingfirst{p}{i}$ ($2 \le i < p$) は $p$ で割り切れる
  • $p \ge 5$ のとき $\stirlingfirst{p}{2}$ は $p^2$ で割り切れる

などといった性質を用いれば、もう少し $f_{p,e}(x)$ の次数を落とせるのだと思う。

ソースコード (Python)

この実装で、$$\binom{199^{197}}{193^{191}} \pmod{181^{179}}$$ を 3 秒程度で計算できる。