概要

粒子の「容器」と「電荷の正負」を合わせて \(2N\) 個の状態とみなし、1 ステップの更新を行列で表します。
あとは、その遷移行列を \(T\) 乗することで、\(T\) ステップ後の状態を高速に求めます。

考察

各容器には正の粒子と負の粒子があります。
そのため、単に容器番号だけを見るのではなく、次のように状態を分けて考えます。

• 状態 \(i\): 容器 \(i\) にある正の粒子
• 状態 \(N+i\): 容器 \(i\) にある負の粒子

つまり、全部で \(2N\) 個の状態があります。

例えば、容器 \(i\) から容器 \(j\) へ粒子が移動するとき、

• 反転ルールがない場合
  • 正の粒子は正のまま
  • 負の粒子は負のまま
• 反転ルールがある場合
  • 正の粒子は負になる
  • 負の粒子は正になる

となります。

この移動はすべて線形です。
つまり、現在の状態ベクトルを \(v\)、1 ステップの遷移を表す行列を \(A\) とすると、1 ステップ後の状態は

\( vA \)

で表せます。

したがって、\(T\) ステップ後の状態は

\( vA^T \)

です。

---

素朴に \(T\) 回シミュレーションすると、\(T\) は最大で \(10^{18}\) なので到底間に合いません。
また、粒子量を有理数としてそのまま管理すると、分母・分子が非常に大きくなります。

そこで、値はすべて \(\bmod\ 10^9+7\) で扱います。
確率

\( P_{i,j} = \frac{R_{i,j}}{\sum_k R_{i,k}} \)

は、分母の逆元を使って

\( R_{i,j} \times \left(\sum_k R_{i,k}\right)^{-1} \)

として計算できます。

\(10^9+7\) は素数なので、分母 \(a\) の逆元はフェルマーの小定理により

\( a^{10^9+5} \bmod 10^9+7 \)

で求められます。

アルゴリズム

まず、サイズ \(2N \times 2N\) の遷移行列 \(A\) を作ります。

インデックスは次のように対応させます。

• 容器 \(i\) の正の粒子: \(i\)
• 容器 \(i\) の負の粒子: \(N+i\)

ただし、実装では \(0\) 始まりにしているため、容器 \(i\) は \(0 \leq i < N\) とします。

---

容器 \(i\) から容器 \(j\) への遷移確率を

\( p = \frac{R_{i,j}}{\sum_k R_{i,k}} \)

とします。

この値を \(\bmod\ 10^9+7\) 上で計算します。

反転ルールがない場合

正は正へ、負は負へ移動します。

つまり、

• 容器 \(i\) の正 \(\rightarrow\) 容器 \(j\) の正
• 容器 \(i\) の負 \(\rightarrow\) 容器 \(j\) の負

なので、行列には次のように加えます。

\( A[i][j] += p \)

\( A[N+i][N+j] += p \)

反転ルールがある場合

正は負へ、負は正へ移動します。

つまり、

• 容器 \(i\) の正 \(\rightarrow\) 容器 \(j\) の負
• 容器 \(i\) の負 \(\rightarrow\) 容器 \(j\) の正

なので、行列には次のように加えます。

\( A[i][N+j] += p \)

\( A[N+i][j] += p \)

---

初期状態では、容器 \(S\) に正の粒子が \(1\) だけあるので、状態ベクトル \(v\) は

\( v[S] = 1 \)

で、それ以外は \(0\) です。

あとは二分累乗で \(A^T\) を計算しながら、必要なときだけベクトルに掛けます。

具体的には、\(T\) の各ビットを見て、

• そのビットが立っていれば \(v = vA\)
• 毎回 \(A = A^2\)

とします。

最終的に、容器 \(G\) の正の粒子量は \(v[G]\) です。

計算量

状態数を \(K = 2N\) とします。

• 行列同士の積: \(O(K^3)\)
• 行列累乗の回数: \(O(\log T)\)
• ベクトルと行列の積: \(O(K^2)\)

したがって全体の計算量は

• 時間計算量: \(O((2N)^3 \log T) = O(N^3 \log T)\)
• 空間計算量: \(O((2N)^2) = O(N^2)\)

\(N \leq 30\) なので、行列サイズは最大でも \(60 \times 60\) であり、十分高速に計算できます。

実装のポイント

• 遷移行列は「現在の状態」を行、「次の状態」を列として作っています。

  • そのため、状態ベクトルは行ベクトルとして扱い、更新は \(v = vA\) です。

• 確率の分母である行和の逆元を、繰り返し二乗法で求めます。

• \(T\) は最大 \(10^{18}\) なので、long long 型で扱います。

• 反転ルールは有向なので、flip[u][v] のみを true にします。逆向きは別扱いです。

  ソースコード

#include <bits/stdc++.h>
using namespace std;

static const long long MOD = 1000000007LL;

using Matrix = vector<vector<long long>>;

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

vector<long long> multiply_vec(const vector<long long>& v, const Matrix& A) {
    int n = (int)v.size();
    vector<long long> res(n, 0);
    for (int i = 0; i < n; i++) {
        if (v[i] == 0) continue;
        long long vi = v[i];
        for (int j = 0; j < n; j++) {
            res[j] = (res[j] + vi * A[i][j]) % MOD;
        }
    }
    return res;
}

long long mod_pow(long long a, long long e) {
    long long r = 1;
    while (e > 0) {
        if (e & 1) r = r * a % MOD;
        a = a * a % MOD;
        e >>= 1;
    }
    return r;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int N, M;
    long long T;
    cin >> N >> M >> T;

    int S, G;
    cin >> S >> G;
    --S;
    --G;

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

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

    int SZ = 2 * N;
    Matrix A(SZ, vector<long long>(SZ, 0));

    for (int i = 0; i < N; i++) {
        long long inv_sum = mod_pow(row_sum[i] % MOD, MOD - 2);
        for (int j = 0; j < N; j++) {
            long long p = (R[i][j] % MOD) * inv_sum % MOD;
            if (p == 0) continue;

            int pos_i = i;
            int neg_i = N + i;
            int pos_j = j;
            int neg_j = N + j;

            if (flip[i][j]) {
                A[pos_i][neg_j] = (A[pos_i][neg_j] + p) % MOD;
                A[neg_i][pos_j] = (A[neg_i][pos_j] + p) % MOD;
            } else {
                A[pos_i][pos_j] = (A[pos_i][pos_j] + p) % MOD;
                A[neg_i][neg_j] = (A[neg_i][neg_j] + p) % MOD;
            }
        }
    }

    vector<long long> v(SZ, 0);
    v[S] = 1;

    while (T > 0) {
        if (T & 1) v = multiply_vec(v, A);
        A = multiply(A, A);
        T >>= 1;
    }

    cout << v[G] << '\n';
    return 0;
}

---

この解説は gpt-5.5-high によって生成されました。