LOJ #2462. 「2018 集训队互测 Day 1」完美的集合

题目链接

LOJ #2462. 「2018 集训队互测 Day 1」完美的集合

题目大意

有一棵 \(N\) 个点的带权树，树上每个节点有重量 \(w_i\) 和价值 \(v_i\) ，在满足节点重量之和 \(\leq M\) 的集合中，称那些价值之和最大的集合为完美集合。

现在要从所有完美集合中选出 \(K\) 个，要求这 \(K\) 个集合的并集中存在一个点 \(x\)，满足对于这些集合中任意一个点 \(y\)，都有 \(dist(x,y)\times v_y\leq Max\)。求满足条件的方案数，答案对 \(11920928955078125\) 取模。

\(N\leq 60\)，\(M\leq 10000\)，\(K,w_i,v_i\leq 10^9\)，\(Max\leq 10^{18}\)，树边权 \(\leq 10^4\)

Part 1

首先考虑求出完美集合的价值，直接做树上背包的话是 \(O(nm^2)\) 的，瓶颈在于需要合并两个区域的重量之和，要想得到较优的时间复杂度，可以设计一个状态，使得每次只需要决策单点的重量。于是我们先枚举一个点作为根，强制连通块包含该点，然后在树的 \(dfs\) 序上做 \(dp\)，设 \(f_{i,j}\) 表示当前考虑到了 \(dfs\) 序上的第 \(i\) 个位置，已选取点重量总和为 \(j\) 时的最大价值和，记 \(x\) 为第 \(i\) 位对应的节点，则有：

• \(f_{i+siz_x,j}\leftarrow f_{i,j}\)
• \(f_{i+1,j+w_x}\leftarrow f_{i,j}+v_j\)

若当前点不选，则子树可不能选，从而直接跳到 \(i+siz_x\) 的位置继续 \(dp\)，否则往下决策。

时间复杂度 \(O(N^2M)\)。

Part 2

对于 \(K\) 个已经选好的完美集合，可以发现满足条件的点 \(x\) 组成了一个连通块（或空集），考虑容斥。对于每一个点，求出当前点可以作为 \(x\) 的选取 \(K\) 个集合方案数，对于每条边，求出两端点同时可以作为 \(x\) 的方案数，可以发现，每个选取方案都会被合法 \(x\) 的连通块中每个点和边都算一次，而在树上，点的数量 \(=\) 边的数量 \(+1\)，从而我们拿点的方案数减去边的方案数，这样每个方案就恰好被计算一次了。

设 \(n\) 为当前限制下可以选取的完美集合数量，显然方案数即为 \(\binom{n}{K}\)，于是考虑计算 \(n\) 。这一点和之前的做法是极其相似的，以钦定的点为根（边就任取一端点），\(f_{i,j}\) 改为存二元组 \((num,val)\)，表示最大价值和为 \(val\)，达到 \(val\) 的集合数量为 \(num\)，然后转移和前面基本相同，不过若 \(dp\) 当前点不满足 \(dist(x,y)\times v_y\leq Max\)，就强制不选，考虑边时另外一端点强制要选即可。

Part 3

现在得到了 \(n\)，我们要求出 \(\binom{n}{K}\) 对 \(11920928955078125\) 取模的结果，可以发现这个质数 \(=5^{23}\)，而 \(n\) 是 \(2^N=2^{60}\) 级别的，\(K\) 是 \(10^9\) 级别的，显然无法通过一般的方式计算。

将 \(\binom{n}{K}\) 换成 \(n!/((n-K)!K!)\)，注意到模数一个小质数的大幂次，于是分开计算 \(n!\)，并将 \(n!\) 拆解为 \(5\) 的指数和非 \(5\) 的倍数部分，非倍数部分存在逆元，为 \(x^{\phi(5^{23})}\)，指数部分直接相减（结果必然为正）然后乘到答案中即可。

\[n!=(\prod_{i=1}^n[i\bot5]i)\cdot(\lfloor\frac{n}{5}\rfloor!\cdot5^{\lfloor\frac{n}{5}\rfloor}) \]

\(\lfloor\frac{n}{5}\rfloor!\) 可以递归求解，考虑求解 \(\prod_{i=1}^n[i\bot5]i\) 。

设多项式 \(f_n(x)=\prod_{i=1}^n[i\bot 5](x+i)\)，我们要求的是 \(f_n(0)\) 。注意到 \(f_{10k}(x)=f_{5k}(x)\cdot f_{5k}(x+5k)\)，如果这不明显，让 \(k\) 是 \(5\) 的倍数，就有 \(f_{2k}(x)=f_{k}(x)\cdot f_k(x+k)\)，这是一个倍增的关系，容易联想到像快速幂那样，使用二进制拆解求出任意一个 \(5\) 的倍数 \(n\) 的 \(f_n(x)\)。

不仅如此，由于我们要求的是 \(f_n(0)\)，当 \(x=0\) 时，\(f_k(x+k)\) 的非常数项全是 \(5\) 的倍数，\(x^{23}\) 往后在模意义下便全部为 \(0\)，所以计算时只保留 \(f_k(x+k)\) 的前 \(23\) 项便可。 而且 \(f_k(x+k)=\sum a_i(x+k)^i\)，我们要求的 \(f_n(0)\) 就是常数项，\(i\geq 23\) 时 \((x+k)^i\) 的常数项变为 \(0\)，所以原来的函数 \(f(x)\) 亦可以只保留前 \(23\) 位，这样我们在做多项式乘法时，暴力相乘就是很快的。

具体来说，我们先预处理出 \(f_{5\times 2^i}(x)\) 这些多项式，计算 \(f_n(x)\) 时若 \(n\) 不是 \(5\) 的倍数，就直接算 \(f_{5\lfloor\frac{n}{5}\rfloor}(x)\)，然后剩余部分直接乘到算出的 \(f(0)\) 里，然后我们对 \(\frac{n}{5}\) 二进制分解，类似快速幂的步骤，先有一个初始多项式 \(f(x)=1\)，遇到 \(\frac{n}{5}\) 中为 \(1\) 的位时将 \(f(x)\) 更新为 \(f_{5\times 2^i}(x)\cdot f(x+5\times 2^i)\)，最终 \(f(x)\) 的常数项即为所求答案。

记 \(P=5^{23}\)，预处理是 \(O(\log n\log^2 P)=O(N\log^2P)\) 的，计算一次 \(f_n(0)\) 也是 \(O(N\log^2P)\) 的，求解 \(\lfloor\frac{n}{5}\rfloor!\) 要递归 \(O(\log n)=O(N)\) 层，所以这一部分复杂度是 \(O(N^2\log^2P)\) 的。

 

综上，总时间复杂度为 \(O(N^2M+N^3\log^2P)\)，\(1.5e8\) 左右，实际上要快一些。

总结

一道题用到了 \(3\) 个比较重要的技巧：

• 对于树上连通块的 \(dp\) 可以使用对 \(dfs\) 序 \(dp\) 的方式进行优化。
• 容斥问题中，若一个方案被一个树形连通块所有位置都计算到，可以使用点边容斥。
• 计算较大的二项式之类的系数时，可以考虑从多项式和倍增的角度去分析。

Code

LOJ 上的 Hack 数据是完美集合价值为 \(0\) 的情形，此时求完美集合数量时要去掉空集的情况。

#include<iostream>
#include<cstring>
#include<vector>
#define mem(a,b) memset(a, b, sizeof(a))
#define rep(i,a,b) for(int i = (a); i <= (b); i++)
#define per(i,b,a) for(int i = (b); i >= (a); i--)
#define N 62
#define V 10100
#define P 23
#define ll long long
#define lll __int128
#define mod 11920928955078125 // 5^23
#define phi 9536743164062500  // mod - mod/5
#define PLL pair<lll, lll>
#define fr first
#define sc second
using namespace std;

int n, M, K; ll Max;
int w[N], v[N];

int head[N], to[2*N], nxt[2*N];
int siz[N], p[N], dis[N][N];
int cnt, num;

ll Max_V;
PLL f[N][V];
lll C[P+2][P+2];
vector<lll> Pow[N+5];

void init(){ mem(head, -1), cnt = -1; }
void add_e(int a, int b, bool id){
    nxt[++cnt] = head[a], head[a] = cnt, to[cnt] = b;
    if(id) add_e(b, a, 0);
}

void dfs(int x, int fa){
    siz[x] = 1, p[++num] = x;
    for(int i = head[x]; ~i; i = nxt[i]){
        if(to[i] == fa) continue;
        dfs(to[i], x), siz[x] += siz[to[i]];
    }
}

void Trans(PLL &a, PLL b){
    if(a.sc < b.sc) a.fr = b.fr;
    else if(a.sc == b.sc) a.fr += b.fr;
    a.sc = max(a.sc, b.sc);
}
void dp(int root, int ex){
    mem(f, 0);
    f[1][0] = {1, 0};
    rep(i,1,n) rep(j,0,M) if(f[i][j].fr){
        int x = p[i];
        if(x != ex) Trans(f[i+siz[x]][j], f[i][j]);
        if(j+w[x] <= M && (ll)dis[root][x] * v[x] <= Max && (ll)dis[ex][x] * v[x] <= Max)
            Trans(f[i+1][j+w[x]], {f[i][j].fr, f[i][j].sc + v[x]});
    }
}

vector<lll> mul(vector<lll> &A, vector<lll> B){
    if(A.empty() || B.empty()) return {};
    int n = A.size(), m = B.size();
    vector<lll> ret(min(n + m - 1, 23));
    rep(i,0,n-1) rep(j,0,m-1) if(i+j < 23)
        (ret[i+j] += A[i] * B[j]) %= mod;
    return ret;
}

vector<lll> shift(vector<lll> &A, lll k){
    int n = A.size();
    vector<lll> ret(n);
    rep(i,0,n-1){
        lll pow = A[i];
        rep(j,0,i) (ret[i-j] += pow * C[i][j]) %= mod, (pow *= k) %= mod;
    }
    return ret;
}
 
lll qpow(lll a, ll b){
    lll ret = 1;
    for(; b; b >>= 1){ if(b&1) (ret *= a) %= mod; (a *= a) %= mod; }
    return ret;
}

void prework(){
    C[0][0] = 1;
    rep(i,1,P){
        C[i][0] = 1;
        rep(j,1,i) C[i][j] = (C[i-1][j] + C[i-1][j-1]) % mod;
    }
    Pow[0] = {1};
    rep(i,1,4) Pow[0] = mul(Pow[0], {i, 1});
    rep(i,1,N) Pow[i] = mul(Pow[i-1], shift(Pow[i-1], 5ll<<(i-1)));
}

PLL fact(lll n){
    if(n == 0) return {1, 0};
    if(n%5){
        PLL p = fact(n-1);
        return {p.fr * n % mod, p.sc};
    }

    lll exp = n/5, prod = 1;
    PLL p = fact(n/5);
    (prod *= p.fr) %= mod, (exp += p.sc) %= mod;

    vector<lll> poly = {1};
    per(i,N,0) if((n/5)>>i&1) poly = mul(Pow[i], shift(poly, 5ll<<i));
    (prod *= poly[0]) %= mod;
    return {prod, exp};
}

lll cal(lll n, lll m){
    if(m < 0 || n < m) return 0;
    PLL p = fact(n);
    lll val = p.fr, exp = p.sc;
    p = fact(n-m), (val *= qpow(p.fr, phi-1)) %= mod, exp -= p.sc;
    p = fact(m), (val *= qpow(p.fr, phi-1)) %= mod, exp -= p.sc;
    return (val * qpow(5, exp)) % mod;
}

int main(){
    cin>>n>>M>>K>>Max;
    rep(i,1,n) cin>>w[i];
    rep(i,1,n) cin>>v[i];

    init(), mem(dis, 0x3f);
    int a, b, c;
    rep(i,1,n-1){
        cin>>a>>b>>c, add_e(a, b, 1);
        dis[a][b] = c, dis[b][a] = c;
    }
    rep(i,1,n) dis[i][i] = 0;
    rep(k,1,n) rep(i,1,n) rep(j,1,n) 
        dis[i][j] = min(dis[i][j], dis[i][k] + dis[k][j]);
    rep(i,1,n) dis[i][0] = dis[0][i] = 0;

    rep(i,1,n){
        num = 0, dfs(i, 0);
        dp(0, 0);
        rep(j,0,M) Max_V = max(Max_V, (ll)f[n+1][j].sc);
    }

    prework();
    lll ans = 0;
    rep(i,1,n+(cnt+1)/2){
        num = 0;
        if(i <= n) dfs(i, 0), dp(i, 0);
        else dfs(to[(i-n-1)*2], 0), dp(to[(i-n-1)*2], to[(i-n-1)*2+1]);
        lll num = 0;
        rep(j,0,M) if(f[n+1][j].sc == Max_V) num += f[n+1][j].fr - (Max_V == 0 && j == 0);
        (ans += i <= n ? cal(num, K) : mod - cal(num, K)) %= mod;
    }
    cout<< (ll)ans <<endl;
    return 0;
}