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 \rightarrow 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?'\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();
}