【演算法筆記】Lagrange 插值法
- 本文總計約 5000 字，閱讀大約需要 30 分鐘。
- 警告！警告！警告！本文有大量的 \(\LaTeX\) 公式渲染，可能會導致載入異常緩慢！
前言
又一次沒有前言，以後要不就不寫前言了吧，怪麻煩的 QwQ。
題目引入
根據代數基本定理，我們知道 \((n+1)\) 個橫座標互不相同的點可以唯一確定一個最高次項為 \(n\) 的多項式。例如：
點 \((1,2)\) 和 \((2,3)\) 就可以唯一確定多項式 \(f(x)=x+1\)；點 \((1,1),(2,4),(3,9)\) 可以唯一確定多項式 \(g(x)=x^2\)。
那麼，給定任意 \(n\) 個點，你能計算出其所唯一確定的多項式 \(f(x)\)
暴力怎麼做
當然，我們依舊是要想一想，暴力怎麼做？
假如這 \(n\) 個點分別為 \((x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)\)，那麼我們可以定義其確定的多項式為 \(f(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}\)。
分別代入這 \(n\) 個點，得到一個方程：
\[\begin{cases} y_1=a_0+a_1x_1+a_2x_1^2+\cdots+a_{n-1}x_1^{n-1},\\ y_2=a_0+a_1x_2+a_2x_2^2+\cdots+a_{n-1}x_2^{n-1},\\ \cdots\\ y_n=a_0+a_1x_n+a_2x_n^2+\cdots+a_{n-1}x_n^{n-1}.\\ \end{cases} \]當然，這個方程，我們就可以用 Gauss 消元法來解決，複雜度為 \(\mathcal{O}(n^3)\)
『哇啊，那這個複雜度太高了吧？能過 \(10^4\) 的資料嗎？』
顯然是不能的，所以，我們要做時間複雜度的優化。
Lagrange 插值
為了優化時間複雜度，著名的數學家 \(\mathcal{J-L.Lagrange}\) 給出了一個更加高效的演算法，即 Lagrange 插值法。
我們知道，對任意多項式 \(f(x)\)，若 \(x=a\)，那麼 \(f(x)\equiv f(a)\pmod{(x-a)}\)。
因此，我們分別代入 \((x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)\)，就有了：
\[\begin{cases} f(x)\equiv y_1\pmod{(x-x_1)},\\ f(x)\equiv y_2\pmod{(x-x_2)},\\ \cdots\\ f(x)\equiv y_n\pmod{(x-x_n)}.\\ \end{cases} \]我們使用中國剩餘定理，得到：
因此，有：
\[f(x)=\sum_{i=1}^n y_i\prod_{i\neq j}\dfrac{x-x_j}{x_i-x_j}. \]上述多項式就是 Lagrange 插值法的公式，通過這個公式，我們可以通過 \(n\) 個點值求出其唯一確定的多項式。
通過這個公式，求多項式的時間複雜度為 \(\Theta(n^2)\)。
程式碼
這個演算法在洛谷上的模板為 P4781【模板】拉格朗日插值，因為我們只需要求 \(f(k)\) 的值，所以直接把 \(x=k\) 代入到上式中即可知：
\[f(k)=\sum_{i=1}^n y_i\prod_{j\neq i}\dfrac{k-x_j}{x_i-x_j}. \]於是就可以直接模擬上述公式了，程式碼如下：
#include <cstdio>
using namespace std;
typedef long long ll;
const ll MOD = 998244353;
const int MAXN = 5001;
int n;
ll x[MAXN], y[MAXN], k, ans;
ll inv(ll val) { //Fermat 小定理求逆元
val = (val % MOD + MOD) % MOD;
ll p = MOD - 2, res = 1;
while(p) {
if(p & 1) res = (res * val) % MOD;
val = (val * val) % MOD; p >>= 1;
}
return res;
}
int main(void) {
scanf("%d%lld", &n, &k);
for(int i = 1; i <= n; ++i)
scanf("%lld%lld", &x[i], &y[i]);
for(int i = 1; i <= n; ++i) { //插值主體
ll tmp = y[i];
for(int j = 1; j <= n; ++j) {
if(i == j) continue;
tmp = (tmp * ((k - x[j]) % MOD + MOD) % MOD * inv(x[i] - x[j])) % MOD;
}
ans = (ans + tmp) % MOD;
}
printf("%lld", ans);
return 0;
}
//by CaO
重心 Lagrange 插值
Lagrange 插值固然效率高，但我們發現其有一個缺點：不可維護。
所謂不可維護，就是如果我們在已知的點的基礎上新增一些未知的點，那麼我們就需要整個重新套用一次公式。
那麼，單次重構的時間複雜度依舊為 \(\Theta(n^2)\)。
『如果我要進行 \(m\) 次重構，那麼總時間複雜度就是 \(\Theta(mn^2)\) 了，如果 \(n\le 10^3\)，\(m\le 10^3\)，那麼時間複雜度就不夠了！』
所以，單次重構的時間複雜度太高，那麼我們就要想辦法降低重構的複雜度。
我們手玩一下上面的那個式子：
\[f(x)=\sum_{i=1}^n y_i\prod_{i\neq j}\dfrac{x-x_j}{x_i-x_j}. \]我們可以把求積符號拆分：
\[\begin{aligned} f(x) = \sum_{i=1}^n \left(y_i\prod_{i\neq j}(x-x_j)\prod_{i\neq j}\dfrac{1}{x_i-x_j}\right)\\ \end{aligned} \]在求和號裡面同時乘 \((x-x_i)\) 併除以 \((x-x_i)\)。就有：
\[f(x) = \sum_{i=1}^n \left(y_i\prod_{j=1}^n (x-x_j)\cdot\dfrac{1}{x-x_i}\prod_{i\neq j}\dfrac{y_i}{x_i-x_j}\right). \]注意到 \(\prod_{j=1}^n x-x_j\) 是一個常數，提出來，則有：
\[\begin{aligned} f(x) &= \prod_{j=1}^n(x-x_j)\left(\sum_{i=1}^n y_i\dfrac{1}{x-x_i}\prod_{i\neq j}\dfrac{1}{x_i-x_j}\right)\\ &= \prod_{i=1}^n(x-x_i)\left(\sum_{i=1}^n \dfrac{y_i}{\prod_{i\neq j}(x_i-x_j)}\cdot \dfrac{1}{x-x_i}\right)\\ \end{aligned} \]定義：
\[A(x)=\prod_{i=1}^n (x-x_i),\\ W_i= \dfrac{y_i}{\prod_{i\neq j}(x_i-x_j)}. \]原式即變為：
\[f(x)=A(x)\sum_{i=1}^n\dfrac{W_i}{x-x_i}. \]我們會發現，每次新插值時，維護 \(A(x)\) 的複雜度為 \(\mathcal{O}(1)\) 的，維護 \(W(1),W(2),\cdots,W(n+1)\) 的總複雜度為 \(\mathcal{O}(n)\) 的。所以單次重構的時間複雜度為 \(\mathcal{O}(n)\)。
這種演算法被稱為重心 Lagrange 插值法，程式碼留給讀者作為練習（其實是我太菜了 QwQ）。
例題
本題目列表會持續更新。