線形漸化式を持つ数列の連続する項を高速に計算する

概要

$d$ 項間の線形漸化式を持つ数列の第 $k$ 項から $k + d$ 項までを $\operatorname{O} \left( d \log \left( d \right) \log \left( k \right) \right)$ の時間計算量で計算するアルゴリズムを説明する。繰り返し二乗法と多項式の余り付き除算を用いる方法と比べると、定数倍が良いと思われる。

問題設定

$P \left( x \right), Q \left( x \right) \in \mathbb{K} \left\lbrack x \right\rbrack$ が $\left\lbrack x ^ 0 \right\rbrack Q \left( x \right) \neq 0, \operatorname{deg} \left( Q \right) = d, \operatorname{deg} \left( P \right) \lt d$ を満たすとする。

この記事内で多項式の除算をしたら、形式的ローラン級数体 $\mathbb{K} \left( \left( x \right) \right)$ 上で考えているものとする。これは形式的冪級数環 $\mathbb{K} \left\lbrack \left\lbrack x \right\rbrack \right\rbrack$ と概ね同じだが、 $c x ^ {- 2}$ のように負の指数の項も有限個持ってよいとするものである。

アルゴリズムの大枠

maspy による記事 [多項式・形式的べき級数]（３）線形漸化式と形式的べき級数 | maspyのHP の議論から、以下の事実が分かる。

$P \left( x \right) / Q \left( x \right)$ の $x ^ k$ 以降の係数を並べた数列の母関数は、 $d$ 次未満の多項式 $S$ を用いて $S \left( x \right) / Q \left( x \right)$ と表される。
$S \left( x \right) = x ^ {- k} P \left( x \right) \pmod {Q \left( x \right)}$ である。

これを用いて、以下のように計算する。

$1 / Q \left( x \right)$ を母関数とする数列の第 $\left\lbrack k, k + d \right)$ 項を計算する (後述)。
これに $Q \left( x \right)$ を掛けることで、 $P \left( x \right) = 1$ の場合の $S \left( x \right)$ が求まる。これは $x ^ {- k} \pmod {Q \left( x \right)}$ である。
$S \left( x \right) = x ^ {- k} P \left( x \right) \pmod {Q \left( x \right)}$ より $S \left( x \right)$ が分かる。
FPS の逆元を計算するアルゴリズムを用いて $S \left( x \right) / Q \left( x \right)$ を計算すれば、元の数列の第 $k$ 項以降を好きなだけ計算できる。

$1 / Q \left( x \right)$ の $x ^ k$ から $x ^ {k + d - 1}$ までの係数の計算

$n = k + d - 1$ とすると、 $\left( x ^ n Q \left( x \right) \right) ^ {- 1}$ の $x ^ {- d + 1}$ から $x ^ 0$ までの係数を求めればよい。

$n = 0$ の場合

$\left\lbrack x ^ 0 \right\rbrack Q \left( x \right) \neq 0$ より $\left\lbrack x ^ 0 \right\rbrack \left( x ^ n Q \left( x \right) \right) ^ {- 1} = \left( \left\lbrack x ^ 0 \right\rbrack Q \left( x \right) \right) ^ {- 1}$ であり、それ以外は $0$ となる。

$n \gt 0, n = 2m$ の場合

$V \left( x ^ 2 \right) = Q \left( x \right) Q \left( - x \right)$ と置き、再帰的に $\left( x ^ m V \left( x \right) \right) ^ {- 1}$ の $x ^ {- d + 1}$ から $x ^ 0$ までの係数を求める。すると $\left( x ^ {2m} V \left( x ^ 2 \right) \right) ^ {- 1}$ の $x ^ {- 2d + 2}$ から $x ^ {0}$ までの係数が求まる。奇数次の係数は必ず $0$ になるから、 $x ^ {- 2d + 1}$ から $x ^ {0}$ までの係数が求まったことになる。ここに $Q \left( - x \right)$ を掛けることで、 $Q \left( - x \right) \left( x ^ {2m} Q \left( x \right) Q \left( - x \right) \right) ^ {- 1} = \left( x ^ n Q \left( x \right) \right) ^ {- 1}$ の $x ^ {- d + 1}$ から $x ^ 0$ までの係数が求まる。

$n \gt 0, n = 2m + 1$ の場合

$V \left( x ^ 2 \right) = Q \left( x \right) Q \left( - x \right)$ と置き、再帰的に $\left( x ^ m V \left( x \right) \right) ^ {- 1}$ の $x ^ {- d + 1}$ から $x ^ 0$ までの係数を求める。すると $\left( x ^ {2m + 1} V \left( x ^ 2 \right) \right) ^ {- 1}$ の $x ^ {- 2d + 1}$ から $x ^ {- 1}$ までの係数が求まる。偶数次の係数は必ず $0$ になるから、 $x ^ {- 2d + 1}$ から $x ^ {0}$ までの係数が求まったことになる。ここに $Q \left( - x \right)$ を掛けることで、 $Q \left( - x \right) \left( x ^ {2m + 1} Q \left( x \right) Q \left( - x \right) \right) ^ {- 1} = \left( x ^ n Q \left( x \right) \right) ^ {- 1}$ の $x ^ {- d + 1}$ から $x ^ 0$ までの係数が求まる。

これで、時間計算量が $\operatorname{\Theta} \left( d \log \left( d \right) \log \left( k + d \right) \right)$ のアルゴリズムが得られた。 $n \leq d$ になった時点で再帰を打ち切り FPS の逆元を計算するアルゴリズムを適用するようにすれば、 $\operatorname{\Theta} \left( d \log \left( d \right) \log \left( k / d + 2 \right) \right) \subseteq \operatorname{O} \left( d \log \left( d \right) \log \left( k \right) \right)$ とすることもできる。

また、詳細は省略するが、アルゴリズム中の様々な畳み込みは DFT の性質を利用して定数倍高速化できる。例えば、 $Q \left( x \right)$ の DFT が分かっているとき $Q \left( - x \right)$ の DFT を $\operatorname{\Theta}(d)$ の時間計算量で計算することができる。

noshi91のメモ

データ構造のある風景

線形漸化式を持つ数列の連続する項を高速に計算する

概要