拉格朗日插值优化 $dp$
拉格朗日插值优化 \(dp\)
题单
拉格朗日插值
对于一个 \(n+1\) 次多项式 \(f(x)\),若它经过 \(n\) 个点 \((x_1,y_1)\sim(x_n,y_n)\)。
设其基本多项式:
你会发现这个多项式经过 \((x_i,y_i)\) 和 \(\forall j\ne i\) \((x_j,0)\),我们把所有这些加起来就会得到经过 \(n\) 个点的多项式 \(f(x)\):
for(int i=1;i<=n;++i){
int s1=y[i],s2=1;
for(int j=1;j<=n;++j){
if(j==i) continue;
s1=s1*(k-x[j])%mod;
s2=s2*(x[i]-x[j])%mod;
}
ans+=s1*get_c(s2)%mod;
ans=(ans+mod)%mod;
}
\(\text{useful small conclusion}\)
- 若 \(p(x)\) 是一个关于 \(x\) 的 \(n\) 次函数,则 \(p(x)-p(x-1)\) 是关于 \(x\) 的 \(n-1\) 次函数,即对一个多项式作差分会导致它次数减一,而作前缀和会导致它次数加一。
- 一个 \(n\) 次多项式乘一个 \(m\) 次多项式会得到一个 \(n+m\) 次多项式。
CF995F Cowmpany Cowmpensation
设 \(f_{x,i}\) 表示以 \(x\) 为根的子树用值域 \([1,i]\) 染色的方案数,则有:
直接 \(dp\) 复杂度是 \(O(n\cdot d)\) 的,
那么我们就需要一个重要的结论来进行一个优化。
重要结论
设 \(x\) 的子树大小为 \(sz_x\),则 \(f_{x,i}\) 是关于 \(i\) 的 \(sz_x\) 次函数。
我们考虑采用归纳法证明这一点:
- 对于叶子节点,\(f_{x,i}=1\),成立。
- 对于非叶子节点,考虑将多项式差分一下:
若 \(y\) 满足结论,则 \(f(x,i)-f(x,i-1)\) 的次数为 \(\sum_{y\in son(x)} sz_y=sz_x-1\),还原回去就是 \(sz_x\)。
所以我们可以求出过 \((x,f(1,x))\) 的多项式方程并用拉格朗日插值插出 \(x=d\) 处的值。
可以使用连续插值但是没有必要。
总体复杂度是 \(O(n^2)\) 的
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N=3e3+5,mod=1e9+7;
int n,d;
vector <int> G[N];
int f[N][N],x[N],y[N];
inline int ksm(int x,int idx){
int res=1;
for(;idx;idx>>=1,x=x*x%mod)
if(idx&1)
res=res*x%mod;
return res;
}
inline int inv(int x){
return ksm(x,mod-2);
}
inline int lag_range(){
int res=0;
for(int i=0;i<=n;++i){
int t=y[i],div=1;
for(int j=0;j<=n;++j){
if(i!=j){
t=t*(d-x[j]+mod)%mod;
div=div*(x[i]-x[j]+mod)%mod;
}
}
res=(res+t*inv(div)%mod)%mod;
}
return res;
}
inline void dfs(int x){
for(int i=1;i<=n;++i)
f[x][i]=1;
for(int y:G[x]){
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;
}
signed main(){
cin>>n>>d;
for(int i=0;i<=n;++i) x[i]=i;
for(int i=2,fa;i<=n;++i){
cin>>fa;
G[fa].push_back(i);
}
dfs(1);
for(int i=1;i<=n;++i)
y[i]=f[1][i];
if(d<=n) cout<<f[1][d]<<endl;
else{
cout<<lag_range()<<endl;
}
}
P4463 calc
我们设 \(f_{i,j}\) 表示前 \(i\) 个数使用的数值域为 \([1,j]\) 的不同合法序列的值的和。
由于必须乱序,这里我们可以只用 \(dp\) 序列中数递增的情况,最后总答案就会乘一个 \(n!\),也就变成了 \(f_{n,k}\cdot n!\)。
\(dp\) 式子就长这样:
显然这样 \(dp\) 是会超时的,我们需要考虑优化。
我们设 \(g_{i,j}=f_{i,j}-f_{i,j-1}\) (差分),则有:
若 \(g_{n,x}\) 可以用一个关于 \(x\) 的多项式来表达,则 \(f_{n,k}\) 也可以,因为对 \(g_{n,x}\) 求前缀和相当于给这个多项式的次数加一
考虑用数学归纳法算出 \(g_{n,x}\) 是一个次数为几的多项式。
- 当 \(n=0\) 时,有 \(g_{n,x}\) 次数为 \(0\)。
- 每次 \(dp\) 会求一次 \(g\) 的前缀和并给 \(g_{n,x}\) 乘上一个 \(x\),次数相当于加了 \(2\)。
- 所以说 \(g_{n,x}\) 就是一个关于 \(x\) 的 \(2\cdot n\) 次多项式,而 \(f_{n,x}\) 就是一个关于 \(x\) 的 \(2\cdot n+1\) 次多项式。
然后使用拉格朗日插值搞出 \(f_{n,k}\) 就完事了,时间复杂度 \(O(4\cdot n^2)\)
#include<bits/stdc++.h>
using namespace std;
const int N=505;
using ll=long long;
int n,m,k,p;
ll f[N][N*2],x[N*2],y[N*2];
inline int ksm(ll x,int idx){
ll res=1;
for(;idx;x=x*x%p,idx>>=1)
if(idx&1) res=res*x%p;
return res;
}
inline int inv(int x){return ksm(x,p-2);}
inline ll lag_range(int k){
if(k<=m)
return y[k];
ll res=0;
for(int i=1;i<=m;++i){
ll mul=y[i],div=1;
for(int j=1;j<=m;++j){
if(i==j) continue;
mul=mul*(k-x[j]+p)%p;
div=div*(x[i]-x[j]+p)%p;
}
res=(res+mul*inv(div)%p)%p;
}
return res;
}
signed main(){
cin>>k>>n>>p;
m=2*n+1;
ll frac=1;
for(int i=1;i<=n;++i)
frac=frac*i%p;
for(int i=0;i<=m;++i) f[0][i]=1;
for(int i=1;i<=n;++i)
for(int j=1;j<=m;++j)
f[i][j]=(f[i-1][j-1]*j%p+f[i][j-1])%p;
for(int i=1;i<=m;++i)
y[i]=f[n][i],x[i]=i;
printf("%d\n",lag_range(k)*frac%p);
}
\(\text{tmp Tiny Summary}\)
拉格朗日插值优化 \(dp\) 一般是建立在有关值域方面的 \(dp\) 问题,面对这类问题我们通常先推导出暴力 \(dp\) 式子,然后再想办法将 \(dp\) 式子表示成关于某一个数的次数为 \(k\) 的多项式,最后通过求出 \(k+1\) 个点来插值。
若次数为 \(k\) 的话最少要使用 \(k+1\) 个点,但是也可以使用更多个点 (准没错),所以推不出来 \(k\) 却又猜出是拉格朗日插值的题可以多跑几个点,也就是可以根据数据范围来计算求几个点插值。
而且这种方式保险且能保证正确性。