読者です 読者をやめる 読者になる 読者になる

階乗 mod 素数

$n! \pmod{p}$ を計算量 $O(p^{1/2} \log{p})$ で求める.

計算式

多項式 $f(x) := \prod_{i=1}^n (x+i)$ を定義したとき,$n! = f(0)$ である.

$v := \lfloor{\sqrt{n}}\rfloor$ とおいて,$g(x) := \prod_{i=1}^{v} (x+i)$ とおくと,

$$
n! = \left( \prod_{i=0}^{v-1} g(vi) \right) \cdot \prod_{i=v^2+1}^{n} i
$$
とかける.後者の積は愚直に $O(p^{1/2})$ で計算できるので無視することにすると,$g(0), g(v), \ldots, g(v(v-1))$ の値を高速に求められればよい.

手法

1. $O(p^{1/2} (\log{p})^2)$ で求める方法.

multipoint evaluation を用いて $g(0), \ldots, g(v(v+1))$ を同時に評価すればよい.ただ,(定数倍の重い)多項式除算を駆使するアルゴリズムのため,$O(p)$ の計算より劇的に早い というわけでもない.

高速よりな実装した場合でも,$p = 10^9+7$ で 0.30 秒前後かかる.

2. $O(p^{1/2} \log{p})$ で求める方法.

$g(x)$ の根 $-1, -2, \ldots, -v$ は等差数列をなしていて,かつ $g(x)$ で評価したい値 $0, v, \ldots, (v-1)v$ も等差数列をなしている.この場合,うまく Lagrange interpolation を用いると,計算量から $\log p$ を落とすことができる.

$d$ 次多項式 $g_d(x) := \prod_{i=1}^{d} (x+i)$ に対して,$g_d(0), g_d(v), \ldots, g_d(dv)$ の値を計算できているとする.このとき,$O(d \log d)$ の計算量で $g_{2d}(0), g_{2d}(v), ..., g_{2d}(2dv)$ を計算できれば,示した計算量を達成できる.

まず,$g_{2d}(x) = g_{d}(x) \cdot g_d(x+d)$ であるので,評価値の組 $G_d(i)$ を $G_d(i) := \left(g_{d}(i), g_{d}(v+i), \ldots, g_{d}(dv+i)\right)$ と定義するとき,$G_d(0)$ から $G_d(d), G_d(dv), G_d(dv+d)$ の 3 つの組を補間できればよい.これらの計算は,Lagrange Interpolation と FFT(NTT) により $O(d \log d)$ で行える [1].

ある $d$ 次多項式 $h(x)$ について $h(0), h(1), \ldots, h(d)$ の値が与えられているとき, $h(m+k)$ の値は Lagrange Interpolation により($m+k-j$ が逆元を持つ範囲で),
$$
\begin{align*}
h(m+k) &= \sum_{i=0}^{d} h(i) \prod _{j=0, i\ne j}^{d} \frac{m+k-j}{i-j} \\
&= \left(\prod_{j=0}^{d} (m+k-j)\right) \left( \sum_{i=0}^{d} \frac{h(i)}{i! (d-i)!(-1)^{d-i}} \cdot \frac{1}{m+k-i} \right)
\end{align*}
$$
と計算できる.ここで,右括弧内が畳み込みの形になっていることに気づくと,$h(m), h(m+1), ..., h(m+d)$ は $O(d \log d)$ で計算できる(middle product を用いるとよい).また,この補間を $G_d(0)$ から $G_d(a)$ に適用する際には,$m := av^{-1}$とすればよい.

以上により,$N! \pmod{p}$ が $O(p^{1/2} \log{p})$ で計算できる.

C++ での実装では,$p = 10^9 + 7, N = (p-1)/2$ でも 0.060 秒程度で計算できる.$p \approx 10^{12}$ でも 3 秒程度で計算できるだろう.

参考文献

[1] Alin Bostan, Pierrick Gaudry, and Eric Schost, Linear recurrences with polynomial coefficients and application to integer factorization and Cartier-Manin operator.