【笔记】拉格朗日插值还原系数表达式
大家都知道拉格朗日插值的公式(已知 \(n\) 个点 \((x_i,y_i)\),求唯一确定的经过这 \(n\) 个点的\(n-1\) 次多项式):
\[f(x)=\sum_{i=1}^ny_i\left(\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}\right)
\]
但是洛谷模板只让求了一个点的点值,就不用把系数表达式算出来。本蒻稽不会 \(O(n\log n)\) 的多项式快速插值又想在省选骗分,只好学一下 \(O(n^2)\) 的插值方法。
上面公式分成两部分,一部分是分子 \(y_i \prod_{j\ne i} \frac 1{x_i-x_j}\),\(O(n)\) 可以计算,但分母部分 \(\prod_{j\ne i}(x-x_j)\) 需要 \(O(n^2)\)。
优化方案是先把 \(F(x)=\prod_{i=1}^n (x-x_i)\) 算出来,然后用 \(F(x)/(x-x_j)\) 来得到 \(f_j(x)\).
设 \(F(x)_i\) 为 \(F(x)\) 第 \(i\) 项的系数,\(f_j(x)_i\) 同理,有:
\[F(x)_0=-x_jf_j(x)_0
\]
\[F(x)_i=f_j(x)_{i-1}-x_jf_j(x)_{i},\quad (i>1)
\]
于是得到:
\[f_j(x)_0=-F(x)_0/x_j
\]
\[f_j(x)_i=(f_j(x)_{i-1}-F(x)_i)/x_j,\quad (i>1)
\]
但是注意这样手动模拟的前提是 \(x_j\ne 0\)。\(x_j=0\) 的情况要特判掉(\(f_j(x)=F(x)/x\)).
算出 \(f_j(x)\) 后,把之前的分子部分乘到系数上。最后 \(f(x)=\sum_{j=1}^n f_j(x)\) 就做完了。
然后把代码写出来就好了(洛谷模板):
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 2005, D = 998244353;
int fpm(int a, int b = D - 2, int c = D) {
int res = 1 % c;
while(b) {
if(b & 1) res = (ll)res * a % c;
b >>= 1, a = (ll)a * a % c;
}
return res;
}
inline int _(int a) { return a + (a >> 31 & D); }
inline void Add(int &a, ll b) { a = (a + b) % D; }
inline void Sub(int &a, ll b) { a = _((a - b) % D); }
inline void Div(int &a, int b) { a = (ll)a * fpm(b) % D; }
inline void Mul(int &a, int b) { a = (ll)a * b % D; }
int n, k, x[N], y[N], f[N];
void trans(int x[], int y[], int f[]) {
static int c[N], F[N], t[N], ix[N];
for(int i = 1; i <= n; i++) c[i] = 1, f[i] = 0;
for(int i = 1; i <= n; i++) {
ix[i] = fpm(x[i]);
for(int j = 1; j <= n; j++)
if(j != i) Mul(c[i], _(x[i] - x[j]));
c[i] = (ll)fpm(_(c[i])) * y[i] % D;
}
for(int i = 0; i <= n; i++) F[i] = 0;
F[0] = 1;
for(int i = 1; i <= n; i++)
for(int j = i; j >= 0; j--)
F[j] = _((F[j - 1] - (ll)x[i] * F[j]) % D);
for(int i = 1; i <= n; i++) {
if(!x[i]) for(int j = 1; j <= n; j++) t[j - 1] = F[j];
else {
t[0] = _((ll)-F[0] * ix[i] % D);
for(int j = 1; j < n; j++)
t[j] = _((ll)(t[j - 1] - F[j]) * ix[i] % D);
}
for(int j = 0; j < n; j++) Add(f[j], (ll)t[j] * c[i] % D);
}
}
int main() {
ios::sync_with_stdio(false); cin.tie(nullptr);
cin >> n >> k;
for(int i = 1; i <= n; i++) cin >> x[i] >> y[i];
trans(x, y, f);
int tmp = 1, ans = f[0];
for(int i = 1; i < n; i++)
Mul(tmp, k), Add(ans, (ll)tmp * f[i]);
cout << ans << '\n';
return 0;
}