拉格朗日插值法
在数值分析中，拉格朗日插值法是以法国 $18$ 世纪数学家约瑟夫·拉格朗日命名的一种多项式插值方法。如果对实践中的某个物理量进行观测，在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。上面这样的多项式就称为拉格朗日（插值）多项式。
设 $n-1$ 次多项式 $f(x)$ 的函数图像过 $n$ 个已知点 $(x_1,y_1),(x_2,y_2),(x_3,y_3),\cdots,(x_n,y_n)$，且对于 $\forall i\neq j,x_i\neq x_j$，那么就可以确定 $f(x)$ 为一个 $n-1$ 次多项式是满足条件的。
因为可以设 $f(x)=a_{n-1}x^{n-1}+a_{n-2}x^{n-2}+\cdots+a_1x+a_0$。
那么就可以根据这 $n$ 个点列出来 $n$ 个不等价的方程，解出来即可确定唯一 $f(x)$。
但是使用高斯消元法解方程的时间复杂度为 $\mathcal O\left(n^3\right)$ 的，有没有更高效的方法呢？
那就需要使用拉格朗日（Lagrange）插值法。
拉格朗日插值法给出了 $f(x)$ 的表达式：
\[f(x)=\sum_{i=1}^ny_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j}\]正确性
考虑到拉格朗日插值法正确当且仅当 $\forall x_i,f(x_i)=y_i$。
将 $x_i$ 代入：
\[f(x_i)=\sum_{j=1}^ny_j\prod_{k\neq j}\frac{x_i-x_k}{x_j-x_k}\]若 $j\neq i$，则在 $\displaystyle y_j\prod\limits_{k\neq j}\dfrac{x_i-x_k}{x_j-x_k}$ 中必然存在 $k=i$，使得 $x_i-x_k=0$，从而使整个乘积为 $0$。
因此：
\[f(x_i)=y_i\prod_{j\neq i}\frac{x_i-x_j}{x_i-x_j}=y_i\prod_{j\neq i}1=y_i\]故，拉格朗日插值法是正确的。
横坐标为连续整数
这种特殊情况可以做到 $\mathcal O(n)$。
设 $n-1$ 多项式 $f(x)$，且已知 $f(1)=y_1,f(2)=y_2,f(3)=y_3,\cdots,f(n)=y_n$。
则有：
\[\begin{aligned} f(x)&=\sum_{i=1}^ny_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j}\\ &=\sum_{i=1}^ny_i\prod_{j\neq i}\frac{x-j}{i-j}\\ &=\sum_{i=1}^ny_i\frac{\prod\limits_{j\neq i}(x-j)}{\prod\limits_{j\neq i}(i-j)}\\ &=\sum_{i=1}^ny_i\frac{\frac{\prod\limits_{j=1}^n(x-j)}{x-i}}{\prod\limits_{j\neq i}(i-j)} \end{aligned}\]对于 $\displaystyle\prod\limits_{j=1}^n(x-j)$ 显然是可以 $\mathcal O(n)$ 处理出来的，除以 $x-i$ 即可得到分子。
而对于分母 $\displaystyle\prod\limits_{j\neq i}(i-j)$，有：
\[\prod_{j\neq i}(i-j)=(-1)^{n-i}(i-1)!(n-i)!\]证明
可以先写一下例子：
$$ \begin{aligned} \prod_{j\neq1}(1-j)&=(1-2)(1-3)(1-4)\cdots(1-n)\\ &=(-1)^{n-1}(2-1)(3-1)(4-1)\cdots(n-1)\\ &=(-1)^{n-1}\times1\times2\times3\times\cdots(n-1)\\ &=(-1)^{n-1}(n-1)!\\ \prod_{j\neq2}(2-j)&=(2-1)(2-3)(2-4)\cdots(2-n)\\ &=(-1)^{n-2}(2-1)(3-2)(4-2)\cdots(n-2)\\ &=(-1)^{n-2}\times1\times1\times2\times\cdots(n-2)\\ &=(-1)^{n-2}(n-2)!\\ \end{aligned} $$
不难发现，$i-j$ 可以转化为 $j-i$ 计算再给整体乘 $-1$。则 $i-j<0$ 有 $n-i$ 个 $j$，即整体乘 $(-1)^{n-i}$。
因此就可以拆成阶乘，$i<j$ 部分是 $(n-i)!$，而 $i>j$ 部分是 $(i-1)!$，于是便得到上式。
所以：
\[\begin{aligned} f(x)&=\sum_{i=1}^ny_i\frac{\frac{\prod\limits_{j=1}^n(x-j)}{x-i}}{(-1)^{n-i}(i-1)!(n-i)!}\\ &=\sum_{i=1}^n(-1)^{n-i}y_i\frac{\prod\limits_{j=1}^n(x-j)}{(i-1)!(n-i)!(x-i)}\\ \end{aligned}\]只需要预处理 $1\sim n$ 的阶乘及其逆元、$x-i$ 的逆元与 $\displaystyle\prod\limits_{j=1}^n(x-j)$ 即可做到单次计算复杂度 $\mathcal O(n)$。如果预处理 $x-j$ 前后缀积， 则可以避免求 $x-i$ 的逆元。
这常用于求解自然数幂次和，即 $\displaystyle\sum\limits_{i=1}^ni^k$。当 $n$ 较小时可以使用线性筛做到 $\mathcal O\left(\dfrac{n}{\ln n}\log k\right)$ 求解。（若 $n,k$ 同阶，则复杂度为 $\mathcal O(n)$）
$n$ 较大时，则可以使用拉格朗日插值求解。可以证明，$\displaystyle\sum\limits_{i=1}^ni^k$ 一定是一个关于 $n$ 的 $k+1$ 次多项式，因此线性筛找 $k+2$ 个值，插值即可。时间复杂度 $\mathcal O(k)$。
参考代码
这是经典题目CF622F The Sum of the k-th Powers的参考代码。使用维护前后缀积的写法。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
constexpr const int N=1e9,K=1e6,P=1e9+7;
int qpow(int base,int n){
int ans=1;
while(n){
if(n&1){
ans=1ll*ans*base%P;
}
base=1ll*base*base%P;
n>>=1;
}
return ans;
}
int y[K+2+1],fact[K+2+1],invFact[K+2+1];
void pre(int k){
fact[0]=1;
for(int i=1;i<=k+2;i++){
fact[i]=1ll*fact[i-1]*i%P;
}
invFact[k+2]=qpow(fact[k+2],P-2);
for(int i=k+1;i>=0;i--){
invFact[i]=invFact[i+1]*(i+1ll)%P;
}
static int vis[K+2+1],prime[K+2+1],size;
y[1]=1;
for(int i=2;i<=k+2;i++){
if(!vis[i]){
vis[i]=i;
prime[++size]=i;
y[i]=qpow(i,k);
}
for(int j=1;j<=size&&i*prime[j]<=k+2;j++){
vis[i*prime[j]]=prime[j];
y[i*prime[j]]=1ll*y[i]*y[prime[j]]%P;
if(i%prime[j]==0){
break;
}
}
}
for(int i=1;i<=k+2;i++){
y[i]=(y[i-1]+y[i])%P;
}
}
int f(int n,int k){
if(n<=k+2){
return y[n];
}
int ans=0;
static int pre[K+2+1],suf[K+2+1+1];
pre[0]=suf[k+3]=1;
for(int i=1;i<=k+2;i++){
pre[i]=1ll*pre[i-1]*(n-i)%P;
}
for(int i=k+2;1<=i;i--){
suf[i]=1ll*suf[i+1]*(n-i)%P;
}
for(int i=1;i<=k+2;i++){
int w=1;
if(k+2-i&1){
w=-1;
}
ans=(ans+w*1ll*y[i]*pre[i-1]%P*suf[i+1]%P*invFact[i-1]%P*invFact[k+2-i]%P)%P;
}
if(ans<0){
ans+=P;
}
return ans;
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
int n,k;
cin>>n>>k;
pre(k);
cout<<f(n,k)<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
例题 AC 代码
普通的拉格朗日插值即可。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
constexpr const int N=2e3,P=998244353;
struct node{
int x,y;
}a[N+1];
int n,k;
int qpow(int base,int n){
int ans=1;
while(n){
if(n&1){
ans=1ll*ans*base%P;
}
base=1ll*base*base%P;
n>>=1;
}
return ans;
}
int f(int k){
int ans=0;
for(int i=1;i<=n;i++){
int pl=a[i].y;
for(int j=1;j<=n;j++){
if(i==j){
continue;
}
pl=1ll*pl*(k-a[j].x)%P*qpow(a[i].x-a[j].x,P-2)%P;
}
ans=(1ll*pl+ans)%P;
}
if(ans<0){
ans+=P;
}
return ans;
}
main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
cin>>n>>k;
for(int i=1;i<=n;i++){
cin>>a[i].x>>a[i].y;
}
cout<<f(k)<<'\n';
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
