「2018 集训队互测 Day 1」完美的集合——树上连通块背包DP、大数阶乘取模

「2018 集训队互测 Day 1」完美的集合

题解

1.连通块背包

我们先考虑怎么求出 $S$ 的最大价值。显然有一个套路性的点分治做法，每次求包含分治中心的连通块的最大价值，用儿子继承父亲DP值的树上连通块背包即可做到 $O(nm)$（不要用论文讲的DFS序的做法，它已经被连通块背包淘汰了，因为太麻烦，尤其是在这一题），这部分复杂度是 $O(nm\log n)$，具体可以看代码。

然后我们可以发现，对于一个合法的解，$k$ 个连通块的交集一定也是个连通块。并且观察式子 $dis(x,y)\ast v_y\le MAX$，可以发现，能够对 $y$ 进行测试的所有点 $x$ 一定是一个包含 $y$ 的连通块，所以可以对 $k$ 个连通块内的点进行测试的点也构成连通块，这个连通块被包含在 $k$ 个连通块的交集以内。

所以我们对这个连通块用“点数-边数”技巧来计数，我们统计 “包含点 $x$，并且块内所有节点可以被 $x$ 测试的完美连通块数量”，设这个数量为 $F(x)$，对每条边 $(u,v)$ 统计 “包含点 $u,v$，并且块内所有节点同时可以被 $u,v$ 测试的完美连通块数量”，设这个数量为 $F(u,v)$，那么答案就是
$$\sum _{x=1}^n\left(\genfrac{}{}{0ex}{}{F(x)}{k}\right)-\sum _{(u,v)\in E}\left(\genfrac{}{}{0ex}{}{F(u,v)}{k}\right)$$

计算 $F$ 也用同样的连通块背包即可，只需要同时维护最大价值和最大值的方案数。进行连通块背包的同时要判断当前节点能否被根测试，不能则要直接返回空的DP值，因为这里已经断掉了，不能贡献连通块。

由于 $n$ 很小，只有 60，所以完美集合的数量是可以用long long存下的。这个部分的总复杂度是 $O(n^{2}m)$。

2.组合数取模

接下来考虑计算组合数。模数 $11920928955078125=5^{23}$，我们只需要把每一个阶乘的质因子 $5$ 提取出来，剩下的可以直接算逆元然后用__int128乘起来。

接下来考虑怎么快速地把一个阶乘用 $a\ast 5^b$ 的形式算出来。

我们设多项式
$$f_n(x)=\prod _{i\in [1,n],5\nmid i}(x+i)$$（我不知道为什么非得设成这个样子，但是这是解决问题的关键）
那么有如下表达式
$$n!=5^{\lfloor \frac{n}{5}\rfloor }\times (\lfloor \frac{n}{5}\rfloor !)\times f_n(0)$$
第二部分可以递归计算，于是问题变成如何快速求解 $f_n(0)$。

我们设 $n=10a+b,b<10$，那么有 $f_{10a}(0)=f_{5a}(0)\times f_{5a}(5a)$。我们可以考虑用 $f_{5a}(x)$ 求出 $f_{10a}(x)$。由于我们多项式里代的所有值都是 5 的倍数（包括 0），所以多项式只需要保留前 23 项即可。

设 $f_{5a}(x)[x^i]$ 表示 $f_{5a}(x)$ 的 $i$ 次项系数，那么在我们求出 $f_{5a}(x)$ 过后，可以很快得到 $f_{5a}(x+d)$ 的各项系数：
$$f_{5a}(x+d)[x^i]=\sum _{j\ge i}{f}_{5a}(x)[x^j]\cdot d^{j-i}\cdot \left(\genfrac{}{}{0ex}{}{j}{j-i}\right)$$
再把两个多项式暴力卷积即可。

我们通过倍增求出 $f_{10a}(x)$ 后，只需要用类似背包的方式加上 $(10a+1)\sim n$ 的值即可。这个部分的总复杂度大约是 $O(n23^2\log k)$。

代码

#include<bits/stdc++.h>//JZM yyds!!
#define ll long long
#define uns unsigned
#define IF (it->first)
#define IS (it->second)
#define END putchar('\n')
using namespace std;
const int MAXN=10005;
const ll INF=1e18;
inline ll read(){
	ll x=0;bool f=1;char s=getchar();
	while((s<'0'||s>'9')&&s>0){if(s=='-')f^=1;s=getchar();}
	while(s>='0'&&s<='9')x=(x<<1)+(x<<3)+(s^48),s=getchar();
	return f?x:-x;
}
int ptf[50],lpt;
inline void print(ll x,char c='\n'){
	if(x<0)putchar('-'),x=-x;
	ptf[lpt=1]=x%10;
	while(x>9)x/=10,ptf[++lpt]=x%10;
	while(lpt)putchar(ptf[lpt--]^48);
	if(c>0)putchar(c);
}
inline ll lowbit(ll x){return x&-x;}

const ll MOD=11920928955078125ll;
const __int128 E=1;
//math tools---
inline void ad(ll&a,ll b){a+=b;if(a>=MOD)a-=MOD;}
inline ll exgcd(ll a,ll b,ll&x,ll&y){
	if(!b){x=1,y=0;return a;}
	ll d=exgcd(b,a%b,y,x);
	y-=a/b*x;return d;
}
inline ll INV(ll a){//a must ⊥ MOD
	ll x;exgcd(a,MOD,x,*new(ll));
	return (x%MOD+MOD)%MOD;
}
inline ll ksm(ll a,ll b,ll mo){
	ll res=1;
	for(;b;b>>=1,a=E*a*a%mo)if(b&1)res=E*res*a%mo;
	return res;
}
ll gc[105][105];
inline int initgc(){
	gc[0][0]=1;
	for(int i=1;i<=100;i++){
		gc[i][0]=1;
		for(int j=1;j<=i;j++)
			gc[i][j]=(gc[i-1][j]+gc[i-1][j-1])%MOD;
	}return 0;
}
int USELESS=initgc();
#define pli pair<ll,int>
#define FI first
#define SE second
ll st[205];
int le;
const int kp=25;
unordered_map<ll,pli>mp;
inline void addnum(vector<ll>&f,ll i){
	f.push_back(0);
	for(int j=f.size()-1;j>0;j--)f[j]=(f[j-1]+E*f[j]*i)%MOD;
	f[0]=E*f[0]*i%MOD;
	if(f.size()>kp)f.resize(kp);
}
inline pli fac(ll n){
	if(mp.find(n)!=mp.end())return mp[n];//记忆化，仅为了只求一遍k的阶乘
	if(n<=10){
		pli res=pli(1,0);
		for(int i=2;i<=n;i++)res.FI*=i;
		while(res.FI%5==0)res.FI/=5,res.SE++;
		return mp[n]=res;
	}
	le=0;
	st[le=1]=n;
	while(st[le]>9)st[le+1]=st[le]/10*5,le++;//倍增预处理
	vector<ll>f,g,h;
	f.push_back(1),f.push_back(1);
	for(int i=2;i<=st[le];i++)if(i%5)addnum(f,i);
	for(int o=le-1;o>0;o--){//非递归倍增
		ll x=st[o],y=st[o+1];
		g=vector<ll>(f.size(),0);
		for(int lim=f.size()-1,i=lim;~i;i--)
			for(int j=lim;j>=i;j--)
				ad(g[i],E*f[j]*ksm(y,j-i,MOD)%MOD*gc[j][i]%MOD);
		h=vector<ll>(f.size()+g.size()-1,0);
		for(int lim=f.size()-1,i=lim;~i;i--)
			for(int j=lim;~j;j--)ad(h[i+j],E*f[i]*g[j]%MOD);
		f=h,y<<=1;
		while(y<x)if((++y)%5)addnum(f,y);
		if(f.size()>kp)f.resize(kp);
	}
	pli res=fac(n/5);
	res.SE+=n/5,res.FI=E*res.FI*f[0]%MOD;
	return mp[n]=res;
}
inline ll C(ll n,ll m){
	if(m>n||m<0)return 0;
	pli a=fac(n),b=fac(m),c=fac(n-m);
	return E*a.FI*INV(b.FI)%MOD*INV(c.FI)%MOD*ksm(5,a.SE-b.SE-c.SE,MOD)%MOD;
}
//---
struct edge{
	int v,to,w;edge(){}
	edge(int V,int T,int W){v=V,to=T,w=W;}
}e[505];
int EN=1,G[505];
inline void addedge(int u,int v,int w){
	e[++EN]=edge(v,G[u],w),G[u]=EN;
	e[++EN]=edge(u,G[v],w),G[v]=EN;
}
int n,m;
ll k,LIM,ans,wt[65],val[65];
struct itn{
	ll f,g;itn(){}
	itn(ll F,ll G){f=F,g=G;}
	inline itn operator+(const ll&b)const{
		return itn(g?f+b:0,g);
	}
}f[65][MAXN];
inline itn MAX(itn a,itn b){
	if(!a.g)a.f=0;
	if(!b.g)b.f=0;
	if(a.f^b.f)return a.f>b.f?a:b;
	else return itn(a.f,a.g+b.g);
}

bool vis[MAXN];
int ct,mn;
inline int pdfs(int x,int fa,int S){
	int siz=1,cp=0;
	for(int i=G[x];i;i=e[i].to){
		int v=e[i].v,sv;
		if(v==fa||vis[v])continue;
		sv=pdfs(v,x,S),siz+=sv,cp=max(cp,sv);
	}cp=max(cp,S-siz);
	if(cp<mn)mn=cp,ct=x;
	return siz;
}
inline void dfs(int x,int fa,ll ds,ll mx){//连通块背包
	for(int i=0;i<=m;i++)f[x][i]=itn(0,0);
	if(ds*val[x]>mx)return;		//不合法的点直接返回
	for(int i=wt[x];i<=m;i++)f[x][i]=f[fa][i-wt[x]]+val[x];//继承父节点的DP值，单点加入
	for(int i=G[x];i;i=e[i].to){
		int v=e[i].v;
		if(v==fa||vis[v])continue;
		dfs(v,x,ds+e[i].w,mx);
		for(int j=0;j<=m;j++)f[x][j]=MAX(f[x][j],f[v][j]);//把包含儿子子树的情况算上
	}
}
itn maxval;
inline void solve(int x,int S){//点分治卡常
	if(S<0)S=pdfs(x,0,0);
	ct=x,mn=1e9,pdfs(x,0,S),x=ct;
	f[0][0]=itn(0,1),dfs(x,0,0,INF);
	for(int i=0;i<=m;i++)maxval=MAX(maxval,f[x][i]);
	vis[x]=1;
	for(int i=G[x];i;i=e[i].to){
		int v=e[i].v;
		if(vis[v])continue;
		solve(v,-1);
	}
}
signed main()
{
	n=read(),m=read(),k=read(),LIM=read();
	for(int i=1;i<=n;i++)wt[i]=read();
	for(int i=1;i<=n;i++)val[i]=read();
	for(int i=1,u,v;i<n;i++)
		u=read(),v=read(),addedge(u,v,read());
	solve(1,n);
	for(int i=1;i<=n;i++)vis[i]=0;
	for(int x=1;x<=n;x++){
		f[0][0]=itn(0,1),dfs(x,0,0,LIM);
		itn dp=itn(0,0);
		for(int i=0;i<=m;i++)dp=MAX(dp,f[x][i]);
		if(dp.f==maxval.f)ad(ans,C(dp.g,k));
	}
	for(int id=2;id<=EN;id+=2){
		int x=e[id].v,y=e[id^1].v;
		if(wt[x]+wt[y]>m)continue;
		for(int i=0;i<=m;i++)f[y][i]=itn(0,0);
		f[y][0]=itn(0,1),dfs(x,y,e[id].w,LIM);
		dfs(y,x,e[id].w,LIM);
		itn dp=itn(0,0);
		for(int i=0;i<=m;i++)dp=MAX(dp,f[y][i]);
		if(dp.f==maxval.f)ad(ans,MOD-C(dp.g,k));
	}
	print(ans);
	return 0;
}