はじめに
線形漸化式の $N$ 項目の計算はフィボナッチ数列のように行列の累乗で計算すればよいが、$d+1$ 項間線形漸化式の場合に素朴には $O(\mathsf{MM}(d) \log N)$ 時間かかる($\mathsf{MM}(d)$ は $d\times d$ 行列の行列積の計算量)。
これを $O(\mathsf{M}(d) \log N)$ に高速化する方法がある($\mathsf{M}(d)$ は $d$ 次多項式同士の積の計算量)。母関数を用いることでとても簡単に理解できるので書き記しておく。
私は RNG - Editorial CodeChef Discussion から以下のアルゴリズムを考えついた。
従来の手法はこのリンク先の記事の最後に書いてあるものであるが、ニュートン法等を用いた多項式剰余の計算もしくは多項式環のモンゴメリ乗算のテクニックを必要とする。
ここで紹介するアルゴリズムは理解も実装も非常に単純で、従来の手法と比べても定数倍高速である。
2020年9月5日追記: このアルゴリズムを論文で発表した https://arxiv.org/abs/2008.08822。
論文では様々な応用についても書いており、多項式剰余べきの計算($x^N \mod \Gamma(x)$ の計算)もほぼ同じ計算量で計算ができることを示している。
素朴な繰り返し二乗法よりも定数倍高速であり、暗号分野等における実用化を期待している。
線形漸化式と母関数
体 $\mathbb{F}$ 上の数列 $a_0,\,a_1,\,\dotsc,$ が $d+1$ 項間線形漸化式を満たすとは、ある $(c_i\in\mathbb{F})_{i=1,\dotsc,d}$ が存在して
$$
a_n = \sum_{i=1}^d c_i a_{n-i}
$$
が任意の $n\ge d$ について成り立つことを言う。
数列の母関数は形式的冪級数として
$$
G(x) := \sum_{i\ge 0} a_i x^i
$$
と定義される。
母関数 $G(x)$ に対して $Q(x) := 1-\sum_{i=1}^d c_i x^i$ を掛けると、$x^d$ 以上の項はすべて消えるので、高々 $d-1$ 次の多項式 $P(x)$ となる。
よって
$$
G(x) = \frac{P(x)}{Q(x)}
$$
となる。また、逆に高々 $d-1$ 次の多項式 $P(x)$ と$d$ 次の多項式 $Q(x)$ が $Q(0)=1$ を満たし、$G(x) = P(x)/Q(x)$ となるとき、
$G(x)$ が $d+1$ 項間線形漸化式を満たす数列の母関数となることは $G(x)Q(x)$ の $d$ 次以上が0になることから簡単に確認できる。
アルゴリズム
$P(x)$ を高々 $d-1$ 次の多項式、$Q(x)$ を $d$ 次の多項式とし、$Q(0)=1$ と仮定する。今から
$$
[x^N]\; \frac{P(x)}{Q(x)}
$$
を計算するアルゴリズムを考えよう(形式的冪級数 $R(x)$ について $[x^N]\;R(x)$ は $R(x)$ の $x^N$ の係数である)。
$N=0$ のときは $P(0)$ が解である。
分母と分子に $Q(-x)$ を掛けると
\begin{align*}
[x^N]\; \frac{P(x)Q(-x)}{Q(x)Q(-x)}
\end{align*}
が得られる。ここで、$Q(x)Q(-x)$ は偶多項式なので、$V(x^2)=Q(x)Q(-x)$ を満たす多項式 $V(x)$ が存在する。
また、$1/V(x^2)$ も偶形式的冪級数である。
分子を $P(x)Q(-x)=U_{\mathrm{e}}(x^2)+xU_{\mathrm{o}}(x^2)$ というように、偶数部分 $U_{\mathrm{e}}(x^2)$ と奇数部分 $xU_{\mathrm{o}}(x^2)$ に分ける。
すると、 $N$ が偶数のときは
\begin{align*}
[x^N]\; \frac{P(x)Q(-x)}{V(x^2)}
=
[x^N]\; \frac{U_{\mathrm{e}}(x^2)}{V(x^2)}
=
[x^{N/2}]\; \frac{U_{\mathrm{e}}(x)}{V(x)}
\end{align*}
同様に $N$ が奇数のときは
\begin{align*}
[x^N]\; \frac{P(x)Q(-x)}{V(x^2)}
=
[x^N]\; \frac{xU_{\mathrm{o}}(x^2)}{V(x^2)}
=
[x^{(N-1)/2}]\; \frac{U_{\mathrm{o}}(x)}{V(x)}
\end{align*}
よって次の漸化式が得られる。
$$
[x^N]\; \frac{P(x)}{Q(x)}
=
\begin{cases}
P(0),&\text{if $N = 0$}\\
[x^{\frac{N}{2}}]\; \frac{U_{\mathrm{e}}(x)}{V(x)},&\text{if $N$ is even}\\
[x^{\frac{N-1}{2}}]\; \frac{U_{\mathrm{o}}(x)}{V(x)},&\text{otherwise}
\end{cases}
$$
が得られる。
多項式 $U_{\mathrm{e}}$, $U_{\mathrm{o}}$ は高々 $d-1$ 次であり、多項式 $V(x)$ は $d$ 次であり、$V(0)=1$ を満たす。
そのため、上記の漸化式を $\lceil \log (N+1)\rceil$ 回繰り返せば解が得られる。
一回の漸化式の適用には次数が高々 $d$ の多項式乗算を2回すればよい。
そのため計算量(演算回数)は $2\mathsf{M}(d)\lceil\log(N+1)\rceil$ である。
アルゴリズムの実装上の工夫
高速フーリエ変換が使える体 $\mathbb{F}$ の場合、$\mathsf{M}(d) = O(d\log d)$ となる。
時間間引き、周波数間引きの高速フーリエ変換の式を考慮することにより、
漸化式の1回の適用の中で計算する高速フーリエ変換(計4回)のサイズを半分にできる。
その結果、計算量は$(2/3)\mathsf{M}(d)\lceil\log(N+1)\rceil$ となる。
実装例はこちら。ただし、周波数領域ではビットリバーサル置換がかかっているので少し読み辛い。
明文化されているのを見たことがないが、FFTダブリング(長さ $k$ のDFTを長さ $2k$ のDFTに変換する手続き)をする際には、周波数領域でビットリバーサル置換をかかっていると
DFTの偶数番目が前半に連続して並んでいるので都合がよい。