【BZOJ4637】期望 Kruskal+矩阵树定理
2024-08-25 21:19:37
【BZOJ4637】期望
Description
在米国有一所大学,名叫万国歌剧与信息大学(UniversalOperaandInformaticasUniversity)。简称UOI大学。UOI大学的建筑与道路分布很有趣,每对建筑之间有且仅有一条直接或者间接的路径相连,更加明确的说,就是成树形分布。其实在设计时,对于大学的N个建筑,总共有M条道路可以修建,每条道路都有一个距离值Disti和一个美学值Valuei。一个设计方案的距离值和美学值定义为该设计方案内包含的道路的距离值与美学值之和。投资方的要求只有设计方案的距离值最小。大神出于对树的喜爱所以决定设计方案必须是一棵树。因为要参加UOI,所以当时大神就急急忙忙地随机选择了一个合法的方案。但其实存在很多合法的方案,假设每种设计方案取的概率是均等的,那么设计方案的美学值期望是多少?
Input
第一行两个整数,N和M,意义如上所述。
第二行到第M+1行,每行4个整数,Xi,Yi,Disti,Valuei,分别表示这条道路连接的两个建筑的编号,距离值以及美学值。
输入保证至少有一种合法方案。
100%的数据保证N<=10000M<=200000
100%的数据保证距离值相同的道路数小于30,同时不保证没有重边。
Output
一行一个整数,即满足总道路长度最小的情况下,设计方案的美学值期望。要求保留5位小数
Sample Input
2 1
1 2 3 4
1 2 3 4
Sample Output
4.00000
题解:傻题细节多啊~
我们先进行Kruskal求距离值的最小生成树,如果有多条边权值相同,则我们将它们放到一起处理。我们再把加入后会被分到同一个连通块中的边放到一起,并把连通块离散化缩成点。因为期望是可加的,所以我们可以枚举其中的每条边,这条边出现的概率就是总的缩点后的图的生成树数目 分之 保证这条边在内时剩余图的生成树数目。拿矩阵树定理算一下就好了,用long double即可过。
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <vector>
#include <cmath>
using namespace std;
const int maxn=10010;
const int maxm=200010;
typedef long double db;
const db eps=1e-6;
struct edge
{
int a,b,c,d;
}p[maxm];
int n,m,tot,tp,now;
int f[maxn],bel[maxn],g[65],vis[maxn],st[40],pa[40],pb[40],qa[40],qb[40],ref[65];
int qs[65][40],qt[65],ps[65][40],pt[65];
db v[40][40],ans;
bool cmp(const edge &a,const edge &b)
{
return a.c<b.c;
}
inline int find(int x)
{
return (f[x]==x)?x:(f[x]=find(f[x]));
}
inline int gind(int x)
{
return (g[x]==x)?x:(g[x]=gind(g[x]));
}
inline int rd()
{
int ret=0,f=1; char gc=getchar();
while(gc<'0'||gc>'9') {if(gc=='-') f=-f; gc=getchar();}
while(gc>='0'&&gc<='9') ret=ret*10+(gc^'0'),gc=getchar();
return ret*f;
}
inline double gauss(int N)
{
db ret=1;
int i,j,k;
for(i=1;i<N;i++)
{
for(j=k=i;j<N;j++) if(fabs(v[j][i])>fabs(v[k][i])) k=j;
if(k!=i) for(ret=-ret,j=1;j<N;j++) swap(v[i][j],v[k][j]);
if(fabs(v[i][i])<1e-6) return 0;
for(j=1;j<N;j++) if(j!=i&&fabs(v[j][i])>1e-6)
{
db tmp=v[j][i]/v[i][i];
for(k=1;k<N;k++) v[j][k]-=v[i][k]*tmp;
}
ret=ret*v[i][i];
}
return ret;
}
inline void calc(int x)
{
int i,j,a,b;
for(i=1;i<=pt[x];i++) ref[ps[x][i]]=i;
for(i=1;i<=qt[x];i++) qa[i]=ref[pa[qs[x][i]]],qb[i]=ref[pb[qs[x][i]]];
memset(v,0,sizeof(v));
for(i=1;i<=qt[x];i++) a=qa[i],b=qb[i],v[a][a]++,v[b][b]++,v[a][b]--,v[b][a]--;
double tmp=gauss(pt[x]);
for(i=1;i<=qt[x];i++)
{
memset(v,0,sizeof(v));
if(qa[i]>qb[i]) swap(qa[i],qb[i]);
for(j=1;j<=qt[x];j++) if(i!=j)
{
a=qa[j],b=qb[j];
if(a==qb[i]) a=qa[i];
if(a>qb[i]) a--;
if(b==qb[i]) b=qa[i];
if(b>qb[i]) b--;
v[a][a]++,v[b][b]++,v[a][b]--,v[b][a]--;
}
ans+=gauss(pt[x]-1)*p[st[qs[x][i]]].d/tmp;
}
qt[x]=pt[x]=0;
}
inline void solve()
{
int i,a,b;
now++,tot=0;
for(i=1;i<=tp;i++)
{
if(vis[find(p[st[i]].a)]!=now) vis[f[p[st[i]].a]]=now,bel[f[p[st[i]].a]]=++tot;
if(vis[find(p[st[i]].b)]!=now) vis[f[p[st[i]].b]]=now,bel[f[p[st[i]].b]]=++tot;
pa[i]=bel[f[p[st[i]].a]],pb[i]=bel[f[p[st[i]].b]];
}
for(i=1;i<=tot;i++) g[i]=i;
for(i=1;i<=tp;i++)
{
a=gind(pa[i]),b=gind(pb[i]);
if(a!=b) g[a]=b;
}
for(i=1;i<=tp;i++) a=gind(pa[i]),qs[a][++qt[a]]=i;
for(i=1;i<=tot;i++) a=gind(i),ps[a][++pt[a]]=i;
for(i=1;i<=tot;i++) if(gind(i)==i) calc(i);
for(i=1;i<=tp;i++)
{
a=find(p[st[i]].a),b=find(p[st[i]].b);
if(a!=b) f[a]=b;
}
tp=0;
}
int main()
{
//freopen("bz4637.in","r",stdin);
n=rd(),m=rd();
int i,a,b,pre=0;
for(i=1;i<=m;i++) p[i].a=rd(),p[i].b=rd(),p[i].c=rd(),p[i].d=rd();
sort(p+1,p+m+1,cmp);
for(i=1;i<=n;i++) f[i]=i;
for(i=1;i<=m;i++)
{
if(p[i].c>pre&&pre) solve();
a=find(p[i].a),b=find(p[i].b);
if(a==b) continue;
st[++tp]=i,pre=p[i].c;
}
solve();
printf("%.5Lf",ans);
return 0;
}//3 3 1 2 1 4 1 3 1 6 2 3 1 8
最新文章
- BZOJ 2093: [Poi2010]Frog
- The 2015 China Collegiate Programming Contest D.Pick The Sticks hdu 5543
- javaSE第二十二天
- 基于css3的3D立方体旋转特效
- 在VC6.0下如何调用Delphi5.0开发的进程内COM
- ionic2 安装(一)
- vue拦截器实现统一token,并兼容IE9验证
- 老婆大人 split,slice,splice,replace的用法
- 关于Python中迭代器的作用
- 《linux就该这么学》第十二节课:第10章,Apache网站服务
- Eclipse常用快捷键(用到想到随时更新)
- docker常用命令整理-在容器中使用service命令
- Topcoder SRM590 Fox And City
- BZOJ 5194--[Usaco2018 Feb]Snow Boots(STL)
- ASP.NET MVC:无法向会话状态服务器发出会话状态请求
- vuex2.0+两个小例子
- aix-syslog
- leetcode130
- Node.js缓存
- Linux 的计划任务
热门文章
- 【iOS开发】关于显示一连串图片组成动画效果UIImageView的使用
- Go语言之高级篇beego框架之model设计构造查询
- boost.property_tree解析xml的帮助类以及中文解析问题的解决(转)
- 这些年我在技术路上做过最虚伪愚蠢的事情,就是在CSDN上刷屏赚分
- sql语句语法
- linux nfs共享文件
- [Aaronyang] 写给自己的WPF4.5 笔记9[复杂数据处理三步曲,数据展示ListView泪奔2/3]
- WEBAPI 的简单示例
- 申请IPV6地址配置IPV6域名
- 菜鸟教程之工具使用(七)——从GIt上导出Maven项目