Overview

This problem asks us to find the amount of positive particles in a specific container after \(T\) steps in a system where particles transition between \(N\) containers, and the charge (positive/negative) is inverted along certain paths. Since \(T\) is extremely large, we represent the state transitions as a matrix and use matrix exponentiation to perform the calculation efficiently.

Analysis

First, let’s consider a naive simulation of the \(T\) transitions. In each step, for each of the \(N\) containers, we must sum the transitions from all other containers. This requires \(O(N^2)\) operations per step. Since \(T\) can be up to \(10^{18}\), an \(O(N^2 T)\) approach will not run within the time limit.

For problems involving “repeated state transitions,” we can represent the transitions as matrix multiplications. By using matrix exponentiation, we can reduce the time complexity to \(O(\text{(number of states)}^3 \log T)\).

Since each container has two states—the “amount of positive particles” and the “amount of negative particles”—the total number of states is \(2N\). Let the state vector \(\vec{x}\) be a vector of length \(2N\), where \(\vec{x}_{2i}\) represents the amount of positive particles in container \(i\), and \(\vec{x}_{2i+1}\) represents the amount of negative particles in container \(i\).

A single-step transition can be represented using a \(2N \times 2N\) transition matrix \(A\) as \(\vec{x}_{new} = A \vec{x}_{old}\). Let the transition probability (or ratio) from container \(i\) to container \(j\) be \(P_{i,j} = \frac{R_{i,j}}{\sum_{j'} R_{i,j'}}\). - If there is no charge inversion: positive particles in container \(i\) transition to positive particles in container \(j\), and negative particles transition to negative particles. That is, \(A_{2j, 2i} = P_{i,j}\) and \(A_{2j+1, 2i+1} = P_{i,j}\). - If there is charge inversion: positive particles transition to negative particles, and negative particles transition to positive particles. That is, \(A_{2j, 2i+1} = P_{i,j}\) and \(A_{2j+1, 2i} = P_{i,j}\).

In the initial state, there is \(1\) unit of positive particles in container \(S\), so \(\vec{x}^{(0)}_{2S} = 1\) and all other elements are \(0\). The state after \(T\) steps is given by \(\vec{x}^{(T)} = A^T \vec{x}^{(0)}\). The value we want to find is the amount of positive particles in container \(G\), which is \(\vec{x}^{(T)}_{2G}\). This is equal to the \((2G, 2S)\) entry of the exponentiated transition matrix \(A^T\).

Algorithm

1. Calculate the transition ratio \(P_{i,j}\) for each transition. As specified in the problem statement, we treat the answer as a rational number and output it modulo \(10^9+7\). Thus, division is performed by multiplying by the modular inverse (modinv) using Fermat’s Little Theorem.
2. Construct the \(2N \times 2N\) transition matrix \(A\). Following the analysis above, place \(P_{i,j}\) in the appropriate positions depending on whether the inversion rule is applied to \((i, j)\).
3. Compute \(A^T\) using the matrix exponentiation algorithm (binary exponentiation). All calculations should be performed modulo \(10^9+7\).
4. The \((2G, 2S)\) entry of the resulting matrix \(A^T\) is the answer.

Complexity

• Time Complexity: \(O(N^3 \log T)\)
  • Since the matrix size is \(2N\), each matrix multiplication takes \(O((2N)^3)\) time, and binary exponentiation performs \(\log T\) multiplications.
• Space Complexity: \(O(N^2)\)
  • Since we store a few matrices of size \(2N \times 2N\).

Key Implementation Points

• Pay close attention to the indices of the transition matrix. If we define the state vector as \(\vec{x}_{new} = A \vec{x}_{old}\), the row index of matrix \(A\) corresponds to the destination (\(j\)), and the column index corresponds to the source (\(i\)). In the code, row1 and row2 represent the destination, while col1 and col2 represent the source.

• By adding if (A.a[i][k] == 0) continue; in the inner loop of the matrix multiplication, we can skip some computations since the transition matrix is sparse. This improves the constant factor.

  Source Code

#include <iostream>
#include <vector>

using namespace std;

const int MOD = 1e9 + 7;

struct Matrix {
    int n;
    vector<vector<long long>> a;
    Matrix(int n) : n(n), a(n, vector<long long>(n, 0)) {}
};

Matrix multiply(const Matrix &A, const Matrix &B) {
    int n = A.n;
    Matrix C(n);
    for (int i = 0; i < n; ++i) {
        for (int k = 0; k < n; ++k) {
            if (A.a[i][k] == 0) continue;
            for (int j = 0; j < n; ++j) {
                C.a[i][j] = (C.a[i][j] + A.a[i][k] * B.a[k][j]) % MOD;
            }
        }
    }
    return C;
}

Matrix mat_pow(Matrix A, long long p) {
    int n = A.n;
    Matrix B(n);
    for (int i = 0; i < n; ++i) B.a[i][i] = 1;
    while (p > 0) {
        if (p & 1) B = multiply(B, A);
        A = multiply(A, A);
        p >>= 1;
    }
    return B;
}

long long modinv(long long a) {
    long long b = MOD, u = 1, v = 0;
    while (b) {
        long long t = a / b;
        a -= t * b; swap(a, b);
        u -= t * v; swap(u, v);
    }
    u %= MOD;
    if (u < 0) u += MOD;
    return u;
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);
    int N, M;
    long long T;
    if (!(cin >> N >> M >> T)) return 0;
    int S, G;
    cin >> S >> G;
    --S; --G;

    vector<vector<long long>> R(N, vector<long long>(N));
    vector<long long> sumR(N, 0);
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; ++j) {
            cin >> R[i][j];
            sumR[i] += R[i][j];
        }
    }

    vector<vector<bool>> F(N, vector<bool>(N, false));
    for (int k = 0; k < M; ++k) {
        int u, v;
        cin >> u >> v;
        --u; --v;
        F[u][v] = true;
    }

    int sz = 2 * N;
    Matrix A(sz);
    for (int i = 0; i < N; ++i) {
        long long inv_sum = modinv(sumR[i]);
        for (int j = 0; j < N; ++j) {
            long long P = R[i][j] * inv_sum % MOD;
            int row1 = 2 * j;
            int row2 = 2 * j + 1;
            int col1 = 2 * i;
            int col2 = 2 * i + 1;
            if (F[i][j]) {
                A.a[row1][col1] = 0;
                A.a[row1][col2] = P;
                A.a[row2][col1] = P;
                A.a[row2][col2] = 0;
            } else {
                A.a[row1][col1] = P;
                A.a[row1][col2] = 0;
                A.a[row2][col1] = 0;
                A.a[row2][col2] = P;
            }
        }
    }

    Matrix B = mat_pow(A, T);

    // x^(0) has 1 at 2*S, 0 elsewhere.
    // So the answer is B.a[2*G][2*S]
    cout << B.a[2 * G][2 * S] << endl;

    return 0;
}

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