線形漸化的数列のN項目の計算

この記事はデータ構造とアルゴリズム Advent Calendar 2020 18日目の記事です。

はじめに

与えられた線形漸化式を満たす数列の $N$ 項目を計算する問題は計算機科学において最も基本的な問題の一つである。この問題に対するアルゴリズムとして Fiduccia のアルゴリズム(1985) が知られている。Fiduccia のアルゴリズムは極めてシンプルかつ高速なアルゴリズムであり35年間最速のアルゴリズムであった。2021年に Bostan と Mori によって演算回数が Fiduccia のアルゴリズムのおよそ 2/3 であるアルゴリズムが発表された。本稿では Fiduccia のアルゴリズムと Bostan–Mori のアルゴリズムの解説をする。
また、Bostan–Mori のアルゴリズムを応用することで、多項式 $\Gamma(x)$ について $x^N\bmod\Gamma(x)$ を繰り返し二乗法の 2/3 倍の演算回数で計算することができる。このアルゴリズムについては「X^N mod P(X)の計算」で紹介している。
なお参考文献は一番下に掲載している。

線形漸化的数列

数列 $(a_n)_{n\ge 0}$ に対してある正の整数 $d$ が存在して線形漸化式

a_n = c_1 a_{n-1} + c_2 a_{n-2} + \dotsb c_d a_{n-d},\qquad \text{for}\quad n\ge d

を満たされるとき、数列 $(a_n )_{n \ge 0}$ は線形漸化的数列であるという。

線形漸化的数列は初項 $a_0, \cdots, a_{d-1}$ と線形漸化式の係数 $c_1, \cdots c_d$ を用いて定義することができる ($c_d \ne 0$ と仮定する)。ここで、正の整数 $d$ を線形漸化式の位数(order)と呼ぶ。線形漸化的数列が満たす線形漸化式は複数有り得るので、位数は線形漸化的数列ではなく線形漸化式に対して定義される。

数列に含まれる「数」は一般的に可換環 $\mathbf{R}$ の元とする。
線形漸化的数列 $N$ 項目の計算とは以下で定義される問題である。

• 入力: $N\in\mathbb{N}$, $a_0, \cdots, a_{d-1} \in \mathbf{R}$, $c_1,\dotsc, c_d\in\mathbf{R}$.
• 出力: 入力で定義される線形漸化的数列の $N$ 項目 $a_N$.

アルゴリズムの計算量

この問題を解くアルゴリズム計算量について以下の二つの尺度がある。

1. 算術計算量(arithmetic complexity): $\mathbf{R}$ 上の演算回数
2. ビット計算量(bit complexity): チューリングマシンで計算する場合の計算量

可換環 $\mathbf{R}$ が有限集合であるとき、$\mathbf{R}$ 上の演算は $O(1)$ のビット計算量となる。一方で $\mathbf{R}=\mathbb{Z}$ などの無限集合の場合は大きな数になる程 $\mathbf{R}$ 上の演算にも時間が掛かる。この場合、算術計算量とビット計算量は大きく異なる。もちろん、ビット計算量は $\mathbf{R}$ 上の演算以外の計算もカウントするという意味で大きく違う。本稿では算術計算量が小さいアルゴリズムについて考える。本稿で紹介するアルゴリズムは算術計算量が低いというだけでなく、$\mathbf{R}$ 上の演算のビット計算量が $O(1)$ であれば、そのビット計算量は算術計算量と同じオーダーになる。

$\mathbf{R}$ 上の $d$ 次多項式同士の乗算の算術計算量を $\mathsf{M}(d)$ と表す。通常の筆算より $\mathsf{M}(d)=O(d^2)$ であることが分かる。Karatsuba法により $\mathsf{M}(d)=O(d^{1.59})$ であることが分かる。また、Toom–Cook法により任意の $\epsilon >0$ について $\mathsf{M}(d)=O(d^{1+\epsilon})$ となることが知られている。現在知られている漸近的に最速のアルゴリズムは Schönhage と Strassen によるもので $\mathsf{M}(d) = O(d\log d\log\log d)$ が得られる。さらに $\mathbf{R}$ が特殊な可換環であるとき、高速フーリエ変換(FFT)のアルゴリズムを用いることで $\mathsf{M}(d)=O(d\log d)$ となる。

Fiduccia のアルゴリズムと Bostan–Mori のアルゴリズム

算術計算量が $O(\mathsf{M}(d)\log N)$ であるアルゴリズムとして

• Fiduccia のアルゴリズム (1985年): $\mathbf{3} \,\mathsf{M}(d)\lceil\log (N+1)\rceil + O(\mathsf{M}(d))$
• Bostan–Mori のアルゴリズム (2021年) $\mathbf{2} \,\mathsf{M}(d)\lceil\log (N+1)\rceil + O(\mathsf{M}(d))$

がある。一般的な可換環のもとで、Bostan–Mori のアルゴリズムは Fiduccia のアルゴリズムと比べて 2/3 の算術計算量となる。36年越しに計算量が(定数倍ではあるが)改善されたことになる。また Bostan–Mori のアルゴリズムはFFTと相性が良く、特殊な体のもとでは改善はより大きなものとなる。

また、任意に与えられた $d$ 次多項式 $\Gamma(x)\in\mathbf{R}[x]$ と $N\in\mathbb{N}$ について $x^N \bmod\Gamma(x)$ を計算する問題についても Bostan–Moriのアルゴリズムを適用することができ、素朴な $\mathbf{R}[x]/\Gamma(x)$ 上の繰り返し二乗法と比べて算術計算量が 2/3 になる。この問題には多項式の因数分解や楕円曲線上の点の個数の数え上げなどたくさんの応用があるため、実用的にはこちらの結果の方がインパクトが大きいと考えられる。このアルゴリズムについては「X^N mod P(X)の計算」で紹介している。

以下では Fiduccia のアルゴリズムと Bostan–Mori のアルゴリズムを説明する。

Fiduccia のアルゴリズム

フィボナッチ数列の N 項目

フィボナッチ数列 $(F_n\in\mathbb{Z})_{n\ge 0}$ は

\begin{align*}
F_0 &= 0,\qquad F_1 = 1\\
F_n &= F_{n-1} + F_{n-2}\qquad\text{for}\quad n\ge 2
\end{align*}

と定義される数列である。最初の11項は 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55 である。
この定義に従って $a_n$ を順次計算すると、$a_N$ は $N-1$ 回の加算で計算できる。より少ない算術計算量で計算するために以下の等式を考える。

\begin{bmatrix}
F_n\\
F_{n+1}
\end{bmatrix}
=
\begin{bmatrix}
0&1\\
1&1
\end{bmatrix}
\begin{bmatrix}
F_{n-1}\\
F_{n}
\end{bmatrix}.

この等式を繰り返し適用することにより以下を得る。

\begin{align*}
\begin{bmatrix}
F_n\\
F_{n+1}
\end{bmatrix}
&=
\begin{bmatrix}
0&1\\
1&1
\end{bmatrix}^n
\begin{bmatrix}
F_{0}\\
F_{1}
\end{bmatrix}\\
&=
\begin{bmatrix}
0&1\\
1&1
\end{bmatrix}^n
\begin{bmatrix}
0\\1
\end{bmatrix}.
\end{align*}

ここで$2\times 2$行列の $N$乗を繰り返し二乗法により計算することで $O(\log N)$ 回の加算と乗算で $F_N$ が計算できる。

行列の N 乗

前節のアルゴリズムを一般的な線形漸化的数列に適用すると

\begin{align*}
v &:=
\begin{bmatrix}
a_{0}\\
a_{1}\\
\vdots\\
a_{d-2}\\
a_{d-1}\\
\end{bmatrix},\qquad
e_0 :=
\begin{bmatrix}
1\\
0\\
\vdots\\
0\\
0\\
\end{bmatrix}\\
C &:=
\begin{bmatrix}
0&1&0&\cdots&0&0\\
0&0&1&\cdots&0&0\\
\vdots&\vdots&\ddots&\ddots&\vdots&\vdots\\
0&0&0&\ddots&1&0\\
0&0&0&\cdots&0&1\\
c_d&c_{d-1}&c_{d-2}&\cdots&c_2&c_1\\
\end{bmatrix}\\
a_N
&= \langle e_0, C^N v\rangle
\end{align*}

という等式が得られる。ここで $\langle\cdot,\cdot\rangle$ はベクトルの内積である。$d\times d$ 行列 $N$ 乗の計算に帰着されたので、$a_N$ を計算するアルゴリズムの算術計算量は $O(\mathsf{MM}(d)\log N)$ である。ここで、$\mathsf{MM}(d)$ は $d\times d$ 行列同士の積の算術計算量である。行列積の定義に従って計算すると $\mathsf{MM}(d)=O(d^3)$ である。Strassen のアルゴリズムなど、より算術計算量の低いアルゴリズムが知られているが、行列の要素数が $d^2$ であることから $\mathsf{MM}(d)=\Omega(d^2)$ である。

Cayley–Hamilton の定理

$\mathbf{R}$ 上の行列の $N$ 乗を $\mathbf{R}$ 上の多項式の $N$ 乗に置き換えることが一般的に可能である。
「可換環 $\mathbf{R}$ 上の行列 $A$ の特性多項式を $\Gamma(\lambda)$ とすると、$\Gamma(A)=0$ が成り立つ」というのが Cayley–Hamilton の定理である。
$x^N$ を $\Gamma(x)$ で割った商と余りをそれぞれ $q(x)$ と $r(x)$ とおくと、$x^N=q(x)\Gamma(x)+r(x)$ が成り立つ。よって $A^N = r(A)$ である。
従って一般的に $A^N$ を計算する方法として

1. $A$ の特性多項式 $\Gamma(\lambda)$ を計算する
2. $r(x) := x^N \bmod \Gamma(x)$ を計算する
3. $r(A)$ を計算し出力する

というアルゴリズムが考えられる。

※これは大変見通しのよいアルゴリズムであるが、一般の行列の $N$ 乗を計算する際には先にフロベニウス標準形を導出する方がより効率的である。

Fiduccia のアルゴリズム

計算したいものは $a_N=\langle e_0, C^Nv\rangle = \langle (C^T)^Ne_0, v\rangle$ である。ここで $C$ の転置 $C^T$ はコンパニオン行列と呼ばれる行列であり、簡単な計算より特性多項式が

\Gamma(x) = x^d - c_1 x^{d-1} - c_2 x^{d-2} - \dotsb - c_d

であることが分かる。そのため Cayley–Hamilton の定理に基づくアルゴリズムを $(C^T)^N$ に適用した場合、最初のステップである特性多項式の計算が必要ない。また、最終的に計算したいのは $(C^T)^Ne_0$ つまり $(C^T)^N$ の一列目であることから、最後のステップである $r(C^T)$ の計算についても簡略化できる。 $i=0,\cdots,d-1$ について $(C^T)^i$ の一列目は第 $i$ 成分が 1 でそれ以外が 0 になる。このことから $r(x) = r_0 + r_1 x + \dotsb r_{d-1}x^{d-1}$ に $C^T$ を代入した $r(C^T)$ の一列目は $[r_0\ r_1\ \dotsm\ r_{d-1}]^T$ である。よって $a_N = \sum_{i=0}^{d-1} a_i r_i$ となる。以下に Fiduccia のアルゴリズムを示す。

入力: $N\in\mathbb{N}$, $a_0, \cdots, a_{d-1} \in \mathbf{R}$, $c_1,\dotsc, c_d\in\mathbf{R}$
出力: 入力で定義される線形漸化的数列の $N$ 項目 $a_N$.
1. $\sum_{i=0}^{d-1} r_i x^i := x^N \bmod \Gamma(x)$ を計算する
2. $\sum_{i=0}^{d-1} a_i r_i$ を計算し出力する

また、Cayley–Hamilton の定理を持ち出さなくても、$C^T$ をベクトルに掛けた結果を考えると $\mathbf{R}[x]/\Gamma(x)$ で $x$ 倍することに対応していることが確認できる。

多項式環におけるモジュラー乗算

$\mathbf{R}[x]$ の次数 $d$ の元同士の乗算の算術計算量を $\mathsf{M}(d)$ とおいたのであった。次数 $d$ の多項式 $\Gamma(x)\in\mathbf{R}[x]$ について、ニュートン法に基づいて $\mathbf{R}[x]/\Gamma(x)$ 上の積を算術計算量 $6\mathsf{M}(d)+O(d)$ で計算するアルゴリズムが知られている。Fiduccia のアルゴリズムでは法となる多項式 $\Gamma(x)$ は固定されているため、算術計算量 $3\mathsf{M}(d)+O(d)$ の事前計算をしておくことで、$\mathbf{R}[x]/\Gamma(x)$ 上の積を一回あたり算術計算量 $3\mathsf{M}(d)+O(d)$ で計算することができる。ここでは詳しく説明しない。

二種類の繰り返し二乗法

さて $x^N\bmod\Gamma(x)$ の計算をするために $\mathbf{R}[x]/\Gamma(x)$ 上の繰り返し二乗法を使う訳だが、二種類の繰り返し二乗法のどちらを用いるかで計算量が定数倍変わる。
一般的にモノイド $M$ の元 $A$ に対する繰り返し二乗法は

A^N =\begin{cases}
e,& \text{if } N = 0\\
A(A^2)^{(N-1)/2}, & \text{if } N \text{ is odd}\\
(A^2)^{N/2}, & \text{otherwise}.
\end{cases}

という漸化式に基づくものと

A^N =\begin{cases}
e,& \text{if } N = 0\\
A(A^{(N-1)/2})^2, & \text{if } N \text{ is odd}\\
(A^{N/2})^2, & \text{otherwise}.
\end{cases}

という漸化式に基づくものがある。ここで $e$ はモノイド $M$ の単位元とする。前者のアルゴリズムは $N$ の二進数表現の各ビットを LSB (least significant bit) から順番に見ていくアルゴリズムであるので LSB-first アルゴリズムと呼ぶことにする。同様に後者のアルゴリズムは MSB (most significant bit) から順番に見ているため MSB-first アルゴリズムと呼ぶ。どちらのアルゴリズムも $N$ の各桁ごとに最大2回の乗算が発生する。このうち一回は一般的な元の二乗の計算であり、もう一回は $A$ を乗ずる計算である。ここで MSB-first アルゴリズムでは $A$ が再帰呼び出しの中で変化しないため、「$A$ を乗ずる」という計算の計算量が特別低い場合に高速化される。$\mathbf{R}[x]/\Gamma(x)$ 上で $x$ を乗ずるというのは算術計算量 $O(d)$ で計算できる。ここで $d$ は $\Gamma(x)$ の次数である。よって $x^N\bmod\Gamma(x)$ を計算する MSB-first アルゴリズムの算術計算量は $(3\mathsf{M}(d)+O(d))\lceil\log(N+1)\rceil + O(\mathsf{M}(d))$ である。Fiduccia のアルゴリズムの算術計算量も同様である。

Bostan–Mori のアルゴリズム

形式的べき級数

Bostan–Mori のアルゴリズムは数列の形式的べき級数表現(母関数)に基づく。$\mathbf{R}$ 上の形式的べき級数 $A(x)\in\mathbf{R}[[x]]$は

A(x) = \sum_{n\ge 0} a_n x^n = a_0 + a_1 x + a_2 x^2 + \dotsb

と表される。これは数列 $(a_n)$ と同一視して差し支えない。数列 $(a_n)$ から定まる形式的べき級数 $A(x)\in\mathbf{R}[[x]]$ を数列 $(a_n)_{n\ge 0}$ の母関数と呼ぶ。 例えばフィボナッチ数列の母関数は

x+x^2+2x^3+3x^4+5x^5+8x^6+13x^7+21x^8+34x^9+55x^{10}+\dotsb

である。

形式的べき級数には和と積が以下のように定義される。

\begin{align*}
\sum_{n\ge 0} a_n x^n + \sum_{n\ge 0} b_n x^n = \sum_{n\ge 0} (a_n + b_n) x^n\\
\left(\sum_{n\ge 0} a_n x^n\right)\cdot\left(\sum_{n\ge 0} b_n x^n\right) = \sum_{n\ge 0} \left(\sum_{m=0}^n a_mb_{n-m}\right) x^n.
\end{align*}

ここで、各演算の結果の $x^n$ の係数は有限和で表わされることに注意しよう。$\mathbf{R}$ 上の無限和や無限積は定義されない。

さらに

\begin{align*}
\frac1{1-x} := \sum_{n\ge 0} x^n
\end{align*}

と定義する。乗算の定義より $(1-x)\frac1{1-x} = 1$ となることが確認できる。 つまり $\frac1{1-x}$ は $1-x$ の乗法の逆元であり、この記法が適切であることが分かる。
同様に一般の形式的べき級数 $B(x)$ について $1 / (1 - xB(x)) := \sum_{n\ge 0} x^nB(x)^n$ と定義する。ここで $\sum_{n\ge 0} x^nB(x)^n$ の $x^k$ の係数は $\sum_{n=0}^k x^nB(x)^n$ の $x^k$ の係数と等しく、$B(x)$ の0次から $k$ 次の係数に対する有限回の演算で表されることが分かる。よって、これもまた well-defined である。

これらのことから定数項が 1 であるような形式的べき級数について乗法の逆元が与えられたわけであるが、同様に定数項が $\mathbf{R}$ の単元(乗法の逆元を持つ元)の場合にも自明に定義できる。さらに、「形式的べき級数 $A(x)$ が乗法の逆元を持つ $\iff$ $A(x)$ の定数項が単元」が成り立つ(自明)。

形式的べき級数 $A(x)$ における $x^N$ の係数を $[x^N]\, A(x)$ と書く。また、$A(0)$ で $A(x)$ の定数項を表す。

線形漸化的数列の母関数

フィボナッチ数列の母関数に $1-x-x^2$ を掛けると

\begin{align*}
&(1-x-x^2)(x+x^2+2x^3+3x^4+5x^5+\dotsb)\\
&=
x+(1-1)x^2+(2-1-1)x^3+(3-2-1)x^4 +(5-3-2)x^5+\dotsb\\
& = x
\end{align*}

となり、多項式(次数が有限の形式的べき級数)となる。このことから、フィボナッチ数列の母関数は $x/(1-x-x^2)$ であることが分かる。同様の命題が一般の線形漸化的数列で成り立つ。

数列 $(a_n)_{n\ge 0}$ とその母関数 $G_a(x)$ について以下が同値である。

1. 数列 $(a_n)_{n\ge 0}$ が位数 $d$ の線形漸化式を満たす
2. ある次数 $d$ の多項式 $Q(x)$ で $Q(0)=1$ を満たすものと次数 $d-1$ 以下の多項式 $P(x)$ が存在し、$G_{a}(x) = P(x)/Q(x)$

証明
1$\rightarrow$2
$Q(x) := 1-\sum_{i=1}^d c_i x^i$ と定義する。すると $Q(x)G_a(x)$ の $n\ge d$ 次の係数は

\begin{align*}
a_n - \sum_{i=1}^d c_i a_{n-i} = 0
\end{align*}

である。つまり $Q(x)G_{a}(x)$ は高々 $d-1$ 次の多項式 $P(x)$ になる。ここで $Q(0)=1$ であることから、$G_{a}(x)=P(x)/Q(x)$ を得る。

2$\rightarrow$1
$Q(x)G_a(x)$ が高々 $d-1$ 次の多項式になる。つまり、$n\ge d$ について $[x^n]\, Q(x)G_a(x)$ が 0 になる。このことから $(a_n)$ が位数 $d$ の線形漸化式を満たすことが分かる。

Bostan–Mori のアルゴリズム

目標は

a_N = [x^N]\, \frac{P(x)}{Q(x)}

を計算することである。ここで $Q(x)$ は $d$ 次多項式、$P(x)$ は高々 $d-1$ 次の多項式とする。また $Q(0)$ は $\mathbf{R}$ の単元であるとする。線形漸化的数列の文脈では $Q(0)=1$ としてよいのだが、一般的に単元としても以下のアルゴリズムには影響がない。$N=0$ のとき、解は $P(0)/Q(0)$ である。以下 $N\ge 1$ の場合を考える。
分母と分子に $Q(-x)$ を掛けることで

a_N = [x^N]\, \frac{P(x)Q(-x)}{Q(x)Q(-x)}

が得られる。ここで、$Q(x)Q(-x)$ は偶多項式であるので、ある $d$ 次多項式 $V(x)$ が存在して、$V(x^2) = Q(x)Q(-x)$ が成り立つ。
分子を偶数次部と奇数次部に分けることで $P(x)Q(-x)=U_\mathrm{e}(x^2) + x U_\mathrm{o}(x^2)$ と表す。

a_N = [x^N]\, \frac{U_\mathrm{e}(x^2) + x U_\mathrm{o}(x^2)}{V(x^2)}

また $1/V(x^2)$ は偶形式的べき級数である。よって

\begin{align*}
a_N &= \begin{cases}
[x^N]\, \frac{U_\mathrm{e}(x^2)}{V(x^2)}& \text{if $N$ is even}\\
[x^N]\, \frac{x U_\mathrm{o}(x^2)}{V(x^2)}& \text{otherwise.}
\end{cases}\\
&= \begin{cases}
[x^{N/2}]\, \frac{U_\mathrm{e}(x)}{V(x)}& \text{if $N$ is even}\\
[x^{(N-1)/2}]\, \frac{U_\mathrm{o}(x)}{V(x)}& \text{otherwise.}
\end{cases}
\end{align*}

が得られる。ここで $V(x)$ は $d$ 次で $V(0)$ は単元、$U_{\mathrm{e}}(x),\,U_{\mathrm{o}}(x)$ は高々 $d-1$ 次であることから、同じ変形を $N=0$ になるまで繰り返し使うことができる。これらの議論から Bostan–Mori のアルゴリズムが得られた。

R.<x> = QQ[]

def even(P):
  return R(P.list()[0::2])

def odd(P):
  return R(P.list()[1::2])

def bostan_mori(N, P, Q):
  while N > 0:
    U = P * Q(-x)
    if N % 2 == 0: P = even(U)
    else: P = odd(U)
    Q = even(Q * Q(-x))
    N //= 2
  return P(0) / Q(0)

print(bostan_mori(100, x, 1-x-x^2))

$N$ の各桁毎に2回の $d$ 次多項式の乗算が発生するので、Bostan–Mori のアルゴリズムの算術計算量は $2\mathsf{M}(d)\lceil \log(N+1)\rceil$ である。もとの線形漸化的数列の $N$ 項目を計算する問題を解く際には $P(x)$ を計算するために一回多項式乗算が必要なので、さらに $\mathsf{M}(d)$ 必要である。

フィボナッチ数列の場合

フィボナッチ数列の母関数は

\frac{x}{1-x-x^2}

である。 これに対して Bostan–Mori のアルゴリズムを適用すると

\begin{align*}
F_N=[x^N]\, \frac{x}{1-x-x^2}
&=
[x^N]\, \frac{x(1+x-x^2)}{(1-x-x^2)(1+x-x^2)}\\
&=
[x^N]\, \frac{x+x^2-x^3}{1-3x^2+x^4}\\
&=
\begin{cases}
[x^{N/2}]\, \frac{x}{1-3x+x^2}& \text{if $N$ is even}\\
[x^{(N-1)/2}]\, \frac{1-x}{1-3x+x^2}&\text{otherwise}
\end{cases}
\end{align*}

が得られる。 また、一般的に

\begin{align*}
[x^N]\, \frac{a+bx}{1-cx+x^2}
&=
[x^N]\, \frac{(a+bx)(1+cx+x^2)}{(1-cx+x^2)(1+cx+x^2)}\\
&=
[x^N]\, \frac{a+(b+ac)x+(a+bc)x^2+bx^3}{1-(c^2-2)x^2+x^4}\\
&=
\begin{cases}
[x^{N/2}]\, \frac{a+(a+bc)x}{1-(c^2-2)x+x^2}& \text{if $N$ is even}\\
[x^{(N-1)/2}]\, \frac{b+ac+bx}{1-(c^2-2)x+x^2}&\text{otherwise}
\end{cases}
\end{align*}

が成り立つ。よって次のようなアルゴリズムが得られる。

def fib(n):
  if n == 0: return 0
  if n % 2 == 0: [a, b] = [0, 1]
  else: [a, b] = [1, -1]
  c = 3
  n //= 2
  while n != 0:
    if n % 2 == 0: b = a + b * c
    else: a = b + a * c
    c = c * c - 2
    n //= 2
  return a

最後の無駄な乗算はビット計算量の意味で無視できないのでこれを取り除き、 さらに $[x^{N-1}]\, 1/(1-x-x^2)$ から始めることにすると

def fib(n):
  if n <= 2: return int(n != 0)
  if n % 2 == 0: [a, b] = [1, 0]
  else: [a, b] = [1, -1]
  c = 3
  n = (n - 1) // 2
  while n != 1:
    if n % 2 == 0: b = a + b * c
    else: a = b + a * c
    c = c * c - 2
    n //= 2
  return b + a * c

が得られる。このアルゴリズムは現在実際に GMP で使われているアルゴリズム https://gmplib.org/manual/Fibonacci-Numbers-Algorithm と同様に $N$ の桁毎に 2 回の整数乗算を含む。フィボナッチ数は指数関数的に大きくなるため、最後の数回の乗算がプログラムの実行時間のほとんどを占めており、算術計算量は実行時間にはあまり関係がない。上記のアルゴリズムを GMP で素朴に実装したとしても GMP の組み込み関数と比べて遜色ない実行速度になる。

フィボナッチ数とリュカ数

このフィボナッチ数を計算するアルゴリズムが何をやっているのかもう少し深く理解してみよう。簡単のため $n$ が2のべきであると仮定すると、アルゴリズムは以下のようになる。

def fib(n):
  if n <= 2: return int(n != 0)
  a = 1
  c = 3
  n = (n - 1) // 2
  while n != 1:
    a = a * c
    c = c * c - 2
    n //= 2
  return a * c

ここで $c$ は 3 で初期化された後 $c \leftarrow c^2-2$ という更新式に従って更新される。これは2べき番目のリュカ数である。
リュカ数列は $L_0=2,\,L_1=1$ と $n\ge 2$ について $L_n=L_{n-1}+L_{n-2}$ で定義される。$\phi_+:=(1+\sqrt{5})/2,\,\phi_-:=(1-\sqrt{5})/2$ と定義すると、Binet の公式

\begin{align*}
F_n &= \frac{\phi_+^n-{\phi^n_{-}}}{\phi_+-\phi_-}&
L_n &= \phi_+^n+{\phi^n_-}
\end{align*}

が成り立つ。 よって

\begin{align*}
F_{2n} / F_n &= \frac{\phi_+^{2n}-{\phi^{2n}_{-}}}{\phi_+^n-{\phi^n_{-}}} = \phi_+^n+{\phi^n_-} = L_n\\
L_{2n} &= \phi_+^{2n}+{\phi^{2n}_-} = (\phi_+^n+{\phi^n_-})^2 - 2 (\phi_+\phi_-)^n = L_n^2 - 2 (-1)^n
\end{align*}

が成り立つ。一つ目の式より $F_{2^k} = \prod_{i=0}^{k-1} L_{2^i}$ が得られる。 また、二つ目の式より $L_{2^k} = L_{2^{k-1}}^2 -2$ が $k\ge 2$ について得られる。 つまり、$L_2=3$ から始めて、2べき番目のリュカ数を計算し、それらの積を取ることで $2^k$ 番目のフィボナッチ数を計算するアルゴリズムであると理解できる。
この $2^k$ 番目のフィボナッチ数に対する Bostan–Mori のアルゴリズムは Cull と Holloway によって ``Product of Lucas numbers'' と呼ばれているアルゴリズムに一致する。一方で $n$ が2べきとは限らない場合、Cull と Holloway はより複雑なアルゴリズムを考えている。また、Knuth は

\begin{align*}
F_{2n} &= F_n L_n&
L_{2n} &= L_n^2 - 2(-1)^2\\
F_{n+1} &= (F_n+L_n)/2&
L_{n+1} &= (5F_n+L_n)/2
\end{align*}

という漸化式に基づく MSB-first なアルゴリズムを提案している。
一般の $n$ に対する Bostan–Mori のアルゴリズムも フィボナッチ数とリュカ数に関する等式で理解できる。Binet の公式より

\begin{align*}
F_{n+m} &= F_nL_m + (-1)^nF_{m-n}
\end{align*}

が得られる。特に $m=2^k$, $k\ge {\lceil \log (n+1)\rceil}$ のとき、

\begin{align*}
F_{n+2^k} &= F_nL_{2^k} + (-1)^nF_{2^k-n}\\
(-1)^nF_{2^{k+1}-n} &= (-1)^nF_{2^{k}-n}L_{2^{k}} + F_n
\end{align*}

が得られる。よって、上記のアルゴリズムは $a = F_n,\, b = (-1)^n F_{2^k-n}$ を $N$ のLSBから順番に読んで更新していくアルゴリズムだと理解できる。

Fiduccia のアルゴリズムの方が Bostan–Mori のアルゴリズムより高速な場合

$\Gamma(x)$ による剰余の計算が特別高速にできる場合、Fiduccia のアルゴリズムの方が高速になることがある。例えば $\Gamma(x)$ が $O(1)$ 個しか非零係数を含まない場合や $\Gamma(x)=x^d-\sum_{n=0}^{d-1}x^n$ の場合(一般化フィボナッチ数列に相当)がある。これらの場合、剰余の計算は算術計算量 $O(d)$ でできるので、Fiduccia のアルゴリズムの算術計算量は $(\mathsf{M}(d)+O(d))\lceil \log(N+1)\rceil$ となる。Bostan–Mori のアルゴリズムでこれらの場合に高速化できるかどうかは未解決である。

MSB-first の Bostan–Mori のアルゴリズム

Bostan–Mori のアルゴリズムは $N$ をLSBから順番に見ていくアルゴリズムである。繰り返し二乗法のように Bostan–Mori のアルゴリズムにも MSB-first のバージョンが存在する。これを用いることにより $d$ 次モニック多項式 $\Gamma(x)$ について $x^n \bmod \Gamma(x)$ を算術計算量 $2\mathsf{M}(d)\lceil\log(N+1)\rceil + \mathsf{M}(d)$ で計算できる。このアルゴリズムについては「X^N mod P(X)の計算」で紹介している。実用的には MSB-first のアルゴリズムの方がインパクトが大きいと考えられる。
例えば $N$ 番目と $N+1$ 番目のフィボナッチ数を計算する MSB-first の Bostan–Mori アルゴリズムは以下のようになる。

def fib_msb(n):
  def aux(n, c):
    if n == 0: return [0, 1]
    [a, b] = aux(n//2, c*c-2)
    if n % 2 == 0: return [a * c, a + b]
    else: return [a + b, b * c]
  if n == 0: return [0, 1]
  [a, b] = aux(n//2, 3)
  if n % 2 == 0: return [a, b - a]
  else: return [b - a, b]

Graeffe 法

$Q(x)Q(-x)$ により偶多項式を作り出す方法は多項式の根を見つけるアルゴリズムである Graeffe 法で用いられている。

本来の Graeffe 法の目的は多項式の根の大きさを分離することであったが、Bostan–Mori のアルゴリズムでは純粋に代数的な目的で Graeffe プロセスを使っている。

また、Schönhage は Graeffe プロセスを用いて、多項式の $\mathbf{R}[[x]]$ における逆元を $N$ 次まで計算するアルゴリズムを考案した。

\begin{align*}
\frac1{Q(x)} \bmod x^{N+1} &= Q(-x) \frac1{Q(x)Q(-x)} \bmod x^{N+1} \\
&= Q(-x) \frac1{V(x^2)} \bmod x^{N+1}\\
&= Q(-x) S(x^2) \bmod x^{N+1} 
\end{align*}

が成り立つ。ここで、 $S(x) := 1/V(x)\bmod x^{\lfloor N/2\rfloor+1}$ である。
よって以下のアルゴリズムが得られた。

R.<x> = QQ[]

def even(P):
  return R(P.list()[0::2])

def inv(N, Q):
  if N == 0: return R(1 / Q(0))
  V = even(Q.multiplication_trunc(Q(-x), N+1))
  S = inv(N // 2, V)
  return Q(-x).multiplication_trunc(S(x^2), N+1)

print(inv(100, 1-x-x^2))

ただし、ニュートン法を用いたアルゴリズムの方が定数倍効率的であるため、このアルゴリズムは通常用いられない。
また、詳細な説明はしないが、 $1/Q(x)$ を $N>d$ 次まで求めたいときは、 $d$ 次までニュートン法や Schönhage の方法で求めた後に、 それを $d$ 次ずつ伸ばしていく方がより算術計算量が低い。

歴史的経緯

行列積を用いたアルゴリズムは Miller と Spencer Brown によって示された。

Krishnamurthy は実多項式 $\Gamma(x)$ の根を求める問題を対応するコンパニオン行列 $C$ の固有値を求める問題とみなし、$C$ の最大固有値をべき乗法でもとめる方法を考えた。そして、Cayley–Hamilton の定理より $r(x) = x^N \bmod \Gamma(x)$ を計算し、$r(C) = C^N$ から $C$ の最大固有値、つまり $\Gamma(x)$ の最大の根をもとめることを考えた。

Fiduccia のアルゴリズムは 1982年に Allerton Conference on Communication, Control and Computing という国際会議で発表され、1985年に SIAM Journal on Computing で論文が発表された。Allerton の論文の情報を得ることは難しく、インターネット上では唯一日本語のページ https://jglobal.jst.go.jp/detail?JGLOBAL_ID=200902049211196468 が見つかるだけである。論文タイトルで検索してもこのページしかヒットしない。また、Knuth によると Richard P. Brent が遅くとも 1981年には同様の結果を得ていたようである。

Bostan–Mori のアルゴリズムは私が 2018年の8月頃に Codechef の Editorial https://discuss.codechef.com/t/rng-editorial から思いついた。当初は論文にするつもりは無かったが、2020年の6月頃に INRIA の Alin Bostan 氏にこのアルゴリズムの新規性と有用性について尋ねたところ共同研究することになり共著で論文を書いた。私一人ではこのような素晴しいイントロは書けなかったし、結果の応用先についてもたくさん教えてもらった。MSB-first のアルゴリズムは彼との議論の中で生まれたものである。

ここからおまけ

代数的べき級数

有理べき級数より一般的な形式的べき級数のクラスに代数的べき級数がある。形式的べき級数 $A(x)$ が代数的であるとは、ある非零の二変数多項式 $E(x,y)\in\mathbf{R}[x,y]$ が存在し $E(x, A(x))=0$ が成り立つことである。特に $E$ の $y$ に関する次数が 1 であるとき、$A(x)$ は有理べき級数であることが分かる。例えば $a_n=\binom{2n}{n}$ の母関数 $A(x)$ は $(1-4x)A(x)^2 - 1 = 0$ を満たすので代数的である。

一般的に代数的べき級数よりも広い形式的べき級数のクラスである D-finite (p-recursive な数列の母関数) な形式的べき級数について、算術計算量 $O(\mathsf{M}(\sqrt{N})\log N)$ で $N$ 項目を計算する baby-step giant-step アルゴリズムが知られている。代数的べき級数に限ってもこれより算術計算量の低いアルゴリズムは知られていない。もし、一般的な可換環もしくは体で $\binom{2N}{N}$ を算術計算量 $o(\sqrt{N})$ で計算するアルゴリズムを発見したら大発見である。

Furstenberg の定理

二変数形式的べき級数から一変数形式的べき級数への写像 $\mathcal{D}$ を $\mathcal{D}(\sum_{n,m\ge 0}a_{n,m}x^ny^m) := \sum_{n\ge 0} a_{n,n}x^n$ と定義する。
Furstenberg は有限体上の代数的べき級数について以下が同値であることを示した。

1. $A(x)$ が代数的
2. ある二変数多項式 $P(x,y)$, $Q(x,y)$ が存在し、$A(x)=\mathcal{D}(P(x,y)/Q(x,y))$

例えば

\binom{2N}{N} = [x^Ny^N]\, \frac1{1-x-y}

という表現がある(これは有限体に限らず任意の可換環で成り立つ)。

Furstenberg は多項式 $P(x,y)$ と $Q(x,y)$ を得る具体的な方法も示した。代数的べき級数 $A(x)$ が $A(0)=0$ を満たし、$E_y(0, 0)\ne 0$ であると仮定する($E_y$ は $E$ の $y$ に関する偏微分である)。このとき

A(x) = \mathcal{D}\left(\frac{y E_y(xy, y)}{\frac1y E(xy,y)}\right)

が成り立つ。ここで $1/(\frac1y E(xy, y))$ が well-defined であることを確認しておこう。$A(0)=0$ より $E(0, 0)=0$ である。よって $E(xy, y)$ は $y$ で割ることができる。また、$E_y(0, 0)\ne 0$ より $E(x, y)$ の $y$ の係数の定数項は非零である。つまり $\frac1y E(xy,y)$ の定数項は非零であり、上記の形式的べき級数は well-defined である。 上記の表現の証明は Furstenberg の論文を参照(そこまで難しくない)。

例として $A(x) = \sum_{n\ge 0}\binom{2n}{n} x^n$ について考える。 $A(0)=1$ なので、上記の定理をそのまま適用することはできない。$B(x):=A(x)-1$とすると、$B(0)=0$ であり、$(1-4x)(1+B(x))^2-1=0$ を満たす。$E(x, y):=(1-4x)(1+y)^2-1$ について、$E_y(x, y)=2(1-4x)(1+y)$ である。そのため $E_y(0,0)=2$ であり上記の定理を適用できる。 すると

\begin{align*}
B(x) &= \mathcal{D}\left(\frac{y 2(1-4xy)(1+y)}{\frac1y ((1-4xy)(1+y)^2-1)}\right)\\
&= \mathcal{D}\left(\frac{2 y (1-4xy)(1+y)}{2-4x+y-8xy-4xy^2}\right)
\end{align*}

という表現が得られる。右辺の有理べき級数を係数が有理数であるとして(有限体ではないが)展開すると

y + \boldsymbol{2 x y} + \frac{1}{2} y^{2} + 4 x^{2} y - \frac{1}{4} y^{3} + 8 x^{3} y + \boldsymbol{6 x^{2} y^{2}} - \frac{1}{2} x y^{3} + \frac{1}{8} y^{4} + 16 x^{4} y + 24 x^{3} y^{2} + \frac{1}{2} x y^{4} - \frac{1}{16} y^{5} + 32 x^{5} y + 72 x^{4} y^{2} + \boldsymbol{20 x^{3} y^{3}} - x^{2} y^{4} - \frac{3}{8} x y^{5} + \frac{1}{32} y^{6} + 64 x^{6} y + 192 x^{5} y^{2} + 116 x^{4} y^{3} + \frac{3}{4} x^{2} y^{5} + \frac{1}{4} x y^{6} - \frac{1}{64} y^{7} + 128 x^{7} y + 480 x^{6} y^{2} + 456 x^{5} y^{3} + \boldsymbol{70 x^{4} y^{4}} - \frac{5}{2} x^{3} y^{5} - \frac{3}{8} x^{2} y^{6} - \frac{5}{32} x y^{7} + \frac{1}{128} y^{8} + 256 x^{8} y + 1152 x^{7} y^{2} + 1504 x^{6} y^{3} + 520 x^{5} y^{4} + \frac{3}{2} x^{3} y^{6} + \frac{1}{8} x^{2} y^{7} + \frac{3}{32} x y^{8} - \frac{1}{256} y^{9} + 512 x^{9} y + 2688 x^{8} y^{2} + 4480 x^{7} y^{3} + 2496 x^{6} y^{4} + \boldsymbol{252 x^{5} y^{5}} - 7 x^{4} y^{6} - \frac{1}{2} x^{3} y^{7} - \frac{7}{128} x y^{9} + \frac{1}{512} y^{10} + O(x, y)^{11}

となり、対角成分が $\binom{2n}{n}$ になっていることが確認できる。整数でない係数が現れるのが気にくわないが $x\mapsto x/2,\,y\mapsto 2y$ と置き換えると係数はすべて整数になる。
$B(x) = \mathcal{D}(1/(1-x-y) -1)$ という表現の方が分母と分子の次数が低いので好ましいが、Furstenberg の方法でこれを得ることはできないようだ。

有限体上の多変数有理べき級数

有限体 $\mathbb{F}_p$ 上の代数的べき級数の $N$ 項目を低い算術計算量で計算する Bostan–Christol–Dumas のアルゴリズムを紹介する(より算術計算量が低いアルゴリズムも知られている)。
目標は

[x^Ny^N]\, \frac{P(x,y)}{Q(x,y)}

を計算することである。ここで一変数のときのように Bostan–Mori のアルゴリズムを適用しても次数が保存されないという問題がある。$Q(x,y)^{p-1}$を分母と分子に掛けることで

[x^Ny^N]\, \frac{P(x,y)Q(x,y)^{p-1}}{Q(x,y)^p}

が得られる。フロベニウス写像の性質より、これは

[x^Ny^N]\, \frac{P(x,y)Q(x,y)^{p-1}}{Q(x^p,y^p)}

に等しい。よって、$k=N\mod p$ の値に従って、分子の $ip + k$ 次だけ ($i=0,\dotsc,d-1$)を残すと Bostan–Mori のアルゴリズムと同様に算術計算量 $O(\mathsf{M}(p^2d_xd_y)\log N)$ のアルゴリズムが得られる。ここで $d_x$ と $d_y$ はそれぞれ $Q(x,y)$ の $x$ と $y$ に関する次数である。このアルゴリズムにおいて分母 $Q(x,y)$ は不変であることから、前計算をして分子の更新のための行列をあらかじめ作っておくと、算術計算量は $O(d_x^2d_y^2\log N + p^2d_xd_y)$ となる。さらなる工夫で算術計算量を減らすことができる。

参考文献

1. C. M. Fiduccia, An efficient formula for linear recurrences, SIAM Journal on Computing, 14(1), pp. 106–112, 1985.
2. A. Bostan and R. Mori, A simple and fast algorithm for computing the $N$-th term of a linearly recurrent sequence, In Proceedings of SIAM Symposium on Simplicity in Algorithms, pp. 118–132, 2021.
3. J. C. P. Miller and D. J. Spencer Brown, An algorithm for evaluation of remote terms in a linear recurrence sequence, The Computer Journal, 9(2), pp. 188–190, 1966.
4. P. Cull and J. L. Holloway, Computing Fibonacci numbers quickly, Information Processing Letters, 32(3), pp. 143–149, 1989.
5. D. E. Knuth, Art of Computer Programming, Volume 2: Seminumerical Algorithms, 3rd edition, Addison-Wesley Professional, 1997.
6. A. Schönhage, Variations on computing reciprocals of power series, Information Processing Letters, 74(1-2), 41-46, 2000.
7. E. V. Krishnamurthy, Solving an algebraic equation by determining high powers of an associated matrix using the Cayley–Hamilton theorem, The Quarterly Journal of Mechanics and Applied Mathematics, 13(4), pp. 508–512, 1960. <!-- 1. A. Bostan, P. Gaudry, and E. Schost, Linear recurrences with polynomial coefficients and application to integer factorization and Cartier–Manin operator, SIAM Journal on Computing, 36(6), pp. 1777–1806, 2007.-->
8. D. V. Chudnovsky and G. V. Chudnovsky, Approximations and complex multiplication according to Ramanujan, In Ramanujan revisited (Urbanahampaign, Ill., 1987), pp. 375–472, Academic Press, Boston, MA, 1988.
9. H. Furstenberg, Algebraic functions over finite fields, Journal of Algebra, 7(2), pp. 271–277, 1967.
10. A. Bostan, G. Christol, and P. Dumas, Fast computation of the $N$th term of an algebraic series over a finite prime field, In Proceedings of the ACM on International Symposium on Symbolic and Algebraic Computation, pp. 119–126, 2016.
11. A. Bostan, X. Caruso, G. Christol, and P. Dumas, Fast coefficient computation for algebraic power series in positive characteristic, The Open Book Series, 2(1), pp. 119–135, 2019.