概要

本問題は、接続された容器間を移動する粒子の拡散をシミュレーションし、指定されたステップ数 \(T\) 経過した後の、特定の容器における「正の電荷を持つ粒子の量」を求める問題です。 ステップ数 \(T\) が最大で \(10^{18}\) と非常に大きいため、愚直なシミュレーションではなく、行列累乗（ダブリング）を用いて高速に計算します。

考察

1. 状態の表現

各容器には「正の電荷を持つ粒子」と「負の電荷を持つ粒子」の2つの状態が存在します。 容器の数が \(N\) 個あるため、システム全体の状態数は \(2N\) 個と考えることができます。

状態を以下のように \(2N\) 次元の列ベクトル \(X\) で表します。 * \(0\) 以上 \(N-1\) 番目の要素：容器 \(1\) から \(N\) の正の粒子の量 * \(N\) 以上 \(2N-1\) 番目の要素：容器 \(1\) から \(N\) の負の粒子の量

初期状態（ステップ \(0\)）では、容器 \(S\) に正の粒子が \(1\) だけ存在し、他はすべて \(0\) です。これを 0-indexed で表すと、ベクトル \(X_0\) の \(S-1\) 番目の要素が \(1\) で、それ以外は \(0\) となります。

2. 遷移の定式化（行列による表現）

1ステップの遷移は線形変換として表すことができます。 容器 \(i\) から容器 \(j\) への遷移割合 \(P_{i,j}\) は、\(\sum_{j'} R_{i,j'}\) の逆元を \(10^9+7\) の下で計算することで、以下のように求められます。 $\(P_{i,j} = R_{i,j} \times \left( \sum_{j'=1}^{N} R_{i,j'} \right)^{-1} \pmod{10^9+7}\)$

この遷移割合と反転ルール \(F_{i,j}\)（反転するなら \(1\)、しないなら \(0\)）を用いて、ステップ \(t\) から \(t+1\) への遷移を表す \(2N \times 2N\) の遷移行列 \(A\) を構築します。 行列 \(A\) の各成分 \(A_{row, col}\) は、「状態 \(col\) から状態 \(row\) への遷移割合」を表します。

• 反転ルールがない場合 (\(F_{i,j} = 0\)): 電荷は変化しないため、正は正へ、負は負へ遷移します。

  • \(A_{j, i} = P_{i,j}\) （正 \(\to\) 正）
  • \(A_{N+j, N+i} = P_{i,j}\) （負 \(\to\) 負）

• 反転ルールがある場合 (\(F_{i,j} = 1\)): 電荷が反転するため、正は負へ、負は正へ遷移します。

  • \(A_{N+j, i} = P_{i,j}\) （正 \(\to\) 負）
  • \(A_{j, N+i} = P_{i,j}\) （負 \(\to\) 正）

3. 行列累乗による高速化

ステップ \(t\) における状態ベクトルを \(X_t\) とすると、1ステップの遷移は \(X_{t+1} = A X_t\) と表されます。 したがって、\(T\) ステップ後の状態ベクトル \(X_T\) は、行列の累積を用いて以下のように表せます。 $\(X_T = A^T X_0\)$

初期状態 \(X_0\) は \(S-1\) 番目の要素のみが \(1\) であるため、\(A^T X_0\) の結果は、行列 \(A^T\) の \(S-1\) 列目そのものになります。 最終的に求めたい値は、容器 \(G\) における正の粒子の量（\(X_T\) の \(G-1\) 番目の要素）なので、これは行列 \(A^T\) の \((G-1, S-1)\) 成分、すなわち A_T[G-1][S-1] と一致します。

\(T \le 10^{18}\) ですが、繰り返し二乗法（ダブリング）を用いることで、\(A^T\) を \(O(N^3 \log T)\) の計算量で求めることができます。

アルゴリズム

1. 逆元の計算: 各行 \(i\) について、行和 \(S_i = \sum_{j} R_{i,j}\) を計算し、フェルマーの小定理を用いて逆元 \(S_i^{-1} \equiv S_i^{MOD-2} \pmod{MOD}\) を求めます。これを用いて各遷移確率 \(P_{i,j}\) を計算します。
2. 遷移行列 \(A\) の構築: サイズ \(2N \times 2N\) の正方行列 \(A\) を用意し、反転ルールの有無に応じて確率 \(P_{i,j}\) を適切に配置します。
3. 行列の累乗計算: 繰り返し二乗法を用いて、遷移行列 \(A\) の \(T\) 乗である \(A^T\) を計算します。
4. 答えの出力: \(A^T\) の \(G-1\) 行 \(S-1\) 列目の要素を出力します。

計算量

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

  • 行列のサイズは \(2N \times 2N\) です。
  • 行列同士の積の計算には \(O((2N)^3) = O(N^3)\) の時間がかかります。
  • 繰り返し二乗法により、行列の積を約 \(\log_2 T\) 回行います。
  • \(N \le 30\)、\(\log_2 T \le 60\) であるため、全体の演算回数は約 \(60^3 \times 60 \approx 1.3 \times 10^7\) 回となり、実行時間制限に十分間に合います。

• 空間計算量: \(O(N^2)\)

  • サイズ \(2N \times 2N\) の行列を数個保持するため、メモリ使用量は極めてわずかです。

実装のポイント

• 行列積の高速化: Pythonなどのインタプリタ言語では、多重ループの処理が遅くなりやすいため、行列積の計算時に「遷移元の値が \(0\) である場合は内側のループをスキップする（枝刈り）」処理を入れることで、実行時間を大幅に短縮できます。

• 1-indexed から 0-indexed への変換: 入力される容器の番号 \(S, G\) や、反転ルールの \(u_k, v_k\) は 1-indexed（1から始まる番号）なので、配列のインデックスとして扱うためにあらかじめ \(1\) を引いて 0-indexed に変換しておきます。

  ソースコード

import sys

def solve():
    input = sys.stdin.read
    data = input().split()
    if not data:
        return
    
    N = int(data[0])
    M = int(data[1])
    T = int(data[2])
    S = int(data[3])
    G = int(data[4])
    
    idx = 5
    R = []
    for i in range(N):
        row = []
        for j in range(N):
            row.append(int(data[idx]))
            idx += 1
        R.append(row)
        
    F = [[0] * N for _ in range(N)]
    for _ in range(M):
        u = int(data[idx]) - 1
        v = int(data[idx+1]) - 1
        idx += 2
        F[u][v] = 1
        
    MOD = 10**9 + 7
    N2 = 2 * N
    
    A = [[0] * N2 for _ in range(N2)]
    for i in range(N):
        sum_ri = sum(R[i])
        inv_sum = pow(sum_ri, MOD - 2, MOD)
        for j in range(N):
            p = (R[i][j] * inv_sum) % MOD
            if p == 0:
                continue
            if F[i][j]:
                A[j][N + i] = p
                A[N + j][i] = p
            else:
                A[j][i] = p
                A[N + j][N + i] = p
                
    def mat_mul(X, Y, N2=N2, MOD=MOD):
        Z = [[0] * N2 for _ in range(N2)]
        for i in range(N2):
            row_X = X[i]
            row_Z = Z[i]
            for k in range(N2):
                x_ik = row_X[k]
                if x_ik == 0:
                    continue
                row_Y = Y[k]
                for j in range(N2):
                    row_Z[j] = (row_Z[j] + x_ik * row_Y[j]) % MOD
        return Z

    def mat_pow(X, p):
        res = [[0] * N2 for _ in range(N2)]
        for i in range(N2):
            res[i][i] = 1
        base = X
        while p > 0:
            if p & 1:
                res = mat_mul(res, base)
            base = mat_mul(base, base)
            p >>= 1
        return res

    A_T = mat_pow(A, T)
    ans = A_T[G - 1][S - 1]
    print(ans)

if __name__ == '__main__':
    solve()

---

この解説は gemini-3.5-flash-high によって生成されました。