拉格朗日插值学习笔记
前言
我们先引入一个问题，给你平面上一些点，让你求穿过这些点的曲线或者说给你一个次数不超过 \(n\) 次的多项式的 \(n+1\) 个函数值，让你求这个函数的表达式。
我们有一个 \(O(n^3)\) 的高斯消元解法，具体来说就是设多项式表达式为 \(y = a_0x^0+a_1x^1+a_2x^2+...\) ，我们把每一项的系数看成一个未知数，把 \(n+1\) 个函数值带入显然会得到 \(n+1\) 个方程组成的线性方程组，然后用高斯消元求解得到 \(n+1\) 个未知数就是每一项的系数。
二 简介：
我们在这里引入一个比高斯消元更方便的方法：拉格朗日插值，我们先给出公式 （后面会有详细证明）。
\(\large\ f(k) = \displaystyle\sum_{i=0}^{n} y[i] \prod_{i\neq j} {k-x[j]\over x[i]-x[j]}\)
其中 \(k\) 是自变量， \(x[i],y[i]\) 是给出的函数值。
拉格朗日插值有两种：朴素法和重心拉格朗日插值法。在求某个多项式的表达式的时候，朴素法的时间复杂度和高斯消元没有什么区别都是 \(O(n^3)\), 但是在求多项式在某一位置的取值的时候复杂度为 \(O(n^2)\) 的。 重心拉格朗日插值在这两种情况下都是 \(O(n^2)\) 。
三：拉格朗日插值
(1) 朴素法:
我们把每个点带入可以得到 \(n+1\) 个方程的线性方程组，即：
\(\begin{cases}y_0 = a_0x_0^0 + a_1x_0^1+a_xx_0^2 + .....+a_nx_0^n \\y_1 = a_0x_1^0 + a_1x_1^1+a_2x_1^2+.....+a_nx_1^n\\y_2 = a_0x_2^0 + a_1x_2^1+a_2x_2^2+.....+a_nx_2^n\\.........\\y_n = a_nx_n^0+a_1x_n^1+a_2x_n^2+.....+a_nx_n^n\end{cases}\)
考虑构造 \(n\) 个函数，其中 第 \(i\) 个函数，\(fi(k)\) 在 \(k = x_i\) 处取值为 \(1\), 在其他地方 \(k = x_j(i\neq j)\)取值为 \(0\)。
由此可以得到 \(f(k) = \displaystyle\sum_{i=1}^{n} y[i] \times fi(k)\) 。
设 \(fi(k) = \displaystyle\prod_{i\neq j} {k-x[j]\over x[i]-x[j]}\) ，当 \(k = x_i\) 时，分子分母相等, 其取值自然为 \(1\), 当 \(k = x_j (i\neq j)\) 时，后面的连乘的分子中必定有一个为 \(0\), 其取值自然为 \(0\)，因此我们构造出来的函数 \(fi(k)\) 符合条件。
把 \(fi(k)\) 代入到 \(f(k)\) 中就可以得到我们上文中所提到的式子。
根据小学知识， \(n+1\) 个点可以确定出唯一一个 \(n\) 次多项式，所以我们这样构造是正确的。
当求多项式在某一位置的具体取值时，根据拉格朗日插值式子直接算就可以。
当求多项式的表达式的时候，我们需要把 \(k\) 当成未知数，后面的连乘式直接展开的复杂度为 \(O(n^2)\), 所以总的复杂度为 \(O(n^3)\).
(2) \(x[i]\) 取值连续时的优化
有的时候题目中给我们的 \(x[i]\) 是连续的，这时候我们可以尝试继续推一下式子，降低一下复杂度。
首先把 \(x[i]\) 替换成 \(i\), 可以得到: \(f(k) = \displaystyle\sum_{i=0}^{n}y[i]\prod_{i\neq j} {k-j\over i-j}\) 。
预处理一下前缀积和后缀积，即 \(pre[i] = \displaystyle\prod_{j < i} k-j\), \(suf[i] = \displaystyle\prod_{j>i} k-j\) ，然后你就会发现分子其实就是 \(pre[i]\times suf(i)\)。之后预处理一个 \(jz[i]\) 表示 \(i!\), 那么分母就可以写成 \(jz[i]\times jz[n-i]\times (-1)^{n-i}\) 。
然后把这些都带入可以得到一个新式子即：
\(\large\ f(k) = \displaystyle\sum_{i=0}^{n}y[i]\prod_{i\neq j} {pre[i]\times suf(i)\over {jz[i]\times jz[n-i] \times (-1)^{n-i}}}\)
在求多项式在某一位置的具体取值的时候，复杂度就可以将为 \(O(n)\) 的。
但在求多项式的表达式的时候，我们还是需要把连乘式展开，复杂度还是 \(O(n^3)\)
(3) 重心拉格朗日插值法
设 \(g(k) = \displaystyle\prod_{i=0} k-x[i]\), 那么连乘式的分子可以写为：\(g(k)\over k-x[i]\) 。
带入可得：\(f(k) = \displaystyle g(k)\sum_{i=0}^{n}{y[i]\over k-x[i]} \prod_{i\neq j} {1\over x[i]-x[j]}\)
设 \(t(i) = \displaystyle\prod_{i\neq j} {1\over x[i]-x[j]}\) ，化简一下可以得到：
\(\large f(k) = g(k)\displaystyle\sum_{i=0}^{n} {y[i]t[i]\over k-x[i]}\)
对于求多项式的表达式的时候，可以枚举 \(i\), 因为除的只有两项，可以 \(O(n)\) 求出 \(g(k)\over k-x[i]\), 总的复杂度 \(O(n)\)。
当插值点的个数增加一个时，将，每个 \(t[i]\) 都除以 \(x[i]-x[n+1]\)，就可以 \(O(1)\) 得到 \(y[i]t[i]\over k-x[i]\), \(O(n)\) 的复杂度得到新的 \(\displaystyle\sum_{i=0}^{n} {y[i]t[i]\over k-x[i]}\), 最后在和新的 \(g(k)\) 相乘，总的复杂度为 \(O(n)\) 。
所以求多项式的表达式时，我们的复杂度就可以将为 \(O(n^2)\) 的啦。
模板代码
#include<iostream>
#include<cstdio>
#include<algorithm>
using namespace std;
#define int long long
const int p = 998244353;
const int N = 5e3+10;
int n,k,ans,x[N],y[N];
inline int read()
{
int s = 0, w = 1; char ch = getchar();
while(ch < '0' || ch > '9'){if(ch == '-') w = -1; ch = getchar();}
while(ch >= '0' && ch <= '9'){s = s * 10 + ch - '0'; ch = getchar();}
return s * w;
}
int ksm(int a,int b)
{
if(a < 0) a = (a + p) % p;
int res = 1;
for(; b; b >>= 1)
{
if(b & 1) res = res * a % p;
a = a * a % p;
}
return res;
}
signed main()
{
n = read(); k = read();
for(int i = 1; i <= n; i++) x[i] = read(), y[i] = read();
for(int i = 1; i <= n; i++)
{
int tmp = 1;
for(int j = 1; j <= n; j++)
{
if(i != j) tmp = tmp * (k - x[j] + p) % p * ksm(x[i]-x[j],p-2) % p;
}
ans = (ans + tmp * y[i] % p) % p;
}
printf("%lld\n",ans);
return 0;
}
