CF446D DZY Loves Games

https://www.luogu.com.cn/problem/CF446D

首先肯定是转换为求出$g[i][j]$表示从第$i$个陷阱到第$j$个陷阱，不经过其他陷阱（相当于是走一步）的概率；
1单独处理一下
然后把这个矩阵$k-2$次方即可(走k-1步)

然后现在来考虑怎么求出$g$

把所有关键点的出边全部删掉
记$f[i][j]$表示原图的第$i$个点(不是关键点)，到第$j$个关键点的概率
这个跑$n$次高斯消元（对于每个结束点跑一次）就可以得出来以每个关键点结束的概率

然后$g[i][j]$就枚举$i$的出边$v$然后通过f来求出

这样是$O(n^4)$，不太行

容易发现，对于每个结束点的矩阵（方程组），只有最后一列（常数）是不同的
所以可以吧常数项全部丢到后面，一次消元，一次性全部求出来

具体实现看代码吧
code:

#include<bits/stdc++.h>
#define N 605
using namespace std;
int tot;
struct MT {
    double a[105][105];
    void clear() {
        for(int i = 0; i <= tot; i ++)
            for(int j = 0; j <= tot; j ++) a[i][j] = 0;
    }
    void I() {
        for(int i = 1; i <= tot; i ++) a[i][i] = 1;
    }
} ans, b;
MT mul(MT a, MT b) {
    MT c; c.clear();
    for(int k = 1; k <= tot; k ++)
        for(int i = 1; i <= tot; i ++)
            for(int j = 1; j <= tot; j ++)
                c.a[i][j] += a.a[i][k] * b.a[k][j];
    return c;
}
MT qpow(MT x, int y) {
    MT ret; ret.clear(); ret.I();
    for(; y; y >>= 1, x = mul(x, x)) if(y & 1) ret = mul(ret, x);
    return ret;
}
double a[N][N];
int o[N], p[N], d[N], gs[N][N], n, m, k;
void Guass(int n) {
    for(int i = 1; i <= n; i ++) {
        int j = i;
        for(int k = i + 1; k <= n; k ++) if(abs(a[j][i]) < abs(a[k][i])) j = k;
        if(j != i) swap(a[j], a[i]);
        double d = a[i][i];
        for(int j = i; j <= n + tot; j ++) a[i][j] /= d;
        for(int j = 1; j <= n; j ++) {
            if(j == i) continue;
            double t = a[j][i] / a[i][i];
            for(int k = i; k <= n + tot; k ++)
                a[j][k] -= t * a[i][k];
        }
    }
}
int main() {
    scanf("%d%d%d", &n, &m, &k);
    for(int i = 1; i <= n; i ++) {
        scanf("%d", &o[i]);
        if(o[i]) p[++ tot] = i;
    }
    for(int i = 1; i <= m; i ++) {
        int u, v;
        scanf("%d%d", &u, &v);
        d[u] ++, d[v] ++; gs[u][v] ++, gs[v][u] ++;
    }
    for(int i = 1; i <= n; i ++) {
        if(!o[i]) {
            a[i][i] = -1;
            for(int j = 1; j <= n; j ++) a[i][j] += 1.0 * gs[i][j] / d[i];
        } else {
            a[i][i] = 1;
            int pos = 0;
            for(int j = 1; j <= tot; j ++) if(p[j] == i) pos = j;
            a[i][n + pos] = 1;
        }
    }
//    for(int i = 1; i <= n; i ++) {
//        for(int j = 1; j <= n + tot; j ++) printf("%.2lf ", a[i][j]); printf("\n");
//    }
    Guass(n);
    for(int i = 1; i <= tot; i ++) ans.a[1][i] = a[1][n + i];
    for(int i = 1; i <= tot; i ++)
        for(int j = 1; j <= tot; j ++) {
            for(int k = 1; k <= n; k ++) b.a[i][j] += gs[p[i]][k] * a[k][j + n];
            b.a[i][j] /= d[p[i]];
        }
    ans = mul(ans, qpow(b, k - 2));
    printf("%.9lf", ans.a[1][tot]);
    return 0;
}