显然,如果选择的 \(k\) 个“合法集合”固定了,那么可以放置装置的点如果存在,那么必然形成一个连通块,也就是说,答案等于所有合法方案中,可以放置装置的点形成的连通块个数之和。而根据点减边的套路,这等价于,枚举每个点,计算有多少种方案满足可以在其放置装置,再枚举每条边,计算有多少种方案满足这条边两个端点都可以放置转置,二者相减。
以点为例,这样那个 \(\text{dis}(x,y)·v_y\le \text{Max}\) 相当于一部分点不能选择,而剩下的点中选择需要形成一个以 \(x\) 为根的连通块,这时候有一个经典的处理树上连通块问题的 DP 模型就是先将 DP 值传给儿子,DFS 一圈以后再将儿子的 DP 值与父亲合并,本质上来说这相当于在欧拉序上进行背包。这样我们可以在 \(O(nm)\) 的时间内计算有多少个合法集合满足 \(x\) 点可以放置装置,\(\dbinom{cnt}{k}\) 就是选出 \(k\) 个这样的集合的方案数。
然后你瞟一眼数据范围发现 \(k\) 数据范围达到 \(10^9\),并且模数很怪,分解质因数一下发现竟然是 \(5^{23}\)。因此我们还要再思考一下怎么计算这个大组合数。显然这等价于计算 \(n!\) 中 \(5\) 的次数与所有与 \(5\) 互质的数的乘积。前者是好求的,考虑怎么计算后者,我们记 \(F_n(x)\) 表示 \(1\sim 5n\) 中所有与 \(5\) 互质的 \(i\) 的 \(x+i\) 的乘积形成的多项式。这样我们对后面零头部分暴力算,剩余部分等价于 \(F_{\lfloor\frac{n}{5}\rfloor}(0)\)。考虑倍增,显然 \(F_{2n}(x)=F_{n}(x)·F_{n}(x+5n)\),后面的那部分可以用二项式定理拆开,由于 \(5n\) 一定是 \(5\) 的倍数,因此二项式定理的时候,\(x^{23}\) 以上的项是没有用的,我们只用存前 \(23\) 项的系数即可。
时间复杂度 \(n^2m+n·(23^2·\log V)\)
typedef __int128_t i128;
const ll MOD=11920928955078125ll;
ll mul(i128 x,i128 y){return x*y%MOD;}
namespace Binomial{
ll pw5[25],c[25][25];
struct poly{
ll v[25];
poly(){memset(v,0,sizeof(v));}
friend poly operator *(const poly &X,const poly &Y){
poly res;
for(int i=0;i<=23;i++)for(int j=0;j+i<=23;j++)
res.v[i+j]=(res.v[i+j]+mul(X.v[i],Y.v[j]))%MOD;
return res;
}
}f[65];//F_{5*2^i}(x)
poly shift(poly x,ll v){
static ll pw[25];poly res;
for(int i=(pw[0]=1);i<=23;i++)pw[i]=mul(pw[i-1],v);
for(int i=0;i<=23;i++)for(int j=0;j<=i;j++)
res.v[j]=(res.v[j]+mul(mul(c[i][j],pw[i-j]),x.v[i]))%MOD;
return res;
}
ll calc_pw5(ll x){if(!x)return 0;return x/5+calc_pw5(x/5);}
void exgcd(ll x,ll y,ll &a,ll &b){
if(!y)return a=1,b=0,void();
exgcd(y,x%y,a,b);ll tmp=a;a=b;b=tmp-x/y*b;
}
ll inv(ll x){ll a,b;exgcd(x,MOD,a,b);return (a+MOD)%MOD;}
ll calc_prd(ll x){
poly prd;prd.v[0]=1;ll sum=0;
for(int i=60;~i;i--)if((x/5)>>i&1)prd=prd*shift(f[i],sum),sum+=5ll<<i;
ll res=prd.v[0];
for(ll i=x/5*5+1;i<=x;i++)if(i%5)res=mul(res,i);
return res;
}
ll calc(ll x){if(!x)return 1;return mul(calc_prd(x),calc(x/5));}
ll binom(ll n,ll k){
if(n<0||k<0||n<k)return 0;
ll pw=calc_pw5(n)-calc_pw5(k)-calc_pw5(n-k);
if(pw>=23)return 0;
else{
ll A=calc(n),B=calc(k),C=calc(n-k);
return mul(mul(A,inv(B)),mul(inv(C),pw5[pw]));
}
}
void init(){
for(int i=(pw5[0]=1);i<=23;i++)pw5[i]=pw5[i-1]*5;
for(int i=0;i<=23;i++){
c[i][0]=1;
for(int j=1;j<=i;j++)c[i][j]=c[i-1][j]+c[i-1][j-1];
}
f[0].v[0]=24;f[0].v[1]=50;f[0].v[2]=35;f[0].v[3]=10;f[0].v[4]=1;
for(int i=1;i<=60;i++)f[i]=f[i-1]*shift(f[i-1],5ll<<i-1);
}
}
ll binom(ll n,ll k){return Binomial::binom(n,k);}
const int MAXN=60;
const int MAXM=1e4;
int n,m,k,w[MAXN+5],v[MAXN+5],hd[MAXN+5],to[MAXN*2+5],nxt[MAXN*2+5],val[MAXN*2+5],ec;ll mx;
void adde(int u,int v,int w){to[++ec]=v;val[ec]=w;nxt[ec]=hd[u];hd[u]=ec;}
int dis[MAXN+5][MAXN+5];
void dfs0(int x,int f,int rt){
for(int e=hd[x];e;e=nxt[e]){
int y=to[e],z=val[e];if(y==f)continue;
dis[rt][y]=dis[rt][x]+z;dfs0(y,x,rt);
}
}
struct dat{
ll mx,cnt;
dat(){mx=-1e18;cnt=0;}
dat(ll _mx,ll _cnt){mx=_mx;cnt=_cnt;}
friend dat operator +(const dat &X,const dat &Y){
dat ret;ret.mx=max(X.mx,Y.mx);ret.cnt=0;
if(X.mx==ret.mx)ret.cnt+=X.cnt;
if(Y.mx==ret.mx)ret.cnt+=Y.cnt;
return ret;
}
}dp[MAXN+5][MAXM+5];
bool ban[MAXN+5];
void dfs(int x,int f){
for(int e=hd[x];e;e=nxt[e]){
int y=to[e];if(y==f||ban[y])continue;
for(int i=0;i+w[y]<=m;i++)dp[y][i+w[y]]=dat(dp[x][i].mx+v[y],dp[x][i].cnt);
dfs(y,x);for(int i=0;i<=m;i++)dp[x][i]=dp[x][i]+dp[y][i];
}
}
ll mx_val;
ll calc_val(int x){
for(int i=1;i<=n;i++)for(int j=0;j<=m;j++)dp[i][j]=dat();
dp[x][w[x]]=dat(v[x],1);dfs(x,0);dat res;
for(int i=0;i<=m;i++)res=res+dp[x][i];
return res.mx;
}
ll calc_vert(int x){
memset(ban,0,sizeof(ban));
for(int i=1;i<=n;i++)if(1ll*dis[i][x]*v[i]>mx)ban[i]=1;
for(int i=1;i<=n;i++)for(int j=0;j<=m;j++)dp[i][j]=dat();
dp[x][w[x]]=dat(v[x],1);dfs(x,0);dat res;
for(int i=0;i<=m;i++)res=res+dp[x][i];
return (res.mx==mx_val)?binom(res.cnt,k):0;
}
ll calc_edge(int x,int y){
memset(ban,0,sizeof(ban));
for(int i=1;i<=n;i++)if(1ll*dis[i][x]*v[i]>mx||1ll*dis[i][y]*v[i]>mx)ban[i]=1;
if(ban[x]||ban[y])return 0;
for(int i=1;i<=n;i++)for(int j=0;j<=m;j++)dp[i][j]=dat();
dp[x][w[x]]=dat(v[x],1);dp[y][w[y]]=dat(v[y],1);dfs(x,y);dfs(y,x);
for(int i=1;i<=m;i++)dp[y][i]=dp[y][i]+dp[y][i-1];
dat res;for(int i=0;i<=m;i++)res=res+dat(dp[x][i].mx+dp[y][m-i].mx,1ll*dp[x][i].cnt*dp[y][m-i].cnt);
return (res.mx==mx_val)?binom(res.cnt,k):0;
}
int main(){
Binomial::init();scanf("%d%d%d%lld",&n,&m,&k,&mx);
for(int i=1;i<=n;i++)scanf("%d",&w[i]);
for(int i=1;i<=n;i++)scanf("%d",&v[i]);
for(int i=1,u,v,w;i<n;i++)scanf("%d%d%d",&u,&v,&w),adde(u,v,w),adde(v,u,w);
for(int i=1;i<=n;i++)dfs0(i,0,i);
for(int i=1;i<=n;i++)chkmax(mx_val,calc_val(i));
ll res=0;
for(int i=1;i<=n;i++)if(w[i]<=m)res=(res+calc_vert(i))%MOD;
for(int i=1;i<=n;i++)for(int e=hd[i];e;e=nxt[e]){
int j=to[e];
if(i<j&&w[i]+w[j]<=m)res=(res-calc_edge(i,j)+MOD)%MOD;
}
printf("%lld\n",res);
return 0;
}
标签:洛谷,23,ll,poly,P9248,2018,mul,ban,return
From: https://www.cnblogs.com/tzcwk/p/luogu-P9248.html