1. 程式人生 > >洛谷3317 SDOI2014重建(高斯消元+期望)

洛谷3317 SDOI2014重建(高斯消元+期望)

qwq
一開始想了個錯的做法。

直接開始說比較正確的做法吧。
首先我們考慮題目的 a n s ans 該怎麼去求
我們令 x x

表示原圖中的某一條邊

a n s = x

t r e e p x
x   n o t t r e e ( 1 p x ) ans = \sum \prod_{x\in tree} p_x \prod_{x\ not\in tree} (1-p_x)

qwq而根據矩陣樹定理,我們可以求出來所有生成樹的邊權乘積的和,也就是前一部分。

現在我們考慮應該怎麼優化第二部分。
qwq
我們經過推理能發現,我們可以用總的除去在生成樹裡面的求出來不在生成樹裡面的。

也就是說
x   n o t t r e e ( 1 p x ) = ( 1 p i ) x t r e e ( 1 p j ) \prod_{x\ not \in tree} (1-p_x)= \frac{\prod (1-p_i)}{\prod_{x\in tree} (1-p_j)}

我們帶回原柿,然後把 p i \prod p_i 提出來

a n s = p x × x t r e e p x 1 p x ans = \prod p_x \times \sum \prod_{x \in tree} \frac{p_x}{1-p_x}

那麼現在,對於後面那一項,我們只需要把所有的邊都設成權值是 x t r e e p x 1 p x \prod_{x \in tree} \frac{p_x}{1-p_x}
然後每個 d [ i ] d[i] 表示與他連線的所有邊權的和。

直接跑矩陣樹定理就能求出來 s u m sum 啦,然後直接用一開始求的 p x \prod p_x ,一減就OK了

但是這裡有一個需要注意的地方就是當 p x p_x 等於 1 1 的時候,我們應該將他的權值設成 1 e p s 1-eps

因為當 p p 等於1的時候, 1 1 p > i n f \frac{1}{1-p} -> inf

然後有因為 1 e p s > i n f \frac{1}{eps}->inf

所以 p = 1 e p s p=1-eps

然後弄完權值直接跑矩陣樹定理就好

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<queue>
#include<map>
#include<set>
#define mk make_pair
#define ll long long
#include<ctime>
using namespace std;
inline int read()
{
  int x=0,f=1;char ch=getchar();
  while (!isdigit(ch)) {if (ch=='-') f=-1;ch=getchar();}
  while (isdigit(ch)) {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
  return x*f;
}
const int maxn = 110;
const double eps = 1e-6;
double a[maxn][maxn];
double d[maxn];
int n;
double ans=1;
void gauss()
{
	int k=1;
	for (int i=1;i<=n;i++)
	{
		int now = k;
		while(now<=n && fabs(a[now][i])<=eps)  now++;
		if (now==n+1) continue;
		for (int j=1;j<=n+1;j++) swap(a[now][j],a[k][j]);
		for (int j=1;j<=n;j++)
		{
			if (j!=k)
			{
				double t = a[j][i]/a[k][i];
				for (int p=1;p<=n+1;p++) a[j][p]-=t*a[k][p];
			}
		}
		++k;
	}
	for (int i=1;i<=n;i++)
	  ans=ans*a[i][i];
}
double ymh=1;
int main()
{
  n=read();
  for (int i=1;i<=n;i++)
    for (int j=1;j<=n;j++)
    {
    	double x;
    	scanf("%lf",&x);
    	if (x==1) x = 1-eps;
		if (i<j) ymh=ymh*(1-x);
    	x=x/(1-x);
    	a[i][j]=-x;
    	d[i]+=x;
    	//d[j]+=x;
	}
  for (int i=1;i<=n;i++) a[i][i]=d[i];
  gauss();
  printf("%.5lf",ans*ymh);
  return 0;
}