【學習筆記】拉格朗日插值
學習多項式的第一步。
參考資料：
1.拉格朗日插值法的簡介
問題：
解法1：高斯消元
顯然\(\deg \geqslant n\)的多項式有無窮個，因為根據高斯消元，一定會出現自由元。
直接把這\(n\)個點列出方程組，用高斯消元求解。
求出多項式後再求出\(f(k)\)即可。
時間複雜度\(O(n^3).\)
解法2：拉格朗日插值法
給出一個關於點\((x_i,y_i)\)的基函式：
\[g(k)=\prod_{j\not=i}^{1\leqslant j\leqslant n} \dfrac{k-x_j}{x_i-x_j} \]容易發現：\(\forall j\not=i,g(x_j)=0.\)因為累乘中總有\(k=x_j\),使得\(k-x_j=x_j-x_j=0,g(x_j)=0.\)
對於\(j=i,g(x_j)=1.\)因為累乘中每一項均為\(\dfrac{x_i-x_j}{x_i-x_j}=1,g(x_j)=1.\)
於是我們的多項式就可以表示為：
\[f(k)=\sum_{i=1}^{n}y_i\times\prod_{j\not=i}^{1\leqslant j\leqslant n} \dfrac{k-x_j}{x_i-x_j} \]這樣\(\forall(x_i,y_i),f(x_i)=y_i\),也可以據此求出\(f(k).\)
大概因為要求逆元，所以時間複雜度為\(O(n^2+n\log n)=O(n^2).\)
\(\tt{code:}\)
#include <bits/stdc++.h> #define ll long long using namespace std; int n; const ll mod = 998244353; ll x[2011], y[2011], k, ans; ll ksm(ll s1, ll s2) { if(!s2) return 1; if(s2 & 1) return s1 * ksm(s1, s2 - 1) % mod; ll ret = ksm(s1, s2 / 2); return ret * ret % mod; } int main() { 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 sum = 1, mul = 1; for(int j = 1; j <= n; j++) { if(i == j) continue; sum *= (k - x[j]); sum %= mod; sum += mod; sum %= mod; mul *= (x[i] - x[j]); mul %= mod; mul += mod; mul %= mod; } sum *= ksm(mul, mod-2); sum %= mod; ans += (sum * y[i]) % mod; ans %= mod; } printf("%lld\n", ans); return 0; }
2.拉格朗日插值法的拓展
問題：
同上，加入\(x_i\)的值連續的限制。
解法：
假設\(x_i\in[1,n],\)則公式化為：
\[f(k)=\sum_{i=1}^{n}y_i\times\prod_{j\not=i}^{1\leqslant j\leqslant n} \dfrac{k-j}{i-j} \]可以維護\(num=\prod_{1\leqslant j\leqslant n}k-j,\)於是分子部分化為\(\dfrac{num}{k-i},\)因為這麼做要求逆元，帶個log，很不優，所以可以預處理字首和\(pre_i\)和字尾和\(suf_i.\)
對於分母，可以化為\(1\times 2\times···\times i-1\times (-1)\times (-2)\times ···\times(i-n)\)的形式。
實質上是一個階乘，|分母|\(=fac_{i-1}\times fac_{n-i},\)然後再處理一下正負即可。
這樣再加上一些\(O(n)\)的預處理\(pre,suf,fac,inv...\)，我們只需要\(O(n)\)即可完成。
3.拉格朗日插值法的具體應用&例題
【例1】CF622F The Sum of the k-th Powers
拉格朗日插值的經典模型，即求\(\sum_{i=1}^n i^k\).
這是個\(k+1\)次的多項式，感性證明一下：
考慮\(n\)次多項式\(g(x)=a_nx^n+a_{n-1}x^{n-1}+...+a_0\).
它的一階差分\(\Delta g(x)=g(x+1)-g(x)=a_n(x+1)^n-a_nx^n+...\).
通過二項式定理，我們發現這時的\(n\)次項會被消掉，即：
每做一次差分，多項式的度便會減一。
令\(f(x)=\sum_{i=1}^ni^k,\)則\(\Delta f(x)=f(x+1)-f(x)=(x+1)^k,\)為\(k\)次多項式。
所以\(f(x)\)為\(k+1\)次多項式。
至此，我們只需確定\(k+2\)個點，做一個插值，就可以求出\(f(n)\)了。
因為\(i\)連續，所以可以用到2.2中的拉格朗日插值拓展做到\(O(k\log k).\)
\(i^k\)也可以線性求出，這裡不多作說明。
\(\tt{code:}\)
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N = 1e6;
const ll mod = 1e9 + 7;
ll n, k, fac[N+11], pre[N+11], suf[N+11], ans;
ll ksm(ll s1, ll s2) {
if(!s2) return 1ll;
if(s2 % 2) return s1 * ksm(s1, s2-1) % mod;
ll ret = ksm(s1, s2/2);
return ret * ret % mod;
}
int main() {
scanf("%lld%lld", &n, &k);
ll sum = 0;
fac[0] = 1;
for(int i = 1; i <= k+2; i++)
fac[i] = fac[i-1] * (ll)i % mod;
pre[0] = suf[k+3] = 1;
for(int i = 1; i <= k+2; i++)
pre[i] = pre[i-1] * (n - i) % mod;
for(int i = k+2; i >= 1; i--)
suf[i] = suf[i+1] * (n - i) % mod;
for(int i = 1; i <= k+2; i++) {
sum += ksm(i, k); sum %= mod;
ll num = pre[i-1] * suf[i+1] % mod;
ll fm = fac[i-1] * fac[k+2-i] % mod;
num *= sum; num %= mod;
num *= ksm(fm, mod-2); num %= mod;
if((k+2-i) % 2 == 0) ans += num;
else ans -= num;
ans %= mod; ans = (ans + mod) % mod;
}
printf("%lld\n", ans);
return 0;
}
【例2】luogu P4593 [TJOI2018]教科書般的褻瀆
問題轉化：先求出\(k,\)即“褻瀆”的次數。忽略大於\(n\)的\(a_i\)（雖然不知道有沒有），容易發現一次“褻瀆”可以使一段血量連續的怪物全部死掉，所以\(k=m+0?1:(\exists a_i=1).\)
又因為怪物的血量互不相同，所以我們可以把求的值轉化為\(\sum_{i=1}^n i^k,\)於是就和【例1】做法相似了。
程式碼就不放了。