拉格朗日插值学习笔记

简介

在数值分析中，拉格朗日插值法是以法国18世纪数学家约瑟夫·拉格朗日命名的一种多项式插值方法。如果对实践中的某个物理量进行观测，在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。上面这样的多项式就称为拉格朗日（插值）多项式。

拉格朗日插值法

由小学知识可知 $n+1$ 个点 $(x_i,y_i)$ 可以唯一地确定一个 $n$ 次多项式 $y = f(x)$。

假设我们现在需要知道这个多项式在 $k$ 点的取值，一个最显然的思路就是利用高斯消元直接求出多项式的系数，时间复杂度为 $O(n^3)$。可以用到复杂度更优的拉格朗日插值法，即：

$$f(k)=\sum_{i=0}^{n} y_i\prod_{i \ne j} \dfrac{k-x[j]}{x[i]-x[j]}$$

如果把 $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$，原先的插值公式可以改写为：

$$f(x)=\sum_{i=0}^{n} y_i\prod_{j \ne i} \dfrac{x-j}{i-j}$$

对于公式中累乘的部分，可以将分母分子分开单独考虑，其中分子为：

$$\prod_{j \ne i} (x-j) =\prod_{i=0}^{n} \dfrac{x-j}{x-i}$$

而分母可以拆成两段阶层来算，需要注意正负符号：

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

于是公式就可以进一步改写为：

$$f(x)=\sum_{i=0}^{n} y_i\dfrac{\prod_{i=0}^{n} (x-j) }{(x-i)\times (-1)^{n-i} \times i! \times (n-i)!}$$

预处理一下就可以做到 $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_{u,i}=f_{u,i-1}+\prod_{v \in son_u} f_{v,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_{u,x}-f_{u,x-1}=\prod_{v \in son_u} f_{v,x}$$

由于 $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;
}