拉格朗日插值

给出 oi-wiki 中插值的定义:

插值是一种通过已知的、离散的数据点推算一定范围内的新数据点的方法。插值法常用于函数拟合中。

拉格朗日插值是一种多项式插值。

给定 nn 个点 (xi,yi)(x_i,y_i),若满足 ∀i≠j,xi≠xj\forall i\not=j,x_i\not =x_j,那么这 nn 个点可以唯一确定一个 n−1n-1 次多项式 y=f(x)y=f(x),求这个多项式 f(x)f(x),并给定参数 kk,求 f(k) mod 998244353f(k) \bmod 998244353 的值。

我们考虑 ∀i∈[1,n]\forall i\in [1,n] 构造 fi(x)f_i(x) 满足:

$$f_i(x_j)=\begin{cases}y_j,j=i\\ 0,j\not=i \end{cases}$$

此时 f(x)=∑i=1nfi(x)f(x)=\sum_{i=1}^nf_i(x),满足 ∀i∈[1,n],f(xi)=yi\forall i\in [1,n],f(x_i)=y_i。

我们尝试构造 fi(x)f_i(x):

$$f_i(x)=\frac{y_i}{\prod_{j\not=i}(x_i-x_j)}\times \prod_{j\not=i}(x-x_j)\\ f_i(x)=y_i\times \frac{\prod_{j\not=i}(x-x_j)}{\prod_{j\not=i}(x_i-x_j)}\\ =y_i\times \prod_{j\not=i}\frac{x-x_j}{x_i-x_j}$$

所以 $f(x)=\displaystyle\sum_{i=1}^n y_i\times \prod_{j\not=i}\dfrac{x-x_j}{x_i-x_j}$。

附上模板题代码:

#include<bits/stdc++.h>
#define FL(i,a,b) for(int i=(a);i<=(b);i++)
#define FR(i,a,b) for(int i=(a);i>=(b);i--)
#define ll long long
#define ull unsigned long long
#define ld long double
#define PII pair<int,int>
using namespace std;
const int MAXN = 2e3 + 10;
const int mod = 998244353;

int n,K,ans=0;
int x[MAXN],y[MAXN];

int qpow(int a,int b){
	int res=1;
	while(b){
		if(b&1) res=1ll*res*a%mod;
		a=1ll*a*a%mod;
		b>>=1;
	}
	return res;
}

int main(){
	scanf("%d%d",&n,&K);
	FL(i,1,n) scanf("%d%d",&x[i],&y[i]);
	FL(i,1,n){
		ll p=y[i],q=1;
		FL(j,1,n){
			if(i!=j){
				p=1ll*p*(K-x[j])%mod;
				q=1ll*q*(x[i]-x[j])%mod;
			}
		}
		ans=(ans+1ll*p*qpow(q,mod-2)%mod)%mod;
	}
	printf("%d\n",(ans+mod)%mod);
}

CF1874E Jellyfish and Hack

感觉这道题的 trick 还是很新奇的。

考虑 dp,设 dpx,ydp_{x,y} 表示长度为 xx 且 fun(P)=yfun(P)=y 的排列 PP 的数量。

枚举 ii,把序列划分为长度 i−1i-1 和 x−ix-i 的两个子序列(即 A1A_1 是 AA 序列中第 ii 大的数)。

发现合法的 yy 的量级是 O(n2)\mathcal{O}(n^2) 的。

转移方程如下:

$dp_{x,y}=\sum_{i=1}^x\binom{x-1}{i-1}\sum_{j=0}^{y-x}dp_{i-1,j}\times dp_{x-i,y-x-j}$。

朴素转移是 O(n6)\mathcal{O}(n^6) 的,考虑优化。

我们令 $F_n(t)=\displaystyle \sum_{i=0}^{\infty}dp_{n,i}t^i$,则有:

$$F_n(t)=\sum_{v=0}^{\infty}\sum_{i=1}^n\binom{n-1}{i-1}\sum_{j=0}^{v-n}dp_{i-1,j}\times dp_{n-i,v-n-j}\times t^v\\ F_n(t)=\sum_{i=1}^n\binom{n-1}{i-1}\sum_{v=0}^{\infty}\sum_{j=0}^{v-n}dp_{i-1,j}\times dp_{n-i,v-n-j}\times t^v\\ 令 k=v-n,\\ F_n(t)=\sum_{i=1}^n\binom{n-1}{i-1}\sum_{v=0}^{\infty}\sum_{j=0}^{k}dp_{i-1,j}\times dp_{n-i,k-j}\times t^{n+k}\\ F_n(t)=\sum_{i=1}^n\binom{n-1}{i-1}\sum_{k=0}^{\infty}(\sum_{j=0}^{k}dp_{i-1,j}\times dp_{n-i,k-j}\times t^k)t^n\\ F_n(t)=\sum_{i=1}^n\binom{n-1}{i-1}F_{i-1}(t)F_{n-i}(t)t^n$$

以及多项式乘法本质上就是序列的卷积。

然后这个多项式的 deg⁡\deg 是 O(n2)\mathcal{O}(n^2) 的(多项式的度(deg⁡\deg)是指该多项式中最高次项的次数)

我们考虑带入 n2+εn^2 + \varepsilon 个值来求出这个多项式的系数。

我们枚举 n,tn,t,预处理 tt 的幂次,可以 O(n4)\mathcal{O}(n^4) 算出 Fn(t)(t∈[1,lim+ε])F_{n}(t)(t\in [1,lim+\varepsilon ])。

我们回忆拉格朗日插值的公式,然后对于每个 tt 代入 Fn(t)F_n(t),用拉格朗日插值可以 O(n2)\mathcal{O}(n^2) 算出 FnF_n 的各项系数,即 dpn,i(i∈[1,n2+ε])dp_{n,i}(i\in [1,n^2+\varepsilon])。

$F_n(x)=\displaystyle\sum_{i=1}^{n^2+\varepsilon} y_i\times \prod_{j\not=i}\dfrac{x-x_j}{x_i-x_j}$。

我们令 $a_i=y_i\times \displaystyle\prod_{j\not=i} \frac{1}{x_i-x_j},L(x)=\prod_{j=1}^n(x-x_j)$,则 $F_n(x)=\displaystyle\sum_{i=1}^{n^2+\varepsilon} a_i\times \frac{L(x)}{x-x_i}$。

我们计算 L(x)L(x) 中 ii 次项的系数 sis_i,然后对每个数 ii,做多项式除法(除以 x−xix-x_i,然后因为一定是能整除的所以可以正着做)求出 ii 次项的系数 bib_i,于是 Fn(x)F_n(x) 第 jj 项的系数(即 dpn,jdp_{n,j})加上 ai×bja_i\times b_j。

最后答案即为 $\displaystyle \sum_{i=lim}^{n^2+\varepsilon}dp_{n,i}$。

#include<bits/stdc++.h>
#define FL(i,a,b) for(int i=(a);i<=(b);i++)
#define FR(i,a,b) for(int i=(a);i>=(b);i--)
#define ll long long
#define ull unsigned long long
#define ld long double
#define PII pair<int,int>
using namespace std;
const int MAXN = 2e2 + 10;
const int mod = 1e9 + 7;
const int eplison = 5;

int n,lim,Mx,ans=0;
int C[MAXN][MAXN],pw[MAXN*MAXN][MAXN];
int x[MAXN*MAXN],y[MAXN*MAXN];
int F[MAXN][MAXN*MAXN],cnt[MAXN*MAXN];
int a[MAXN*MAXN],s[MAXN*MAXN],b[MAXN*MAXN];

int qpow(int a,int b){
	int res=1;
	while(b){
		if(b&1) res=1ll*res*a%mod;
		a=1ll*a*a%mod;
		b>>=1;
	}
	return res;
}

void Solve(){
	FL(i,1,Mx) x[i]=i,y[i]=F[n][i];
	FL(i,1,Mx){
		int Mul=1;
		FL(j,1,Mx)
			if(i!=j) Mul=1ll*Mul*(x[i]-x[j]+mod)%mod;
		a[i]=1ll*y[i]*qpow(Mul,mod-2)%mod;
	}
	s[0]=1;
	FL(i,1,Mx){
		FR(j,i,1)
			s[j]=(s[j-1]-1ll*x[i]*s[j]%mod+mod)%mod;
		s[0]=1ll*(mod-s[0])*x[i]%mod;
	}
	FL(i,1,Mx){
		int inv=qpow(mod-x[i],mod-2);
		b[0]=1ll*s[0]*inv%mod;
		cnt[0]=(cnt[0]+1ll*a[i]*b[0]%mod)%mod;
		FL(j,1,Mx){
			b[j]=1ll*(s[j]-b[j-1]+mod)*inv%mod;
			cnt[j]=(cnt[j]+1ll*a[i]*b[j]%mod)%mod;
		}
	}
}

int main(){
	scanf("%d%d",&n,&lim),Mx=n*(n+1)/2;
	if(lim>Mx) puts("0"),exit(0);
	Mx+=eplison;
	FL(i,1,Mx){
		pw[i][0]=1;
		FL(j,1,n) pw[i][j]=1ll*pw[i][j-1]*i%mod;
	}
	FL(i,0,MAXN-1){
		C[i][0]=C[i][i]=1;
		FL(j,1,i-1) C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
	}
	FL(i,1,Mx) F[1][i]=i,F[0][i]=1;
	FL(i,2,n){
		FL(j,1,Mx){
			FL(k,1,i)
				F[i][j]=(F[i][j]+1ll*C[i-1][k-1]*F[k-1][j]%mod*F[i-k][j]%mod)%mod;
			F[i][j]=1ll*F[i][j]*pw[j][i]%mod;
		}
	}
	Solve();
	FL(i,lim,Mx) ans=(ans+cnt[i])%mod;
	printf("%d\n",ans);
} 
省集题目(无题号)

给定 nn 个节点的树,每个节点有权值 xi∈[Li,Ri]x_i\in[L_i,R_i],终点在树的一个叶子节点,但是对你未知。

你需要在走到终点的时候恰好选择 kk 个节点,且选择的节点的 xix_i 按选择顺序单调不增。

但是对于每个节点 uu,只有在你对其进行是否选择的决策后你才能得知终点在 uu 的哪个儿子的子树中。

所以你会选择一个节点集合 SS,在经过集合中的节点时必然选择这个节点,在经过集合外的节点时必然不选择这个节点。

并且保证每个叶子节点作为终点的情况都能满足条件。

问对于所有可能的权值序列 xx,选取节点集合的方案数的总和,对 10045358091004535809 取模。

n≤200,k≤20,Ri≤109n\le 200,k\le 20,R_i\le 10^9。

我们称在集合 SS 中的节点为关键点,要求根到每个叶子的路径都满足恰好经过 kk 个关键点,且关键点权值按路径顺序单调不增。

我们考虑朴素树形 dp,令 fu,i,jf_{u,i,j} 表示在 uu 子树内给每个节点赋值并选点,满足 uu 到每个叶子的路径上恰好选 ii 个点,且所有关键点权值 ≤j\le j 的方案数。

我们令 gu,i,j=∏v∈son(u)fv,i,jg_{u,i,j}=\prod_{v\in son(u)}f_{v,i,j},考虑 uu 是否选入集合:

  • 如果不选 uu,那么 xux_u 不影响 jj 的值,可以随便填,贡献为 (Ru−Lu+1)⋅gu,i,j(R_u-L_u+1)\cdot g_{u,i,j}。
  • 否则令 xu=y∈[Lu,Ru]x_u=y\in [L_u,R_u],贡献为 ∑y=Lumin⁡(Ru,j)gu,i−1,y\sum_{y=L_u}^{\min(R_u,j)}g_{u,i-1,y}。

因此 fu,i,jf_{u,i,j} 的转移式子是:

$$f_{u,i,j}=(R_u-L_u+1)\cdot g_{u,i,j}+\sum_{y=L_u}^{\min(R_u,j)}g_{u,i-1,y}$$

时间复杂度 O(nkV2)\mathcal{O}(nkV^2),可以前缀和优化到 O(nkV)\mathcal{O}(nkV),其中 V=109V=10^9。

发现 jj 的取值范围是 10910^9 的,直接枚举显然不可行,所以考虑离散化,把 Li,Ri+1L_i,R_i+1 放进离散化数组里得到 tottot 个不同的数,中间有 tot−1tot-1 段。

发现同一段中的数,和各个节点的取值范围的包含关系不变。

考虑对每段分别处理,设当前在处理第 XX 段,也就是 [bX,bX+1)[b_X,b_{X+1}),发现 ff 是关于 p=j−bXp=j-b_X 是一个次数不超过 nn 的多项式 F(p)F(p)(这个可以从叶子开始推,每个点的选择方案是带 pp 的式子,显然乘起来是多项式),于是我们求出其在 p=0∼n+1p=0\sim n+1 的点值然后用拉格朗日插值还原多项式并求出整段的答案。

具体地,我们固定 XX,维护 Fu,i(p)=fu,i,bX+p(p∈[0,n+1])F_{u,i}(p)=f_{u,i,b_X+p}(p\in [0,n+1]),定义 hisu,ihis_{u,i} 表示在处理当前段 XX 之前,所有已经处理过且属于 uu 的权值范围的 yy 的 gu,i,yg_{u,i,y} 的和。

对于节点 uu:

  • 如果 xux_u 的取值范围覆盖这一段,那么对于上界 j=bX+pj=b_X+p,选 uu 的贡献就是 Su,i(p)=hisu,i+∑q=0pgu,i,bX+qS_{u,i}(p)=his_{u,i}+\sum_{q=0}^pg_{u,i,b_X+q}。

  • 否则选 uu 的贡献只有以前贡献和 hisu,ihis_{u,i}。

转移后,我们需要把这段的新值加入 hisu,ihis_{u,i}:

我们已经得到 p=0∼n+1p=0\sim n+1 时 SS 的值,用拉格朗日插值可以求出 SS,然后带入 p=bX+1−bX−1p=b_{X+1}-b_X-1,此时多项式等于 $S_{u,i}(p)=his_{u,i}+\sum_{q=0}^{b_{X+1}-b_X-1}g_{u,i,b_X+q}$,此时对应的 j=bX+1−1j=b_{X+1}-1,正好是我们想求的 hisu,ihis_{u,i} 的新值。

然后最后用同样的方法还原答案,即令 X=tot−1X=tot-1,求 F1,k(btot−btot−1−1)=f1,k,btot−1F_{1,k}(b_{tot}-b_{tot-1}-1)=f_{1,k,b_{tot}-1} 的值,此时对应的 j=btot−1j=b_{tot}-1,也就是还没有选择任何点且没有任何取值限制时的答案。

时间复杂度 O(n3K)\mathcal{O}(n^3K)。

#include<bits/stdc++.h>
#define FL(i,a,b) for(int i=(a);i<=(b);i++)
#define FR(i,a,b) for(int i=(a);i>=(b);i--)
#define ll long long
#define ull unsigned long long
#define ld long double
#define PII pair<int,int>
using namespace std;
const int MAXN = 2e2 + 10;
const int MAXK = 20 + 5;
const int mod = 1004535809;

int n,K,X;
int L[MAXN],R[MAXN],sum[MAXN];
int f[MAXN][MAXK][MAXN],g[MAXK][MAXN];
int his[MAXN][MAXK];
vector<int>G[MAXN];
int b[MAXN<<1],tot=0;

int qpow(int a,int b){
	int res=1;
	while(b){
		if(b&1) res=1ll*res*a%mod;
		a=1ll*a*a%mod;
		b>>=1;
	}
	return res;
}

int suf[MAXN],pre[MAXN],exc[MAXN];

void init(){
	FL(i,0,n+1){
		int Mul=1;
		FL(j,0,n+1)
			if(j!=i) Mul=1ll*Mul*(i-j+mod)%mod;
		exc[i]=qpow(Mul,mod-2);
	}
}

int Get(int *f,int ps){
	int R=b[ps+1]-b[ps]-1;
	suf[n+1]=(R-(n+1)+mod)%mod;
	FR(i,n,0) suf[i]=1ll*suf[i+1]*(R-i+mod)%mod;
	pre[0]=R;
	FL(i,1,n+1) pre[i]=1ll*pre[i-1]*(R-i+mod)%mod;
	int Mul=1,res=0;
	FL(i,0,n+1){
		res=(res+1ll*f[i]*exc[i]%mod*(i?pre[i-1]:1)%mod*(i+1<=n+1?suf[i+1]:1)%mod)%mod;
		Mul=1ll*Mul*(R-i+mod)%mod; 
	}
	return res;
}

void dfs1(int u,int fth){
	for(int v:G[u]){
		if(v==fth) continue;
		dfs1(v,u),sum[u]++;
	}
}
void dfs2(int u,int fth){
	for(int v:G[u]){
		if(v==fth) continue;
		dfs2(v,u);	
	}
	FL(j,0,K) FL(k,0,n+1) g[j][k]=(!sum[u]&&j?0:1);
	for(int v:G[u]){
		if(v==fth) continue;
		FL(j,0,K) FL(k,0,n+1) g[j][k]=1ll*g[j][k]*f[v][j][k]%mod;
	}
	int len=b[R[u]+1]-b[L[u]];
	FL(j,0,K) FL(k,0,n+1) f[u][j][k]=1ll*len*g[j][k]%mod;
	FL(j,0,K-1){
		int prv=his[u][j];
		if(L[u]<=X&&X<=R[u]){
			int sum=prv;
			FL(k,0,n+1){
				sum=(sum+g[j][k])%mod;
				f[u][j+1][k]=(f[u][j+1][k]+sum)%mod,g[j][k]=sum;
			}
			his[u][j]=Get(g[j],X);
		}
		else FL(k,0,n+1) f[u][j+1][k]=(f[u][j+1][k]+prv)%mod;
	}
}

int main(){
	scanf("%d%d",&n,&K);
	FL(i,1,n) scanf("%d",&L[i]);
	FL(i,1,n) scanf("%d",&R[i]);
	FL(i,1,n-1){
		int u,v;
		scanf("%d%d",&u,&v);
		G[u].push_back(v);
		G[v].push_back(u);
	}
	FL(i,1,n) b[++tot]=L[i],b[++tot]=R[i]+1;
	sort(b+1,b+tot+1);
	tot=unique(b+1,b+tot+1)-b-1;
	FL(i,1,n)
		L[i]=lower_bound(b+1,b+tot+1,L[i])-b,
		R[i]=lower_bound(b+1,b+tot+1,R[i]+1)-b-1;
	init();
	dfs1(1,0); 
	FL(i,1,tot-1) X=i,dfs2(1,0);
	printf("%d\n",Get(f[1][K],tot-1));
}