拉格朗日插值学习笔记
简介
在数值分析中，拉格朗日插值法是以法国18世纪数学家约瑟夫·拉格朗日命名的一种多项式插值方法。如果对实践中的某个物理量进行观测，在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。上面这样的多项式就称为拉格朗日（插值）多项式。
拉格朗日插值法
由小学知识可知 \(n+1\) 个点 \((x_i,y_i)\) 可以唯一地确定一个 \(n\) 次多项式 \(y = f(x)\)。
假设我们现在需要知道这个多项式在 \(k\) 点的取值，一个最显然的思路就是利用高斯消元直接求出多项式的系数，时间复杂度为 \(O(n^3)\)。可以用到复杂度更优的拉格朗日插值法，即：
如果把 \(x_i\) 带入公式中，除了第 \(i\) 项外的每一项分子中都有 \(x_i-x_i\)，那么最终只会留下 \(y_i\)，因此正确性可以保证。单次求值的时间复杂度为 \(O(n^2)\)。
模板题code：
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int N=2023,mod=998244353;
int n,k,x[N],y[N],ans;
void add(int &x,int y){x+=y;if(x>mod) x-=mod;}
int mul(int a,int b){int res=1;while(b) ((b&1)&&(res=1ll*res*a%mod,0)),a=1ll*a*a%mod,b>>=1;return res;}
int main()
{
scanf("%d%d",&n,&k);for(int i=1;i<=n;i++) scanf("%d%d",&x[i],&y[i]);
for(int i=1;i<=n;i++)
{
int tmp=y[i];
for(int j=1;j<=n;j++)
{
if(j==i) continue;
tmp=1ll*tmp*(k-x[j])%mod*mul((x[i]-x[j]),mod-2)%mod;
}
add(ans,(tmp+mod)%mod);
}
printf("%d\n",ans);
return 0;
}
横坐标是连续整数的拉格朗日插值
如果已知的 \(n+1\) 个点的横坐标是连续的，那么可以将单词求值的复杂度优化到 \(O(n)\)。
由于横坐标连续，假设这 \(n+1\) 个点的横坐标分别为 \(0,1,2\ldots n\)，原先的插值公式可以改写为：
对于公式中累乘的部分，可以将分母分子分开单独考虑，其中分子为：
而分母可以拆成两段阶层来算，需要注意正负符号：
于是公式就可以进一步改写为：
预处理一下就可以做到 \(O(n)\)。
教科书般的亵渎
小豆喜欢玩游戏，现在他在玩一个游戏遇到这样的场面，每个怪的血量为 \(a_i\)，且每个怪物血量均不相同，小豆手里有无限张“亵渎”。亵渎的效果是对所有的怪造成 \(1\) 点伤害，如果有怪死亡，则再次施放该法术。我们认为血量为 \(0\) 怪物死亡。
小豆使用一张“亵渎”会获得一定的分数，分数计算如下，在使用一张“亵渎”之后，每一个被亵渎造成伤害的怪会产生 \(x^k\)，其中 \(x\) 是造成伤害前怪的血量为 \(x\) 和需要杀死所有怪物所需的“亵渎”的张数 \(k\)。
\(m\leq50\)，\(n\leq10^{13},1\leq a_i \leq n\) 且 \(a_i\) 各不相同。
思路：
本题的关键在于求出多项式 \(\sum_{i=1}^{n} i^k\) 的值。
对于函数 \(f(n)=\sum_{i=1}^{n} i^k\)，可以采用归纳法证明其一定为一个 \(k+1\) 次多项式。
那么只需要取连续的 \(k+2\) 个点，就可以在 \(O(n)\) 时间内求出 \(f(n)\) 的值。
code：
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int N=66,mod=1e9+7;
int f[N],pre[N],n,m,k,a[N];
int mul(int a,int b){int res=1;while(b) ((b&1)&&(res=1ll*res*a%mod,0)),a=1ll*a*a%mod,b>>=1;return res;}
void Add(int &a,int b){a+=b;if(a>=mod) a-=mod;}
void Sub(int &a,int b){a-=b;if(a<0) a+=mod;}
int calc(int x)
{
if(x<=k+2) return f[x];
int y=0,sum=1,ans=0;
for(int i=1;i<=k+2;i++) sum=1ll*sum*(x-i)%mod;
for(int i=1;i<=k+2;i++)
{
Add(y,mul(i,k));int val=mul((1ll*(x-i)*pre[i-1]%mod*pre[(k+2)-i]%mod),mod-2);
if(k+2-i&1) Sub(ans,1ll*y*val%mod*sum%mod);else Add(ans,1ll*val*y%mod*sum%mod);
}
// printf("%d\n",ans);
return ans;
}
void solve()
{
scanf("%d%d",&n,&m);k=m+1;int ans=0;for(int i=1;i<=m;i++) scanf("%d",&a[i]),Sub(ans,mul(a[i],k));
for(int i=1;i<=k+2;i++) f[i]=(f[i-1]+mul(i,k))%mod;sort(a+1,a+m+1);Add(ans,calc(n));
for(int i=1;i<=m;i++) Add(ans,calc(n-a[i]));
for(int i=1;i<=m;i++) for(int j=i-1;j>=1;j--) Sub(ans,mul(a[i]-a[j],k));
printf("%d\n",ans);
}
int main()
{
pre[0]=1;for(int i=1;i<N;i++) pre[i]=1ll*pre[i-1]*i%mod;
int T;scanf("%d",&T);while(T--) solve();
return 0;
}
CF995F Cowmpany Cowmpensation
一棵 \(n\) 个节点的树，给每个节点分配工资（\([1, D]\)），子节点不能超过父亲节点的工资，问有多少种分配方案。
$ 1 \le n \le 3000 $ , $ 1 \le D \le 10^9 $
思路：
由于 \(n\) 的最大值很小，而 \(D\) 的最大值很大，可以考虑容斥。
除此之外，回到暴力的 dp 转移方程来，令 \(f_{i,j}\) 表示以 \(i\) 为根的子树用 \([1,i]\) 染色的方案数，可以得到状态转移方程：
那么最终的答案就是 \(f_{1,d}\)。这种做法的复杂度为 \(O(nd)\)，考虑用拉格朗日插值优化。
设 \(u\) 的子树大小为 \(s_u\)。那么，可以猜测 \(f_{u,x}\) 是关于 \(x\) 的 \(s_u\) 次函数。
证明可以用到归纳法：
-
当 \(u\) 为叶子节点时，\(f(u,x)=1\)，成立。
-
当 \(v \in son_u\) 时成立，考虑 \(u\)：
由于 \(f_{v,x}\) 的次数为 \(s_v\)，则 \(f_{u,x}-f_{u,x-1}\) 的次数为 \(\sum_{v \in son_u} s_v=s_u -1\)。再还原查分，那么 \(f_{u,x}\) 就是关于 \(x\) 的 \(s_u-1+1=s_u\) 次函数。
因此，只需要预处理出 \(f_{u,1},f_{u,2} \ldots f_(u,s_u+1)\) 的值，就可以用插值求出 \(f_{u,d}\) 的值。
code：
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int N=3010,mod=1e9+7;
int n,d,h[N],idx,f[N][N],fac[N];
struct edge{int v,nex;}e[N];
void add(int u,int v){e[++idx]=edge{v,h[u]};h[u]=idx;}
void Add(int &a,int b){a+=b;if(a>=mod) a-=mod;}
void Sub(int &a,int b){a-=b;if(a<0) a+=mod;}
void Mul(int &a,int b){a=1ll*a*b%mod;}
int mul(int a,int b){int res=1;while(b) ((b&1)&&(res=1ll*res*a%mod,0)),a=1ll*a*a%mod,b>>=1;return res;}
void dfs(int u)
{
for(int i=1;i<=n+1;i++) f[u][i]=1;
for(int i=h[u];i;i=e[i].nex)
{
int v=e[i].v;dfs(v);
for(int j=1;j<=n+1;j++) Mul(f[u][j],f[v][j]);
}
for(int i=1;i<=n+1;i++) f[u][i]=(f[u][i-1]+f[u][i])%mod;
}
int main()
{
scanf("%d%d",&n,&d);for(int v=2,u;v<=n;v++) scanf("%d",&u),add(u,v);dfs(1);
if(d<=n+1) printf("%d\n",f[1][d]);
else
{
fac[0]=1;for(int i=1;i<=n+1;i++) fac[i]=1ll*fac[i-1]*i%mod;int sum=1,y,ans=0;
for(int i=1;i<=n+1;i++) Mul(sum,d-i);
for(int i=1;i<=n+1;i++)
{
y=f[1][i];int val=1ll*y*sum%mod*mul(1ll*fac[i-1]*fac[n+1-i]%mod*(d-i)%mod,mod-2)%mod;
if(n+1-i&1) Sub(ans,val);else Add(ans,val);
}
printf("%d\n",ans);
}
return 0;
}
填树
给定一棵 \(n\) 个点的树，每个的权值 \(w_i\) 为 \(0\) 或为 \([l_i,r_i]\) 之间的整数，选定一棵树上的路径，要求满足 \(\max_{w_i}-\min_{u_i}=K\)，求这样的方案数以及权值之和。
\(1 \leq n \leq 200\)，\(1 \leq l_i \leq r_i \leq {10}^9\)，\(1 \leq K \leq {10}^9\)。
思路：
考虑钦定点权的范围为 \([w,w+K]\)，那么只需要树形dp就可以得出这一部分的答案。考虑将点权依照 \(l_i,r_i,w,w+K\) 分割成 \(O(n)\) 个区间。
根据上一题的经验，令 \(f_{u,x}\) 为点权在 \([x,x+K]\) 范围内时的方案数。那么 \(\sum f_{u,x}\) 就是一个关于 \(x\) 的 \(n+1\) 次函数、同理，权值和就是一个关于 \(x\) 的 \(n+2\) 次多项式，那么就可以利用拉格朗日插值在 \(O(n)\) 时间内求出这部分区间的答案，总的时间复杂度即为 \(O(n^3)\)。
最后还需要注意，在钦定范围为 \([w,w+K]\) 时，点权的实际取值范围可能在 \([w+1,w+K]\)，那么只需要多做一遍插值减去这部分的贡献即可。
code：
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int N=1010,mod=1e9+7;
int n,K,num[N],l[N],r[N],f[N],g[N],d1[N],d2[N],pre[N],suf[N],Y1[N],Y2[N],sig[2],ans1,ans2,h[N],idx,infac[N];
struct edge{int v,nex;}e[N];
void add(int u,int v){e[++idx]=edge{v,h[u]};h[u]=idx;}
void Add(int &a,int b){a+=b;if(a>=mod) a-=mod;if(a<0) a+=mod;}
void dfs(int u,int fa,int L,int R)
{
f[u]=g[u]=d1[u]=d2[u]=0;
int ql=max(L,l[u]),qr=min(R,r[u]);
int cntu=max(qr-ql+1,0),sumu=(ql<=qr)?1ll*(ql+qr)*(qr-ql+1)/2%mod:0;
f[u]=d1[u]=cntu,g[u]=d2[u]=sumu;
for(int i=h[u];i;i=e[i].nex)
{
int v=e[i].v;if(v==fa) continue;dfs(v,u,L,R);
(f[u]+=f[v])%=mod;(g[u]+=g[v])%=mod;
g[u]=(g[u]+1ll*d2[u]*d1[v]+1ll*d1[u]*d2[v])%mod;
d2[u]=(d2[u]+1ll*d2[v]*cntu+1ll*d1[v]*sumu)%mod;
f[u]=(f[u]+1ll*d1[u]*d1[v])%mod;d1[u]=(d1[u]+1LL*d1[v]*cntu)%mod;
}
}
void solve(int w)
{
int m=0;num[++m]=w+1,num[++m]=1e9+1;
for(int i=1;i<=n;i++)
{
if(l[i]>w) num[++m]=l[i];
if(l[i]>w) num[++m]=r[i];
if(l[i]>K) num[++m]=l[i]-K+w;
if(r[i]>K) num[++m]=r[i]-K+w;
}
sort(num+1,num+m+1);m=unique(num+1,num+m+1)-num-1;
for(int i=1;i<m;i++)
{
int k=min(num[i+1]-num[i],n+3),x=num[i+1]-num[i];pre[0]=suf[k+1]=1;
for(int j=1;j<=k;j++) dfs(1,0,num[i]+j-1,num[i]+j+K-w-1),Y1[j]=f[1],Y2[j]=g[1];
for(int j=2;j<=k;j++) Add(Y1[j],Y1[j-1]),Add(Y2[j],Y2[j-1]);
for(int j=1;j<=k;j++) pre[j]=1ll*pre[j-1]*(x-j+mod)%mod;
for(int j=k;j>=1;j--) suf[j]=1ll*suf[j+1]*(x-j+mod)%mod;
for(int j=1;j<=k;j++)
{
int res=1ll*pre[j-1]*suf[j+1]%mod*infac[j-1]%mod*infac[k-j]%mod*sig[k-j&1]%mod;
Add(ans1,1ll*sig[w]*res%mod*Y1[j]%mod);Add(ans2,1ll*sig[w]*res%mod*Y2[j]%mod);
}
}
}
int main()
{
// freopen("tree3.in","r",stdin);
scanf("%d%d",&n,&K);for(int i=1;i<=n;i++) scanf("%d%d",&l[i],&r[i]);
for(int u,v,i=1;i<n;i++) scanf("%d%d",&u,&v),add(u,v),add(v,u);infac[0]=infac[1]=sig[0]=1;sig[1]=mod-1;
for(int i=2;i<N;i++) infac[i]=1ll*(mod-mod/i)*infac[mod%i]%mod;
for(int i=2;i<N;i++) infac[i]=1ll*infac[i-1]*infac[i]%mod;
solve(0),solve(1);printf("%d\n%d\n",ans1,ans2);
return 0;
}
浙公网安备 33010602011771号