[洛谷P3343][ZJOI2015]地震后的幻想乡

这道题我用了一个特别神奇的做法，然后被卡精度卡成90分了，然后我不要face地用了__float128过了此题。

一看就是一道积分题。

对于任意实数$x\in [0,1]$，我们可以求出最小生成树的最大边权不超过x的概率，如果我们把x看作未知数，那么对于这张图我们可以得到一个关于x的式子$F(x)$

那么现在原问题可以分解为两个小问题：

如何求出$F(x)$

已知$F(x)$怎么计算答案

先来解决第一个。考虑如果x是一个确定的实数而不是未知数怎么做。我们设$f(sta)$为包含sta状态的节点的最小生成树最大边不超过x的概率。那么$f(sta)$相当于是计算用不超过x的边能使图联通的概率。

我们可以用总概率1减去不合法的概率来进行转移。

如何计算不合法的概率？考虑枚举编号最小的节点所在联通块有哪些点，然后子集枚举进行转移。为了节省时间，我们不妨只处理所有包含1号点的sta即可。那么有：

$$f(sta)=1-\displaystyle\sum_{sta1\subseteq sta,1\in sta1} f(sta1)\times (1-x)^{(\displaystyle\sum_{i=1}^n\displaystyle\sum_{j=1}^n [i\in sta1][j \in sta-sta1])}$$

这时我们发现如果把x作为未知数，那么算出来的一定是一个多项式，我们用$f(sta)$维护多项式，然后用多项式运算代替一部分实数运算即可。

由于多项式最高次项是$n^2$级别的，所以这部分时间复杂度是$O(3^n n^4)$的，算了一下$n=10$的时候大约是5亿，应该可以使用FFT优化，但是仔细一想发现没有这个必要，因为这个东西完全跑不满，常数绝对不会超过$\frac{1}{10}$，所以放心暴力算即可。

那么得到了$F(x)$，我们可以直接计算答案了，答案就是：

$$\int_{x=0}^1 F(x)-F(x-dx)$$

设$F(x)=\displaystyle\sum_{i=0}^k a_ix^i$

那么可得答案为：

$$\int_{x=0}^1 \displaystyle\sum_{i=0}^k a_ix^i-a_i(x-dx)^i$$

$$=\int_{x=0}^1 \displaystyle\sum_{i=0}^k a_i\displaystyle\sum_{j=1}^i C_i^j dx^j(-1)^{j-1}$$

$$=\int_{x=0}^1 \displaystyle\sum_{i=0}^k ia_idx$$

$$=\displaystyle\sum_{i=0}^k \frac{i}{i+1}a_i$$

于是就做完了。

但是这个方法直接写会炸精度，得到90分的高分（虽说如果是正式比赛这个分数也够了），然后我face不要地用了__float128过了（事实上用高精度也是可以的，如果你愿意去写实数的高精度那也行，反正我这个菜鸡是不想写）

附上我用了__float128的代码：

#include<bits/stdc++.h>
using namespace std;
typedef __float128 ld;
struct Poly{
    vector<ld> a;
    Poly(ld x0=0,ld x1=0){
        a.clear();
        if(x0!=0) a.push_back(x0);
        if(x1!=0) a.push_back(x1);
    }
    unsigned size(){ return a.size(); }
    void resize(int n){ a.resize(n); }
    inline ld& operator[](int x){ return a[x]; }
    inline Poly operator+(Poly b){
        Poly c; c.a.resize(max((unsigned)a.size(),b.size()));
        for(int i=0;i<c.size();++i)
            c[i]=(i<a.size()?a[i]:0)+(i<b.size()?b[i]:0);
        while(c.size()>0&&c.a.back()==0) c.a.pop_back();
        return c;
    }
    inline Poly operator-(){
        Poly c; c.a.resize(a.size());
        for(int i=0;i<a.size();++i) c[i]=-a[i];
        return c;
    }
    inline Poly operator*(Poly b){
        if(a.size()==0||b.size()==0) return Poly(0,0);
        Poly c; c.a.resize(a.size()+b.size()-1);
        for(int i=0;i<a.size();++i)
            for(int j=0;j<b.size();++j)
                c[i+j]+=a[i]*b[j];
        while(c.size()>0&&c.a.back()==0) c.a.pop_back();
        return c;
    }
    inline long double solve(){
        ld ret=0;
        for(int i=0;i<a.size();++i)
            ret+=a[i]*i/(i+1);
        return ret;
    }
};
Poly f[1<<10|5],g[1<<10|5],p[105];
bool e[15][15],mark[1<<10|5];
struct DSF{
    int fa[15];
    void init(int n){ for(int i=1;i<=n;++i) fa[i]=i; }
    int find_set(int x){ return fa[x]=(fa[x]==x?x:find_set(fa[x])); }
    void merge(int x,int y){ fa[find_set(x)]=find_set(y); }
}dsf;
bool getbit(int sta,int i){ return sta&(1<<i-1); }
int main(){
    int n,m; scanf("%d%d",&n,&m);
    for(int i=1;i<=m;++i){
        int a,b; scanf("%d%d",&a,&b);
        e[a][b]=1; e[b][a]=1;
    }
    for(int sta=0;sta<1<<n;++sta){
        dsf.init(n); mark[sta]=1;
        for(int i=1;i<=n;++i)
            if(getbit(sta,i))
                for(int j=1;j<=n;++j)
                    if(getbit(sta,j))
                        if(e[i][j]) dsf.merge(i,j);
        int last=0;
        for(int i=1;i<=n;++i)
            if(getbit(sta,i)){
                if(last==0) last=dsf.find_set(i);
                if(last!=dsf.find_set(i)) mark[sta]=0;
            }
    }
    p[0]=Poly(1,0);
    for(int i=1;i<=n*n;++i)
        p[i]=p[i-1]*Poly(1,-1);
    f[1]=Poly(1,0);
    for(int sta=2;sta<1<<n;++sta){
        if(!getbit(sta,1)||!mark[sta]) continue;
        f[sta]=Poly(1,0);
        for(int sta1=(sta-1)&sta;sta1>0;sta1=(sta1-1)&sta){
            if(!getbit(sta1,1)) continue;
            int cnt=0;
            for(int i=1;i<=n;++i)
                if(getbit(sta1,i))
                    for(int j=1;j<=n;++j)
                        if(e[i][j]&&!getbit(sta1,j)&&getbit(sta,j))
                            ++cnt;
            f[sta]=-f[sta1]*p[cnt]+f[sta];
        }
    }
    printf("%.6Lf",f[(1<<n)-1].solve());
    return 0;
}