【学习笔记】拉格朗日插值

前言

我们要知道根据 k+1 个点可以确定一个最高 k 次的多项式，首先要是想真的知道他是谁，高斯消元欢迎你。但是我们很多情况下无法承受高斯消元特别高的复杂度，而且我们往往只需要知道这个函数在某一点的值就好了，所以拉格朗日插值应运而生。

拉格朗日插值一般就是你发现某个东西像一个多项式，而且这个多项式最高 k 次，题目中还给定了这个多项式在 k+1 处所取得的值，那么就可以愉快的拉插过了，注意这里的点是指的类似 (x_i,y_i) 这种坐标，但是其实题目不可能那么明显，但本质还是这个东西

常见的拉格朗日插值形式

注意下文中的 k 均指，我们要求的那个位置。本文均是以 x 下标为 [1,n] 进行讨论。

任意条件下

f(k) = \sum_{i=1}^ny_i\prod_{j \not=i} \dfrac{k - x_j}{x_i-x_j}

虽然不会证但是会发现这是正确的。

时间复杂度 O(n^2)

取值连续，为 [1,n]

那么就可以直接将 x_i ，变为 i，因为我们的 x 这一维本来就是记录他们的位置。

那么就是考虑如何快速计算以下的式子：

\prod_{j\not=i}\dfrac{k - j}{i-j}

我们记以下的两个式子，分别为 k-j 的前缀积和后缀积（高级的名字）：

pre_i = \prod_{j=1}^ik-j \\ suf_i = \prod_{j=i}^n k-j

注意 pre_0 = 1,suf_{n+1} = 1

考虑分母如何化简： 我们发现，分母其实是什么啊？首先要知道 j\in[1,i-1]\bigcup[i+1,n]，这一点看到没有人提，但是当时的我来说确实是没有悟出来。那么下面这一项其实就相当于

\dfrac{1}{(i-1)!} \times \dfrac{1}{(n-i)!} \times (-1)^{n-i}

看着是不是非常的神奇。

首先因为 j\not=i ，我们统计的值是 i-j，所以我们可以考虑将 j 以 i 为界分开，分成两块，一块小于 i，另一块大于 i。

我们考虑对于小于 i 的那一块如何处理：很明显嘛，这个时候就有 j\in[1,i-1]，用 i 去减这个范围，那么就是 [1,i-1] ，所以这一块就是 (i-1)!

考虑对于大于 i 的那一块怎么处理：这个时候既有 j\in[i+1，n]，用 i 去减这个范围，能得到 [i-n,-1]，我们会发现这也不好搞，那么我们就考虑全部乘以一个负一，那么就能得到这个范围 [1,n-i] 会发现这就是 (n-i)! 所以这块就是 (n-i)! 但是不要我们给这 n-i 个数都乘了 -1 ，所以最后自然还要乘回来，所以最后还要乘以 (-1)^{n-i}

分母化完了，下面就要来解决上面留下来的分子的问题了： 这个其实非常好理解，因为 j\not=i，所以通过前后缀积不乘 k - i 这一项就好了

那么分子分母都化完了，我们的式子就可以直接化了（请注意下标）：

f(k) = \sum_{i=1}^ny_i\times\dfrac{pre_{i-1}\times suf_{i+1}}{(i-1)!\times(n-i)!\times(-1)^{n-i}}

这还是不是很好看，我们就记 fac_i = i!，继续化：

f(k) = \sum_{i=1}^ny_i\times\dfrac{pre_{i-1}\times suf_{i+1}}{fac_{i-1}\times fac_{n-i}\times(-1)^{n-i}}

上面是 $x\in[1,n]$的式子，下面给出 x\in[0,n] 的式子，推导方式同上，注意此时 pre,suf 数组也要调整一下范围

f(k) = \sum_{i=0}^ny_i\times\dfrac{pre_{i-1}\times suf_{i+1}}{fac_i\times fac_{n-i}\times(-1)^{n-i}}

注意：因为我们大多数题目要取模，即要对阶乘数组求逆元，那么最好把阶乘数组预处理为逆元，就是求出 fac_n 的逆元，然后从后向前一点点乘就好了 时间复杂度 O(n) 下面是一道模板题的完整代码：

点击查看代码

#include<bits/stdc++.h>
using namespace std;
const long long MOD = 998244353;
const long long MAXN = 1e6+5;
long long pre[MAXN],suf[MAXN],fac[MAXN],y[MAXN],n;
long long power(long long a,long long b){  //求阶乘
long long res = 1;
while(b){
if(b & 1){
res = res * a % MOD;
}
a = a * a % MOD;
b>>=1;
}
return res;
}
long long inv(long long a){  //求逆元
return power(a,MOD-2) % MOD;
}
void prework(long long k){
pre[0] = k;
for(long long i=1; i<=n; i++){   //预处理前缀积
pre[i] = pre[i-1] * ((k - i + MOD) % MOD) % MOD;
}
suf[n+1] = 1;
for(long long i=n; i>=0; i--){   //预处理后缀积
suf[i] = suf[i+1] * ((k - i + MOD) % MOD) % MOD;
}
fac[0] = 1;
for(long long i=1; i<=n; i++){  //预处理阶乘数组
fac[i] = fac[i-1] * i % MOD;
}
fac
 = inv(fac
);  //预处理逆元
for(int i=n-1; i>=0; i--){
fac[i] = (fac[i+1] * (i+1))%MOD;
//逆元可以理解为 负一次方
}
}
long long get_val(long long k){
prework(k);
long long ans = 0;
for(long long i=1; i<=n; i++){  //按照公式乘的
long long res = y[i];
res = res * ((pre[i-1]* suf[i+1]) % MOD) % MOD;
res = res * (fac[i - 1] * fac[n-i] % MOD)%MOD;
res = res * ((n-i) & 1 ? -1:1) % MOD;
res = (res + MOD) % MOD;
ans = (ans + res) % MOD;
}
return (ans + MOD)%MOD;
}
int main(){
long long k;
scanf("%lld%lld",&n,&k);
for(long long i=1; i<=n; i++){
scanf("%lld",&y[i]);
}
printf("%lld\n",get_val(k));
return 0;
}

重心拉格朗日插值

重心拉格朗日插值解决的也是上文提到的这种问题，但是多了一种操作：还会加入某些点，就是再告诉你这个多项式在某些地方的取值，但是这个时候我们每加入一个点跑一遍拉格朗日插值很明显复杂度爆炸了，所以重心拉格朗日插值就诞生了

我们考虑记：

g = \prod_{i=1}^n (k-x_i) \\ t_i = \dfrac{y_i}{\prod_{j\not=i}x_i - x_j}

先把我们最原始的式子摆在这里：

f(k) = \sum_{i=1}^ny_i\prod_{j \not=i} \dfrac{k - x_j}{x_i-x_j}

考虑把 g 带进去优化：

f(k) = g\sum_{i=1}^n(\dfrac{y_i}{k - x_i}\prod_{j\not=i}\dfrac{1}{x_i - x_j})

会发现一个很神奇的东西，y_i 带了个分母，因为当我们将 g 乘进去后会发现，后面的分子不含有 k - x_i 这一项，因为限制条件 j\not=i ，所以应该把 g 中多出来的这一项通过除法消掉

把 t_i 带进去优化：

f(k) = g\sum_{i=1}^n\dfrac{t_i}{k-x_i}

我们考虑每带入一个新值，O(n) 更新 t_i，因为我们不用每一次都重新算一遍 t_i ，直接在原来的上面加上新加入的值的影响就好了，每查询一个值 O(n) 扫一遍

注意我们的 n 也就是点的个数会随时变化

整体的复杂度为 O(n^2)。

经典例题

K次方求和

题目链接：luogu CF622F 题意：给定 n,k 求下面这个式子的值

\sum_{i=1}^ni^k

这个东西是一个 k+1 次多项式，所以我们就找 k+2 个点带进去拉格朗日插值就能求出 n 对应的值，考虑优化就可以找连续的 k 个点，这样复杂度就为 O(k) 代码：

点击查看代码

#include<bits/stdc++.h>
using namespace std;
const long long MOD = 1e9+7;
const long long MAXN = 1e6+5;
long long suf[MAXN],pre[MAXN],y[MAXN],fac[MAXN];
long long power(long long a,long long b){
long long res = 1;
while(b){
if(b & 1){
res = (res * a) % MOD;
}
b>>=1;
a = (a * a) % MOD;
}
return res;
}
long long inv(long long a){
return power(a,MOD-2);
}
void prework(long long k,long long n){  //代表要求在 k 的值，现在有 n 个点
pre[0] = 1;
for(long long i=1; i<=n; i++){
pre[i] = (pre[i-1] * ((k - i + MOD) % MOD)) % MOD;
}
suf[n+1] = 1;
for(long long i=n; i>=0; i--){
suf[i] = (suf[i+1] * ((k - i + MOD) % MOD)) % MOD;
}
fac[0] = 1;
for(long long i=1; i<=n; i++){
fac[i] = (fac[i-1] * i) % MOD;
}
fac
 = inv(fac
);
for(long long i=n-1; i>=0; i--){
fac[i] = (fac[i+1] * (i+1)) % MOD;
}
}
long long get_val(long long k,long long n){
prework(k,n);
long long ans = 0;
for(long long i=1; i<=n; i++){
long long res = y[i];
res = res * (pre[i-1] * suf[i+1] % MOD)%MOD;
res = res * (fac[i-1] * fac[n-i] * (n - i & 1? -1 : 1) % MOD) % MOD;
res = (res + MOD) % MOD;
ans = (ans + res) % MOD;
}
return (ans + MOD)%MOD;
}
int main(){
//	freopen("in.txt","r",stdin);
//	freopen("out.txt","w",stdout);
long long n,k;
cin>>n>>k;
for(long long i=1; i<=k+2; i++){
y[i] = (y[i-1] + power(i,k)) % MOD;
}
printf("%lld\n",get_val(n,k+2));
return 0;
}