【题解】LOJ2462完美的集合(树DP 魔改Lucas)
【题解】LOJ2462完美的集合(树DP 魔改Lucas)
省选模拟考这个？？？？？？？？？？？？？？？？？？？
题目大意:
有一棵树，每个点有两个属性，一个是重量\(w_i\)一个是价值\(v_i\)。我们称一个点集\(S\)合法当且仅当
- 该集合是一个联通块\(\qquad (1)\)
- 该集合的所有点的重量和\(\le m\)，输入中给定\(m\),\(\qquad (2)\)
- 该集合的所有点的价值和是(全局)所有可能的价值和中最大的\(\qquad (3)\)
我们称一个点集的集合\(B=\{S_i\}\)合法当且仅当
- \(|B|=k\)，输入中给定\(k\),\(\qquad (4)\)
- 存在一个点\(x\)，对于所有\(S\in B\)有\(x \in S\)。\(\qquad (5)\)
- 对于上述点\(x\)，存在一个\(x\)对于所有\(S\in B\)，对于所有\(y\in S\)有\(\mathrm{dis}(x,y)\times v_y\le Max\)。输入中给定\(Max\)。\(\qquad (6)\)
请输出不同的\(B\)的个数，答案对\(5^{23}\)取模。
\(n\le 60,m\le 10000,k,w_i,v_i\le 10^9,Max\le 10^{18}\)
分两个部分解决问题因为这道题是二合一。
考虑找到所有包含\(x\)的\(S\)们。这些\(S\)的并集在树上构成了一个联通块，由于我们确定了一个\(x\)，所以每个点\(y\)是否在\((6)\)中合法是确定的，因此我们从原树上扣出一个和\(x\)联通的联通块(记为树\(T_x\))，在这个上面找到所有的\(S\)。
沿用这个博客t2的一些方法https://www.cnblogs.com/winlere/protected/p/11788856.html
外校的同学可能打不开，这里摘抄过来
定义一下二元组的运算:
\(e_1=(x_1,y_1),e_2=(x_2,y_2)\)，设\(u=\max\{x_1,x_2\}\)
\[e_1+e_2=(u,[x_1=u]y_1+[x_2=u]y_2) \]值得注意的是这个东西满足结合律和交换律
我们枚举一个\(x\)并且令其为根，设\(dp(i,j)=(a,b)\)，表示\(i\)是\(S\)中最浅的点，\(S\)的重量和是\(j\)，且价值和是\(a\)，这样的\(S\)的方案数是\(b\)。(因为x是根，也就是当前最浅的点，所以\(dp(x,j)\)的含义就是树上所有包含\(x\)的合法的\(S\)的情况了)
二元组的运算法则和链接里面那道题是一样的。
所以我们到此满足了\((1),(2),(3),(5),(6)\)的限制，条件\((3)\)关于全局的限制可以把所有\(x\)枚举完之后得到最大值再统计答案。
考虑如何转移\(dp\)，直接树上\(DP\)的复杂度是\(O(nm^2)\)的。这是因为在儿子转移父亲的时候做了一个完全背包，但是我们实际上是01背包，然后题解给出了一个办法可以做成01背包
状态改为\(dp(i,j)\)表示考虑前\(i\)个dfs序,选择的联通块的价值为j，转移顺序改变一下，按照\(T_x\)的dfs序的倒序转移，转移是(注意这里的加法是上面定义的!):
- 选择该dfs序的点\(dp(i,j)=dp(i,j)+ dp(i+1,j-w_{now})\)，其中\(now\)表示当前点的编号。(记得更新二元组的x)
- 不选择该dfs序的点\(dp(i,j)=dp(i,j)+dp(i+siz[now],j)\)，\(now\)的意思同上。\(+siz[now]\)表示跳过一整个子树(因为这个点我们不选，又由于\(S\)是要和\(x\)(根)联通，所以要跳过整颗子树。这样子DP对于一个点，要么选，要么不选它的整颗子树，可以保证最终\(dp\)到根上的方案都是合法的，相当于从最近的考虑过了的点转移过来)
这样我们就做了一个01背包了，转移只要for一个j，不需要for第二个j。复杂度\(O(nm)\)可以接受
现在我们可以得到\(h_x\)表示包含\(x\)的\(S\)的方案数。问\(B\)的方案，只要满足\((4)\)那么直接就是\(h_x\choose k\)了。
但是这里有些问题，对于一个\(B\)，可能有多个\(x\)使得它合法。设这个\(x\)的集合为\(|X|\)，一个方案我们总共算了\(|X|\)次。怎么去重?
然后题解给出一个性质
对于一个\(B\)，使得这个\(B\)合法的\(x\)的集合\(|X|\)在树上构成一个联通块
假设一条链\(1--2--3\)，\(1,3\)可以但是\(2\)不行，这种情况不可能存在
因为:
- 对于\(1\)的子树,既然\(3\)行，\(2\)为啥不行，距离还短些。
- 对于\(2\)的子树,既然\(1,3\)行为啥\(2\)不行，距离还短些。
- 对于\(3\)的子树,既然\(1\)行，\(2\)为啥不行，距离还短些。
因此\(X\)是树上一个联通块
然后由于一个联通块，点的个数=边的个数+1，因此去重就有思路了，设\(h_{e=(a,b)}\)表示必须\(ab\)点都是可以作为\(x\)的所有合法的\(S\)的个数，求\(h_e\)的方法和\(h_x\)一样，扣出来的树是\(T_a\cap T_b\)，按dfs倒序转移的时候判断下\(now==i\)，如果是这样就强制转移"选择"的情况。
那么答案就是
这样对于一个\(B\)，我们算了\(|X|-(|X|-1)=1\)次。
到这里复杂度\(O(n^2m)\)
二合一的第一部分完结，现在问题就是求一个组合数。。。。
显然\(h\le 2^{60}\)拿long long存下，现在要求一个\(n,m\)都很大的组合数，模\(5^{23}\)。
考虑exLucas。这里只需要求\(p^i\)下的组合数，然而\(exLucas\)是\(O(p^i)\)的，搞不得。
然而考虑这里的瓶颈是啥，其实是求\(g(n)\)的时候我们是暴力求循环节\(O(p^i)\)，然后给出一个不暴力的做法....
搞个生成函数
所以\(g(n)=f_n(0)\)。
然后考虑\(f\)的倍增...
对于\(f_{5k}(x)\)，\(x^{>23}\)次都是无意义的，因为系数都有一个\(5^{23}\)。而最终我们只需要求\(f_{n}(0)\)(也就是常数项)，所以对于任何\(f_{t}(x)\)，我们只需要保留前\(24\)项。
求那么\(f_{10k}\)可以递归到\(f_{5k}\)，问题规模缩小了一半。问题是\(f_{5k}(x)\)不一定能递归下去，其实只要递归到\(\le 5k\)最近的\(10\)的倍数即可，再暴力乘上最多\(9\)项形如\((x+i)\)的式子。处理\(f_n(x)\)也是同样的办法。
因为多项式长度是常数(24)，所以复杂度是\(T(n)=T(n/2)+\text{不大不小的常数}=O(\text{不大不小的常数}\log n)\)
一个坑点是,"最大价值和",没有x的合法限制,所以要单独DP出来....
什么叫二合一啊(战术仰头)
//@winlere
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
using namespace std; typedef long long ll;
typedef vector<ll> poly;
int n;
ll m,k,Max_val;
const ll mod=11920928955078125;
template<class T> T MOD(const T&a){return a>=mod?a-mod:a;}
template<class T> T MOD(const T&a,const T&b){return (__int128)a*b%mod;}
ll ksm(ll ba,ll p){
ll ret=1;
while(p){
if(p&1) ret=MOD(ret,ba);
ba=MOD(ba,ba);
p>>=1;
}
return ret;
}
const ll phi=ksm(5,22)*4;
const int SS=24;
poly operator * (poly a,poly b){
if(a.empty()||b.empty()) return poly(24,0);
poly ret(a.size()+b.size()-1,0);
for(int t=0,ed=a.size();t<ed;++t)
for(int i=0,ed2=b.size();i<ed2;++i)
ret[t+i]=MOD(ret[t+i]+MOD(a[t],b[i]));
return ret.resize(min<int>(SS,ret.size())),ret;
}
namespace C{
ll cnt(ll n){
ll ret=0;
while(n) ret+=n/5,n/=5;
return ret;
}
poly PP[101];
ll ch[SS+1][SS+1];
poly fac(ll n){
if(n<=100) return PP[n];
ll k=n/10*10,w=1;
poly f=fac(k/2),ret(f.size(),0);
for(int t=0,ed=f.size();t<ed;++t,w=MOD<ll>(w,k/2))
for(int i=t;i<ed;++i)
ret[i-t]=MOD(ret[i-t]+MOD(MOD(w,ch[i][t]),f[i]));
ret=ret*f;
for(ll g=k+1;g<=n;++g)
if(g%5) ret=ret*(poly){MOD(g),1};
return ret;
}
ll jc(ll n){
ll ret=1;
while(n) ret=MOD(ret,fac(n)[0]),n/=5;
return ret;
}
void init(){
PP[0]={1,0};
for(int t=1;t<=100;++t)
if(t%5) PP[t]=PP[t-1]*(poly){t,1};
else PP[t]=PP[t-1];
for(int t=0;t<=SS;++t)
for(int i=ch[t][0]=1;i<=t;++i)
ch[t][i]=MOD(ch[t-1][i]+ch[t-1][i-1]);
}
ll c(ll n){
if(n<k) return 0;
if(k==1) return n;
ll ret=ksm(5,cnt(n)-cnt(k)-cnt(n-k));
if(!ret) return ret;
ret=MOD(ret,MOD(jc(n),ksm(MOD(jc(k),jc(n-k)),phi-1)));
return ret;
}
}
namespace DDPP{
const int maxn=60+5;
vector<pair<int,int>> e[maxn];
void add(int fr,int to,int w){
e[fr].push_back({to,w});
e[to].push_back({fr,w});
}
int siz[maxn],dfn[maxn],arc[maxn],r[maxn],w[maxn],v[maxn],time;
ll dis[maxn];
struct DATA{
ll cnt,val;
DATA(ll a=0,ll b=0):cnt(a),val(b){}
DATA operator + (DATA x){return {(val>=x.val)*cnt+(x.val>=val)*x.cnt,(val>=x.val)*val+(val<x.val)*x.val};}
DATA operator + (ll x){return {cnt,cnt?val+x:0};}
}dp[maxn][10001];
void dfs(int now,int last){
siz[now]=0;
if(Max_val!=-1&&dis[now]*v[now]>Max_val) return;
dfn[now]=++time; arc[time]=now; siz[now]=1;
for(auto t:e[now])
if(t.first^last)
dis[t.first]=dis[now]+t.second,dfs(t.first,now),siz[now]+=siz[t.first];
}
DATA solve(int rt,int must){
if(must==-1) return {0,0};
memset(arc,0,sizeof arc);
memset(dp,0,sizeof dp);
time=0;
ll W=must?lower_bound(e[rt].begin(),e[rt].end(),(pair<int,int>){must,0})->second:0;
if(Max_val!=-1&&(W*v[rt]>Max_val||W*v[must]>Max_val)) return {0,0};
siz[rt]=1; dis[rt]=W; dfn[rt]=time=1; arc[1]=rt;
if(must) dis[must]=W,dfs(must,rt),siz[rt]+=siz[must];
for(auto t:e[rt])
if(t.first!=must)
dis[t.first]=W+t.second,dfs(t.first,rt),siz[rt]+=siz[t.first];
dp[0][0]={1,0};
for(int t=time;t>0;--t){
int cur=arc[t],last=arc[t+1],sub=arc[t+siz[cur]];
if(cur==must||cur==rt)
for(int i=w[cur];i<=m;++i)
dp[cur][i]=dp[last][i-w[cur]]+v[cur];
else{
for(int i=w[cur];i<=m;++i)
dp[cur][i]=(dp[last][i-w[cur]]+v[cur])+dp[sub][i];
for(int i=0;i<w[cur];++i)
dp[cur][i]=dp[sub][i];
}
}
DATA ret(0,0);
for(int t=w[rt];t<=m;++t)
ret=ret+dp[rt][t];
return ret;
}
void dfs0(int now,int last){
r[now]=last;
for(auto t:e[now])
if(t.first^last)
dfs0(t.first,now);
}
void init(){
dfs0(1,-1);
for(int t=1;t<=n;++t) sort(e[t].begin(),e[t].end());
}
}
int main(){
#ifndef ONLINE_JUDGE
freopen("yukinoshita_yukino.in","r",stdin);
freopen("yukinoshita_yukino.out","w",stdout);
#endif
cin.tie(0); cout.tie(0); ios::sync_with_stdio(0);
cin>>n>>m>>k>>Max_val;
for(int t=1;t<=n;++t) cin>>DDPP::w[t];
for(int t=1;t<=n;++t) cin>>DDPP::v[t];
for(int t=1,a,b,c;t<n;++t)
cin>>a>>b>>c,DDPP::add(a,b,c);
C::init(); DDPP::init();
ll ret=0,g=0,sav=Max_val;
Max_val=-1;
for(int t=1;t<=n;++t){
auto t1=DDPP::solve(t,0);
if(t1.val>g) g=t1.val;
}
Max_val=sav;
for(int t=1;t<=n;++t){
auto t1=DDPP::solve(t,0),t2=DDPP::solve(t,DDPP::r[t]);
if(t1.val>g) ret=0,g=t1.val;
if(t1.val==g) ret=MOD(ret+C::c(t1.cnt));
if(t2.val==g) ret=MOD(ret-C::c(t2.cnt)+mod);
}
cout<<ret<<endl;
return 0;
}
