BZOJ 3640: JC的小苹果 [概率DP 高斯消元 矩阵求逆]

3640: JC的小苹果

题意：求1到n点权和(le k)的概率

---

sengxian orz的题解好详细啊

容易想到(f[i][j])表示走到i点权为j的概率

按点权分层，可以DP

但是对于(val[i]=0)的点，就不是DAG了，必须使用高斯消元

每层消元一次？复杂度(O(SN^3))，boom！！！

发现每次的系数矩阵一样啊

[Ax=b ightarrow x=A^{-1}b ]

我们求出(A)矩阵的逆，然后直接让常数向量乘逆就可以了，因为常数矩阵是向量，一次的复杂度(O(N^2))

然后就可以(O(SN^2))了

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;
typedef unsigned long long ll;
const int N=155, M=1e4+5;
const double eps=1e-8;
inline int read() {
	char c=getchar(); int x=0, f=1;
	while(c<'0' || c>'9') {if(c=='-')f=-1; c=getchar();}
	while(c>='0' && c<='9') {x=x*10+c-'0'; c=getchar();}
	return x*f;
}

int n, m, s, val[N], u, v;
struct edge{int v, ne;}e[M];
int cnt=1, h[N], de[N];
inline void ins(int u, int v) { e[++cnt]=(edge){v, h[u]}; h[u]=cnt; }
struct Matrix {
	double a[N][N];
	Matrix() {memset(a, 0, sizeof(a));}
	double* operator [](int x) {return a[x];}
	void im() {for(int i=1; i<=n; i++) a[i][i]=1;}
	void print() {for(int i=1; i<=n; i++) for(int j=1; j<=n; j++) printf("%lf%c",a[i][j], j==n?'
':' ');puts("");}
}a, c;
Matrix inverse(Matrix a) {
	Matrix c; c.im();
	for(int i=1; i<=n; i++) {
		int r=i;
		for(int j=i; j<=n; j++) if(abs(a[j][i])>abs(a[r][i])) r=j;
		if(r!=i) for(int j=1; j<=n; j++) swap(a[i][j], a[r][j]), swap(c[i][j], c[r][j]);

		double t = a[i][i];
		for(int j=1; j<=n; j++) a[i][j]/=t, c[i][j]/=t;

		for(int k=1; k<=n; k++) if(k!=i) {
			double t = a[k][i];
			for(int j=1; j<=n; j++) a[k][j] -= t*a[i][j], c[k][j] -= t*c[i][j];
		}
	}
	return c;
}

double b[N], f[N][M];
void dp() {
	double ans=0;
	for(int now=s; now>0; now--){
		memset(b, 0, sizeof(b));
		if(now==s) b[1]=1;
		for(int u=1; u<=n; u++) if(val[u] && now+val[u]<=s) 
			for(int i=h[u];i;i=e[i].ne) if(e[i].v != n) 
				b[u] += f[e[i].v][now+val[u]] / de[e[i].v];

		for(int i=1; i<=n; i++) 
			for(int j=1; j<=n; j++) f[i][now] += a[i][j] * b[j];

		ans += f[n][now];
	}
	printf("%.8lf", ans);
}
int main() {
	freopen("in","r",stdin);
	n=read(); m=read(); s=read();
	for(int i=1; i<=n; i++) val[i]=read();
	for(int i=1; i<=m; i++) {
		u=read(); v=read();
		ins(u, v); de[u]++;
		if(u!=v) ins(v, u), de[v]++;
	}
	for(int u=1; u<=n; u++) {
		a[u][u]=1;
		if(!val[u]) for(int i=h[u];i;i=e[i].ne) if(e[i].v != n) a[u][e[i].v] -= 1.0/de[e[i].v];
	}
	a = inverse(a);
	dp();
}