bzoj 1778 [Usaco2010 Hol]Dotp 驱逐猪猡(高斯消元)
2016-03-31 18:27
441 查看
【题意】
炸弹从1开始运动,每次有P/Q的概率爆炸,否则等概率沿边移动,问在每个城市爆炸的概率。
【思路】
设M表示移动一次后i->j的概率。Mk为移动k次后的概率,则有:
Mk=M^k
设S={ 1,0,0,0,… }
设pi为移动i步后到对应点爆炸的概率矩阵,则有:
p0=P/Q * S
p1=P/Q * S * M1
…
p+oo=P/Q * S * Mn
则答案为:sigma{ pi },0<=i<+oo
即:
Ans=P/Q * S * sigma{ M^i } ,0<=i<+oo
根据等差数列的求和公式:
Ans= P/Q * S * (I-M^(+oo))/(I-M)
= P/Q * S * I/(I-M)
=> Ans(I-M) = P/Q * S
其中I为单位矩阵,I[j][j]=1,其余为0。
于是可以得到n个线性方程,可以使用高斯消元法求解。
注意Ans*(I-M),所以第i个方程的第j项系数为-(1-P/Q)*Mji。
【代码】
炸弹从1开始运动,每次有P/Q的概率爆炸,否则等概率沿边移动,问在每个城市爆炸的概率。
【思路】
设M表示移动一次后i->j的概率。Mk为移动k次后的概率,则有:
Mk=M^k
设S={ 1,0,0,0,… }
设pi为移动i步后到对应点爆炸的概率矩阵,则有:
p0=P/Q * S
p1=P/Q * S * M1
…
p+oo=P/Q * S * Mn
则答案为:sigma{ pi },0<=i<+oo
即:
Ans=P/Q * S * sigma{ M^i } ,0<=i<+oo
根据等差数列的求和公式:
Ans= P/Q * S * (I-M^(+oo))/(I-M)
= P/Q * S * I/(I-M)
=> Ans(I-M) = P/Q * S
其中I为单位矩阵,I[j][j]=1,其余为0。
于是可以得到n个线性方程,可以使用高斯消元法求解。
注意Ans*(I-M),所以第i个方程的第j项系数为-(1-P/Q)*Mji。
【代码】
#include<set> #include<cmath> #include<queue> #include<vector> #include<cstdio> #include<cstring> #include<iostream> #include<algorithm> #define trav(u,i) for(int i=front[u];i;i=e[i].nxt) #define FOR(a,b,c) for(int a=(b);a<=(c);a++) using namespace std; typedef long long ll; const int N = 400; const int M = 1e5+10; ll read() { char c=getchar(); ll f=1,x=0; while(!isdigit(c)) { if(c=='-') f=-1; c=getchar(); } while(isdigit(c)) x=x*10+c-'0',c=getchar(); return x*f; } struct Edge { int v,nxt; }e[M]; int en=1,front ; void adde(int u,int v) { e[++en]=(Edge){v,front[u]}; front[u]=en; } int n,m,du ; double P,a ; void gause() { for(int i=1;i<=n;i++) { int r=i; for(int j=i+1;j<=n;j++) if(fabs(a[j][i])>fabs(a[r][i])) r=i; for(int j=1;j<=n+1;j++) swap(a[i][j],a[r][j]); for(int j=n+1;j>=i;j--) for(int k=i+1;k<=n;k++) a[k][j]-=a[k][i]/a[i][i]*a[i][j]; } for(int i=n;i;i--) { for(int j=i+1;j<=n;j++) a[i][n+1]-=a[i][j]*a[j][n+1]; a[i][n+1]/=a[i][i]; } } int main() { int x,y; n=read(),m=read(),x=read(),y=read(); P=(double)x/y; FOR(i,1,m) { x=read(),y=read(); adde(x,y),adde(y,x); du[x]++,du[y]++; } FOR(u,1,n) { a[u][u]=1; trav(u,i) { int v=e[i].v; a[u][v]-=(1.0-P)/du[v]; } } a[1][n+1]=P; gause(); FOR(i,1,n) printf("%.9lf\n",a[i][n+1]); return 0; }
相关文章推荐
- SCTP协议详解与实例
- thinkphp 中 ajax的使用
- PHPCMS_数据库配置
- 4.1 Zend_Config
- PHP中的连贯操作
- php检查漏洞防护补丁-防护XSS,SQL,文件包含等多种高危漏洞
- ThinkPHP处理海量数据分表机制详细代码
- php输出json格式数据的例子
- ZendFramework-1.12.17
- PHP C扩展初探
- yii框架开启debug和gii
- matplotlib-绘制精美的图表
- php 获取今日、昨日、上周、本月的起始时间戳和结束时间戳的方法
- yii报错400
- sphinx使用步骤
- php中常用的运算符
- php programmer should know?
- update the gedit for php programming(searched not try yet)
- 127报错解决方法
- PHP缓存Xcache安装