题目链接
题解
设矩阵\(F\)为从\(i\)出发到\(j\)停止的概率(对应\(f_{i,j}\))，矩阵\(G\)为从\(i\)出发到\(j\)无数次的概率之和(对应\(g_{i,j}\))，概率矩阵为P(对应\(p_{i,j}\))。
对于矩阵\(F\)容易得到：
\[f_{i,j}=g_{i,j}\times p_{j,j} \]对于矩阵\(G\)可得：
\[g_{i,j}=\sum\limits_{k\neq j}{g_{i,k}\times p_{k,j}}+[i=j] \]\([i=j]\)意思是从\(i\)出发有个第0次到达\(i\)的概率为1，所以要额外加1。
观察式子，可以发现和矩阵乘法非常像，只是多了一个\(k\neq j\)的条件。只需在原本矩阵乘积的基础上减去\(k=j\)的情况即可。
\[g_{i,j}=\sum\limits_{k=1}^n{g_{i,k}\times p_{k,j}}+[i=j]-g_{i,j}\times p_{j,j} \\ g_{i,j}+g_{i,j}\times p_{j,j}-[i=j]=\sum\limits_{k=1}^n{g_{i,k}\times p_{k,j}} \]令矩阵\(D\)为
\[d_{i,j}=\left\{ \begin{aligned} p_{i,j} & & i=j \\ 0 & & i \neq j \\ \end{aligned} \right. \]可得(\(E\)为单位矩阵)
\[F=G\times D \\ G+G\times D-E = G \times P \]化简得
\[G=(E+D-P)^{-1} \\ F=G\times D \]直接套矩阵求逆的板子即可
#include <bits/stdc++.h>
#define endl '\n'
#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 << " \n"[j == n];
}
}
}
}