はじめに

線形漸化式の $n$ 項目の計算はフィボナッチ数列のように行列の累乗で計算すればよいが、$k+1$ 項間線形漸化式の場合に素朴には $O(k^3 \log n)$ 時間かかる。 これを $O(k\log k \log n)$ に高速化する方法がある。多項式を用いることでとても簡単に理解できるので書き記しておく。 私は RNG - Editorial CodeChef Discussion から以下のアルゴリズムを考えついた。 通常はこのリンク先の記事の最後に書いてある手法が紹介されるが、それは多項式環上のモンゴメリ乗算など非直感的なテクニックを必要とする。 ここで紹介するアルゴリズムは $k$ 次多項式の乗算が高速フーリエ変換で $O(k\log k)$ 時間で計算できることだけ分かっていれば理解できる。

線形漸化式と母関数

体 $\mathbb{F}$ 上の数列 $a_0,\,a_1,\,\dotsc,$ が $k+1$ 項間線形漸化式を満たすとは、ある $(c_i\in\mathbb{F})_{i=1,\dotsc,k}$ が存在して $$ a_n = \sum_{i=1}^k c_i a_{n-i} $$ が任意の $n\ge k$ について成り立つことを言う。 数列の母関数は形式的冪級数として $$ G(x) := \sum_{i\ge 0} a_i x^i $$ と定義される。 母関数 $G(x)$ に対して $Q(x) := 1-\sum_{i=1}^k c_i x^i$ を掛けると、$x^k$ 以上の項はすべて消えるので、高々 $k-1$ 次の多項式 $P(x)$ となる。 よって $$ G(x) = \frac{P(x)}{Q(x)} $$ となる。また、逆に高々 $k-1$ 次の多項式 $P(x)$ と$k$ 次の多項式 $Q(x)$ が $Q(0)=1$ を満たし、$G(x) = P(x)/Q(x)$ となるとき、 $G(x)$ が $k+1$ 項間線形漸化式を満たす数列の母関数となることは $G(x)Q(x)$ の $k$ 次以上が0になることから簡単に確認できる。 形式的羃級数 $G(x)$ の偶数次数(または奇数次数)だけからなる形式的羃級数は $\mathcal{E}(G(x)) := (G(x) + G(-x))/2$ (または $\mathcal{O}(G(x)) := (G(x) - G(-x))/2$) である。

アルゴリズム

$P(x)$ を高々 $k-1$ 次の多項式、$Q(x)$ を $k$ 次の多項式とし、$Q(0)=1$ と仮定する。今から $$ [x^n] \frac{P(x)}{Q(x)} $$ を計算するアルゴリズムを考えよう(形式的冪級数 $R(x)$ について $[x^n]R(x)$ は $R(x)$ の $x^n$ の係数である)。 $n=0$ のときは $P(0)$ が解である。 $n$ が偶数のときは \begin{align*} [x^n] \frac{P(x)}{Q(x)} &= [x^n] \left(\frac{P(x)}{Q(x)} + \frac{P(-x)}{Q(-x)}\right)\frac12\\ &= [x^n] \frac{\mathcal{E}(P(x)Q(-x))}{Q(x)Q(-x)} \end{align*} 同様に $n$ が奇数のときは \begin{align*} [x^n] \frac{P(x)}{Q(x)} &= [x^n] \left(\frac{P(x)}{Q(x)} - \frac{P(-x)}{Q(-x)}\right)\frac12\\ &= [x^n] \frac{\mathcal{O}(P(x)Q(-x))}{Q(x)Q(-x)} \end{align*} よって次の漸化式が得られる。 $$ [x^n] \frac{P(x)}{Q(x)} = \begin{cases} [x^n] \frac{\mathcal{E}(P(x)Q(-x))}{Q(x)Q(-x)},&\text{if $n$ is even}\\ [x^n] \frac{\mathcal{O}(P(x)Q(-x))}{Q(x)Q(-x)},&\text{otherwise} \end{cases} $$ ここで分母 $Q(x)Q(-x)$ は偶多項式である。 そのため、$Q'(x^2) = Q(x)Q(-x)$となる $k$ 次多項式 $Q'(x)$ が存在する。 さらに $x^2$ を $x$ に置き換えることによって $$ [x^n] \frac{P(x)}{Q(x)} = \begin{cases} [x^{\frac{n}{2}}] \frac{P_e'(x)}{Q'(x)},&\text{if $n$ is even}\\ [x^{\frac{n-1}{2}}] \frac{P_o'(x)}{Q'(x)},&\text{otherwise} \end{cases} $$ が得られる。 ここで $P_e'(x^2) := \mathcal{E}(P(x)Q(-x))$, $xP'_o(x^2) := \mathcal{O}(P(x)Q(-x))$ である。 多項式 $P'_e$, $P'_o$ は高々 $k-1$ 次であり、多項式 $Q'(x)$ は $k$ 次であり、$Q'(0)=1$ を満たす。 そのため、上記の漸化式を $\lfloor \log n\rfloor+1$ 回繰り返せば解が得られる。 一回の漸化式の適用に次数が高々 $k$ の多項式乗算を2回すればよい。これは高速フーリエ変換を使うと $O(k\log k)$ 時間でできる。 よってアルゴリズムの時間計算量は $O(k\log k\log n)$ である。

アルゴリズムの実装上の工夫

時間間引き、周波数間引きの高速フーリエ変換の式を考慮することにより、 漸化式の1回の適用の中で計算する高速フーリエ変換(計4回)のサイズを半分にできる。 例えば このコード。ただし、周波数領域ではビットリバーサル置換がかかっているので少し読み辛い。