Summary

This problem asks us to efficiently simulate a diffusion process where each particle has a “location” and a “charge” up to \(T\) steps.
Since \(T\) is very large (up to \(10^{18}\)), we solve this using matrix exponentiation.

Analysis

For each container, if we directly track the amount of positive particles \(x_i^+\) and negative particles \(x_i^-\), the number of states would be \(2N\).

However, the problem becomes simpler if we consider the following two quantities instead:

• Total quantity: \(a_i = x_i^+ + x_i^-\)
• Difference quantity: \(b_i = x_i^+ - x_i^-\)

In this case, the amount of positive particles can be restored later as:

\[ x_i^+ = \frac{a_i + b_i}{2} \]

About the total quantity

Even if the charge is inverted, the total quantity of particles does not change.

Therefore, if particles move from container \(i\) to container \(j\) with a transition probability \(P_{i,j}\), the total quantity is simply updated as:

\[ a'_j = \sum_i P_{i,j} a_i \]

About the difference quantity

The difference quantity \(b_i = x_i^+ - x_i^-\) flips its sign when the charge is inverted.

For example, if the charge is inverted during the transition from container \(i\) to \(j\):

• Positive particles become negative
• Negative particles become positive

Thus, the difference after arrival is:

\[ x_i^- P_{i,j} - x_i^+ P_{i,j} = - (x_i^+ - x_i^-) P_{i,j} = -b_i P_{i,j} \]

If the charge is not inverted, it is simply:

\[ b_i P_{i,j} \]

Therefore, the difference quantity can be written as:

\[ b'_j = \sum_i c_{i,j} P_{i,j} b_i \]

where

\[ c_{i,j} = \begin{cases} 1 & \text{if not inverted} \\ -1 & \text{if inverted} \end{cases} \]

Why a naive approach is too slow

Simulating step-by-step requires calculating all transitions at each step, which takes \(O(TN^2)\) time.

Since \(T\) is up to \(10^{18}\), this is far too slow to pass within the time limit.

Therefore, we represent the single-step update as a matrix and calculate \(T\) steps all at once using matrix exponentiation.

Algorithm

We represent the state as a column vector.

Let us define the matrix \(A\) for the total quantity as follows:

\[ A_{j,i} = P_{i,j} \]

This represents the transition probability from container \(i\) to container \(j\).

Then, the total quantity vector \(a\) can be updated as:

\[ a' = A a \]

Next, we define the matrix \(B\) for the difference quantity:

\[ B_{j,i} = \begin{cases} P_{i,j} & \text{if not inverted} \\ -P_{i,j} & \text{if inverted} \end{cases} \]

Then, the difference quantity vector \(b\) can be updated as:

\[ b' = B b \]

In the initial state, container \(S\) has exactly \(1\) positive particle, so we have:

\[ a_S = 1,\quad b_S = 1 \]

and all other elements are \(0\).

After \(T\) steps, we have:

\[ a^{(T)} = A^T a^{(0)} \]

\[ b^{(T)} = B^T b^{(0)} \]

Finally, the amount of positive particles in container \(G\) can be obtained by:

\[ x_G^+ = \frac{a_G^{(T)} + b_G^{(T)}}{2} \]

Handling Fractions

The transition probability is:

\[ P_{i,j} = \frac{R_{i,j}}{\sum_k R_{i,k}} \]

Since we need to output the answer modulo \(10^9+7\), division is handled using the modular multiplicative inverse.

That is, \(P_{i,j}\) is computed in the program as:

\[ R_{i,j} \times \left(\sum_k R_{i,k}\right)^{-1} \pmod{10^9+7} \]

Since \(10^9+7\) is prime, the modular inverse can be found using Fermat’s Little Theorem as:

\[ x^{-1} \equiv x^{MOD-2} \pmod{MOD} \]

We also multiply by the modular inverse of \(2\) for the final division by \(2\).

\[ \frac{1}{2} \equiv 500000004 \pmod{10^9+7} \]

Complexity

• Time Complexity: \(O(N^3 \log T)\)
• Space Complexity: \(O(N^2)\)

Although we maintain two matrices (one for the total quantity and one for the difference quantity), this only affects the runtime by a constant factor, so the overall time complexity remains \(O(N^3 \log T)\).

Implementation Details

Be careful with the directions of the matrix indices.

The transition probability from container \(i\) to container \(j\) is stored in the matrix as:

\[ A_{j,i} \]

This is because we update the column vector \(v\) as:

\[ v' = A v \]

In the code, this is implemented as follows:

total_mat[j][i] = val

In the matrix for the difference quantity, we invert the sign only when the inversion rule applies.

if flip[i][j]:
    diff_mat[j][i] = MOD - val
else:
    diff_mat[j][i] = val

A negative value \(-val\) is represented as MOD - val under modulo arithmetic.

Finally, we restore the positive particle quantity from the total and difference quantities as follows:

ans = (total_vec[G] + diff_vec[G]) * INV2 % MOD

Source Code

import sys

MOD = 1000000007
INV2 = (MOD + 1) // 2

def mat_mul(A, B, n):
    res = []
    mod = MOD
    for i in range(n):
        row_res = [0] * n
        Ai = A[i]
        for k in range(n):
            a = Ai[k]
            if a:
                Bk = B[k]
                for j in range(n):
                    row_res[j] += a * Bk[j]
        for j in range(n):
            row_res[j] %= mod
        res.append(row_res)
    return res

def mat_vec_mul(A, v, n):
    res = [0] * n
    mod = MOD
    for i in range(n):
        s = 0
        Ai = A[i]
        for j in range(n):
            s += Ai[j] * v[j]
        res[i] = s % mod
    return res

def main():
    data = list(map(int, sys.stdin.buffer.read().split()))
    idx = 0

    N = data[idx]
    M = data[idx + 1]
    T = data[idx + 2]
    idx += 3

    S = data[idx] - 1
    G = data[idx + 1] - 1
    idx += 2

    R = []
    row_sums = []
    for _ in range(N):
        row = data[idx:idx + N]
        idx += N
        R.append(row)
        row_sums.append(sum(row))

    flip = [[False] * N for _ in range(N)]
    for _ in range(M):
        u = data[idx] - 1
        v = data[idx + 1] - 1
        idx += 2
        flip[u][v] = True

    total_mat = [[0] * N for _ in range(N)]
    diff_mat = [[0] * N for _ in range(N)]

    for i in range(N):
        inv_sum = pow(row_sums[i], MOD - 2, MOD)
        Ri = R[i]
        for j in range(N):
            val = Ri[j] * inv_sum % MOD
            total_mat[j][i] = val
            if flip[i][j]:
                diff_mat[j][i] = MOD - val if val else 0
            else:
                diff_mat[j][i] = val

    total_vec = [0] * N
    diff_vec = [0] * N
    total_vec[S] = 1
    diff_vec[S] = 1

    t = T
    while t:
        if t & 1:
            total_vec = mat_vec_mul(total_mat, total_vec, N)
            diff_vec = mat_vec_mul(diff_mat, diff_vec, N)
        t >>= 1
        if t:
            total_mat = mat_mul(total_mat, total_mat, N)
            diff_mat = mat_mul(diff_mat, diff_mat, N)

    ans = (total_vec[G] + diff_vec[G]) * INV2 % MOD
    print(ans)

if __name__ == "__main__":
    main()

---

This editorial was generated by gpt-5.5-high.