「LOJ2462」完美的集合

题目

小 A 有一棵 \(N\) 个点的带边权的树，树的每个节点有重量 \(w_i\) 和价值 \(v_i\)。

现在小 A 要从中选出若干个节点形成一个集合 \(S\)，满足这些节点重量之和 \(\leq M\) 并且构成一个连通块。小 A 是一个完美主义者，因此他只会选择节点价值之和最大的那些 \(S\)。我们称这样的集合 \(S\) 为完美的集合。

现在小 \(A\) 要从所有完美的集合中选出 \(K\) 个，并对这 \(K\) 个完美的集合分别进行测试。在这 \(K\) 次测试开始前，小 A 首先需要一个点 \(x\) 来放置他的测试装置，这个测试装置的最大功率为 \(\mathrm{Max}\)。

接下来的每次测试，小 A 会对测试对象 \(S\) 中的所有点进行一次能量传输，对一个点 \(y\) 进行能量传输需要的功率为 \(\operatorname{dist}(x,y)\times v_y\)，其中 \(\operatorname{dist}(x,y)\) 表示点 \(x,y\) 在树上的最短路长度。因此，如果 \(S\) 中存在一个点 \(y\)，满足 \(\operatorname{dist}(x,y)\times v_y>\mathrm{Max}\)，测试就会失败。同时，为了保证能量传输的稳定性，测试装置所在的点 \(x\) 需要在集合 \(S\) 中，否则测试也会失败。

现在小 A 想知道，有多少种从所有完美的集合选出 \(K\) 个的方法，使得他能找到一个放置测试装置的点，来完成他的测试呢？

你只需要输出方案数对 \(11920928955078125\) 取模的结果。

所有数据满足 \(1\le N\le 60,1\le M\le 10^4, 1\le C_i\le 10^4, 1\le K,w_i\le 10^9,0\le v_i\le 10^9,1\le \mathrm{Max}\le 10^{18}\)。

分析

“选 \(K\) 个”的要求是依托答辩，先不管它，考虑 \(K=1\) 的情况。

注意到，对于定点 \(y\)，满足 \(\operatorname{dist}(x,y)\times v_y>\mathrm{Max}\) 的 \(x\) 构成树上的一个连通块。于是，对于任意一个集合 \(S\)，可以放置测试装置的点也构成树上的一个连通块。因此，做一个“点-边容斥”即可将问题转化为“使得点 \(x\)（或边 \((x,y)\) 的两个端点）可以放置测试装置的完美的集合有多少个？”。

统计根为某顶点的连通块的信息是容易的，此处就有一个 \(O(nm)\) 的算法。于是，先用一个 \(O(n^2m)\)（或 \(O(nm\log n)\)）的算法计算出完美的集合的结点价值之和，再花 \(O(n^2m)\) 的时间做上述容斥，复杂度即为 \(O(n^2m)\)。

加入“选 \(K\) 个”的要求后，我们发现对于一个集簇 \(\mathcal F\)，可以放置测试装置的点还是形成一个树上的连通块，于是“点-边容斥”还是可行的。只不过，如果“使得点 \(x\)（或边 \((x,y)\) 的两个端点）可以放置测试装置的完美的集合”的个数为 \(t\)，原先贡献为 \(\pm t\)，现在贡献为 \(\pm\binom{t}{K}\)，这一点和「十二省联考 2019」希望一模一样。

（正片开始）

问题的主要矛盾来到了“计算 \(\binom{t}{K}\)”，其中 \(t,K\) 都是（相对于通常情况而言）非常大的数。不过，因为 \(t\le 2^n\)，所以上下指标都还能存得下真实值。

此时就不得不研究模数的性质了。进行一个质因数分解，发现 \(M=11920928955078125=5^{23}\)，非常的 smooth。所以，我们考虑使用扩展 Lucas 算法。

扩展 Lucas 算法关键为：对于 \(n\)，求出 \(\prod_{1\le k\le n,5\nmid k}k\)。在模数较小时，我们可以直接利用乘积式的周期性预处理，但现在不行。考虑常规策略——分治计算。特别地，因为结果和模 \(5\) 的余数有关，我们考虑从它开始进行分治。

假如现在要对于 \(n>1\)，求出 \(g(n)=\prod_{1\le k\le n,5\nmid k}k\)。设 \(m=\lceil\log_5 n\rceil-1\) 为满足 \(5^m<n\) 的最大的整数，并设 \(r=\lceil\frac n {5^m}\rceil-1\)。当 \(m\ge 1\) 时，我们可以将乘积式划分为：

\[\prod_{t=0}^{r-1}\left(\prod_{t\times 5^m< k\le (t+1)\times 5^m,5\nmid k}k\right)\cdot \prod_{r\times 5^m< k\le n,5\nmid k}k \]

注意到 \(r\times 5^m+k\) 对于 \(5\) 的整除性和 \(k\) 相同，所以有：

\[\prod_{r\times 5^m< k\le n,5\nmid k}k=\prod_{0< k\le n-r\times 5^m,5\nmid k}(k+r\times 5^m) \]

要想快速地计算这个式子，我们可以考虑计算出 \(f_{n}(x)=\prod_{0<k\le n,5\nmid k}(x+k)\)，这样就将 \(n\) 减小，变成了子问题。对于前面若干个长度为 \(5^m\) 的部分做类似处理，我们要算的就是：

\[f_n(x)=\prod_{t=0}^{r-1}f_{5^m}(x+t\times 5^m)\cdot f_{n-r\times 5^m}(x+r\times 5^m) \]

问题来了：这个多项式可能很长，怎么办？观察到三条性质：

\(g(n)=[x^0]f_n(x)\)，也就是我们只需要保证它的常数项是对的就好。
答案对于 \(5^{23}\) 取模，而复合的一次式的常数项一定是 \(5\) 的倍数。
合并子问题时，基本运算为“复合一次式”和“卷积”，其中只有“复合一次式”会产生高次项到低次项的影响。而计算 \(f(x+c)\) 的过程中，\(x^{t}\) 到 \(x^0\) 的贡献必然带有 \(c^t\) 的因子。即便是多次复合，包含参数的因式也是一个次数为 \(t\) 的齐次式。

以上三条可以导出，次数不低于 \(23\) 的项都可以被舍弃，因为它们不会再对常数项造成影响。这样多项式就变得很轻巧了，可以快速地完成运算。

Remark.

相对间接的处理方法。我们只有一个 \(\pmod {5^{23}}\) 的条件，要想从它出发搞出一个 \(\pmod {x^c}\)，就只能借助系数贡献的性质，也就是考虑系数中 \(5\) 的指数。

特别地，当 \(m=0\) 时，\(5^m\) 不再是 \(5\) 的倍数。不过，这个边界情况下，直接暴力算一下多项式的乘积即可。

具体操作的时候，需要对于每一个 \(m\) 都预处理出 \(f_{5^m}(x)\)，最好还可以预处理一个 \(\prod_{t=0}^{r-1}f_{5^m}(x+t\times 5^m)\) 的前缀积。

复杂度不算了，能过就是了。

代码

#include <bits/stdc++.h>

#define rep( i, a, b ) for( int i = (a) ; i <= (b) ; i ++ )
#define per( i, a, b ) for( int i = (a) ; i >= (b) ; i -- )

typedef long long LL;
typedef  __int128 ExLL;

const LL mod = 11920928955078125;
const int MAXLOG = 30, MAXN = 65, MAXM = 10005;

template<typename _T>
inline void Read( _T &x ) {
    x = 0; char s = getchar(); bool f = false;
    while( s < '0' || '9' < s ) { f = s == '-', s = getchar(); }
    while( '0' <= s && s <= '9' ) { x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar(); }
    if( f ) x = -x;
}

template<typename _T>
inline void Write( _T x ) {
    if( x < 0 ) putchar( '-' ), x = -x;
    if( 9 < x ) Write( x / 10 );
    putchar( x % 10 + '0' );
}

inline LL Sub( LL x, const LL &v ) { return ( x -= v ) < 0 ? x + mod : x; }
inline LL Add( LL x, const LL &v ) { return ( x += v ) >= mod ? x - mod : x; }

inline LL& SubEq( LL &x, const LL &v ) { return ( x -= v ) < 0 ? ( x += mod ) : x; }
inline LL& AddEq( LL &x, const LL &v ) { return ( x += v ) >= mod ? ( x -= mod ) : x; }

namespace PureCalculation {
    struct Poly {
        LL coe[23];

        Poly(): coe{} {}
    };

    Poly bas[MAXLOG], pref[MAXLOG][6];
    LL sml[23][23], pw[23], pw5[MAXLOG];
    LL vpK, factK; int K;

    inline LL Inv( LL base, LL indx = 9536743164062500 - 1 ) {
        LL ret = 1;
        while( indx ) {
            if( indx & 1 ) ret = ( ExLL ) ret * base % mod;
            base = ( ExLL ) base * base % mod, indx >>= 1;
        }
        return ret;
    }

    inline Poly operator * ( const Poly &a, const Poly &b ) {
        Poly ret; ExLL tmp;
        rep( i, 0, 22 ) {
            tmp = 0;
            rep( j, 0, i )
                tmp += ( ExLL ) a.coe[j] * b.coe[i - j];
            ret.coe[i] = tmp % mod;
        }
        return ret;
    }

    inline Poly Shift( const Poly &f, const LL &c ) {
        if( ! c ) return f;
        Poly ret; ExLL tmp;
        pw[0] = 1; 
        rep( i, 1, 22 ) pw[i] = ( ExLL ) pw[i - 1] * c % mod;
        rep( i, 0, 22 ) {
            tmp = 0;
            rep( j, i, 22 )
                tmp += ( ExLL ) f.coe[j] * pw[j - i] * sml[j][i];
            ret.coe[i] = tmp % mod;
        }
        return ret;
    }

    inline Poly PartialFactorial( const int &lvl, const LL &n ) {
        if( lvl == 0 ) return pref[0][n];
        int idx = ( n - 1 ) / pw5[lvl];
        if( ! idx ) return PartialFactorial( lvl - 1, n );
        return Shift( PartialFactorial( lvl - 1, n - idx * pw5[lvl] ), idx * pw5[lvl] ) * pref[lvl][idx];
    }

    inline LL Factorial( const LL &n ) {
        if( n == 0 ) return 1;
        return ( ExLL ) PartialFactorial( 25, n ).coe[0] * Factorial( n / 5 ) % mod;
    }

    inline LL Legendre( LL n ) {
        LL ret = 0;
        for( ; n ; ret += n /= 5 );
        return ret;
    }

    inline LL Binom( const LL &n ) {
        if( n < K ) return 0;
        LL idx = Legendre( n ) - vpK - Legendre( n - K );
        if( idx >= 23 ) return 0;
        LL a = Factorial( n ), c = Factorial( n - K );
        return ( ExLL ) a * Inv( c ) % mod * factK % mod * pw5[idx] % mod;
    }

    inline void Init( const int &k ) {
        K = k, pw5[0] = 1;
        rep( i, 1, 25 ) pw5[i] = pw5[i - 1] * 5;
        rep( i, 0, 22 ) {
            sml[i][0] = 1;
            rep( j, 1, i )
                sml[i][j] = Add( sml[i - 1][j], sml[i - 1][j - 1] );
        }
        bas[0].coe[0] = bas[0].coe[1] = 1;
        pref[0][0].coe[0] = 1;
        rep( i, 1, 4 ) {
            Poly tmp;
            tmp.coe[0] = i, tmp.coe[1] = 1;
            pref[0][i] = pref[0][i - 1] * tmp;
        }
        pref[0][5] = pref[0][4];
        rep( i, 1, 25 ) {
            bas[i] = pref[i - 1][5];
            pref[i][0].coe[0] = 1;
            rep( j, 1, 5 )
                pref[i][j] = pref[i][j - 1] * Shift( bas[i], pw5[i] * ( j - 1 ) );
        }
        vpK = Legendre( K );
        factK = Inv( Factorial( K ) );
    }
}

namespace OnTree {
    struct Edge {
        int to, nxt, w;
    } Graph[MAXN << 1];

    struct Values {
        LL val, cnt;
    
        Values(): val( -1 ), cnt( 0 ) {}
        Values( LL V ): val( V ), cnt( 1 ) {}
        Values( LL V, LL C ): val( V ), cnt( C ) {}

        inline void operator += ( const Values &q ) {
            if( q.val  > val ) val = q.val, cnt = 0;
            if( q.val == val ) cnt += q.cnt;
        }

        inline Values operator + ( const Values &q ) const {
            if( val > q.val ) return *this;
            if( val < q.val ) return q;
            return Values( val, cnt + q.cnt );
        }

        inline Values operator * ( const Values &q ) const {
            return Values( val + q.val, cnt * q.cnt );
        }
    };

    Values dp[MAXN][MAXM];
    LL dist[MAXN][MAXN];

    int seq[MAXN], siz[MAXN], tot = 0;
    int edgFr[MAXN], edgTo[MAXN];
    int head[MAXN], cnt = 1;

    int wei[MAXN], val[MAXN];

    int N, M, K; LL lim;

    inline void AddEdge( const int &from, const int &to, const int &W ) {
        Graph[++ cnt].to = to, Graph[cnt].nxt = head[from];
        Graph[cnt].w = W, head[from] = cnt;
    }

    inline void Input() {
        Read( N ), Read( M ), Read( K ), Read( lim );
        rep( i, 1, N ) Read( wei[i] );
        rep( i, 1, N ) Read( val[i] );
        rep( i, 1, N - 1 ) {
            int u, v, w;
            Read( u ), Read( v ), Read( w );
            AddEdge( u, v, w ), AddEdge( v, u, w );
            edgFr[i] = u, edgTo[i] = v;
        }
    }

    void ProcessDist( const int &u, const int &fa, LL *d ) {
        for( int i = head[u], v ; i ; i = Graph[i].nxt )
            if( ( v = Graph[i].to ) ^ fa )
                d[v] = d[u] + Graph[i].w, ProcessDist( v, u, d );
    }
    
    void DFS( const int &u, const int &fa ) {
        seq[++ tot] = u, siz[u] = 1;
        for( int i = head[u], v ; i ; i = Graph[i].nxt )
            if( ( v = Graph[i].to ) ^ fa )
                DFS( v, u ), siz[u] += siz[v];
    }

    inline void ClearDP() {
        rep( i, 1, N + 1 )
            rep( j, 0, M )
                dp[i][j] = Values();
    }

    inline void Solve() {
        Input();
        rep( i, 1, N )
            ProcessDist( i, 0, dist[i] );
        PureCalculation :: Init( K );

        Values glb;
        rep( i, 1, N ) {
            tot = 0, DFS( i, 0 );
            ClearDP(), dp[1][0] = Values( 0 );
            rep( j, 1, N ) {
                int u = seq[j];
                rep( k, 0, M ) 
                    if( dp[j][k].val >= 0 )
                        dp[j + siz[u]][k] += dp[j][k];
                if( wei[u] <= M ) {
                    Values delt( val[u] );
                    rep( k, 0, M - wei[u] )
                        if( dp[j][k].val >= 0 )
                            dp[j + 1][k + wei[u]] += dp[j][k] * delt;
                }
            }
            rep( k, 0, M ) 
                glb += dp[N + 1][k];
        }

        LL ans = 0;
        rep( i, 1, N ) {
            if( wei[i] > M ) continue;
            tot = 0, DFS( i, 0 );
            ClearDP(), dp[1][0] = Values( 0 );
            rep( j, 1, N ) {
                int u = seq[j];
                if( u != i ) {
                    rep( k, 0, M )
                        if( dp[j][k].val >= 0 )
                            dp[j + siz[u]][k] += dp[j][k];
                }
                if( wei[u] <= M && ( ExLL ) dist[i][u] * val[u] <= lim ) {
                    Values delt( val[u] );
                    rep( k, 0, M - wei[u] )
                        if( dp[j][k].val >= 0 )
                            dp[j + 1][k + wei[u]] += dp[j][k] * delt;
                }
            }
            LL num = 0;
            rep( k, 0, M )
                if( glb.val == dp[N + 1][k].val )
                    num += dp[N + 1][k].cnt;
            AddEq( ans, PureCalculation :: Binom( num ) );
        }
        rep( i, 1, N - 1 ) {
            int x = edgFr[i], y = edgTo[i];
            if( wei[x] + wei[y] > M ) continue;
            tot = 0, DFS( x, 0 );
            ClearDP(), dp[1][0] = Values( 0 );
            rep( j, 1, N ) {
                int u = seq[j];
                if( u != x && u != y ) {
                    rep( k, 0, M )
                        if( dp[j][k].val >= 0 )
                            dp[j + siz[u]][k] += dp[j][k];
                }
                if( wei[u] <= M && ( ExLL ) dist[x][u] * val[u] <= lim &&
                                   ( ExLL ) dist[y][u] * val[u] <= lim ) {
                    Values delt( val[u] );
                    rep( k, 0, M - wei[u] )
                        if( dp[j][k].val >= 0 )
                            dp[j + 1][k + wei[u]] += dp[j][k] * delt;
                }
            }
            LL num = 0;
            rep( k, 0, M )
                if( glb.val == dp[N + 1][k].val )
                    num += dp[N + 1][k].cnt;
            SubEq( ans, PureCalculation :: Binom( num ) );
        }
        Write( ans ), putchar( '\n' );
    }
}

int main() {
    OnTree :: Solve();
    return 0;
}

posted @ 2023-05-19 22:33 crashed 阅读(60) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

crashed

「LOJ2462」完美的集合

题目

分析

代码

公告