E - 粒子の拡散シミュレーション / Particle Diffusion Simulation 解説 by admin
gemini-3.5-flash-high概要
本問題は、接続された \(N\) 個の容器の間で粒子が拡散する様子をシミュレーションし、\(T\) ステップ後に特定の容器 \(G\) に存在する正の電荷を持つ粒子の量を求める問題です。ステップ数 \(T\) が最大 \(10^{18}\) と非常に大きいため、愚直なシミュレーションではなく、行列累乗（繰り返し二乗法）を用いた高速化を行います。
考察
1. 愚直なシミュレーションの限界
各ステップにおける粒子の遷移は、現在の各容器の粒子量（正・負）の線形結合（足し算と定数倍）で表されます。 しかし、ステップ数 \(T\) が最大 \(10^{18}\) であるため、1ステップずつ順番に状態を更新していく愚直なシミュレーションでは実行時間制限に間に合いません（TLEとなります）。
2. 線形性と状態の拡張
各容器 \(i\) について、正の粒子の量 \(x_i^{+}\) と負の粒子の量 \(x_i^{-}\) の \(2\) つの状態を管理する必要があります。容器の数が \(N\) 個なので、全体で \(2N\) 個の状態が存在することになります。
これらの状態を一つの列ベクトル \(\mathbf{x}\) としてまとめます。 $\(\mathbf{x} = \begin{pmatrix} x_1^{+} \\ \vdots \\ x_N^{+} \\ x_1^{-} \\ \vdots \\ x_N^{-} \end{pmatrix}\)$
1ステップの遷移は、このベクトル \(\mathbf{x}\) に対する線形写像、すなわち \((2N) \times (2N)\) の遷移行列 \(A\) を用いた行列とベクトルの積 \(\mathbf{y} = A \mathbf{x}\) として表現できます。
3. 行列累乗による高速化
1ステップの遷移が \(\mathbf{x}^{(1)} = A \mathbf{x}^{(0)}\) と表せるならば、\(T\) ステップ後の状態は次のようになります。 $\(\mathbf{x}^{(T)} = A^T \mathbf{x}^{(0)}\)$
行列の \(T\) 乗（\(A^T\)）は、繰り返し二乗法（バイナリ法）を用いることで \(O((2N)^3 \log T)\) の計算量で求めることができます。これにより、\(T = 10^{18}\) の場合でも高速に計算が可能です。
アルゴリズム
1. 遷移確率の計算とモジュラ逆数
容器 \(i\) から容器 \(j\) への遷移確率 \(P_{i,j}\) は次のように定義されます。 $\(P_{i,j} = \frac{R_{i,j}}{\sum_{j'=1}^{N} R_{i,j'}}\)$
有理数を \(\pmod{10^9+7}\) で計算するため、分母の \(D_i = \sum_{j'=1}^{N} R_{i,j'}\) について、フェルマーの小定理を用いて逆元 \(D_i^{-1} \pmod{10^9+7}\) を計算します。これにより、遷移確率は次のように掛け算で表せます。 $\(P_{i,j} \equiv R_{i,j} \times D_i^{-1} \pmod{10^9+7}\)$
2. 遷移行列 \(A\) の構築
サイズ \((2N) \times (2N)\) の遷移行列 \(A\) を構築します。 インデックス \(0\) から \(N-1\) を「正の電荷」、\(N\) から \(2N-1\) を「負の電荷」に対応させます。
容器 \(i\) から容器 \(j\) への遷移について： * 反転ルールが存在しない場合 (\(F_{i,j} = 0\)): * 正の粒子は正の粒子のまま遷移します： \(A[j][i] \leftarrow A[j][i] + P_{i,j}\) * 負の粒子は負の粒子のまま遷移します： \(A[j+N][i+N] \leftarrow A[j+N][i+N] + P_{i,j}\) * 反転ルールが存在する場合 (\(F_{i,j} = 1\)): * 正の粒子は負の粒子に反転します： \(A[j+N][i] \leftarrow A[j+N][i] + P_{i,j}\) * 負の粒子は正の粒子に反転します： \(A[j][i+N] \leftarrow A[j][i+N] + P_{i,j}\)
3. 行列累乗と回答の出力
初期状態 \(\mathbf{x}^{(0)}\) は、容器 \(S\) に正の粒子が \(1\) だけ存在し、他はすべて \(0\) です。 したがって、求める \(T\) ステップ後の容器 \(G\) の正の粒子の量は、累積遷移行列 \(A^T\) の \((G, S)\) 成分、すなわち \(A^T[G][S]\)（0-indexed では \(A^T[G-1][S-1]\)）に一致します。
計算量
時間計算量: \(O(N^3 \log T)\)
- 遷移行列のサイズは \(2N \times 2N\) です。
- 行列の乗算 1 回につき \(O((2N)^3) = O(N^3)\) の計算量がかかります。
- 繰り返し二乗法による累乗計算で、乗算を \(O(\log T)\) 回行うため、全体の計算量は \(O(N^3 \log T)\) となります。
- \(N \le 30\)、 \(T \le 10^{18}\) のとき、演算回数はおおよそ \(60^3 \times 60 \approx 1.3 \times 10^7\) 回となり、実行時間制限に余裕で間に合います。
空間計算量: \(O(N^2)\)
- \(2N \times 2N\) の行列を保持するために必要なメモリは \(O(N^2)\) です。\(N \le 30\) のため、メモリ使用量は極めてわずかです。
実装のポイント
0-indexed への変換: 入力される容器の番号 \(S, G\) や反転ルールの \(u_k, v_k\) は 1-based なので、実装内では \(1\) を引いて 0-based に変換して処理します。
行列積の高速化（枝刈り）: 行列積 \(C_{i,j} = \sum_k A_{i,k} B_{k,j}\) を計算する際、内側のループで \(A_{i,k} == 0\) の場合に処理をスキップ（
continue）することで、定数倍の高速化を図っています。モジュラ逆数: \(10^9+7\) が素数であるため、
modInverseは繰り返し二乗法を用いて \(O(\log(\text{MOD}))\) で計算します。ソースコード
#include <iostream>
#include <vector>
using namespace std;
const int MOD = 1000000007;
// Function to compute (base^exp) % MOD
long long power(long long base, long long exp) {
long long res = 1;
base %= MOD;
while (exp > 0) {
if (exp % 2 == 1) res = (res * base) % MOD;
base = (base * base) % MOD;
exp /= 2;
}
return res;
}
// Function to compute modular inverse
long long modInverse(long long n) {
return power(n, MOD - 2);
}
typedef vector<vector<long long>> Matrix;
// Matrix multiplication under modulo MOD
Matrix multiply(const Matrix& A, const Matrix& B) {
int n = 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;
for (int j = 0; j < n; ++j) {
C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD;
}
}
}
return C;
}
// Binary exponentiation for Matrix
Matrix matrixPower(Matrix A, long long p) {
int n = A.size();
Matrix res(n, vector<long long>(n, 0));
for (int i = 0; i < n; ++i) res[i][i] = 1;
while (p > 0) {
if (p & 1) res = multiply(res, A);
A = multiply(A, A);
p >>= 1;
}
return res;
}
int main() {
// Optimize standard I/O operations for performance
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; // Convert to 0-based indexing
vector<vector<long long>> R(N, vector<long long>(N));
vector<long long> D(N, 0);
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
cin >> R[i][j];
D[i] = (D[i] + R[i][j]) % MOD;
}
}
vector<vector<int>> F(N, vector<int>(N, 0));
for (int k = 0; k < M; ++k) {
int u, v;
cin >> u >> v;
--u; --v; // Convert to 0-based indexing
F[u][v] = 1;
}
vector<long long> D_inv(N);
for (int i = 0; i < N; ++i) {
D_inv[i] = modInverse(D[i]);
}
int sz = 2 * N;
Matrix A(sz, vector<long long>(sz, 0));
// Construct the transition matrix A
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
long long P = (R[i][j] * D_inv[i]) % MOD;
if (P == 0) continue;
int f = F[i][j];
if (f == 0) {
// No inversion: positive to positive, negative to negative
A[j][i] = (A[j][i] + P) % MOD;
A[j + N][i + N] = (A[j + N][i + N] + P) % MOD;
} else {
// Inversion: positive to negative, negative to positive
A[j][i + N] = (A[j][i + N] + P) % MOD;
A[j + N][i] = (A[j + N][i] + P) % MOD;
}
}
}
// Compute A^T
Matrix AT = matrixPower(A, T);
// The answer is the positive charge in container G, starting from container S
long long ans = AT[G][S];
cout << ans << "\n";
return 0;
}
この解説は gemini-3.5-flash-high によって生成されました。
投稿日時:
最終更新: