Official

E - 粒子の拡散シミュレーション / Particle Diffusion Simulation Editorial by admin Editorial for printing

gemini-3.5-flash-high

Overview

This problem asks us to simulate the diffusion of particles moving between connected containers and find the “amount of positively charged particles” in a specific container after a given number of steps \(T\). Since the number of steps \(T\) can be as large as \(10^{18}\), instead of a naive simulation, we compute this quickly using matrix exponentiation (doubling).

Analysis

1. State Representation

In each container, there are two states: “positively charged particles” and “negatively charged particles”. Since there are \(N\) containers, the total number of states in the system can be considered as \(2N\).

We represent the state as a \(2N\)-dimensional column vector \(X\) as follows: * Elements \(0\) to \(N-1\): The amount of positive particles in containers \(1\) to \(N\). * Elements \(N\) to \(2N-1\): The amount of negative particles in containers \(1\) to \(N\).

In the initial state (step \(0\)), there is exactly \(1\) positive particle in container \(S\), and all other values are \(0\). Expressing this in 0-indexed terms, the \((S-1)\)-th element of the vector \(X_0\) is \(1\), and all other elements are \(0\).

2. Formulating Transitions (Matrix Representation)

A single-step transition can be expressed as a linear transformation. The transition probability \(P_{i,j}\) from container \(i\) to container \(j\) can be found by calculating the modular inverse of \(\sum_{j'} R_{i,j'}\) modulo \(10^9+7\) as follows: $\(P_{i,j} = R_{i,j} \times \left( \sum_{j'=1}^{N} R_{i,j'} \right)^{-1} \pmod{10^9+7}\)$

Using these transition probabilities and the inversion rules \(F_{i,j}\) (\(1\) if inverted, \(0\) otherwise), we construct a \(2N \times 2N\) transition matrix \(A\) representing the transition from step \(t\) to \(t+1\). Each entry \(A_{row, col}\) of the matrix \(A\) represents the “transition probability from state \(col\) to state \(row\)”.

  • When there is no inversion rule (\(F_{i,j} = 0\)): The charge does not change, so positive transitions to positive, and negative transitions to negative.

    • \(A_{j, i} = P_{i,j}\) (positive \(\to\) positive)
    • \(A_{N+j, N+i} = P_{i,j}\) (negative \(\to\) negative)
  • When there is an inversion rule (\(F_{i,j} = 1\)): The charge is inverted, so positive transitions to negative, and negative transitions to positive.

    • \(A_{N+j, i} = P_{i,j}\) (positive \(\to\) negative)
    • \(A_{j, N+i} = P_{i,j}\) (negative \(\to\) positive)

3. Acceleration via Matrix Exponentiation

Letting \(X_t\) be the state vector at step \(t\), a single-step transition is represented as \(X_{t+1} = A X_t\). Therefore, the state vector \(X_T\) after \(T\) steps can be expressed using matrix multiplication as follows: $\(X_T = A^T X_0\)$

Since only the \((S-1)\)-th element of the initial state \(X_0\) is \(1\), the result of \(A^T X_0\) is exactly the \((S-1)\)-th column of the matrix \(A^T\). The final value we want to find is the amount of positive particles in container \(G\) (the \((G-1)\)-th element of \(X_T\)), which corresponds to the \((G-1, S-1)\) entry of the matrix \(A^T\), i.e., A_T[G-1][S-1].

Although \(T \le 10^{18}\), by using binary exponentiation (doubling), we can compute \(A^T\) with a time complexity of \(O(N^3 \log T)\).

Algorithm

  1. Calculate Modular Inverses: For each row \(i\), calculate the row sum \(S_i = \sum_{j} R_{i,j}\), and find the modular inverse \(S_i^{-1} \equiv S_i^{MOD-2} \pmod{MOD}\) using Fermat’s Little Theorem. Use this to calculate each transition probability \(P_{i,j}\).
  2. Construct the Transition Matrix \(A\): Prepare a square matrix \(A\) of size \(2N \times 2N\), and place the probabilities \(P_{i,j}\) appropriately depending on the presence of the inversion rules.
  3. Compute Matrix Exponentiation: Use binary exponentiation to compute \(A^T\), which is the transition matrix \(A\) raised to the power of \(T\).
  4. Output the Answer: Output the element at row \(G-1\), column \(S-1\) of \(A^T\).

Complexity

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

    • The size of the matrix is \(2N \times 2N\).
    • Computing the product of two matrices takes \(O((2N)^3) = O(N^3)\) time.
    • By using binary exponentiation, we perform matrix multiplication approximately \(\log_2 T\) times.
    • Since \(N \le 30\) and \(\log_2 T \le 60\), the total number of operations is approximately \(60^3 \times 60 \approx 1.3 \times 10^7\), which is well within the time limit.
  • Space Complexity: \(O(N^2)\)

    • Since we only store a few matrices of size \(2N \times 2N\), the memory usage is extremely small.

Implementation Points

  • Optimizing Matrix Multiplication: In interpreted languages like Python, nested loops can be slow. You can significantly reduce the execution time by adding a pruning optimization during matrix multiplication: “skip the innermost loop if the value from the source state is \(0\)”.

  • Converting 1-indexed to 0-indexed: Since the input container indices \(S, G\) and the inversion rules \(u_k, v_k\) are 1-indexed, we should subtract \(1\) from them beforehand to convert them to 0-indexed for array indexing.

    Source Code

import sys

def solve():
    input = sys.stdin.read
    data = input().split()
    if not data:
        return
    
    N = int(data[0])
    M = int(data[1])
    T = int(data[2])
    S = int(data[3])
    G = int(data[4])
    
    idx = 5
    R = []
    for i in range(N):
        row = []
        for j in range(N):
            row.append(int(data[idx]))
            idx += 1
        R.append(row)
        
    F = [[0] * N for _ in range(N)]
    for _ in range(M):
        u = int(data[idx]) - 1
        v = int(data[idx+1]) - 1
        idx += 2
        F[u][v] = 1
        
    MOD = 10**9 + 7
    N2 = 2 * N
    
    A = [[0] * N2 for _ in range(N2)]
    for i in range(N):
        sum_ri = sum(R[i])
        inv_sum = pow(sum_ri, MOD - 2, MOD)
        for j in range(N):
            p = (R[i][j] * inv_sum) % MOD
            if p == 0:
                continue
            if F[i][j]:
                A[j][N + i] = p
                A[N + j][i] = p
            else:
                A[j][i] = p
                A[N + j][N + i] = p
                
    def mat_mul(X, Y, N2=N2, MOD=MOD):
        Z = [[0] * N2 for _ in range(N2)]
        for i in range(N2):
            row_X = X[i]
            row_Z = Z[i]
            for k in range(N2):
                x_ik = row_X[k]
                if x_ik == 0:
                    continue
                row_Y = Y[k]
                for j in range(N2):
                    row_Z[j] = (row_Z[j] + x_ik * row_Y[j]) % MOD
        return Z

    def mat_pow(X, p):
        res = [[0] * N2 for _ in range(N2)]
        for i in range(N2):
            res[i][i] = 1
        base = X
        while p > 0:
            if p & 1:
                res = mat_mul(res, base)
            base = mat_mul(base, base)
            p >>= 1
        return res

    A_T = mat_pow(A, T)
    ans = A_T[G - 1][S - 1]
    print(ans)

if __name__ == '__main__':
    solve()

This editorial was generated by gemini-3.5-flash-high.

posted:
last update: