階乗 mod 素数

$n! \bmod 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)$ の計算より劇的に早い というわけでもない。

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{aligned} 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{aligned} $$ と計算できる。ここで、右括弧内が畳み込みの形になっていることに気づくと、$h(m), h(m+1), ..., h(m+d)$ は $O(d \log d)$ で計算できる。また、この補間を $G_d(0)$ から $G_d(a)$ に適用する際には、$m := av^{-1}$とすればよい。

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

追記

数列 $(a_i)_{i=0}^{\infty}$ が p-recursive となっている場合に利用できることが多い方法です。(行列積の部分を上の方法を用いて高速化します。)

  • $\left[\begin{matrix}n! \end{matrix}\right] = \left(\prod\limits_{i=0}^{n-1} \left[\begin{matrix}i + 1 \end{matrix}\right]\right) \left[\begin{matrix}1 \end{matrix}\right]$
  • $\left[\begin{matrix}!n \\ !(n+1) \end{matrix}\right] = \left(\prod\limits_{i=0}^{n-1} \left[\begin{matrix}0 & 1 \\ i+1 & i+1 \end{matrix}\right]\right) \left[\begin{matrix}1 \\ 0 \end{matrix}\right]$ (行列積は右から順に $i = 1, \dots$) (!n: subfactorial)
    • $\left[\begin{matrix}!n \\ (-1)^{n+1} \end{matrix}\right] = \left(\prod\limits_{i=0}^{n-1} \left[\begin{matrix}i+1 & 1 \\ 0 & -1 \end{matrix}\right]\right) \left[\begin{matrix}1 \\ -1 \end{matrix}\right]$
  • $\left[\begin{matrix}(n+1)! \\ \sum_{i=0}^n i! \end{matrix}\right] = \left(\prod\limits_{i=0}^{n} \left[\begin{matrix}i +1 & 0 \\ 1 & 1 \end{matrix}\right]\right) \left[\begin{matrix}1 \\ 0 \end{matrix}\right]$
  • $\left[\begin{matrix} \binom{n}{m+1} \\ \sum_{i=0}^{m} \binom{n}{i} \end{matrix}\right] = \left(\prod\limits_{i=0}^{m} \left[\begin{matrix}(n-i)/(i+1) & 0 \\ 1 & 1 \end{matrix}\right]\right) \left[\begin{matrix}1 \\ 0 \end{matrix}\right] = \frac{1}{(m+1)!} \left(\prod\limits_{i=0}^{m} \left[\begin{matrix}n - i & 0 \\ i + 1 & i +1 \end{matrix}\right]\right) \left[\begin{matrix}1 \\ 0 \end{matrix}\right]$
  • $\left[\begin{matrix} \genfrac{[}{]}{0pt}{}{n+1}{1} \\ \genfrac{[}{]}{0pt}{}{n+1}{2} \\ \genfrac{[}{]}{0pt}{}{n+1}{3} \end{matrix}\right] = \left(\prod\limits_{i=0}^{n-1} \left[\begin{matrix}i+1 & 0 & 0 \\ 1 & i+1 & 0 \\ 0 & 1 & i+1 \end{matrix}\right]\right) \left[\begin{matrix}1 \\ 0 \\ 0\end{matrix}\right]$

(後ろ 2 つの記事は怪しげな導出をしていますが、行列で関係式を作るのが楽だと思います。)

参考文献

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