ブログのとさか

技術的な話をしたりしなかったり

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

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

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

イントロ

m項間漸化式の第nO (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

粗い測定ですが、mが小さいときから大きいときまで、一貫して総和の漸化式+kitamasa法が高速であることが確認できました。

一応gistにもソースコードを上げておきます。ブログの説明とは細かい添字が違ったりするので注意してください。

https://gist.github.com/tosaka2/e1b4e9dcb7892568d16570075d85941a

余談 O (m + n)

c_{i}=m - iである場合、総和の漸化式への変換を使ってm項間漸化式の第n項をO (m + n)で求めることができます。

これに限らず係数がある条件を満たす場合は、総和を経由することにより、第n項をO (m + n)で求められます。*8

このアルゴリズムm \leq 10 ^ 5, n \leq 10 ^ 6のように、 mが大きく、 m nのスケールが近いとき、kitamasa法を使うよりも高速な解法になり得ます。

所感

コンテストでこれが必要になることは無いと思うので、まぁ小ネタです。
息抜きに書くつもりがはてなブログの数式の仕様と格闘していたら思いの外時間がかかってしまいました...

何か間違い等あればご指摘ください。

*1:この書き方だとm+1項間漸化式な気もしますが、蟻本に合わせてこの式をm項間漸化式と呼ぶことにします。

*2:プログラミングコンテストチャレンジブック 第2版 P.180~

*3:コンパニオン行列のべき乗

*4:高速kitamasa法を用いるとO (m log m log n)まで計算量が落ちます

*5:はてなで上手く数式がかけないのでこの書き方になりました...

*6:c'_0, c'_1, \ldots, c'_{m}s_0 , s_1 , \ldots , s_{m}を前処理で計算する必要がありますが、これはどちらもO (m)で計算できます

*7:特に総和の漸化式+kitamasa法では係数が負になるので、元の数列の係数が全て正でも負のmodの処理をしてあげる必要があります

*8:総和に変換→ O(1)で次の項が求められる式に変形(これができることが条件)→s_0からs_{n}まで順に求める→a_{n}=s_{n}-s_{n -1}n項目が求まる