【学习笔记】拉格朗日插值(Lagrange)
前言
本学习笔记同步发布于 My Blog
蒟蒻的数学简直太垃圾了,什么都不会,所以要从最简单的东西开始......要是出现错误,请在下方评论区指出。
正题
拉格朗日插值
众所周知,一个 $n$ 次多项式函数可以被 $n + 1$ 个点唯一确定。
那么现在给定 $n + 1$ 个点:
$$(x_0,y_0),(x_1,y_1)......,(x_n,y_n)$$
显然需要满足条件 :
$$\forall i \ne j,x_i \ne x_j$$
(否则不是函数了)
显然我们可以通过高斯消元做到 $O(n^3)$,但是它太慢了......
所以我们需要拉格朗日插值,它可以做到 $O(n^2)$ 的优秀复杂度。
我们现在想要 $n + 1$ 个函数,第 $i$ 个 $f_i(x)$ 满足:在 $x = x_i$ 时,$f_i(x) = y_i$ ,否则 $f_i(x) = 0$。
则我们需要的函数显然为 :
$$g(x) = \sum_{i = 0}^n f_i(x)$$
现在的关键就是构造所有的 $f_i$。
如果我们只是希望 $\forall x \ne x_i,f_i(x) = 0$ 的话,显然有一个构造的方法:
$$f_i(x) = \prod_{j \ne i} {(x - x_j)}$$
这个应该很显然。
那么现在需要做的,就是让这个函数在 $x = x_i$ 时 $f_i(x) = y_i$。
首先肯定乘上一个 $y_i$。
$$f_i(x) = \prod_{j \ne i} {(x - x_j)} \times y_i$$
现在要做的,就是让 $x = x_i$ 时,$y_i$ 乘上的是个 $1$。
也就是我们需要考虑,当 $x = x_i$ 的时候:
$$\prod_{j \ne i}(x-x_j) \times ? = 1$$
则我们可以乘上一个:
$$\frac{1}{\prod_{j \ne i}(x_i - x_j)}$$
那么 $f_i$ 就出来啦:
$$f_i(x) = \frac{\prod_{j \ne i}(x-x_j)}{\prod_{j \ne i}(x_i - x_j)}\times y_i$$
我们需要的函数 $g(x)$:
$$g(x) = \sum_{i = 0}^n f_i(x)$$
写在一起就是:
$$g(x) = \sum_{i = 0}^n (y_i \times\prod_{j \ne i}{\frac{(x-x_j)}{(x_i - x_j)})}$$
对于特殊点的优化
如果这些点满足 $x_i = i$ ,那么这个函数可以写成:
$$g(x) = \sum_{i = 0}^n (y_i \times \prod_{j \ne i}{\frac{(x-j)}{(i - j)})}$$
考虑优化掉 $\prod_{j \ne i}{\frac{(x-j)}{(i - j)}}$ 的这个 $O(n)$。
假设我们现在要求的是 $g(k)$,处理分子可以设置数组:
$$pre[i] = \prod_{j=0}^i k-j$$
$$suf[i] = \prod_{j=i}^n k-j$$
这个显然可以 $O(n)$ 处理。
那么就变成:
$$g(k) = \sum_{i = 0}^n ( y_i \times pre[i-1]\times suf[i+1] \times \prod_{j \ne i}{\frac{1}{(i - j)})}$$
分母其实就是一个类似阶乘相乘的东西,令:
$$a[i] = a[i - 1] \times i,a[0]=1$$
$$b[i] = b[i - 1] \times i \times -1,b[0]=1$$
显然也是 $O(n)$ 的。
那么就会变成:
$$g(k) = \sum_{i = 0}^n ( y_i \times \frac{pre[i-1]\times suf[i+1]}{a[i] \times b[n-i]})$$
重心拉格朗日插值
终于记起来把最后一点点补上了 同时,如果我们需要支持动态插入一个点,然后询问某一个单点的值时,每次暴力做都需要 $\mathcal{O}(k^2)$ 的时间,不太划算,所以考虑使用 重心拉格朗日插值法。
$$f(x) = \sum_{i = 1}^n (y_i \times\prod_{j \ne i}{\frac{(x-x_j)}{(x_i - x_j)})}$$
看看这个柿子,若我们要求 $f(k)$,则可以令:
$$g=\prod_{i = 1}^n (k-x_i)$$
那么可以化成:
$$f(k) = g\sum_{i = 1}^n (\dfrac{y_i}{k-x_i} \times\prod_{j \ne i}{\frac{1}{x_i - x_j})}$$
令:
$$t_i=\dfrac{y_i}{\prod_{j\ne i} x_i-x_j}$$
那么:
$$f(k)=g\sum_{i=1}^n \dfrac{t_i}{k-x_i}$$
每次询问,可以直接 $\mathcal{O}(n\log n)$ 求解 $g$ 与后面的柿子,$\log$ 是由于需要逆元。插入一个新点,则同样需要 $\mathcal{O}(n\log n)$ 的时间对前面所有的 $t$ 更新并求出新的 $t$。
拉格朗日插值求多项式系数
发现以上都在求点值,但如果多项式只有一个,询问却非常多,不妨求出这个多项式的系数然后每次 $\mathcal{O}(n)$ 算。
$$f(x) = \sum_{i = 1}^n (y_i \times\prod_{j \ne i}{\frac{x-x_j}{x_i - x_j})}$$
继续观察这个柿子 发现除了 $\prod_{j\ne i} x-x_j$ 之外都是常数,不妨考虑直接模拟多项式乘法、除法来直接求出上述的柿子展开的形式,然后把系数乘上就可以了。
板子 P4781 【模板】拉格朗日插值
$\texttt{Code:}$
//Code By CXY07
//Lagrange interpolation formula
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int MAXN = 2e3 + 10;
const int mod = 998244353;
int n,k,ans;
int x[MAXN],y[MAXN];
inline int inv(int x) {
int res = 1,b = mod - 2;
while(b) {
if(b & 1) (res *= x) %= mod;
(x *= x) %= mod;
b >>= 1;
}
return res;
}
signed main () {
scanf("%lld%lld",&n,&k);
for(register int i = 1;i <= n; ++i)
scanf("%lld%lld",&x[i],&y[i]);
for(register int i = 1,a,b;i <= n; ++i) {
a = b = 1;
for(register int j = 1;j <= n; ++j) {
if(i == j) continue;
(a *= (k - x[j] + mod) % mod) %= mod;
}
for(register int j = 1;j <= n; ++j) {
if(i == j) continue;
(b *= (x[i] - x[j] + mod) % mod) %= mod;
}
(ans += a * inv(b) % mod * y[i] % mod) %= mod;
}
printf("%lld\n",ans);
return 0;
}
习题
CF622F The Sum of the k-th Powers
给定 $n,k,n \le 10^9,k \le 10^6$ 求:
$$\sum_{i=1}^n i^k mod (10^9+7)$$
$\texttt{Solution:}$
发现这是一个 $k+1$ 次多项式 不会证
拿 $k+1$ 个点直接插值就好了,这些点要小、好算。
发现是 $O(k^2)$ 好像还有个log 的,因此我们选择取 $x=1$ ~ $k+2$ 的点,就可以用上面的 方法 $O(k)$ 解决 貌似带个log
[TJOI2018]教科书般的亵渎
题意贼模糊,自己去看 $8$。
$\texttt{Solution:}$
首先可以发现 $k=m+1$。
发现只有很少的数字 $(\le 50)$ 没有,本身的 $n$ 又贼大。
发现连续的一段会被一个亵渎全部减完,贡献类似 $\sum_{i=1}^n i^k$ 减去那些还不为 $0$ 的断点的贡献,直接插值即可,复杂度应该是 $O(m^2 log m)$ 的(有个快速幂)。
关键代码(把 $0$ 当作了第一个断点):
for(register int i = 0;i <= m; ++i) {//s[0] = 0
(ans += Lagrange(n - s[i])) %= mod;
for(register int j = i + 1;j <= m; ++j)
(ans -= ksm(s[j] - s[i],k)) %= mod;
}
咕咕咕