概要

\(N\) 個の容器間で粒子が遷移し、特定の経路では電荷（正・負）が反転するシステムにおいて、\(T\) ステップ後に特定の容器にある正の粒子の量を求める問題です。\(T\) が非常に大きいため、状態遷移を行列で表現し、行列累乗を用いて高速に計算します。

考察

まず、素朴に \(T\) 回の遷移をシミュレーションすることを考えます。各ステップで \(N\) 個の容器それぞれに対して、他のすべての容器からの遷移量を足し合わせる必要があります。1ステップあたり \(O(N^2)\) の計算が必要となり、\(T\) が最大 \(10^{18}\) であるため、\(O(N^2 T)\) では実行時間制限に間に合いません。

このような「状態遷移を何度も繰り返す」問題では、遷移を行列の掛け算として表現し、行列累乗を用いることで計算量を \(O(\text{(状態数)}^3 \log T)\) に落とすことができます。

各容器には「正の粒子の量」と「負の粒子の量」の2つの状態があるため、全体の状態数は \(2N\) となります。 状態ベクトル \(\vec{x}\) を長さ \(2N\) のベクトルとし、\(\vec{x}_{2i}\) を容器 \(i\) の正の粒子の量、\(\vec{x}_{2i+1}\) を容器 \(i\) の負の粒子の量とします。

1ステップの遷移は、\(2N \times 2N\) の遷移行列 \(A\) を用いて \(\vec{x}_{new} = A \vec{x}_{old}\) と表せます。 容器 \(i\) から容器 \(j\) への遷移割合を \(P_{i,j} = \frac{R_{i,j}}{\sum_{j'} R_{i,j'}}\) とします。 - 電荷の反転がない場合：容器 \(i\) の正の粒子は容器 \(j\) の正の粒子として、負の粒子は負の粒子として遷移します。 すなわち、\(A_{2j, 2i} = P_{i,j}\), \(A_{2j+1, 2i+1} = P_{i,j}\) - 電荷の反転がある場合：正の粒子は負の粒子として、負の粒子は正の粒子として遷移します。 すなわち、\(A_{2j, 2i+1} = P_{i,j}\), \(A_{2j+1, 2i} = P_{i,j}\)

初期状態では容器 \(S\) に正の粒子が1単位あるため、\(\vec{x}^{(0)}_{2S} = 1\)、それ以外は \(0\) です。 \(T\) ステップ後の状態は \(\vec{x}^{(T)} = A^T \vec{x}^{(0)}\) となります。 求める値は容器 \(G\) の正の粒子の量、つまり \(\vec{x}^{(T)}_{2G}\) であり、これは遷移行列の累乗 \(A^T\) の \((2G, 2S)\) 成分に等しくなります。

アルゴリズム

1. 各遷移の割合 \(P_{i,j}\) を求めます。問題文より答えを有理数として扱い、最終的に \(\bmod 10^9+7\) で出力するため、割り算はフェルマーの小定理を利用したモジュラ逆元（modinv）を掛けることで行います。
2. \(2N \times 2N\) の遷移行列 \(A\) を構築します。上記の考察に従い、反転ルール \((i, j)\) の有無に応じて適切な位置に \(P_{i,j}\) を格納します。
3. 行列累乗のアルゴリズム（繰り返し二乗法）を用いて \(A^T\) を計算します。計算はすべて \(\bmod 10^9+7\) で行います。
4. 計算結果の行列 \(A^T\) の \((2G, 2S)\) 成分が答えとなります。

計算量

• 時間計算量: \(O(N^3 \log T)\)
  • 行列サイズは \(2N\)、行列乗算が \(O((2N)^3)\)、繰り返し二乗法で \(\log T\) 回の乗算を行うため。
• 空間計算量: \(O(N^2)\)
  • \(2N \times 2N\) の行列を数個保持するため。

実装のポイント

• 遷移行列のインデックスに注意が必要です。状態ベクトルを \(\vec{x}_{new} = A \vec{x}_{old}\) とすると、行列 \(A\) の行インデックスが遷移先（\(j\)）、列インデックスが遷移元（\(i\)）に対応します。コード内では row1, row2 が遷移先、col1, col2 が遷移元を表しています。

• 行列乗算の内部ループで if (A.a[i][k] == 0) continue; とすることで、遷移行列はスパース（疎）なため計算を一部省略できます。これにより定数倍が改善されます。

  ソースコード

#include <iostream>
#include <vector>

using namespace std;

const int MOD = 1e9 + 7;

struct Matrix {
    int n;
    vector<vector<long long>> a;
    Matrix(int n) : n(n), a(n, vector<long long>(n, 0)) {}
};

Matrix multiply(const Matrix &A, const Matrix &B) {
    int n = A.n;
    Matrix C(n);
    for (int i = 0; i < n; ++i) {
        for (int k = 0; k < n; ++k) {
            if (A.a[i][k] == 0) continue;
            for (int j = 0; j < n; ++j) {
                C.a[i][j] = (C.a[i][j] + A.a[i][k] * B.a[k][j]) % MOD;
            }
        }
    }
    return C;
}

Matrix mat_pow(Matrix A, long long p) {
    int n = A.n;
    Matrix B(n);
    for (int i = 0; i < n; ++i) B.a[i][i] = 1;
    while (p > 0) {
        if (p & 1) B = multiply(B, A);
        A = multiply(A, A);
        p >>= 1;
    }
    return B;
}

long long modinv(long long a) {
    long long b = MOD, u = 1, v = 0;
    while (b) {
        long long t = a / b;
        a -= t * b; swap(a, b);
        u -= t * v; swap(u, v);
    }
    u %= MOD;
    if (u < 0) u += MOD;
    return u;
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);
    int N, M;
    long long T;
    if (!(cin >> N >> M >> T)) return 0;
    int S, G;
    cin >> S >> G;
    --S; --G;

    vector<vector<long long>> R(N, vector<long long>(N));
    vector<long long> sumR(N, 0);
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
            cin >> R[i][j];
            sumR[i] += R[i][j];
        }
    }

    vector<vector<bool>> F(N, vector<bool>(N, false));
    for (int k = 0; k < M; ++k) {
        int u, v;
        cin >> u >> v;
        --u; --v;
        F[u][v] = true;
    }

    int sz = 2 * N;
    Matrix A(sz);
    for (int i = 0; i < N; ++i) {
        long long inv_sum = modinv(sumR[i]);
        for (int j = 0; j < N; ++j) {
            long long P = R[i][j] * inv_sum % MOD;
            int row1 = 2 * j;
            int row2 = 2 * j + 1;
            int col1 = 2 * i;
            int col2 = 2 * i + 1;
            if (F[i][j]) {
                A.a[row1][col1] = 0;
                A.a[row1][col2] = P;
                A.a[row2][col1] = P;
                A.a[row2][col2] = 0;
            } else {
                A.a[row1][col1] = P;
                A.a[row1][col2] = 0;
                A.a[row2][col1] = 0;
                A.a[row2][col2] = P;
            }
        }
    }

    Matrix B = mat_pow(A, T);

    // x^(0) has 1 at 2*S, 0 elsewhere.
    // So the answer is B.a[2*G][2*S]
    cout << B.a[2 * G][2 * S] << endl;

    return 0;
}

この解説は or-glm5.2-high によって生成されました。