はじめに

この解説は線形漸化的数列のN項目の計算の続きである。記号の定義等は引き継ぐ。ここでは MSB-first な Bostan–Mori のアルゴリズムを紹介する。その応用として $d$ 次多項式 $\Gamma(x)$ について $x^N \bmod \Gamma(x)$ が算術計算量 $2\mathsf{M}(d)\lceil\log(N+1)\rceil + \mathsf{M}(d)$ で計算できること示す。素朴な繰り返し二乗法の算術計算量 $(3\mathsf{M}(d)+O(d))\lceil\log(N+1)\rceil + 3\mathsf{M}(d)$ の 2/3 倍以下の算術計算量となっている。繰り返し二乗法が改善できるなんて驚きではないだろうか？
元の Bostan–Mori アルゴリズムと転置の関係があることを以下では説明しているが、アルゴリズムそのものを理解するためには直接的な導出から読めば十分である。

行列の転置とアルゴリズムの転置

$\mathbf{R}$ 上のベクトル $x$ を入力としてとり、$\mathbf{R}$ 上の行列 $A$ を $x$ に掛ける計算について考える。この計算について算術的なアルゴリズムを考える。ここで考える算術的なアルゴリズムとはベクトルの要素には $\mathbf{R}$ 上の演算のみを通じてアクセスするようなアルゴリズムである。ベクトルの要素を比較した結果で条件分岐をしたり、要素の1ビットだけを読んだりすることはできない。
算術的なアルゴリズムは $\mathbf{R}$ の重み付き有向非巡回グラフを使って表すことができる。入力ベクトルの要素は入次数0の頂点(入力ノード)、出力ベクトルの要素は出次数0の頂点(出力ノード)に対応する。入力ノードと出力ノード以外のノードは入次数も出次数も正であるとする。各有向辺には $\mathbf{R}$ の元が重みとして与えられている。入力ベクトルの要素は有向辺を通って接続するノードに渡される。$\mathbf{R}$ の元は有向辺を通るとき、その重みが掛けられる。入次数も出次数も正のノードでは入力辺によって入力された値の和を出力辺を通して他のノードに渡す。出力ノードは入力辺を通じて入力された値の和を計算結果として保持する。このようにして算術的なアルゴリズムを重み付き有向非巡回グラフで表す。

転置原理 (Tellegen の原理)
$A$ を $m\times n$ 行列ですべて 0 の行や列を含まないものとする。$A$ をベクトルに作用する算術計算量 $C$ のアルゴリズムがあるとき、$A^T$ をベクトルに作用する算術計算量 $C-n+m$ のアルゴリズムが存在する。

証明
$A$ をベクトルに作用する算術計算量 $C$ のアルゴリズムを表す重み付き有向非巡回グラフを $G$ とおく。$G$ の有向辺の向きを逆にしたグラフを $G^T$ と表す。このとき、$G^T$ は $A^T$ をベクトルに作用するアルゴリズムを表している。このことを確認するために $A$ の $i,\,j$ 成分と $G$ の関係について考える。
$A$ の $i,\,j$ 成分は

A_{i,j} = \sum_{p\colon\, i\rightsquigarrow j} \prod_{e\in p} w_e

と表せる。ここで $p$ は $G$ における入力ノード $i$ から出力ノード $j$ へのパス全体から選ばれる。また、$w_e$ は有向辺 $e$ の重みである。
このことから、 $G^T$ が表す行列の $j,\,i$ 成分は $G$ が表す行列の $i,\,j$ 成分と等しいことが分かる。よって、 $G^T$ は $A^T$ をベクトルに作用するアルゴリズムを表している。

これらのアルゴリズムの算術計算量について考える。有向グラフが表すアルゴリズムにおいて $\mathbf{R}$ の乗算回数は有向辺のうち重みが $\pm 1$ でないものの本数である。また $\mathbf{R}$ の加算回数は「入次数 $-1$」を入次数が正のノードについて和をとったものである。「$A$ がすべて 0 の行や列を含まない」という仮定より、入力ノード以外に入次数0のノードは存在しないし、出力ノード以外に出次数0のノードは存在しない。よって加算回数は「有向辺の本数 $-$ 入力ノード以外のノード数」と等しい。よって $G^T$ が表す算術的アルゴリズムは $G$ が表す算術的アルゴリズムと比べて、$\mathbf{R}$ の乗算回数は同じで加算回数は $-n+m$ される。証明終わり。

このようにして、機械的にアルゴリズムの転置を取ることができる。しばしば、この転置を用いて定数倍算術計算量の低いアルゴリズムが設計できる。

Middle product

行列の転置で得られる問題として middle product がある。
$d$ 次多項式 $Q(x)$ について $A(x) \mapsto A(x) Q(x)$ という線形変換を考える。ここで $A(x)$ は高々 $s$ 次の多項式としよう。
多項式 $A(x)=\sum_{n=0}^s a_nx^n$ を $[a_0\,a_1\,\dotsc a_s]^T=:a$ というベクトルで表すことにする。
一方で、多項式 $B(x)=\sum_{n=0}^{s+d} b_nx^n$ を $[b_{s+d}\,b_{s+d-1}\,\dotsc b_0]^T=:b$ というベクトルで表すことにする(低次数の係数が下にくることに注意)。
また、線形変換 $A(x) \mapsto A(x) Q(x)$ を表す行列を $Q$ と書くことにすると、

\langle b, Qa\rangle = [x^{s+d}]\, B(x)Q(x)A(x) = [x^{s+d}]\, (Q(x)B(x))A(x) = \langle c, a\rangle

となる。ここで $c$ は $s+1$ 次元ベクトルで第 $n\in\{0,1,\dotsc,s\}$ 成分に $Q(x)B(x)$ の $s+d-n$ 次の係数を持つものである。
よって $b\mapsto Q^Tb = c$ という線形変換は「$b$ が表す多項式 $B(x)$ に $Q(x)$ を掛けて $d$ 次から $s+d$ 次を取り出したもの」という計算をしていることが分かる。
これを middle product と呼ぶ。転置原理により $Q^T$ を掛ける算術計算量は $Q$ を掛ける算術計算量とほぼ同じはずである。特に $s=d$ のとき、middle product の算術計算量は $\mathsf{M}(d)+d$ となる。

Bostan–Mori のアルゴリズムの転置

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))

この計算の中の while 文の1ステップを考える。$Q(x)Q(-x)$ という計算を含んでいるため、$Q(x)$ については線形ではない。一方で $N$ と $Q(x)$ を固定しているとすると、$P(x)$ の更新については線形の計算になっている。
よって、このアルゴリズムの出力は

\langle e_0, A_{n-1}A_{n-2}\dotsm A_0 p\rangle

と書ける。ただし $n := \lceil\log(N+1)\rceil$ とする。また、$p$ は入力の多項式 $P(x)$ を表すベクトル(低次が上)、$e_0$ は一番上が $1/Q(0)$ でそれ以外は 0 のベクトルとし、$A_i$ はアルゴリズムの $i$ 回目のループにおける $P(x)$ の更新を表す行列とする。ここで

\begin{align*}
Q_0(x)&:=Q(x),& Q_i(x^2)&:=Q_{i-1}(x)Q_{i-1}(-x)
\end{align*}

と定義すると、$A_i$ は「$Q_i(-x)$ を掛けてから偶数番目もしくは奇数番目の係数だけを取りだす」という意味の行列になる。これらの行列の転置 $A_i^T$ を使って

\langle A_0^T A_1^T\dotsm A_{n-1}^T e_0, p\rangle

と計算結果を表すことができる。転置原理より $A_i^T$ を掛けるのは $A_i$ を掛けるのと同じ算術計算量でできるため、もとのアルゴリズムと同じ計算量の新しいアルゴリズムが設計できた(厳密に言うと最後に内積を計算するために算術計算量が $2d-2$ だけ増えているが)。この新しいアルゴリズムは $A_{n-1}^T$ から順番に掛けるので、$N$ を MSB から順番に読んでいく MSB-first なアルゴリズムである。また、初項から決まる $P(x)$ は最後に内積を取る際に初めて登場し、計算の大部分は初項とは独立に計算ができる。これは Fiduccia のアルゴリズムと同様である。

直接的な導出

前節である意味既にアルゴリズムは導出されているのだが、アルゴリズムを機械的に転置したのでよく意味が分からないものになっている。ここでは MSB-first な Bostan–Mori のアルゴリズムを意味の分かる形で直接導出する。
本来の目的は $[x^N]\,P(x)/Q(x)$ を計算することであった。ここで、$P(x)$ は高々 $d-1$ 次の多項式であるため $1/Q(x)$ の $N-d+1$ 次から $N$ 次の係数が分かればいいことになる。

\mathcal{F}_{N,\,d}\left(\sum_{n\ge 0} a_n x^n\right) := \sum_{n=0}^{d-1} a_{N-d+1+n} x^n

と定義する。ただし、$n<0$ のとき $a_n=0$ とする。すると $[x^{d-1}]\, P(x)\mathcal{F}_{N,\,d}(1/Q(x))$ がもとめたい解 $[x^N]\,P(x)/Q(x)$ となる。

よって $\mathcal{F}_{N,\,d}(1/Q(x))$ を計算することを目標とする。

$N=0$ のとき $\mathcal{F}_{N,\,d}(1/Q(x)) = x^{d-1}/Q(0)$ である。また、

\begin{align*}
\mathcal{F}_{N,\,d}\left(\frac1{Q(x)}\right) &= \mathcal{F}_{N,d}\left(\frac{Q(-x)}{Q(x)Q(-x)}\right)\\
&= \mathcal{F}_{N,\,d}\left(Q(-x) x^{N-2d+1}\mathcal{F}_{N,2d}\left(\frac1{Q(x)Q(-x)}\right)\right)\\
&= \mathcal{F}_{2d-1,\,d}\left(Q(-x) \mathcal{F}_{N,2d}\left(\frac1{V(x^2)}\right)\right)\\
\end{align*}

が成り立つ。ここで $V(x^2):=Q(x)Q(-x)$ とする。二つ目の等式では $1/(Q(x)Q(-x))$ の $N-2d$ 次以下の係数を 0 にした。三つめの等式では

\mathcal{F}_{N,\,d}(xA(x))=\mathcal{F}_{N-1,\,d}(A(x))

を使った。

ここで $W(x):=\mathcal{F}_{\lfloor N/2\rfloor,\,d}(1/V(x))$ と定義すると、簡単な計算により、

\mathcal{F}_{N,\,2d}\left(\frac1{V(x^2)}\right)=S(x):=\begin{cases}
xW(x^2),&\text{if $N$ is even}\\
W(x^2),&\text{otherwise}
\end{cases}

であることが分かる。よって

\mathcal{F}_{N,\,d}\left(\frac1{Q(x)}\right) = \mathcal{F}_{2d-1,\,d}(Q(-x)S(x))

が得られる。よって以下のアルゴリズムが得られる。

R.<x> = QQ[]

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

def bostan_mori_msb(N, Q):
  d = Q.degree()
  if N == 0: return x^(d-1) / Q(0)
  V = even(Q * Q(-x))
  W = bostan_mori_msb(N // 2, V)
  if N % 2 == 0: S = x * W(x^2)
  else: S = W(x^2)
  B = S * Q(-x)
  return R(B.list()[d:2*d])

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

このアルゴリズムの算術計算量を考えよう。$Q(x)Q(-x)$ と $S(x)Q(-x)$ の計算で $\mathbf{R}$ 上の演算が発生している。前者は $d$ 次 × $d$ 次で後者は $(2d-1)$ 次 × $d$ 次である。ここで、$S(x)Q(-x)$ についてはすべての係数をもとめる必要はなく、$d$ 次から $(2d-1)$ 次の結果だけ分かればよい。これはまさに middle product であり、LSB-first なアルゴリズムにおける $P(x)Q(-x)$ の計算の転置に相当する。そのため、算術計算量は $\mathsf{M}(d)+d$ となる。上記の sage プログラムでは一般的な多項式乗算を使っているため $(2d-1)$ 次 × $d$ 次の乗算が発生してしまっているが、これを middle product のアルゴリズムに置き換えるとよい。また、middle product の入力 $S(x)$ はその係数の半分は 0 であるため、実際には算術計算量は $\mathsf{M}(d)+d$ ではなく $\mathsf{M}(d)$ になる。

これらのことから MSB-first の Bostan–Mori のアルゴリズムでは1ステップの算術計算量が $2\mathsf{M}(d)$ となることが分かった。よって全体の算術計算量は $2\mathsf{M}(d)\lceil \log(N+1)\rceil$ である。

X^N mod P(X)

Fiduccia のアルゴリズムにおいて線形漸化的数列の $N$ 項目は $a_N=\langle r, v\rangle$ と書けた。ここで $r$ は $x^N\bmod\Gamma(x)$ を表すベクトル(低次が上)で $v$ は初項を表すベクトルであった。
一方、MSB-first な Bostan–Mori のアルゴリズムでは $a_N=\langle u, p\rangle$ と書ける。ここで $u$ は $\mathcal{F}_{N,d}(1/Q(x))$ を表すベクトル(低次が下)で $p$ は $P(x)$ を表すベクトル(低次が上)であった。よって、少し計算をすれば $u$ から $r$ がもとまりそうである。

定理
任意の $d$ 次多項式 $\Gamma(x)$ について、$Q(x) := x^d\Gamma(1/x),\,U(x) := \mathcal{F}_{N,d}(1/Q(x)),\, V(x) := U(x)Q(x)\bmod x^d$ とすると $x^N \bmod \Gamma(x) = x^{d-1}V(1/x)$ が成り立つ。

証明
$x^N$ を $\Gamma(x)$ で割った商と余りを $q(x),\,r(x)$ と置く。すると $x^N = q(x)\Gamma(x)+r(x)$ が成り立つ。
ここで $x$ を $1/x$ に起き替えて $x^N$ を左右に掛ける(ベクトルとして見ると逆順にすることに対応)と

1 = q_{\mathrm{rev}}(x) Q(x) + x^{N-d+1} r_{\mathrm{rev}}(x)

が得られる。ここで $q_{\mathrm{rev}}(x):= x^{N-d}q(1/x),\, r_{\mathrm{rev}}(x):=x^{d-1}r(1/x)$ とする。よって

\frac1{Q(x)} = q_{\mathrm{rev}}(x) + x^{N-d+1} \frac{r_{\mathrm{rev}}(x)}{Q(x)}

が成り立つ。ここで $q_{\mathrm{rev}}(x)$ の次数は高々 $N-d$ なので、$U(x) = \frac{r_{\mathrm{rev}}(x)}{Q(x)} \bmod x^d$ が得られる。
よって $U(x)Q(x) = r_{\mathrm{rev}}(x) \bmod x^d$ が成り立つ。証明終わり。

このことから以下のアルゴリズムが得られる。

def modexp(N, G):
  d = G.degree()
  Q = G.reverse(d)
  U = bostan_mori_msb(N, Q)
  V = U.multiplication_trunc(Q, d)
  return V.reverse(d-1)

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

$\Gamma(x)$ の最上位係数は単元であるとする。また、 modexp の中の d と bostan_mori_msb の中の d が異なるとまずいので、$\Gamma(0)\ne0$ を仮定する。この仮定は本質的ではなく、そもそも $x^k$ が $\Gamma(x)$ を割るとき、$x^N \bmod\Gamma(x) = x^k\left(x^{N-k}\bmod (\Gamma(x)/x^k)\right)$ であるので $\Gamma(0)\ne 0$ を仮定して一般性を失わない。
このアルゴリズムは MSB-first の Bostan–Mori のアルゴリズムの他に $U(x)Q(x)$ の乗算を含む(必要なのは $x^{d-1}$ までの係数であるが)。この計算は高々 $\mathsf{M}(d)$ の算術計算量ですむので、$x^N\bmod\Gamma(x)$ を計算する算術計算量は $2\mathsf{M}(d)\lceil \log(N+1)\rceil + \mathsf{M}(d)$ となる。素朴な繰り返し二乗法と比べて 2/3 倍以下の算術計算量となっている。

線形漸化的数列の N 項目から N+k-1 項目の計算

MSB-first のアルゴリズムを使って線形漸化的数列の $N$ 項目から $N+k-1$ 項目の計算ができる。説明はいずれ追記する。

def slice(N, k, P, Q):
  d = Q.degree()
  W = Q.inverse_series_trunc(d)
  def next(U):
    S = U * Q
    P = -R(S.list()[d:2*d])
    return P.multiplication_trunc(W, d)
  U = bostan_mori_msb(N, Q)
  U += x^d * next(U)
  L = (U * P).list()[d-1:2*d-1]
  A = R(L)
  while len(L) < k:
    A = next(A)
    L += A.list()
  return L[:k]

このアルゴリズムの算術計算量は $2\mathsf{M}(d)\lceil\log(N+1)\rceil + 4 \mathsf{M}(d) + 2\lceil\frac{k}{d}\rceil \mathsf{M}(d)$ である。

応用

$x^N\bmod\Gamma(x)$ の計算は様々な応用を持つ。「線形漸化的数列のN項目の計算」に書いたように行列の $N$ 乗の計算に応用できる。また $\mathbf{R}$ 上の多項式の因数分解や楕円曲線上の点の数え上げのアルゴリズムにも用いられる。Bostan–Mori のアルゴリズムが実際にこれらの問題に適用されることを期待している。

参考文献

C. M. Fiduccia, An efficient formula for linear recurrences, SIAM Journal on Computing, 14(1), pp. 106–112, 1985.
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.
P. Bürgisser, M. Clausen, and M. A. Shokrollahi, Algebraic Complexity Theory, Springer-Verlag Berlin Heidelberg, 1997.
A. Bostan, G. Lecerf, and E. Schost, Tellegen's principle into practice, In Proceedings of the 2003 international symposium on Symbolic and algebraic computation, pp. 37–44, 2003.
G. Hanrot, M. Quercia, and P. Zimmermann, The middle product algorithm I, Applicable Algebra in Engineering, Communication and Computing, 14(6), pp. 415–438, 2004.

X^N mod P(X) の計算