地震后的幻想乡

题意

一幅图的每条边的权值为互不相关的 $[0,1]$内的连续随机变量。求图的最小生成树的最大边的权值的期望，结果保留 6 位小数。

概率密度函数

在研究连续型随机变量的过程中，我们一般用概率密度函数 $\rho_a(x)$描述一个随机变量的取值。

$\rho_a(x)$在 $x$处的取值越大，说明随机变量 $a$有更大的概率取 $x$。但是，由于我们研究的 $a$是一个连续型随机变量，虽说 $a$可能为 $x$，但它的值为 $x$的概率实际上是 $0$，这样也就不太好比较 $a$为某个值的概率的大小。因此为了方便，我们用一个更直观的方式表示 $a$取值的分布。

$$
P(a\in[l,r]) = \int_l^r\rho_a(x)dx
$$

这也恰好符合”$\rho_a(x)$越大，$a$越容易取得 $x$” 这一特性。同时，根据这个定义，我们知道，$\int_{-\infty}^{+\infty}\rho_a(x)dx = 1$。

问题规约

求最小生成树的最大边的权值的期望并不是一件容易的事。但是如下的一个问题实际上比较容易解决：原图是否存在一个最大边不大于 $x$的生成树。

这个问题相当于：每条边有 $x$的概率存在，有 $1-x$的概率不存在，求原图联通的概率。利用状态压缩 dp 可以在 $O(n3^n)$的时间复杂度内解决。

设 $F[S]$为原图内 $S$代表的点集联通的概率。那么我们选定该点集内的一个点为 “观察点”，我们尝试观察出这个点集中与这个点联通的点有哪些。如果并不是所有点都与观察点联通，那么这个点集就是不联通的。恰好，每个不连通的情况都不重复地对应着一个观察点所在的联通块。所以我们可以断言，点集 $S$联通的概率，等于 1 减去观察点所在的每一种联通块存在的概率 (这个联通块不等于 S)。

由于我们只需要求全图联通的概率，我们将观察点设为 1 号节点即可，这样已经可以补充不漏地计算到所有情况，即：
$$
\forall S\ (1\in S)\ \ : \ \ F[S] = 1-\sum_{T\subsetneqq S, 1\in T}F[T]G(S-T,T)
$$
其中 $G(A,B)$表示点集 $A,B$之间的所有边都断开的概率。

那么怎么把原问题规约到这个问题上来呢？如果我们枚举最小生成树的最大边的权值 $x$，那么每条边出现的概率就是 $x$，进而通过上面的状态转移方程求出的每一个状态都会是一个关于 $x$的多项式。而多项式 $F[S](x)$就表示 “S 的最小生成树最大边不超过 $x$” 的概率。

那么，我们又要求 $F[s](x)$的期望，怎么办呢？可以这样：

先对 $F(x)$求导，$F'(x)$表示最小生成树最大边等于 $x$的概率密度函数。为什么呢？反过来考虑：设 $f(x)$表示最小生成树最大边的概率密度函数，则 $F=\int_0^af(x)dx$表示的是最小生成树最大边不超过 $x$的概率。由于 $F(0)=0$，所以求导再定积分的结果就是原函数，这样证明是没有毛病的。

则 $\int_0^1xF'(x)$为最小生成树最大边期望。为什么呢？当我们在计算期望时，当函数取值为 $a$时，它对期望的贡献也为 $a$。因此当最大边为 $x$时，它的贡献就是 $x$，而最大边为 $x$的概率又是通过 $F'(x)$描述的。所以将两个多项式卷积，意味着点值的直接相乘，这样求出的就是期望。

所以总的来说，我们只需要将状态压缩 DP 的 $F[S]$设为一个多项式即可，原先的乘法、加法、减法全部用多项式的运算代替。最后利用上述方法求出期望就行了。

实现上的问题

在运算的过程中，自变量 $x$的系数都是整数，因此我们可以用 long long 存储。但是在最终求答案时，我们需要将多项式积分，这时会出现小数。因此需要将 long long 转化为浮点数。在转化过程中会丢失精度，再将转化后的数除以一个较小的整数，这时会严重影响结果。因此我们需要一个技巧避免精度问题：

对于较大的 $a$和较小的 $b$，如下的方法可以避免精度问题：
$$
\frac{a}{b}=\lfloor\frac{a}{b}\rfloor + \frac{a\mod b}{b}
$$

代码

#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <assert.h>
#define mov(x) (1<<(x))

using namespace std;

typedef long long ll;

struct POLY
{
    ll a[50];
    int d;

    POLY () {memset(a, 0, sizeof(a)); d = 0;}
    POLY deriv ()
    {
        POLY c; c.d = d-1;
        for(int i=0; i<c.d; i++) c.a[i] = a[i+1]*(i+1);
        return c;
    }
    POLY operator + (const POLY& t)
    {
        POLY c; c.d = max(d, t.d);
        for(int i=0; i<c.d; i++) c.a[i] = a[i]+t.a[i];
        return c;
    }
    POLY operator - (const POLY& t)
    {
        POLY c; c.d = max(d, t.d);
        for(int i=0; i<c.d; i++) c.a[i] = a[i]-t.a[i];
        return c;
    }
    POLY operator * (const POLY& t)
    {
        POLY c; c.d = d+t.d-1;
        for(int i=0; i<d; i++)
            for(int j=0; j<t.d; j++)
                c.a[i+j] += a[i]*t.a[j];
        return c;
    }
    double f(double x)
    {
        double y = 0, xi = 1;
        for(int i=0; i<d; i++) y += xi*a[i], xi*=x;
        return y;
    }
    void print()
    {
        for(int i=0; i<d; i++) printf("%lld ", a[i]);
        putchar('\n');
    }
};

POLY F[1024], MI[100];
int N, M, E[12], SE[1024][1024], BN[1024];

void dp()
{
    F[1].d = 1;
    F[1].a[0] = 1;
    for(int u=2; u<mov(N); u++)
    {
        if(!(u & 1)) continue;
        for(int s=(u-1)&u; s; s=(s-1)&u)
            if(s & 1)
                F[u] = F[u] + F[s]*MI[SE[s][u^s]];
        for(int i=0; i<F[u].d; i++) F[u].a[i] *= -1;
        F[u].a[0] ++;
    }
}

void calc()
{
    POLY pos = F[mov(N)-1];
    POLY rho = pos.deriv();
    POLY mul; mul.d = 2, mul.a[0] = 0, mul.a[1] = 1;
    POLY ans = rho*mul;
    long double prt = 0;
    for(int i=0; i<ans.d; i++) prt += (long double)(ans.a[i]%(i+1))/(long double)(i+1);
    while(prt < 0) prt++;
    while(prt >= 1) prt--;
    printf("%.6LF\n", prt);
}

void init()
{
    MI[0].d = 1;
    MI[0].a[0] = 1;
    MI[1].d = 2;
    MI[1].a[0] = 1, MI[1].a[1] = -1;
    for(int i=2; i<=M; i++) MI[i] = MI[i-1]*MI[1];
    for(int i=1; i<mov(N); i++) BN[i] = BN[i>>1] + (i&1);
    for(int s=1; s<mov(N); s++)
    {
        int u = (~s)&(mov(N)-1);
        for(int t=u; t; t=(t-1)&u)
            for(int i=0; i<N; i++)
                if(s & mov(i))
                    SE[s][t] += BN[E[i]&t];
    }
}

void input()
{
    scanf("%d%d", &N, &M);
    for(int i=0; i<M; i++)
    {
        int a, b;
        scanf("%d%d", &a, &b);
        E[a-1] |= mov(b-1);
        E[b-1] |= mov(a-1);
    }
}

int main()
{
    input();
    init();
    dp();
    calc();
    return 0;
}