Overview

This problem asks us to find the amount of positively charged particles in a specific container after \(T\) steps in a simulation where particles with positive and negative charges transition on a directed graph. Since \(T\) is extremely large, we compute this efficiently using matrix exponentiation.

Discussion

The transition of particles can be represented as linear calculations using the amounts of positive and negative particles at each step. Naively simulating this \(T\) times would take \(O(NT)\) time, which results in TLE (Time Limit Exceeded) since \(T \leq 10^{18}\).

Therefore, we consider representing the transitions using a matrix. For each container \(i\), we maintain two states: the “amount of positive particles” and the “amount of negative particles”, managing \(2N\) states in total.

Let’s consider the transition from container \(i\) to container \(j\). If there is an inversion rule (\(F_{i,j} = 1\)), positive particles arrive as negative, and negative particles arrive as positive. If there is no inversion rule (\(F_{i,j} = 0\)), they arrive with their charges unchanged.

By summarizing this into a \(2N \times 2N\) transition matrix \(A\), we can represent a single-step update as multiplying the state vector by the matrix \(A\). The state after \(T\) steps can be obtained by multiplying by the matrix \(A\) raised to the power of \(T\). Using matrix exponentiation, we can compute this efficiently in \(O((2N)^3 \log T)\).

Algorithm

1. State Definition: We represent the amount of positive particles in container \(i\) (\(0 \le i < N\)) at index \(2i\), and the amount of negative particles at index \(2i+1\).
2. Constructing the Transition Matrix: Using the transition probability \(P_{i,j}\) from container \(i\) to container \(j\) and the inversion flag \(F_{i,j}\), we set the elements of the transition matrix \(A\) as follows:
   • \(A_{2j, 2i} = (1 - F_{i,j}) P_{i,j}\)
   • \(A_{2j, 2i+1} = F_{i,j} P_{i,j}\)
   • \(A_{2j+1, 2i} = F_{i,j} P_{i,j}\)
   • \(A_{2j+1, 2i+1} = (1 - F_{i,j}) P_{i,j}\) Note that \(P_{i,j}\) is obtained by dividing \(R_{i,j}\) by the row sum. In modular arithmetic, this is calculated by multiplying by the modular multiplicative inverse of the row sum.
3. Matrix Exponentiation: Since the initial state has “1 positive particle in container \(S\)”, the \(2(S-1)\)-th element of the state vector is \(1\) and all other elements are \(0\). Thus, the desired amount of positive particles in container \(G\) after \(T\) steps is equal to the entry at \((2(G-1), 2(S-1))\) of the matrix \(A^T\).
4. Computing Matrix Exponentiation: We compute \(A^T\) using binary exponentiation (exponentiation by squaring) and output the corresponding entry.

Complexity

• Time Complexity: Let the size of the matrix be \(n = 2N\). Since one matrix multiplication takes \(O(n^3)\) and matrix exponentiation performs \(O(\log T)\) matrix multiplications, the total time complexity is \(O(N^3 \log T)\). This is fast enough since \(N \le 30\).
• Space Complexity: Since the size of the matrix is \(2N \times 2N\), the space complexity is \(O(N^2)\).

Implementation Points

• When performing division under a modulo, we compute it by multiplying by the “modular inverse of the row sum” using Fermat’s Little Theorem. We can utilize pow(s, MOD - 2, MOD).

• In the innermost loop of the matrix multiplication, skipping calculations when the element is \(0\) optimizes the constant factor.

• Care should be taken during implementation to avoid confusing the state indices (e.g., \(2i\), \(2i+1\)).

  Source Code

import sys

def solve():
    input = sys.stdin.readline
    N, M, T = map(int, input().split())
    S, G = map(int, input().split())
    R = [list(map(int, input().split())) for _ in range(N)]
    F = [[0] * N for _ in range(N)]
    for _ in range(M):
        u, v = map(int, input().split())
        F[u-1][v-1] = 1
    
    MOD = 10**9 + 7
    invS = [0] * N
    for i in range(N):
        s = sum(R[i]) % MOD
        invS[i] = pow(s, MOD - 2, MOD)
        
    n = 2 * N
    A = [[0] * n for _ in range(n)]
    for i in range(N):
        for j in range(N):
            p = R[i][j] * invS[i] % MOD
            f = F[i][j]
            # y_j^+ = (1-F) * x_i^+ * p + F * x_i^- * p
            A[2*j][2*i] = (1 - f) * p % MOD
            A[2*j][2*i+1] = f * p % MOD
            # y_j^- = F * x_i^+ * p + (1-F) * x_i^- * p
            A[2*j+1][2*i] = f * p % MOD
            A[2*j+1][2*i+1] = (1 - f) * p % MOD
            
    def mat_mul(A, B, n):
        C = [[0] * n for _ in range(n)]
        for i in range(n):
            Ai = A[i]
            Ci = C[i]
            for k in range(n):
                aik = Ai[k]
                if aik:
                    Bk = B[k]
                    for j in range(n):
                        Ci[j] = (Ci[j] + aik * Bk[j]) % MOD
        return C

    def mat_pow(A, p, n):
        res = [[0] * n for _ in range(n)]
        for i in range(n):
            res[i][i] = 1
        while p > 0:
            if p & 1:
                res = mat_mul(res, A, n)
            A = mat_mul(A, A, n)
            p >>= 1
        return res
    
    res = mat_pow(A, T, n)
    ans = res[2*(G-1)][2*(S-1)]
    print(ans)

if __name__ == '__main__':
    solve()

---

This editorial was generated by or-glm5.2-high.