bzoj4637: 期望

Description

在米国有一所大学，名叫万国歌剧与信息大学（UniversalOperaandInformaticasUniversity）。简称UOI大学。UO

I大学的建筑与道路分布很有趣，每对建筑之间有且仅有一条直接或者间接的路径相连，更加明确的说，就是成树

形分布。其实在设计时，对于大学的N个建筑，总共有M条道路可以修建，每条道路都有一个距离值Disti和一个美

学值Valuei。一个设计方案的距离值和美学值定义为该设计方案内包含的道路的距离值与美学值之和。投资方的要

求只有设计方案的距离值最小。大神出于对树的喜爱所以决定设计方案必须是一棵树。因为要参加UOI，所以当时

大神就急急忙忙地随机选择了一个合法的方案。但其实存在很多合法的方案,假设每种设计方案取的概率是均等的

，那么设计方案的美学值期望是多少？

Input

第一行两个整数，N和M，意义如上所述。

第二行到第M+1行，每行4个整数，Xi，Yi，Disti，Valuei，分别表示这条道路连接的

两个建筑的编号，距离值以及美学值。

输入保证至少有一种合法方案。

100%的数据保证N<=10000M<=200000

100%的数据保证距离值相同的道路数小于30，同时不保证没有重边。

Output

 一行一个整数，即满足总道路长度最小的情况下，设计方案的美学值期望。要求保留5位小数

按边权升序加边同时缩点，忽略缩点产生的自环，对同一权值且 加入此权值的边后在同个联通块内 的一组边，用矩阵树定理计算出生成树个数以及删去每条边后的生成树个数，于是可得一条边在最小生成树中的概率，统计答案。似乎要用long double才能过，但计算行列式时无论用long double还是模意义下的整数计算都能过。

#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long double ld;
#define double ld
int _(){
    int x=0,c=getchar(),f=1;
    while(c<48)c=='-'&&(f=-1),c=getchar();
    while(c>47)x=x*10+c-48,c=getchar();
    return x*f;
}
int es[2007],enx[2007],e0[2007],ep=2,ev[2007];
struct edge{
    int a,b,c,d;
    void init(){
        a=_();b=_();c=_();d=_();
    }
    void adde(){
        es[ep]=b;enx[ep]=e0[a];ev[ep]=d;e0[a]=ep++;
        es[ep]=a;enx[ep]=e0[b];ev[ep]=d;e0[b]=ep++;
    }
}e[200007];
bool operator<(edge a,edge b){
    return a.c<b.c;
}
double ans=0;
int n,m,f[100007],id[100007],idp=0,idt[100007],tk=0,ed[100007],ID[1007],IDP;
int v[33][33],v1,ee[2007][3],eep,os[2007][2],op=0;
int get(int x){
    int a=x,c;
    while(x!=f[x])x=f[x];
    while(x!=f[a])c=f[a],f[a]=x,a=c;
    return x;
}
void gid(int x){
    if(idt[x]!=tk)idt[x]=tk,id[x]=++idp;
}
void dfs(int w){
    if(ed[w]!=tk)ed[w]=tk,ID[w]=++IDP;
    for(int i=e0[w],u;i;i=enx[i]){
        u=es[i];
        if(ed[u]!=tk)dfs(u);
        if(ID[w]<ID[u]){
            ++v[ID[w]][ID[w]];
            ++v[ID[u]][ID[u]];
            --v[ID[w]][ID[u]];
            --v[ID[u]][ID[w]];
            ee[eep][0]=ID[w];
            ee[eep][1]=ID[u];
            ee[eep++][2]=ev[i];
        }
    }
}
const double _0=1e-9;
double solve(int n){
    static double a[33][33];
    double s=1;
    for(int i=1;i<=n;++i)
        for(int j=1;j<=n;++j)
            a[i][j]=v[i][j];
    for(int i=1;i<=n;++i){
        if(fabs(a[i][i])<_0){
            int t=i;
            for(int j=i+1;j<=n;++j)if(fabs(a[j][i])>fabs(a[t][i]))t=j;
            if(fabs(a[t][i])<_0)return 0;
            for(int j=i;j<=n;++j)swap(a[i][j],a[t][j]);
        }
        for(int j=i+1;j<=n;++j)if(a[j][i]){
            double x=a[j][i]/a[i][i];
            for(int k=i;k<=n;++k)a[j][k]-=x*a[i][k];
        }
    }
    for(int i=1;i<=n;++i)s*=a[i][i];
    return s;
}
void chk(int x){
    if(ed[x]==tk)return;
    IDP=0;++tk;eep=0;
    dfs(x);
    double v0=solve(IDP-1);
    for(int i=0;i<eep;++i){
        int x=ee[i][0],y=ee[i][1];
        --v[x][x],--v[y][y],++v[x][y],++v[y][x];
        double v1=v0-solve(IDP-1);
        ans+=v1/v0*ee[i][2];
        ++v[x][x],++v[y][y],--v[x][y],--v[y][x];
    }
    for(int i=1;i<=IDP;++i){
        for(int j=1;j<=IDP;++j)v[i][j]=0;
    }
}
int main(){
    n=_();m=_();
    for(int i=1;i<=n;++i)f[i]=i;
    for(int i=0;i<m;++i)e[i].init();
    std::sort(e,e+m);
    for(int i=0,j=0;i<m;){
        for(++tk,idp=0;j<m&&e[i].c==e[j].c;++j);
        for(int k=i;k<j;++k){
            int x=get(e[k].a),y=get(e[k].b);
            if(x==y){
                e[k].a=-1;
                continue;
            }
            os[op][0]=e[k].a;os[op++][1]=e[k].b;
            gid(x);gid(y);
            e[k].a=id[x];
            e[k].b=id[y];
            e[k].adde();
        }
        while(op)--op,f[get(os[op][0])]=get(os[op][1]);
        for(;i<j;++i)if(~e[i].a){
            chk(e[i].a);
            chk(e[i].b);
        }
        for(int t=1;t<=idp;++t)e0[t]=0;
        ep=2;
    }
    printf("%.5Lf",ans);
    return 0;
}