概要

\(N\) 個の容器の上を粒子が拡散し、各容器に「正電荷」「負電荷」の粒子量が管理されます。\(T\) ステップ後の容器 \(G\) における正電荷粒子の量を、線形変換の繰り返しとして捉え、行列累乗で高速に求める問題です。

考察

状態は線形変換で表せる

この問題の最大のポイントは、1 ステップの更新が線形変換（行列とベクトルの積）で表せることです。

各容器 \(i\) には「正電荷量 \(x_i^+\)」と「負電荷量 \(x_i^-\)」の 2 種類の値があります。これらをすべて並べた長さ \(2N\) のベクトルを状態とみなします。

更新式

\[y_j^{+} = \sum_{i} \Bigl( (1 - F_{i,j}) x_i^{+} P_{i,j} + F_{i,j} x_i^{-} P_{i,j} \Bigr)\]

\[y_j^{-} = \sum_{i} \Bigl( (1 - F_{i,j}) x_i^{-} P_{i,j} + F_{i,j} x_i^{+} P_{i,j} \Bigr)\]

を見ると、\(y\) の各成分は \(x\) の各成分の1 次式（係数倍の和）になっています。つまり、ある \(2N \times 2N\) 行列 \(A\) を用いて

\[\mathbf{y} = A \mathbf{x}\]

と書けます。これを \(T\) 回繰り返すと

\[\mathbf{x}_T = A^T \mathbf{x}_0\]

となります。

素朴なアプローチの問題点

1 ステップずつ愚直に更新すると \(O(T \cdot N^2)\) かかりますが、\(T \leq 10^{18}\) なので絶対に間に合いません。そこで 行列累乗（繰り返し二乗法） を使い、\(A^T\) を \(O(\log T)\) 回の行列積で求めます。

有理数の扱い

\(P_{i,j} = R_{i,j} / (\sum_{j'} R_{i,j'})\) には割り算が含まれます。答えは有理数ですが、出力は \(\bmod\ 10^9+7\) での値です。\(10^9+7\) は素数なので、割り算はモジュラ逆元（フェルマーの小定理による \(a^{-1} \equiv a^{p-2}\)）に置き換えて計算します。行和の逆元をあらかじめ求めておけば、すべての計算を \(\bmod\ 10^9+7\) の整数演算で行えます。

アルゴリズム

状態のインデックス付け

長さ \(2N\) のベクトルにおいて、 - インデックス \(i\)（\(0 \leq i < N\)）：容器 \(i\) の正電荷量 - インデックス \(N+i\)：容器 \(i\) の負電荷量

とします。

遷移行列 \(A\) の構築

各 \((i, j)\) について \(P_{i,j} = R_{i,j} \cdot (\text{行和})^{-1}\) を計算し、 - 反転なし（\(F_{i,j}=0\)）： - \(A[j][i] \mathrel{+}= P_{i,j}\)（正→正） - \(A[N+j][N+i] \mathrel{+}= P_{i,j}\)（負→負） - 反転あり（\(F_{i,j}=1\)）： - \(A[j][N+i] \mathrel{+}= P_{i,j}\)（負→正） - \(A[N+j][i] \mathrel{+}= P_{i,j}\)（正→負）

ここで \(A[行][列]\) は「\(\mathbf{y} = A\mathbf{x}\)」の行列を表します。

答えの取り出し

初期ベクトル \(\mathbf{x}_0\) は「容器 \(S\) の正電荷」が \(1\)、他はすべて \(0\) です。つまり \(\mathbf{x}_0\) はインデックス \(S\) だけが \(1\) の単位ベクトルです。

\(\mathbf{x}_T = A^T \mathbf{x}_0\) の容器 \(G\) の正電荷量はインデックス \(G\) の成分なので、これは

\[\mathbf{x}_T[G] = (A^T)[G][S]\]

つまり \(A^T\) の \((G, S)\) 成分そのものです。行列累乗の結果から 1 要素を取り出すだけで答えが得られます。

計算量

行列のサイズは \(D = 2N \leq 60\) です。

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

\(N \leq 30\)、\(T \leq 10^{18}\) なので、\(60^3 \times 60 \approx 1.3 \times 10^7\) 程度で十分高速です。

実装のポイント

• 割り算は逆元で：行和 \(\sum_{j} R_{i,j}\) の逆元を inv で事前計算し、\(P_{i,j} = R_{i,j} \cdot \text{invrow}[i]\) とします。すべての値を \(\bmod\ 10^9+7\) で保持します。

• \(P_{i,j}=0\) のスキップ：\(R_{i,j}=0\) の遷移は寄与しないので飛ばします（反転ルールが設定されていても影響しません）。

• 行列の向きに注意：状態を「行列 \(\times\) 列ベクトル」として扱う場合、更新式の係数をそのまま \(A[\text{出力先}][\text{入力元}]\) に格納します。最終的に答えは \((G, S)\) 成分です。

• 同時更新：1 ステップの更新は「更新前の値のみ」から計算されますが、行列 \(\times\) ベクトルの形式は自動的にこの同時更新を表現しているため、特別な配慮は不要です。

• \(T\) が大きい：\(T \leq 10^{18}\) なので long long で受け取り、繰り返し二乗法を用います。

  ソースコード

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

const long long MOD = 1000000007;

long long pw(long long a, long long b){
    a%=MOD; if(a<0)a+=MOD;
    long long r=1;
    while(b>0){
        if(b&1) r=r*a%MOD;
        a=a*a%MOD;
        b>>=1;
    }
    return r;
}
long long inv(long long a){ return pw(a, MOD-2); }

int main(){
    int N, M; long long T;
    scanf("%d %d %lld", &N, &M, &T);
    int S, G;
    scanf("%d %d", &S, &G);
    S--; G--;
    vector<vector<long long>> R(N, vector<long long>(N));
    for(int i=0;i<N;i++)for(int j=0;j<N;j++) scanf("%lld",&R[i][j]);
    vector<vector<int>> F(N, vector<int>(N,0));
    for(int k=0;k<M;k++){
        int u,v; scanf("%d %d",&u,&v); u--;v--;
        F[u][v]=1;
    }
    // P[i][j] = R[i][j] * inv(rowsum_i)
    vector<long long> invrow(N);
    for(int i=0;i<N;i++){
        long long s=0;
        for(int j=0;j<N;j++) s+=R[i][j];
        invrow[i]=inv(s%MOD);
    }
    int D=2*N;
    // state: positive i -> i, negative i -> N+i
    vector<vector<long long>> A(D, vector<long long>(D,0));
    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++){
            long long P = (R[i][j]%MOD)*invrow[i]%MOD;
            if(P==0) continue;
            if(F[i][j]==0){
                // no flip
                A[j][i] = (A[j][i] + P)%MOD;      // y_j+ from x_i+
                A[N+j][N+i] = (A[N+j][N+i] + P)%MOD; // y_j- from x_i-
            } else {
                A[j][N+i] = (A[j][N+i] + P)%MOD;  // y_j+ from x_i-
                A[N+j][i] = (A[N+j][i] + P)%MOD;  // y_j- from x_i+
            }
        }
    }
    // matrix power A^T
    auto mul=[&](vector<vector<long long>>&X, vector<vector<long long>>&Y){
        vector<vector<long long>> Z(D, vector<long long>(D,0));
        for(int i=0;i<D;i++){
            for(int k=0;k<D;k++){
                if(X[i][k]==0) continue;
                long long v=X[i][k];
                for(int j=0;j<D;j++){
                    Z[i][j]=(Z[i][j]+v*Y[k][j])%MOD;
                }
            }
        }
        return Z;
    };
    vector<vector<long long>> result(D, vector<long long>(D,0));
    for(int i=0;i<D;i++) result[i][i]=1;
    long long e=T;
    vector<vector<long long>> base=A;
    while(e>0){
        if(e&1) result=mul(result,base);
        base=mul(base,base);
        e>>=1;
    }
    // initial vector: e[S]=1 (positive at S)
    // after T steps: y = result * e, y[G] = result[G][S]
    printf("%lld\n", result[G][S]%MOD);
    return 0;
}

---

この解説は claude4.8opus-high によって生成されました。