階乗 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.