Overview

We consider the combination of the particle’s “container” and its “charge sign (positive/negative)” as \(2N\) states, and represent a 1-step update using a matrix.
Then, by raising this transition matrix to the power of \(T\), we can quickly find the state after \(T\) steps.

Analysis

Each container has positive and negative particles.
Therefore, instead of just looking at the container numbers, we separate the states as follows:

• State \(i\): Positive particle in container \(i\)
• State \(N+i\): Negative particle in container \(i\)

In other words, there are \(2N\) states in total.

For example, when a particle moves from container \(i\) to container \(j\):

• If there is no flip rule:
  • A positive particle remains positive.
  • A negative particle remains negative.
• If there is a flip rule:
  • A positive particle becomes negative.
  • A negative particle becomes positive.

These transitions are all linear.
That is, if we let \(v\) be the current state vector and \(A\) be the matrix representing a 1-step transition, the state after 1 step can be represented as:

\( vA \)

Therefore, the state after \(T\) steps is:

\( vA^T \)

---

A naive simulation of \(T\) steps will not run in time because \(T\) can be up to \(10^{18}\).
Also, if we manage the particle amount directly as rational numbers, the numerators and denominators will become extremely large.

Therefore, we perform all calculations modulo \(10^9+7\).
The probability

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

can be calculated using the modular inverse of the denominator as:

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

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

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

Algorithm

First, we construct a transition matrix \(A\) of size \(2N \times 2N\).

We map the indices as follows:

• Positive particle in container \(i\): \(i\)
• Negative particle in container \(i\): \(N+i\)

Note that in the implementation, we use 0-based indexing, so container \(i\) is \(0 \leq i < N\).

---

Let the transition probability from container \(i\) to container \(j\) be:

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

We calculate this value modulo \(10^9+7\).

If there is no flip rule

Positive moves to positive, and negative moves to negative.

That is:

• Positive in container \(i\) \(\rightarrow\) Positive in container \(j\)
• Negative in container \(i\) \(\rightarrow\) Negative in container \(j\)

Therefore, we add to the matrix as follows:

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

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

If there is a flip rule

Positive moves to negative, and negative moves to positive.

That is:

• Positive in container \(i\) \(\rightarrow\) Negative in container \(j\)
• Negative in container \(i\) \(\rightarrow\) Positive in container \(j\)

Therefore, we add to the matrix as follows:

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

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

---

In the initial state, there is exactly \(1\) positive particle in container \(S\), so the state vector \(v\) has:

\( v[S] = 1 \)

and \(0\) for all other elements.

Then, we compute \(A^T\) using binary exponentiation, multiplying it with the vector only when necessary.

Specifically, we examine each bit of \(T\):

• If the bit is set, we update \(v = vA\).
• We update \(A = A^2\) at each step.

Finally, the amount of positive particles in container \(G\) will be \(v[G]\).

Complexity

Let the number of states be \(K = 2N\).

• Matrix-matrix multiplication: \(O(K^3)\)
• Number of matrix exponentiation steps: \(O(\log T)\)
• Vector-matrix multiplication: \(O(K^2)\)

Therefore, the overall complexity is:

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

Since \(N \leq 30\), the maximum matrix size is \(60 \times 60\), which is small enough to be computed very quickly.

Implementation Details

• The transition matrix is constructed with “current state” as rows and “next state” as columns.

  • Therefore, the state vector is treated as a row vector, and the update is \(v = vA\).

• The modular inverse of the row sum (the denominator of the probability) is computed using binary exponentiation.

• Since \(T\) can be up to \(10^{18}\), we use a 64-bit integer type (such as long long in C++).

• Since the flip rules are directed, we only set flip[u][v] to true. The reverse direction is treated separately.

  Source Code

#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;
}

---

This editorial was generated by gpt-5.5-high.