問題設定

未知の \(k\) 次の多項式 \(f(x)\) を考えます．これについて \(y_0 = f(x_0)\), \(y_1 = f(x_1)\), ..., \(y_k = f(x_k)\) が得られているとします．このとき，任意の \(x\) について \(f(x)\) の値を求めてみましょう． \(k\) 次式は \(k+1\) 点によって一意に定まります*1．

いくつかの点での値から，一般の点での値を復元することを補間と言います．

以下では Lagrange（ラグランジュ）補間と呼ばれる方法について説明します．

考察

以下の \(k\) 次多項式をひらめきます．この \(x_0\), ..., \(x_k\) は上のものです． \[l_j(x) = \prod_{0\le m\le k, m\neq j} \frac{x-x_m}{x_j-x_m} = \frac{x-x_0}{x_j-x_0} \cdots \frac{x-x_{j-1}}{x_j-x_{j-1}} \frac{x-x_{j+1}}{x_j-x_{j+1}}\cdots \frac{x-x_k}{x_j-x_k}\] これについて，\(l_j(x_j) = 1\), \(l_j(x_i) = 0\) (\(i \neq j)\) であることに注意しましょう．

この \(l_j(x)\) を用いた以下の \(k\) 次多項式を考えます． \[L(x) = \sum_{j=0}^k y_j l_j(x)\] \(l_j(x)\) の性質から，\(L(x_i) = y_i\) となってくれます．

さて，\(L\) は \(x_0\) から \(x_k\) の \(k+1\) 点での値が \(f\) と一致していますので，一意性から \(L = f\) であることが言えます．うれしい．

よって，あとは \(L(x)\) を高速に求められればよいです．どうしよう．

もっと考察

\(l(x) = (x-x_0)(x-x_1)\cdots(x-x_k)\) とおき，これの微分を考えてみましょう．積の微分公式から次のようにできます． \[l'(x) = \sum_{i=0}^k \prod_{j=0, j\neq i}^k (x-x_j)\] ですので，\(l'(x_j) = \prod_{i=0, i\neq j}^k (x-x_i)\) となります．これは \(l_j(x)\) の分母と等しいです*2．また，\(\frac{l(x)}{(x-x_j)}\) は \(l_j(x)\) の分子と等しくなります．ゼロ割りのことを考えたくないので，各 \(x_i\) からは目を背けます．

結局，次のように書けることになります． \[L(x) = \sum_{j=0}^k y_j\cdot \frac{l(x)/(x-x_j)}{l'(x_j)}\] \(l(x)\) を分離して，次のようにします． \[L(x) = \left(\sum_{j=0}^k \frac{l'(x_j)^{-1}}{x-x_j}\cdot y_j\right) \cdot l(x)\]

\(l(x) = \prod_{i=0}^k (x-x_i)\) でしたので，与えられた \(x\) について \(l(x)\) を計算するのは \(O(k)\) 時間でできます．\(\sum\) の各項について，\(l'(x_j)^{-1}\) を求める部分以外は \(O(1)\) でできるので，あとは \(l'(x_j)^{-1}\) について考えます．

さて，\(l'\) は次の通りに書けるわけでした． \[l'(x_j)^{-1} = \prod_{i=0, i\neq j}^k (x_j-x_i)^{-1}\] う〜ん，一般のこれを定数時間で求めるのは大変そう．\(O(k)\) 時間で求めるのはできますね．

限定した状況での高速化

いま，\(f\) の関数値を求める \(k+1\) 個の点 \(x_i\) をすきに選べる状況を考えます．

えびちゃん：\(k+1\) 回の質問に答えるよ．聞かれた \(x_i\) に対して \(f(x_i)\) を教えるね．その後えびちゃんが \(x\) を聞くから，\(f(x)\) の値を当ててね．\(k \le 10^6\) で，\(f\) の次数は \(k\) だよ．

あなた：わかった！

あなたは \(l'(x_j)\) を高速に計算できるように各 \(x_i\) を選びたいです．

ばらばらに選ぶと大変そうなので，\(x_i = i\) としてみます．とすると，\(l'(j)^{-1}\) は以下のように書けます． \[l'(j)^{-1} = \prod_{i=0, i\neq j}^k (j-i)^{-1}\]

\[ \begin{aligned} l'(j)^{-1} &= \prod_{i=0, i\neq j}^k (j-i)^{-1} \\ &= ( (j-0)\cdots(j-(j-1)))^{-1} \cdot ( (j-(j+1))\cdots(j-k))^{-1} \\ &= (j\cdots 1)^{-1} \cdot ( (-1)\cdots(-(k-j)))^{-1} \\ &= j!^{-1} \cdot (-1)^{k-j} \cdot (k-j)!^{-1} \end{aligned} \]

ここで，\(j\) と \(k-j\) はそれぞれ \(k\) 以下の非負整数ですので，\(0\le i \le k\) について \(i!^{-1}\) のテーブルを作っておけば \(O(1)\) 時間で求めることができます．

\( (-1)^{k-j}\) についても \(k-j\) の偶奇を判断して \(O(1)\) 時間で求められます．

これにより，全体として \(O(k)\) 時間で求められることになります．

実際には，階乗は非常に大きくなりますので，\(f(x)\) を適当な素数 \(p\) で割ったあまりを答える問題が多いでしょう．このとき，\(i!^{-1}\) は \(p\) を法としたときの \(i!\) の逆元であり，また \(-1\) は \(p-1\) として計算します．

例

では，例として \(k\) 乗和を求める問題を考えてみましょう．与えられた \(k\), \(x\), \(p\) について次の値を求めてください． \[f_k(n) = \left(\sum_{i=1}^n i^k\right) \bmod{p}\] \(1\le k\le 10^6\), \(1\le n\le 10^9\), \(p = 10^9+7\)．

\(f_k(n) - f_k(n-1) = n^k\) となることから，\(f_k(n)\) は \(k+1\) 次式になります*3．実際，\(f_1(n) = n(n+1)/2\) は \(2\) 次式で，\(f_2(n) = n(n+1)(2n+1)/6\) は \(3\) 次式ですね．

このことから，まず \(f_k(0)\) から \(f_k(k+1)\) までの値を求め，補間により \(f_k(n)\) を得るとよいです．