「ZJOI2015」地震后的幻想乡

Posted 2022-03-14 crashed

题目

点这里看题目。

分析

首先，设原图的最小生成树的边集为 \\(T\\)，则容易得到：

\\[\\beginaligned E(\\max_x\\in Te_x) &=\\int_0^1P(t<\\max_x\\in Te_x)\\mathrm dt \\endaligned \\]

而可以发现 \\(P(t<\\max_x\\in Te_x)\\) 恰好为加入了 \\(\\le t\\) 的边之后，图仍然不连通的概率。因此我们可以直接枚举一个边集 \\(E\'\\subseteq E\\)，使得加入了 \\(E\'\\) 后图不连通，期望则可以被表述为：

\\[\\beginaligned \\int_0^1P(t<\\max_x\\in Te_x)\\mathrm dt &=\\sum_E\'\\subseteq E\\int_0^1t^|E\'|(1-t)^m-|E\'|\\mathrm d t\\\\ &=\\sum_E\'\\subseteq E\\int_0^1\\sum_j=0^m-|E\'|(-1)^j\\binomm-|E\'|jt^j+|E\'|\\mathrm d t\\\\ &=\\sum_j=0^m\\sum_k=0^m-j(-1)^j\\binomm-kj\\left(\\int_0^1t^j+k\\mathrm d t\\right)\\left(\\sum_E\'\\subseteq E[|E\'|=k]\\right)\\\\ &=\\sum_j=0^m\\sum_k=0^m-j(-1)^j\\binomm-kj\\frac1j+k+1\\left(\\sum_E\'\\subseteq E[|E\'|=k]\\right) \\endaligned \\]

所以，我们尝试计算出 \\(\\sum_E\'\\subseteq E[|E\'|=k]\\)，整个问题也便迎刃而解了。注意 \\(E\'\\) 还需要保证不连通，这其实颇有一点麻烦——我们反之保证连通。设 \\(f_S,k\\) 表示在 \\(S\\subseteq V\\) 的导出子图中，能够保证 \\(S\\) 内点连通且大小恰好为 \\(k\\) 的边集的数量。设 \\(S\\) 的导出子图的边集为 \\(E(S)\\)。正难则反，我们可以想到容斥不连通的情况：

\\[f_S,k=\\binom|E(S)|k-\\sum_T\\subseteq S,T\\neq\\varnothing[\\min T=\\min S]\\sum_j=0^kf_T,j\\binom|E(S\\setminus T)|k-j \\]

如果直接暴力 DP 是 \\(O(3^nm^2)\\) 的。不过注意到，第二维可以看作是生成函数的指数，我们相当于在进行生成函数运算。而我们又知道，最终 \\(f_V\\) 的生成函数的指数一定是 \\(\\le m\\) 的，因此我们可以预先插入 \\(m+1\\) 个点值，最后再做拉格朗日插值得到每一项系数。

最终我们得到了一个 \\(O(3^nm+m^2)\\) 的算法。

---

实现还需要一点技巧：

为了方便，由于 \\(\\binom4522\\approx4.1\\times 10^12\\) 并不是很大，long long 也还装得下，我们可以选一个略大于 \\(\\binom4522\\) 的质数，在它的模域下做运算，这样就简单了许多。

最后，由于我们需要输出浮点数结果，运算精度也是我们需要考虑的问题。注意到，答案最终一定是在 \\((0,1)\\) 范围内，因此我们可以输出它 \\(\\bmod 1\\) 之后的结果。这样的话，像 \\(\\binomm-kj\\frac1j+k+1\\bmod 1\\) 这样的运算，就可以转化为 \\(\\frac1j+k+1\\left(\\binomm-kj\\bmod (j+k+1)\\right)\\)。我们可以直接将分子对分母取模，最后算的时候再计算 \\(ans \\bmod 1\\) 即可。

如果不是不调整精度过不了完全图，我会卡吗我？

小结：

1. 仍然需要注意常见的运算技巧。这里主要的处理期望的技巧，还是把贡献滚到前缀上，从而转化为概率。

2. 注意一下图上的子集 DP 的方法。一般来说，我们会倾向于通过容斥计算，在无向图的连通性和有向图的强连通中都可以这么做；另一方面，强连通图中也有用耳分解的算法。

3. 注意实现中的小技巧。尤其是处理精度这一块的，平时接触不多更要注意。

   一般来说，我们会把浮点数压缩到一个较小的范围内，从而保证精度。这里进行 \\(\\bmod 1\\) 及后续运算就是为了达成这一点。

   你也可以说，我是发现所有整数部分都被抵消了才想到了这个方法

代码

#include <cmath>
#include <cstdio>

#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;

const LL mod = 4116715363889;
const int MAXN = 105, MAXS = ( 1 << 10 ) + 5;

template<typename _T>
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>
void write( _T x ) 
    if( x < 0 ) putchar( \'-\' ), x = -x;
    if( 9 < x ) write( x / 10 );
    putchar( x % 10 + \'0\' );


int fin[MAXN];

LL f[MAXS], g[MAXS];
int cnt[MAXS];

LL C[MAXN][MAXN];

LL ways[MAXN], coe[MAXN];

int N, M;

inline LL Mul( const LL a, const LL b ) 
    return ( a * b - ( LL ) ( ( long double ) a / mod * b ) * mod + mod ) % mod;


inline LL Qkpow( LL, LL );
inline LL Inv( const LL a )  return Qkpow( a, mod - 2 ); 
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 Qkpow( LL base, LL indx ) 
    LL ret = 1;
    while( indx ) 
        if( indx & 1 ) ret = Mul( ret, base );
        base = Mul( base, base ), indx >>= 1;
    
    return ret;


LL Query( const int x ) 
    for( int S = 0 ; S < ( 1 << N ) ; S ++ ) 
        f[S] = 0;
        per( i, M, 0 ) f[S] = Add( Mul( f[S], x ), C[cnt[S]][i] );
        g[S] = f[S];
        for( int T = ( S - 1 ) & S ; T ; T = ( T - 1 ) & S )
            if( ( T & ( - T ) ) == ( S & ( - S ) ) )
                g[S] = Sub( g[S], Mul( g[T], f[S ^ T] ) );
    
    return g[( 1 << N ) - 1];


inline long double Mod1( const long double x )  
    return x - floorl( x );


int main() 
    read( N ), read( M );
    rep( i, 1, M ) 
        int u, v; read( u ), read( v );
        cnt[( 1 << ( u - 1 ) ) | ( 1 << ( v - 1 ) )] ++;
    
    rep( i, 0, M ) 
        C[i][0] = C[i][i] = 1;
        rep( j, 1, i - 1 ) C[i][j] = Add( C[i - 1][j], C[i - 1][j - 1] );
    
    for( int i = 0 ; i < N ; i ++ )
        for( int S = 0 ; S < ( 1 << N ) ; S ++ )
            if( ! ( S >> i & 1 ) ) cnt[S | ( 1 << i )] += cnt[S];
    coe[0] = 1;
    rep( i, 1, M + 1 )
        per( j, M + 1, 0 ) 
            coe[j] = Mul( coe[j], mod - i );
            if( j ) coe[j] = Add( coe[j], coe[j - 1] );
        
    rep( i, 1, M + 1 ) 
        LL inv = Inv( i ), tmp = 1;
        rep( j, 0, M + 1 ) 
            if( j ) coe[j] = Sub( coe[j], coe[j - 1] );
            coe[j] = Mul( coe[j], mod - inv );
        
        rep( j, 1, M + 1 ) if( i ^ j )
            tmp = Mul( tmp, Sub( i, j ) );
        tmp = Mul( Query( i ), Inv( tmp ) );
        rep( j, 0, M ) ways[j] = Add( ways[j], Mul( tmp, coe[j] ) );
        per( j, M + 1, 0 ) 
            coe[j] = Mul( coe[j], mod - i );
            if( j ) coe[j] = Add( coe[j], coe[j - 1] );
        
    
    long double ans = 0;
    rep( i, 0, M ) 
        rep( j, 0, M - i ) 
            int cur = i + j + 1;
            int res = 1ll * ( C[M - i][j] % cur ) * ( Sub( C[M][i], ways[i] ) % cur ) % cur;
            fin[cur] = j & 1 ? ( ( fin[cur] - res ) % cur + cur ) % cur : ( fin[cur] + res ) % cur;
        
    rep( i, 1, M + 1 ) ans = Mod1( ans + ( long double ) fin[i] / i );
    printf( "%.6Lf\\n", ans );
    return 0;