【題解】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揹包
狀態不變，轉移順序改變一下，按照\(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\)，如果是這樣就強制轉移"選擇"的情況。
那麼答案就是
\[\sum_i {h_i \choose k} -\sum_{e=(a,b)} {h_e \choose k} \]
這樣對於一個\(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)\)，然後給出一個不暴力的做法....
搞個生成函式
\[f_n(x)=\prod_{5|i}^n (x+i) \]
所以\(g(n)=f_n(0)\)。
然後考慮\(f\)的倍增...
\[f_{10k}(x)=f_{5k}(x)f_{5k}(x+5k) \]
對於\(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)\)
什麼叫二合一啊(戰術仰頭)
//@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);
}
return ret;
}
const ll phi=ksm(5,22)*4;
namespace C{
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,ed=b.size();i<ed;++i)
ret[t+i]=MOD(ret[t+i]+MOD(a[t],b[i]));
return ret.resize(min<int>(SS,ret.size())),ret;
}
ll cnt(ll n){
ll ret=0;
while(n) ret+=n/5,n/=5;
return ret;
}
poly PP[101];
poly fac(ll n){
if(n<=10) return PP[n];
ll k=n/10*10,sav=0;
auto f=fac(k/2),ret=f;
for(int i=f.size()-1;~i;--i)
sav=MOD(MOD(sav+f[i]),k/2);
ret[0]=MOD(ret[0]+sav);
ret=ret*f;
for(ll g=k;g<=n;++g)
if(g%5) ret=ret*(poly){MOD(g),1};
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];
}
ll c(ll n){
if(n<k) return 0;
if(n==k) return 1;
ll ret=ksm(5,cnt(n)-cnt(k)-cnt(n-k));
ret=MOD(ret,MOD(fac(n)[0],ksm(MOD(fac(k)[0],fac(n-k)[0]),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(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){
time=0;
int g=must?lower_bound(e[rt].begin(),e[rt].end(),(pair<int,int>){must,0})-e[rt].begin():0;
swap(e[rt][0],e[rt][g]);
int W=must?e[rt].front().second:0;
dfn[rt]=time=1; arc[1]=rt;
for(auto t:e[rt])
dis[t.first]=W,dfs(t.first,rt);
memset(dp,0,sizeof dp);
dp[arc[time]][w[arc[time]]]={1,v[arc[time]]};
dp[arc[time]][0]={1,0};
for(int t=time-1;t>0;--t){
int cur=arc[t],last=arc[t+1],sub=arc[t+siz[cur]];
if(last==must)
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]]+dp[sub][i-w[cur]]+v[cur];
}
DATA ret(0,0);
for(int t=w[rt];t<=m;++t)
ret=ret+dp[rt][t];
swap(e[rt][0],e[rt][g]);
return ret;
}
void init(){
for(int t=1;t<=n;++t) sort(e[t].begin(),e[t].end());
}
}
int main(){
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();
ll ret=0;
for(int t=1;t<=n;++t){
auto t1=DDPP::solve(t,0),t2=DDPP::solve(t,DDPP::r[t]);
if(t1.val>Max_val) ret=0,Max_val=t1.val;
if(t1.val==Max_val) ret=MOD(ret+C::c(t1.cnt));
if(t2.val==Max_val) ret=MOD(ret-C::c(t2.cnt)+mod);
}
cout<<ret<<endl;
return 0;
}