多项式插值
多项式插值
本文介绍了一些常见的插值方法和它的证明，还有一些他们对其他问题的优化。
所谓插值，就是给定一个 \(n+1\) 个点，这些点能确定一个唯一确定最高次项为 \(n\) 的多项式 \(A(x)\)，我们要求 \(A(k)\) 的值。
使用高斯消元法暴力求解多项式的时间复杂度为 \(O(n^3)\) ，但是万能的数学家们通过一些奇怪的转换发明了更为有效的插值法。
拉格朗日插值法
先定义拉格朗日基本多项式：
这个式子有比较好的性质：
如果我们带入 \(x=x_i\) ，我们惊奇的发现上式中的每个分式都为 1 ，所以答案为 1 。
如果我们带入 \(\forall x=x_j,j\ne i\) ，都会有一个分式满足 \(\frac{x_j-x_j}{x_i-x_j}\) ，我们容易发现分子为 0 ，所以答案为 0 。
可以认为：
拉格朗日插值定理是基于拉格朗日基本多项式的，设多项式为 \(L(x)\) ，\(n\) 个点为 \((x_i,y_i)\) ，我们要找到多项式在 \(k\) 点的取值。
有：
展开就是拉格朗日插值定理：
其中等号右边是 \(\sum_{i=1}^ny_i\ell_i(x)\)，它正好满足 \(f(x_i)=y_i\) 。
用幼儿园知识我们可以知道，用 \(n\) 个点可以确定一个 \(n+1\) 次多项式（没有自由元的情况）。
所以我们右边的式子满足 \(f(x_i)=y_i\) ，所以它自然和实际上的多项式是等价的。于是就能快速求值了。
这个过程是 \(O(n^2)\) 的。当然，这个方法是通法，当每个点的值是连续去到的时，我们可以在 \(O(n)\) 的时间复杂度内解决这个问题。
Ex拉格朗日插值法
当 \(x_i\) 连续取到时，我们有更优秀的性质，当 \(x_i=i\) 时，带入可得：
离线处理一些东西：
式子变为：
还是挺显然的，每次都离线处理一下就是 \(O(n)\) 了。
重心拉格朗日插值法
当我们拥有的点不断增多时，我们就需要一个在线的做法。
观察一下原来的式子：
这个式子有一些一直不变的项，设 \(g=\prod_{i=1}^nx-x_j=({\prod_{i\not = j}x-x_j})\times({x-x_i})\) ，原式变为：
再设 \(t_i=\frac{y_i}{\prod_{j=1,j\ne i}^nx_i-x_j}\) ，原式变为：
每次插入点时计算它的 \(t_i\) 即可。
拉格朗日插值法求系数
通过拉格朗日插值法我们能在 \(O(n^2)\) 的时间复杂度之内求出这个多项式的所有系数。
我们将式子转化为 ：
其中 :
\(a_i\) 和 \(g(x)\) 都是能在 \(O(n^2)\) 范围内很好预处理的。所以我们每次暴力模拟即可。
例题 P3643 [APIO2016] 划艇
题意简述：构造一个严格递增序列，有一些点可以不选，第 \(i\) 个数的范围为 \([a_i,b_i]\) ，求方案数。
数据范围：\(N\le 500,1\le a_i \le b_i \le 10^9\) .
Solution
如果你提前做一些拉格朗日插值的板子，你一定见过形如这样的式子：\(\sum_{i=1}^ni^k\) 。
有一个不太好证明的结论，这是一个 \(k+1\) 次多项式，CF622F。
它们大致有一个共同的性质：\(f_0(x)=C,f_i(x)=\sum_{j=1}^xf_{i-1}(j)\)（一开始是一个常函数，每次递推加上某段前缀和）。
对于这道题，我们简单列一下式子，设 \(f(i,j)\) 表示前 \(i\) 个元素中可选可不选，最后一个元素的值为 \(j\) ，状态转移方程为：
然后直接上拉格朗日插值即可。
代码如下：
#include<bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N=1e3+10,mod=1e9+7;
ll power(ll a,ll b){
ll ans=1;
for(;b;b>>=1){
if(b&1)ans=ans*a%mod;
a=a*a%mod;
}
return ans;
}
ll incf[N];
void init(int n=N-1){
ll mul=1;
for(int i=1;i<=n;i++)mul=mul*i%mod;
incf[n]=power(mul,mod-2);
for(int i=n-1;~i;i--)incf[i]=incf[i+1]*(i+1)%mod;
}
ll pre[N],nxt[N];
ll Lagrange(int n,ll x,ll a[]){
if(x<=n)return a[x];
pre[0]=x,nxt[n+1]=1;
ll ret=0;
for(int i=1;i<=n;i++)pre[i]=pre[i-1]*(x-i)%mod;
for(int i=n;~i;i--)nxt[i]=nxt[i+1]*(x-i)%mod;
for(int i=0;i<=n;i++){
ret=(ret+a[i]*(i>0?pre[i-1]*incf[i]%mod:1)%mod*nxt[i+1]%mod*incf[n-i]%mod*((n-i)%2?mod-1:1)%mod+mod)%mod;
}
return ret;
}
int n,L[N],R[N],b[N<<1],top;
ll f[N][N],sf[N][N],t[N],st[N][N],ans;
int main(){
init();
scanf("%d",&n);
for(int i=1;i<=n;i++){
scanf("%d%d",&L[i],&R[i]);
R[i]++;
b[++top]=L[i],b[++top]=R[i];
}
sort(b+1,b+top+1);
top=unique(b+1,b+top+1)-(b+1);
for(int i=1;i<=n;i++){
L[i]=lower_bound(b+1,b+top+1,L[i])-b;
R[i]=lower_bound(b+1,b+top+1,R[i])-b;
}
for(int i=1;i<=n;i++){
for(int j=L[i];j<R[i];j++){
ll sum=1;
for(int k=0;k<=n;k++){
t[k]=(sf[i-1][j-1]+sum)%mod;
sum=(sum+st[j][k])%mod;
st[j][k]=(st[j][k]+t[k])%mod;
}
for(int k=1;k<=n;k++)t[k]=(t[k]+t[k-1])%mod;
f[i][j]=Lagrange(n,b[j+1]-b[j]-1,t);
ans=(ans+f[i][j])%mod;
}
for(int j=1;j<top;j++)
sf[i][j]=((sf[i][j-1]+sf[i-1][j]-sf[i-1][j-1]+f[i][j])%mod+mod)%mod;
}
printf("%lld\n",ans);
return 0;
}
例题 CF995F Cowmpany Cowmpensation
题意简述：树形结构，给每个节点分配工资 \([1,d]\)，子节点不能超过父亲节点的工资，问有多少种分配方案。
数据范围： \(1\le n\le 3000 , 1\le d\le10^9\) 。
Solution
这题感觉比上题水。
和上题类似的是，我们发现 \(n\) 比较小， \(d\) 比较大。我们先把它们同时扔进状态里面。设 \(f(x,i)\) 表示以 \(x\) 为根的子树中，给 \(x\) 分配 \(y\) 的工资，子树的方案数。
我们容易写出转移：
我们先看，最下面一层的答案都为 1 ，每次往上一层我们对一堆东西求和后再乘起来。
手玩一下，从最下面一层往上一层的过程中和 \(d\) 有关，层层向上。所以我们猜测最后答案 \(f(1,d)\) 也是一个多项式，且次数不超过 \(n\) （太弱了不会严谨证明，只能感性理解一下）。
这道题和上道题的区别是，这道题每个点的限制都是一样的。所以我们直接先做一次 dp ，然后，做拉格朗日插值即可。
\(O(n^2)\) 就能过且瓶颈在于 dp ，所以我们就直接用最板的拉格朗日插值即可。
代码如下：
#include<bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N=3e3+10,M=N*2,mod=1e9+7;
ll power(ll a,ll b){
ll ans=1;
for(;b;b>>=1){
if(b&1)ans=ans*a%mod;
a=a*a%mod;
}
return ans;
}
int head[N],ver[M],nxt[M],tot=1;
void add(int x,int y){
ver[++tot]=y,nxt[tot]=head[x],head[x]=tot;
}
int n,m;
ll f[N][N];
void dfs(int x){
for(int i=1;i<=n;i++)f[x][i]=1;
for(int i=head[x];i;i=nxt[i]){
int y=ver[i];
dfs(y);
for(int j=1;j<=n;j++)
f[x][j]=f[x][j]*f[y][j]%mod;
}
for(int i=1;i<=n;i++)
f[x][i]=(f[x][i]+f[x][i-1])%mod;
}
ll x[N],y[N];
ll Lagrange(int n,ll k,ll x[],ll y[]){
ll ret=0;
for(int i=1;i<=n;i++){
ll tmp1=y[i],tmp2=1;
for(int j=1;j<=n;j++){
if(i==j)continue;
tmp1=tmp1*(k-x[j]+mod)%mod;
tmp2=tmp2*(x[i]-x[j]+mod)%mod;
}
ret=(ret+tmp1*power(tmp2,mod-2))%mod;
}
return ret;
}
int main(){
scanf("%d%d",&n,&m);
for(int i=2;i<=n;i++){
int fa;scanf("%d",&fa);
add(fa,i);
}
dfs(1);
for(int i=1;i<=n;i++)x[i]=i;
for(int i=1;i<=n;i++)y[i]=f[1][i];
printf("%d\n",Lagrange(n+1,m,x-1,y-1));
return 0;
}
感觉拉格朗日插值就很骗分，如果发现一个状态转移方程无法直接优化，又发现这个 dp 式子就相当于若干维前缀和，大胆拉格朗日插值做即可。
建议配合题单食用：拉格朗日插值优化 dp。
浙公网安备 33010602011771号