m項間漸化式の第n項までの和を$O(m ^ 2 log n)$で

この記事では $m$ 項間漸化式の第 $n$ 項までの和を $O (m ^ 2 log n)$ で求める方法について説明します。

@mt_caret がnth Fibonacci number in O(logn)という記事を書いていたのを見て、以前ブログに書こうと思っていて完全に忘却していたネタを思い出したので書きました。

イントロ

$m$ 項間漸化式の第 $n$ 項 $O (m ^ 3 log n)$

$m$ 項間の線形漸化式*1
${ \displaystyle a_{n+m} = \sum_{i=0}^{m-1} c_i a_{n+i} }$

は、次のように行列で表現し、繰り返し二乗法を使って $n$ 乗を計算することで $O (m ^ 3 log n)$ で求めることができます。

$$A= \begin{pmatrix} c_{m -1} & \cdots & c_{1} & c_{0} \\ 1 & \cdots & 0 & 0 \\ \vdots & \ddots & \vdots & \vdots \\ 0 & \cdots & 1 & 0 \end{pmatrix} $$

$$ \begin{pmatrix} a_{n+m -1} \\ a_{n+m -2} \\ \vdots \\ a_{n} \\ \end{pmatrix} =A \begin{pmatrix} a_{n+m -2} \\ a_{n+m -3} \\ \vdots \\ a_{n -1} \\ \end{pmatrix} =A ^ n \begin{pmatrix} a_{m -1} \\ a_{m -2} \\ \vdots \\ a_{0} \\ \end{pmatrix} $$

このことは、はじめに紹介した記事や蟻本*2で説明されています。

$m$ 項間漸化式の第 $n$ 項までの和 $O (m ^ 3 log n)$

似たような方法で第 $n$ 項までの和を求めることもできます。

$$ S_i = I + A + \ldots + A ^ {i -1} $$

とすると、ブロック行列を用いて

$$ \begin{pmatrix} A ^ n \\ S_n \\ \end{pmatrix} =\begin{pmatrix} A & 0 \\ I & I \\ \end{pmatrix}\begin{pmatrix} A ^ {n -1} \\ S_{n -1} \\ \end{pmatrix} =\begin{pmatrix} A & 0 \\ I & I \\ \end{pmatrix} ^ n \begin{pmatrix} I \\ 0 \\ \end{pmatrix} $$

と書くことができるので、繰り返し二乗法を適用すると $O (m ^ 3 log n)$ で $S_{n+1}$ が求められます。

後は $S_{n+1}$ に $\begin{pmatrix} a_{m-1}, a_{m-2} , \ldots , a_0 \end{pmatrix} ^ {\mathrm{T}}$ を掛ければ、第 $n$ 項までの和 $a_{0} + a_{1} + \ldots + a_{n}$ が求まります。

$$ S_{n+1} \begin{pmatrix} a_{m -1} \\ a_{m -2} \\ \vdots \\ a_{0} \\ \end{pmatrix} = ( I + A + \ldots + A ^ {i -1} + A ^ i ) \begin{pmatrix} a_{m -1} \\ a_{m -2} \\ \vdots \\ a_{0} \\ \end{pmatrix} $$

$$ = \begin{pmatrix} a_{m -1} \\ a_{m -2} \\ \vdots \\ a_{0} \\ \end{pmatrix} + \begin{pmatrix} a_{m} \\ a_{m -1} \\ \vdots \\ a_{1} \\ \end{pmatrix} + \begin{pmatrix} a_{n+m -1} \\ a_{n+m -2} \\ \vdots \\ a_{n} \\ \end{pmatrix} $$

$$ = \begin{pmatrix} a_{m -1} + a_{m -2} + \ldots + a_{n+m -1} \\ a_{m -2} + a_{m -3} + \ldots + a_{n+m -2} \\ \vdots \\ a_{0} + a_{1} + \ldots + a_{n} \\ \end{pmatrix} $$

このことは蟻本で行列の累乗和を求める方法として説明されています。

kitamasa法による $m$ 項間漸化式の第 $n$ 項の計算 $O (m ^ 2 log n)$

$m$ 項間漸化式の第 $n$ 項はより高速に計算する方法があります。
日本の競技プログラミング界隈ではkitamasa法と呼ばれます。*3

ここでは詳しい説明はしません。
気になる人はm項間漸化式の高速なアルゴリズム - 競技プログラミングをするんだよや高速 Kitamasa 法 - みさわめも等のブログ記事を読んで下さい。

ここで重要な点は、kitamasa法を用いると $m$ 項間漸化式の第 $n$ 項が $O (m ^ 2 log n)$ で計算できるという点です。*4

本題

$m$ 項間漸化式の第 $n$ 項までの和 $O (m ^ 2 log n)$

$m$ 項間漸化式の第 $n$ 項までの和は、行列の累乗和を陽に求めなくても計算できます。

アイデアは簡単で、第 $n$ 項までの総和を

${ \displaystyle s_{n} = \sum_{i=0}^{n} a_{n} }$

として、 $m$ 項間漸化式から、総和の漸化式

${ \displaystyle s_{n+m'} = \sum_{i=0}^{m'} c'_i s_{n+i} }$

を構成することを考えます。

結論だけ言うと、次のような関係になります。*5

$m'=m+1$ つまり $m+1$ 項間漸化式 $s_{n+m+1}$ について、

$c'_0=-c_0$
$c'_{m}=2 c_{m-1}$
$c'_i=c_i - c_{i+1}$

この関係は ${ \displaystyle a_{n+m} = \sum_{i=0}^{m-1} c_i a_{n+i} }$ と $s_{n+1}=s_{n}+a_{n+1}$ から導くことができるので、計算してみてください。

これにより、 $m$ 項間漸化式の第 $n$ 項までの和を求める問題を、 $m+1$ 項間漸化式の第 $n$ 項を求める問題に変換できました。

すでに説明したとおり、 $m$ 項間漸化式の第 $n$ 項はkitamasa法により $O (m ^ 2 log n)$ で求めることができるので、 $m$ 項間漸化式の第 $n$ 項までの和も $O (m ^ 2 log n)$ で求められるようになりました。*6

実験

$m$ 項間漸化式の第 $n$ 項までの和を求める計算を以下の3つの方法で行い、実行速度を比較しました。

前から順番に求めて足していく $O (m n)$
行列の累乗和 $O (m ^ 3 log n)$
総和の漸化式＋kitamasa法 $O (m ^ 2 log n)$

言語はC++です。答えは非常に大きくなりうるので $10 ^ 9 + 7$ でmodを取っています*7。その他の詳細な設定はリンク先で確認してください。

kitamasa法のソースコードはこちらからお借りしました。
累乗和の実装は蟻本を参考にしています。

結果は以下の表のようになりました。

	$$m=3,n=10 ^ 6$$	$$m=10,n=10 ^ 6$$	$$m=100,n=10 ^ 8$$
前から足していく	21.5772ms	75.9941 ms	-（時間がかかりすぎるので省略）
行列の累乗和	0.024ms	0.6107 ms	1764.7 ms
総和の漸化式+kitamasa	0.0153ms	0.0991 ms	12.0769 ms
詳細	https://wandbox.org/permlink/uVmogwgBVcEINwYf	https://wandbox.org/permlink/6HbnojdzEU61JMgP	https://wandbox.org/permlink/GYXVyPbczVkOEen4