21航电5E

题目链接

Problem - 7016

题解

设矩阵(F)为从(i)出发到(j)停止的概率(对应(f_{i,j}))，矩阵(G)为从(i)出发到(j)无数次的概率之和(对应(g_{i,j}))，概率矩阵为P(对应(p_{i,j}))。

对于矩阵(F)容易得到：

[f_{i,j}=g_{i,j} imes p_{j,j} ]

对于矩阵(G)可得：

[g_{i,j}=sumlimits_{k eq j}{g_{i,k} imes p_{k,j}}+[i=j] ]

([i=j])意思是从(i)出发有个第0次到达(i)的概率为1，所以要额外加1。

观察式子，可以发现和矩阵乘法非常像，只是多了一个(k eq j)的条件。只需在原本矩阵乘积的基础上减去(k=j)的情况即可。

[g_{i,j}=sumlimits_{k=1}^n{g_{i,k} imes p_{k,j}}+[i=j]-g_{i,j} imes p_{j,j} \ g_{i,j}+g_{i,j} imes p_{j,j}-[i=j]=sumlimits_{k=1}^n{g_{i,k} imes p_{k,j}} ]

令矩阵(D)为

[d_{i,j}=left{ egin{aligned} p_{i,j} & & i=j \ 0 & & i eq j \ end{aligned} ight. ]

可得((E)为单位矩阵)

[F=G imes D \ G+G imes D-E = G imes P ]

化简得

[G=(E+D-P)^{-1} \ F=G imes D ]

直接套矩阵求逆的板子即可

#include <bits/stdc++.h>

#define endl '
'
#define IOS std::ios::sync_with_stdio(0); cin.tie(0); cout.tie(0)
#define mp make_pair
#define seteps(N) fixed << setprecision(N) 
typedef long long ll;

using namespace std;
/*-----------------------------------------------------------------*/

ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}
#define INF 0x3f3f3f3f

const int N = 3e3 + 10;
const int M = 998244353;
const double eps = 1e-5;

inline ll qpow(ll a, ll b, ll m) {
    ll res = 1;
    while(b) {
        if(b & 1) res = (res * a) % m;
        a = (a * a) % m;
        b = b >> 1;
    }
    return res;
}

ll p[N][N << 1];
ll d[N][N];

bool Gauss(ll a[][N << 1], int n) { //高斯消元求矩阵的逆
    for(int i = 1; i <= n; i++) {
        for(int j = 1; j <= n; j++) {
            a[i][j + n] = 0;
        }
        a[i][i + n] = 1;
    }
    for(int i = 1; i <= n; i++) {
        int r = i;
        for(int j = i + 1; j <= n; j++) {
            if(a[j][i] > a[r][i]) r = j;
        }
        if(r != i) swap(a[i], a[r]);
        if(!a[i][i]) return false;
        ll inv = qpow(a[i][i], M - 2, M);
        for(int j = 1; j <= n; j++) {
            if(j == i) continue;
            ll da = a[j][i] * inv % M; //a[j][i]可能会被更新，所以要先保存
            for(int k = i; k <= (n << 1); k++) {
                ll t = a[i][k] * da % M;
                a[j][k] = ((a[j][k] - t) % M + M) % M;
            }
        }
        for(int j = i; j <= (n << 1); j++) a[i][j] = a[i][j] * inv % M;
    }
    return true;
}


int main() {
    IOS;
    int t;
    cin >> t;
    while(t--) {
        int n;
        cin >> n;
        for(int i = 1; i <= n; i++) {
            ll sum = 0;
            for(int j = 1; j <= n; j++) {
                d[i][j] = 0;
                cin >> p[i][j];
                sum += p[i][j];
            }
            sum = qpow(sum, M - 2, M);
            for(int j = 1; j <= n; j++) {
                p[i][j] = p[i][j] * sum % M;
            }
        }
        for(int i = 1; i <= n; i++) {
            for(int j = 1; j<= n; j++) {
                if(i == j) {
                    d[i][j] = p[i][j];
                    p[i][j] = 1;
                } else {
                    p[i][j] = -p[i][j];
                }
            }
        }
        Gauss(p, n);
        for(int i = 1; i <= n; i++) {
            for(int j = 1; j <= n; j++) {
                ll res = 0;
                res = p[i][j + n] * d[j][j] % M;
                res %= M;
                cout << res << " 
"[j == n];
            }
            
        }
    }
}