Skip to content →

Multipoint Evaluation の アルゴリズム

multipoint evaluation とは、与えられた \(m\) 点での多項式の値を評価する問題です:

Multipoint Evalution :
\(m\) 点 \(x_1, x_2, \ldots, x_m\) が与えられます。\(f(x) = \sum_{i=0}^n a_i x^i\) の各点での値 \(f(x_1), f(x_2), \ldots, f(x_m)\) を求めてください。

愚直に計算すると \(O(nm)\) ですが、工夫すると \(O(m \log(m)^2 + n\log(n))\) にできます。

アルゴリズム

アルゴリズムを構築します。カギになるのは剰余の定理です:

剰余の定理 :
任意の \(a\) に対して \(f(a) = f(x) \bmod (x-a)\) が成立。

剰余の定理より、\(f(x) \bmod (x-x_1), f(x) \bmod (x-x_2), \ldots, f(x) \bmod (x-x_n)\) を列挙する問題に言い換えられます。\(k\) 次多項式の割り算には \(O(k \log (k))\) 時間が掛かるので、そのまま剰余の定理を適用すると \(O(nm \log(n))\) になってしまいます。そこで、doublingを使います。

まず、次の二つの条件を満たす根付き木 “subproduct tree” を構築します。

    1. \(m\) 個の葉があり、それぞれ \(x-x_1, x-x_2, \ldots, x-x_m\) が書いてあります。
    2. 葉以外の頂点には、子に書かれた多項式の積が書いてあります。

\(m = 4\) のときは下図のようになります。\(1\) が書かれた葉を全2分木になるまでつけ足しても答えは変わらず、\(m\) のオーダーも変化しません。従って、\(m=2^d\) だと仮定して一般性を失いません。subproduct tree の構築では \(2^k\)次多項式の掛け算が \(2^{d-k}\) 回現れます。掛け算の計算量が \(O(k2^k)\) なので全体では \(O(\sum_{k=0}^d 2^{d-k}k2^k) = O(2^dd^2) = O(m \log(m)^2)\) です。

次に、subproduct tree に書かれた多項式で \(f(x)\) を割った余りを、各頂点に書いた根付き木 “subremainder tree”を構築します。subremainder tree の葉に書き込まれた余りが、求めたい値です。親に書き込まれた多項式が \(f \bmod g_0g_1\)、自身に書き込まれた多項式が \(f \bmod g_0\) だとすると、\((f \bmod g_0g_1) \bmod g_0 = f \bmod g_0 \) が成り立ちます。この式に従って根から順に BFS で余りを書き込んでいきます。

subremainder tree の計算量は subproduct tree と同様に計算すると \(O(m\log(m)^2 + n\log(n))\) となります(2項目は根の余りを計算する計算量)。全体の計算量も \(O(m\log(m)^2+n\log(n))\) となります。

評価点が等比数列のとき(Chirp Z-transform)

評価する点が等比数列 \(x_1=q^1,x_2=q^2,\ldots,x_m=q^m\) のとき、\(O(n\log(n))\) で multipoint evaluation できます。FFTの一般化です。

\(2ij = (i+j)^2-i^2-j^2\) を用いて、
\begin{align} f(q^j) =& \sum_{i=0}^{n-1} a_iq^{ij} \\
=& q^{-j^2/2} \sum_{i=0}^{n-1} a_iq^{-i^2/2} \cdot q^{(i+j)^2/2}\end{align}となります。\(\sum\) の中が畳み込みの形なので、FFTで \(O((n+m)\log(n+m))\) です。\(q^{1/2}\) が \(f(x)\) の係数体に含まれないときは、添加して拡大する必要があることに注意。

multipoint evaluationの転置

余談ですが、\(f(x) := \sum_{i=0}^n a_i x^i \) を \(x_0, \ldots, x_m\) で求めよという問題は、
$$\begin{pmatrix}f(x_{0}) \\ f(x_{1}) \\ \vdots \\ f(x_m) \end{pmatrix}=\begin{pmatrix}x_{0}^0 & x_{0}^1 & \ldots & x_{0}^n \\ x_{1}^0 & x_{1}^1 & \ldots & x_{1}^n \\ \vdots & \vdots & \ddots & \vdots \\ x_{m}^0 & x_{m}^1 & \ldots & x_{m}^n \end{pmatrix}\begin{pmatrix}a_{0} \\ a_{1} \\ \vdots \\ a_n \end{pmatrix}$$という行列とベクトルの積を求める問題に言い換えられます。Tellegen’s principle より、転置を取った$$\begin{pmatrix}v_{0} \\ v_{1} \\ \vdots \\ v_n \end{pmatrix}=\begin{pmatrix}x_{0}^0 & x_{1}^0 & \ldots & x_{m}^0 \\ x_{0}^1 & x_{1}^1 & \ldots & x_{m}^1 \\ \vdots & \vdots & \ddots & \vdots \\ x_{0}^n & x_{1}^n & \ldots & x_{m}^n \end{pmatrix}\begin{pmatrix}b_{0} \\ b_{1} \\ \vdots \\ b_m \end{pmatrix}$$を(入出力を除き)同じ計算量で求めるアルゴリズムを機械的に導出できます。逆にこちらの問題は \(\sum_{i=1}^n v_it^i = \sum_{i=1}^n \frac{b_i}{1-x_it} \bmod t^{n+1}\) を求めよという問題に言い換えられ、分割統治で通分していくことで multipoint evaluation と同じ計算量 \(O(m\log(m)^2+n\log(n))\) になります。そして、Tellegen’s principle により multipoint evaluation のアルゴリズムが機械的に導出できます。転置を取る前後のどちらの方が解きやすかったでしょうか。

演習問題

問題1:\(N! \bmod p\) (\(p\) : 素数) を \(O(\sqrt{p}\log(p)^2)\) で求めて下さい。

解答:\(M:=\lfloor \sqrt{N} \rfloor\) と置きます。まず、\(f(x):=\prod_{i=1}^{M}(i+x)\) を分割統治法を用いて \(O(M\log(M)^2)\) で求めます。\(f(x)\) と \(N!\) の間に $$N!=f(0)f(M)f(2M) \cdots f(M(M-1)) \prod_{i=M^2+1}^N i$$ が成り立ちます。\(f(0),f(M),f(2M),\ldots,f(M(M-1))\) は multipoint evaluation により \(O(M\log(p)^2)\) で求まります。\(O(N-M^2)=O(M)\) なので、全体で所望の計算量を達成できました。

一般に、p-recursive な数列の \(N\) 項目は multipoint evaluation を使うと高速に求まります。更に工夫すると、\(O(\log p)\) を落とせることが知られています[4]。

参考文献

  1. Borodin, Allan, and Robert Moenck. “Fast modular transforms.” Journal of Computer and System Sciences 8.3 (1974): 366-386.
    • このアルゴリズムが提案された論文
  2. Bostan, Alin, Grégoire Lecerf, and Éric Schost. “Tellegen’s principle into practice.” Proceedings of the 2003 international symposium on Symbolic and algebraic computation. 2003.
    • Tellegens’s principle を用いて multipoint evaluation を定数倍高速化しています。
  3. アルゴリズム/データ構造を語る会第1回tkoさん
  4. Alin Bostan, Pierrick Gaudry, and Eric Schost, Linear recurrences with polynomial coefficients and application to integer factorization and Cartier-Manin operator.
掛け算の doubling はよく見ますが、割り算の doubling とは驚き。

Published in データ構造

Comments

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

CAPTCHA