21航电5E - random walk2（高斯消元）

题目链接

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}\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];
            }
            
        }
    }
}

标签：21,int,sum,random,矩阵,times,neq,高斯消,define
来源： https://www.cnblogs.com/limil/p/15116927.html