Eternal Tours(yukicoder No.1241)を線形代数で殴る

昨日の yukicoder コンテストで面白い問題が出題されました.AtCoder のコンテストで出題予定だった問題と被ったというハプニングもあり,話題になりました.
No.1241 Eternal Tours - yukicoder

想定解法は「鏡像法」と「二次元フーリエ変換」を用いるものでしたが,この記事ではこれらを使わない別解を紹介します.

別解では,行列のクロネッカー和など,競プロ界隈ではおそらくあまり知られていない線形代数の知識を使います.

問題の定式化

問題に合わせて全て 1-indexed で書きます.$H := 2^X-1$,$W := 2^Y-1$ とおきます.

北から $i$ 番目,西から $j$ 番目の区画を区画 $(i, j)$ と呼びます.区画 $(i, j)$ での行動後に区画 $(k, l)$ にいることが可能なとき $m_{(i, j), (k, l)} = 1$ と定め,そうでないとき $m_{(i, j), (k, l)} = 0$ と定めます.

$m_{(i, j), (k, l)}$ を各成分に持つ $HW \times HW$ 行列を $M$ とすると,問題の答えは $M^T$ の第 $((a, b), (c, d))$ 成分に等しいです.以下では行列 $M$ の固有値・固有ベクトルを計算することで,この値を求めます.

($M^T$ は $M$ の転置ではなく $M$ の $T$ 乗です.転置は $M^\top$ と書くのが好みです.)

行列 $M$ の構造

まず行列 $M$ の構造を調べましょう.

正の整数 $n$ に対し,$n \times n$ 単位行列を $I_n$ と書き,$n \times n$ 行列 $C_n$ を\[C_n :=\begin{bmatrix}\quad& 1 &\\[0.2em]1 &\quad& 1 & \\[-0.3em]& 1 &\quad& \ddots\\[-0.3em]& & \ddots&\quad& 1 \\[0.3em]& & & 1 & \quad\\\end{bmatrix}\]と定めます.つまり,行列 $C_n$ の $(i, j)$ 成分は $|i-j| = 1$ のとき $1$ で,そうでないとき $0$ です.

これを用いると,行列 $M$ は\begin{equation}\label{M_block}M =\red{\begin{bmatrix}C_W & I_W \\[-0.3em]I_W & C_W & \ddots \\[-0.3em]& \ddots & \ddots & I_W\\[0.3em]& & I_W & C_W\end{bmatrix}+ I_{HW}}\end{equation}と書けます.右辺第 $1$ 項の行列は,$W \times W$ マスの小さなブロックが $H \times H$ 個並んでいるという感じです.

これで行列 $M$ の構造が分かりやすくなりました!

勘の良い方はお気づきかもしれませんが,行列 $C_n$ は頂点数 $n$ のパスグラフ隣接行列です.また,式 \eqref{M_block} の右辺第 $1$ 項は問題で与えられた $H \times W$ のグリッドグラフの隣接行列で,第 $2$ 項 $I_{HW}$ は「その場に留まる」ことに対応しています.

クロネッカー積

式 \eqref{M_block} の右辺第 $1$ 項には行列 $C_W, I_W$ が繰り返し現れているため,クロネッカー積(テンソル積)が使えそうです.クロネッカー積は,以下のように定義されます:

$m \times n$ 行列 $A = (a_{i,j})$ と $p \times q$ 行列 $B$ に対し,$A$ と $B$ のクロネッカー積 $\red{A \otimes B}$ は\[A \otimes B =\begin{bmatrix}a_{1,1}B & a_{1,2}B & \cdots & a_{1,n}B \\a_{2,1}B & a_{2,2}B & \cdots & a_{2,n}B \\\vdots & \vdots & & \vdots\\a_{m,1}B & a_{m,2}B & \cdots & a_{m,n}B\end{bmatrix}\]で定義される.これは $mp \times nq$ 行列である.

式 \eqref{M_block} から,行列 $M$ はクロネッカー積を用いて\begin{equation}\label{M-kron-sum}M = \red{C_H \otimes I_W + I_H \otimes C_W + I_{HW}}\end{equation}と書けることが分かります.かなりシンプルになりました!

クロネッカー和の固有値

式 \eqref{M-kron-sum} の $C_H \otimes I_W + I_H \otimes C_W$ の部分に注目してください.実はこれにもクロネッカー和という名前がついていて,固有値・固有ベクトルについて嬉しい性質が知られています:

$m \times m$ 行列 $A$ と $n \times n$ 行列 $B$ に対し,$A$ と $B$ のクロネッカー和 $\red{A \oplus B}$ は\[A \oplus B = \red{A \otimes I_n + I_m \otimes B}\]で定義される.これは $mn \times mn$ 行列である.

$A$ の固有値・固有ベクトルの組を $(\lambda_i, \boldsymbol{u}_i)$ $(1 \leq i \leq m)$ とし,$B$ の固有値・固有ベクトルの組を $(\mu_j, \boldsymbol{v}_j)$ $(1 \leq j \leq n)$ とすると,クロネッカー和 $A \oplus B$ の固有値・固有ベクトルの組は\[\red{(\lambda_i + \mu_j, \ \boldsymbol{u}_i \otimes \boldsymbol{v}_j)}\qquad(1 \leq i \leq m,\ 1 \leq j \leq n)\]と書ける.

$\boldsymbol{u}_i \otimes \boldsymbol{v}_j$ は,$\boldsymbol{u}_i$ を $m\times 1$ 行列と,$\boldsymbol{v}_j$ を $n \times 1$ 行列と見なしたときのクロネッカー積です.

上の性質から,行列 $C_H, C_W$ の固有値・固有ベクトルが分かったら,行列 $M = \red{C_H \oplus C_W + I_{HW}}$ の固有値・固有ベクトルも分かります!

$C_H$ の固有値・固有ベクトル

$C_H, C_W$ は同じ形をしているので,$C_H$ についてのみ考えます.

実は $C_H$ は結構有名な行列で,その固有値・固有ベクトルはきれいに書けることが知られています:

行列 $C_H = C_{2^X-1}$ の固有値は\[\lambda_i = \red{2 \cos\Big( \frac{i \pi}{2^X} \Big)} \qquad (1 \leq i < 2^X) \] である. \[ u_{i,j} = \red{\frac{1}{\sqrt{2^{X-1}}} \sin\Big( \frac{ij \pi}{2^X} \Big)} \qquad (1 \leq i,j < 2^X) \] と定めると,$\boldsymbol{u}_i = [u_{i,1},\dots,u_{i,2^X-1}]^\top$ は固有値 $\lambda_i$ に対応する固有ベクトルである. さらに,$\boldsymbol{u}_1,\dots,\boldsymbol{u}_{2^X-1}$ は正規直交基底をなす.

参考: 三重対角行列の特殊形の固有値は綺麗 | 高校数学の美しい物語

問題の答え

以上の議論から,行列 $M = C_H \oplus C_W + I_{HW}$ の固有値が\[\lambda_{(i, j)} = \red{2\cos\Big( \frac{i \pi}{2^X} \Big) + 2\cos\Big( \frac{j \pi}{2^Y} \Big) + 1}\]$(1\leq i < 2^X,\ 1 \leq j < 2^Y)$ であり,これに対応する固有ベクトルの第 $(k, l)$ 成分が \[ v_{(i,j),(k,l)} = \red{\frac{1}{\sqrt{2^{X+Y-2}}} \sin\Big( \frac{ik \pi}{2^X} \Big) \sin\Big( \frac{jl \pi}{2^Y} \Big)} \] と書けることが分かりました!

今求めた $M$ の固有ベクトルたちが正規直交基底であることにも注意すると,問題の答え,すなわち $M^T$ の第 $((a,b), (c,d))$ 成分は,\begin{equation}\label{ans}\redunder{\sum_{i=1}^{2^X-1} \sum_{j=1}^{2^Y-1} \lambda_{(i,j)}^T v_{(i,j), (a,b)} v_{(i,j), (c,d)}}\end{equation}となります.

答えを $\mathrm{mod}\ \ 998244353$ で聞かれているため,実際には $2\cos z = e^{\mathrm{i}z} + e^{-\mathrm{i}z}$ や $1$ の原始 $2^{X+1}$ 乗根などを使って $\mathrm{mod}\ \ 998244353$ で計算する必要があります.

式 \eqref{ans} を素朴に計算すると,時間計算量は $O(2^{X+Y} (X + Y + \log T))$ となり,十分高速です.

余談

実は,行列 $M$ のようなグリッドグラフ由来の行列が,「パスグラフ由来の行列のクロネッカー和」を用いて書けることは(数値線形代数界隈や代数的グラフ理論界隈では)有名です.このことは,例えば波動方程式を差分近似して数値的に解く際などにも用いられます.

AC コード

AtCoder Librarymodint を使っています.

#include <bits/stdc++.h>
#include <atcoder/modint>
using namespace std;
using namespace atcoder;
#define rep2(i, m, n) for (int i = (m); i < (n); ++i)

using mint = modint998244353;

int main() {
  long long X, Y, T, a, b, c, d; cin >> X >> Y >> T >> a >> b >> c >> d;

  mint r = 3; // 998244353 の原始根
  mint xi_x = r.pow(119*(1<<(22-X))); // 1 の 2^{X+1} 乗根
  mint xi_y = r.pow(119*(1<<(22-Y))); // 1 の 2^{Y+1} 乗根

  mint ans = 0;
  rep2(i, 1, 1<<X) rep2(j, 1, 1<<Y) {
    mint lam = 1; // 行列 M の固有値
    lam += xi_x.pow(i) + 1 / xi_x.pow(i);
    lam += xi_y.pow(j) + 1 / xi_y.pow(j);

    mint res = lam.pow(T);
    res *= xi_x.pow(i*a) - 1 / xi_x.pow(i*a);
    res *= xi_y.pow(j*b) - 1 / xi_y.pow(j*b);
    res *= xi_x.pow(i*c) - 1 / xi_x.pow(i*c);
    res *= xi_y.pow(j*d) - 1 / xi_y.pow(j*d);
    ans += res;
  }
  ans /= mint(2).pow(X+Y+2);

  cout << ans.val() << '\n';
}
タイトルとURLをコピーしました