DFTを行列の対角化の視点から捉える

はじめに

行列の \(n\) 乗の計算をする時,対角化すると,固有値を \(n\) 乗するだけになるので計算が楽になります (競プロだと繰り返し二乗法で十分速いのであまり使い所がないですが,手計算でやるときは対角化をすると圧倒的に楽になります).

ところで,数列の \(n\) 乗を計算する時,数列を DFT してから,各要素を \(n\) 乗し,それをIDFT すれば \(O(N\log N)\) で高速に求まります.

これって行列の対角化に似てませんか?実は DFT は行列の対角化として捉えることができるのです.

DFT のおさらい

以下の記事が非常によくまとまっているのでそちらを見るのがおすすめですが,一応こちらにもかんたんに定義や基本的な性質を書いておきます.

qiita.com


数列 \((a_n)_{n=0,\dots,N-1}\) の離散 Fourier 変換 (DFT) \((A_n)_{n=0,\dots,N-1}\) は,\(W_N = \exp\left(-i\frac{2\pi}{N} \right)\) (1の\(N\)乗根) とすると,次の式で定義されます.

$$
A_n = \sum_{k=0}^{N-1} a_k W_N^{k}
$$

以降もとの数列を小文字で,それの DFT を大文字で書きます.また数列の長さは全部 \(N\) とします.

また,通常の畳み込みの代わりに巡回畳み込みを考えることにします.これは,添字を \(\mathrm{mod}\,N\) で考えて,畳み込み後の数列の長さが変わらないようにするだけです.つまり,数列 \(a, b\) の巡回畳み込み \(c\) は次式で定義されます.

$$
c_n = \sum_{k=0}^{N-1} a_{n-k\,\mathrm{mod}\,N} b_k \quad \text{for } n = 0,\dots, N-1
$$

通常の畳み込みを計算したいときは,適当に数列を延長して0を詰めて巡回畳み込みを計算すれば良いです.

DFT の著しい性質は,巡回畳み込みが要素ごとの積に置き換わることです.つまり,\(C_n = A_n B_n \) が成り立ちます.実際,

$$
C_n = \sum_{k=0}^{N-1} c_k W_N^{nk} \\
= \sum_{k=0}^{N-1} \left( \sum_{l=0}^{N-1} a_{k-l} b_l \right) W_N^{nk} \\
= \sum_{l=0}^{N-1} \left( \sum_{k=0}^{N-1} a_{k-l}W_N^{n(k-l)} \right) b_l W_N^{nl} \\
= \sum_{l=0}^{N-1} \left( \sum_{k=0}^{N-1} a_k W_N^{n(k-l)} \right) b_l W_N^{nl} \\
= \left( \sum_{k=0}^{N-1} a_k W_N^{nk} \right) \left( \sum_{l=0}^{N-1} b_l W_N^{nl} \right) \\
= A_n B_n
$$

3行目の変形では,和を取る順番を入れ替えて,因子 \(W_N^{nk}\) を分解しました.4行目の変形では,\(k-l \rightarrow k\) と置き換えました (添字を \(\mathrm{mod}\,N\) で考えており,さらに \(W_N^N=1\) なのでこれができる)

要素ごとの積は \(O(N)\) で求まるので,DFT が高速に計算できれば数列の畳み込みが高速に計算できます.DFT を定義どおりに計算すると \(O(N^2)\) となり,愚直に畳み込みを計算するのと変わりませんが,DFT を分割統治で高速に \(O(N\log N)\) で計算する高速 Fourier 変換 (FFT) というアルゴリズムが知られており,それを用いると畳み込みが \(O(N\log N)\) で求まります.

数列を線形写像として表す

DFTを行列の対角化として捉えたいので,まず数列を行列で表すことから始めます.数列を\(N\) 次元ベクトルとして捉えます.数列 \(x \in \mathbb{R}^N \) に,数列 \(a \in \mathbb{R}^N \) との巡回畳込み \(ax \in \mathbb{R}^N\) を対応させる写像を \(f_a(x)\) とします. これは線形写像です.

\(f_a\) の表現行列 \(M_a\) は以下のように表せます.

$$
M_a = \begin{pmatrix}
a_0 & a_{N-1} & \dots & a_1 \\
a_1 & a_0 & \dots & a_2 \\
\vdots & \vdots & \ddots & \vdots \\
a_{N-1} & a_{N-2} & \dots & a_{0}
\end{pmatrix}
$$

\(M_{ab}=M_aM_b, \quad M_{a^n} = M_a^n\) です.これで,数列の畳み込みが行列の掛け算として表現できました.

表現行列の対角化

\(M_a\) を対角化してみましょう.よく見ると \(M_a\) は巡回行列になっているので,この構造を利用すると,\(M_a\) の固有ベクトルは \(\boldsymbol{v}_n = \left( W_N^{n0}, W_N^{n1}, \dots, W_N^{n(N-1)} \right)^\intercal \text{for } n=0,\dots, N-1 \) となり,\(v_n\) に対応する固有値は \(A_n = \sum_{k=0}^{N-1} a_k W_N^{nk}\) になっていることがわかります.

これはDFTになっていますね?

よって,
$$
W = \begin{pmatrix}
W_N^{00} & W_N^{10} & \dots & W_N^{(N-1)0} \\
W_N^{01} & W_N^{11} & \dots & W_N^{(N-1)1} \\
\vdots & \vdots & \ddots & \vdots \\
W_N^{0(N-1)} & W_N^{1(N-1)} & \dots & W_N^{(N-1)(N-1)}
\end{pmatrix}
$$
とすれば,

$$
M_a = W\mathrm{diag}\left( A_0, A_1, \dots, A_{N-1} \right) W^{-1}
$$

と対角化できるわけです.

畳み込みの高速計算

\(a^n\) の表現行列 \(M_{a^n}\) は \(M_a^n\) として表せました.\(M_a\) が対角化できたので,その \(n\) 乗は固有値を \(n\) 乗するだけで簡単に求まります.

$$
M_{a^n} = M_a^n = W\mathrm{diag}\left( A_0^n, A_1^n, \dots, A_{N-1}^n \right)W^{-1}
$$

これは,DFTして,各要素を \(n\) 乗して,IDFTするという操作に対応しています.これで,DFTを用いて数列の \(n\) 乗を計算するアルゴリズムが,行列を対角化して固有値を \(n\) 乗するというアルゴリズムの特別な場合であることがわかりました.

ところで,対角化に用いた行列 \(W\) は,数列 \(a\) に依存しない形になっています.つまり,\(M_x\) は \(x\) が何であろうと全部 \(W\) で対角化できるということです.これすごいですね.

\(M_a = W\mathrm{diag}\left( A_0, A_1, \dots, A_{N-1} \right)W^{-1},\,M_b = W\mathrm{diag}\left( B_0, B_1, \dots, B_{N-1} \right) W^{-1}\) ですから,

$$
M_{ab} = M_aM_b = W\mathrm{diag}\left( A_0B_0, A_1B_1, \dots, A_{N-1}B_{N-1} \right) W^{-1}
$$

と計算できます.これは,DFT して要素ごとの積をとってIDFTするという操作に対応しています.一般の行列では,同じ行列で対角化できるとは限りませんから,行列積も対角化すると高速に計算できるというのは巡回行列の間で成り立つ特別な性質です.

おわりに

DFTすると畳み込みが積になるというのは,式は追えるんだけどあまり直観的な理解ができないでいました.今回,行列のDFTは対角化に対応しているという解釈を知って,少しは見通しが良くなった感覚があります.まあ直観的な理解に至ったかと言われるとまだ微妙ですが.

Hadamard 変換もこんな感じで解釈できるのでしょうか.

ICPC アジア横浜大会 2021 参加記

国内予選の参加記はこちら

sotanishy.hatenablog.com

チームはsuzukaze_Aobayama.メンバーは相変わらずむげんさん,仮の人君,僕.

目標は,とりあえず全体の半分以上の20位,欲を言えば15位.国内予選での順位(6位)を考えるとかなり低めの目標と思われるかも知れないが,国内予選での異常な順位はたまたま早解きに大成功しただけのまぐれ当たりだと思っており,アジアでは強い人達にちゃんと完数で差をつけられてしまうだろうと思った.まあ,国内は予選通過するために必死だったのに対し,アジアでは何も失うものがないので,順位をあまり気にせずに楽しみたいとも思っていた.

3/6 模擬地区

僕はCIKを解いた.Kが面白かった.全体では9完16位で,悪くはないという印象.

3/12 TUPC

suzukaze_Aobayamaのメンバーは皆運営陣として参加した.うちの大学で有志コンをやるのは初めてなので運営の仕方とかわからず色々と大変だったが,うまく行ったようでよかった.チームはレートに関わらず3人まで許容することにしたので,ICPCのチーム練として参加してくれたチームが多かったように見えた.僕はBIKのwriterとCHの共同writerをしていた.楽しんでいただけましたか?

3/14 チーム練

10時に始めるつもりだったが僕が起床に失敗したため30分遅れてしまった.チームメイト,すまん.......本番は起きれるのだろうか.
2021 ICPC Southeastern Europe Regional Contest を走った.14問中8問解いた.まあ解くべきものは解いたかなあと.仮の人くんがずっと天才をしていた.

昼ごはんを食べるタイミングを逃し,お腹を鳴らしながら問題を解いていた.本番ではコンテスト中に適当なタイミングで昼を食べないといけない.

3/15 開会式,リハーサル

12時に青葉山に集合し,本番で使わせてもらう部屋の準備などをして,そのままそこで13時からの開会式に出た.開会式は,オープニング動画がかっこよかった.

www.youtube.com

東京大学の圧倒的強者感.

リハーサルコンテスト.ABCの3完.リハーサルなので適当に簡単なのが出ると思っていたが,去年のアジアの問題がそのまま流用されていたらしい.去年のアジアこんなに難しかったのか?2時間で3完しかできないってだいぶ厳しい.......

終わったあとは翌日のコンテスト中に食べる食料の買い出しに行った.ICPCからもらったお菓子はゼリー,チョコ,ジュース以外は食べてしまったので補給.本番中はカフェインと糖分を大量に摂取するつもりである.

夜はいくつかAtCoderで問題を解いたり,知識の復習をしたり,ゆるゆりを観たりして,早めに寝た.

3/16 本番

7時頃起床に成功.朝ごはんを食べて山を登り8時半頃に参加場所に到着.起床に失敗したチームはどれくらいいるのかな〜と思いながらZoomに入ると,めちゃくちゃ人がいてびっくりした.まあ,流石に本番くらいは起きるか.......なんかトラブルが生じたらしく,開始が10分遅れた.

9時10分開始.いつもどおりむげんさんがA,仮の人くんがB,僕が後ろのKから見る.Kはちょっと考えてみて,何もとっかかりがつかめなかったので,とりあえず飛ばす.

次にJを開く.Jは,二次元平面上の2点を選んだときに定まる十字型の領域を考えて,その領域に他のすべての点が含まれるような2点の組は何個あるか数えるという問題.どうせ1点固定して,それと組になりうる点を適当にデータ構造で探索するんだろうなあと思う.頑張ればできそうだがめんどくさそうなのでどうしようかなあと悩んでいたところで,Aを通してDをみていたむげんさんが,「これは君たちの問題ですよ」といっているのでちょっとDの問題概要を伝えてもらう.なるほど,僕が結構好きな見た目をしている.そこで,Jをむげんさんに任せて,僕が代わりにDを考えることにする.むげんさんはデータ構造やるだけみたいな問題に異常に強いので託したほうが速いと判断した(あとで振り返るとこれがかなり良かった).

仮の人くんがBに苦戦しているようだったが,順位表を見ると通しているところがあまりなく,ペナが量産されていたので,その旨を伝える.それで冷静になれたらしく嘘をちゃんと直してノーペナで通してくれる.ナイスプレー.国内予選のときも僕がCで苦しんでいたときに,「どこも通ってないのでゆっくりやっていいですよ」と声をかけてくれたのがかなり精神的に助かった記憶があるので,こういうのはチーム戦のいいところ.

Dを見る.数列が与えられるので,そこからなんか木を作って,その木の直径を最大化するという問題.木はす木好きなので,取り組んでみる.40分ほどじーっとサンプルを睨んでいると,少しずついい感じの性質が見えてくる.作れる木の形というのはかなり限られている.そこで,数列を2つに分けて,左右それぞれでできるだけ長いパスを作って真ん中でつなげるというアルゴリズムで解けそうだとわかる.その最長のパスの長さは,左から順に見ていってスタックを掘っていくといい感じに求まりそうだ.さっと実装してみるとサンプルが一発で合う.半信半疑で提出.AC!!!FAだったらしい.気持ちえ〜〜〜〜〜!!

ここまで開始1時間ほど.その時点で3完2位.国内予選と同じ流れだが,流石にアジアになれば早解きだけでは勝てなさそうなので,あまり気にしない.順位表を見ると,通っているのがGとJくらいで,Jはむげんさんが大丈夫そうだと言っているのでGを見る.

Gは,ラベル付き二分木の個数の数え上げ.ある頂点に子があるとき,少なくとも一方の子のラベルは自分よりも大きくなくてはいけない.また,それぞれの頂点について,それが持ちうる子の数の最大値と最小値が指定される.これらの条件をすべて満たす二分木の個数の数を求める.ラベルの大小に関する条件は,頂点を後ろから順番に決めていけばうまく扱えそう.制約が \(O(n^3)\)っぽいので,連結成分の数と,空いている場所の数をそれぞれ持たせてDPをするんだろうなと思う.遷移は,今見ている頂点をある木の根の親にするか,別の頂点の子にして自分は子を持たないか,で良さそう.実装してみて,サンプルがあったので提出.WAが返ってくる.なにか見落としている遷移のパターンが無いかと考えると,自分をある木の根の親にして,さらに別の頂点の子にするパターンがあり得ることに気づく.しかしこれはヤバそう.木の根の親にしたあと,自分の子孫にはつないではいけないので,今持たせている状態だけでは遷移できない.困った.Jを通したむげんさんがGの考察に合流し,一緒に考える.ここで,よくよく考えると,根にしてから子にすると考えるのではなく,子にしてから根にすると考えると,今持たせている状態だけでちゃんと遷移できることに気づく.実装して提出.AC!!!そこの順番入れ替えるだけで解けるようになるんだ〜と感動していた.

ここまでで2時間が経過している.ABDGJの5完.この時点でまだ2位.は?またこの流れですか.このまま2位だったらWorld Finalですか,とうっかり口が滑ってしまう.そういうことは,思っても言わないものですよ.

仮の人くんがFにずっと取り組んでいる.結構考察が進んでいるらしいので任せて,僕はほかを見ることにする.順位表を見ると,他に通されているのがCしかない.とりあえずむげんさんと一緒にCを見るが,めちゃくちゃ嫌な見た目をしている.文字列をデコードする手順と,デコード後の文字列が与えられるので,デコードしてその文字列になるような元の文字列を構築する問題.やばいのは,それを反転しても正しくデコードできなければいけない.さらにそのような文字列のうち最短かつ辞書順最小のものを求めなくてはいけない.まじでやりたくないが,他の問題で通されているものがない以上Cに取り組むのが一番現実的なので,考えてみる.制約的に区間DPかなあと思ってしばらくその方針で考えてみるが,全然わからない.1時間ほどそれで考えていたが,考えれば考えるほどヤバいケースが思いつくので,区間DPは捨てる.次に取った方針が,\(n^2\)状態のDPで,遷移をするときにAとLを全探索するというもので,愚直にやると\(O(n^4)\)で間に合わなさそうだが,頑張って高速化するのかなあと思う.しばらく考えていると天啓が降ってきて,AとLは高々9であることに気づく(それはそう).あまりにも単純だが重要な情報だったので思わず奇声を上げてしまう.なんだかそれで行けそうだったので,実装をしてみる.すでにCを2時間位考えていて,順位表は凍結されていた.バグらせまくったがなんとかサンプルがあったので出してみる.TLE.DPで文字列を陽に持っていたので,最適解の長さだけ持たせて具体的な解はあとで復元するというように変更してみる.これもバグらせまくったがなんとか実装が終わり,終了10分くらい前に投げてみる.WA.え?どうしてだろう,とコードを睨んでいたが全然わからない.とりあえず適当なケースで出力を見てみるか,と思って 01 を入力してみたら,もうそれで全然出力が違った.困った......残り5分でそれを直す方法も思いつかず,時間切れ.

Fは,仮の人くんとむげんさんがずっと取り組んでいて,かなりいいところまで行っていたらしいが,そちらもデバッグが間に合わず(実は,あと1行足したら通るという所まで来ていたらしい.すげー).結局のこり3時間何も解けずに5完で終了.凍結前が6位で,そこから少し下がることが見込まれたが,それでも10位以内には入れそうな雰囲気がある.

コンテスト後

しばらく残って余韻に浸ったり感想戦をしたりしたいところだが,閉会式が始まるまでの空き時間を逃すと帰宅するタイミングがなくなるので速攻帰る.家が遠いので,家についたときにはすでに解説が始まっていたが,とりあえずYes/Noには間に合って良かった.

結果は10位だった.国内予選に続き,想定していたよりずっといい順位を取れた.まさか top 10 に入れるとは思っていなかった.

同じく東北大の Aobayama_dropout も8位と(おめでとうございます!),東北大から top 10 に2チーム入った.top 10 に2チーム,東大と東北大しかないらしい.実は東北大って強豪校ですか?

順位表: https://icpcsec.firebaseapp.com/standings/

freeeから企業賞を頂いた.国内予選でも同じところから企業賞をもらったので,freeeのTシャツがこれで2枚目になる.freee T-shirts for free......

その後Gatherでの懇親会に参加した.Gartic phoneという伝言ゲームで遊んだ.競プロに関することは伝わり方が異常に良くて面白かった.

振り返り

楽しかったというのが1番の感想.序盤,早解きに成功したときは最高の気分だったし,その後も最後まで手を動かしていられた.5時間コンは,終盤不可能しか残らなくてお通夜になってしまうというのがよくあるが,今回は最後までやることが絶えなかった.

suzukaze_Aobayama は,メンバー全員が黄色でかつ得意分野がいい感じにばらけているので,うまく問題を割り振れば中難易度までの早解きには非常に強いと思っている.今回も早解きには大成功し,suzukaze_Aobayama らしい戦い方ができたと思う.

ただやっぱり,橙や赤のメンバーを擁するチームには,完数で負けてしまう.高難易度を時間をかけてでも通す力がこのチームの課題だ.そのためには個人の力が重要になってくる.僕がいつまでも黄色に安住していないで早く橙になるべきなんだよな.来年度の国内までには橙になりたい.

ICPC に関わった全人類に感謝します.特に,チームメイトに感謝したい.仮の人くん,数学と天才パズルを解いてくれるので頼もしかった.来年も数え上げを全部君に投げます.むげんさん,去年のstate_of_the_artから2年間同じチームで戦ったが,今回でICPC引退となる.あなたがいなくなったら,誰が幾何とデータ構造を解けばいいんですか?

新たな3人目のメンバーは誰になるんでしょう.

ICPC 国内予選 2021 参加記

チーム

昨年のチーム state_of_the_art のメンバーのクマノミさんが,外部の大学院に進学されたので,人員補充をしました.東北大でかなりレートを伸ばしてきていた仮の人君に声をかけて,チームを組みました.

むげんさん:典型をたくさん知っている.データ構造に強い.
仮の人君:数学.数え上げは任せた.
ぼく:ぼくの強みは何でしょう...?自己評価は,良く言えばオールラウンダー,悪く言えば器用貧乏

来年はむげんさんがICPCに出れなくなるので,このチームは1年毎に1人入れ替わることになりそう.

チーム名は,去年東北大からアジアに出場したAobayama_dropoutとAobayama_sugarstepをリスペクトして,僕らもAobayama_*にしたいねと話していたのですが,僕がふとsuzukaze_Aobayamaというのを思いついたので言ってみたら思いのほか好評だったのでそれになりました.ぞいぞい.

出たコンテストの振り返り

ICPCチームで出たチームコンを振り返ります.

9/23, 25 ACPC

ACPCのday1, day2に出ました.両日とも,仮の人君が天才プレイを連発してくれました.僕とむげんさんはかなりデータ構造に偏っているので,数学パズルを解いてくれるメンバーが来てくれてとても心強い.なかなか良い結果となりました.

10/23 PG BATTLE

PG BATTLEでは,かつおぶしを解いた僕が大戦犯となり,不甲斐ない成績となりました.やっぱりfull-feedbackしか勝たないんだ.

10/24 模擬国内

去年は,模擬国内で5完して大成功した思い出があります.今年も5完を目標にします.

僕がCを読みます.今年も辛いのがCに来ました.ペナと血反吐を吐きながら1時間ほどかけてAC.中盤がどれだけ解けるかがアジアへの鍵となるので,そこにたっぷり時間を使うためにもこういうのをすばやく解けるようにならなければいけない.

僕がCで苦しんでいる間にチームメイトがさっさとABDを通してくれていました.さらにEも解けそうだと言っているので,僕はFを読みます.大まかな方針は比較的早く立ったのですが,細かいところを詰めるのが下手で実装が大爆発.ちゃんと紙で最後まで詰めてから実装すべきなんだよな.

チームメイトもEが詰めきれず苦戦していたようですが,終盤でちゃんとACしてくれました.僕は最後の数分でようやくサンプルがあったので投げたらWA.絶望.

結果5完,全体47位,学内4位.目標は達成しましたが,5完が遅かったのと,Fが結構通されていたからFも通したかったので,あまり芳しく無い結果となりました.東北大では学内2位がアジアの必要条件なので,だいぶ厳しい.4チームが5完以上って,東北大のレベル上がり過ぎじゃないか?

10/31 KUPC

5時間コンは最高ですね!

僕はCを見ます.実装が下手くそでペナりながら1時間ほどでAC.その間にABDが通っています.Eを仮の人君が考察してくれたので実装をひったくります.しかし投げても投げても通りません!カス!やめた!と叫んで僕はGを見ます.むげんさんが「じゃあぼくがまっさらな気持ちで実装します」というのでお願いします.通ります.え?Gがぜんぜんわからなかったので次に通されていたIを見ます.HLDやるだけ,AC.仮の人君がFを,むげんさんがGを通す.チームメイト,天才しかいない.終盤僕はMを見ていましたが,むずかしかったです.時間切れ.

8完24位.今回はまじで僕が足を引っ張ってばかりだったので多少の悔しさはありますが,それでも24位はかなり良い成績だと思うので良しとします.5時間ありましたが終始やることがあったので楽しかった.

このあと3時間のAGCがありました.AGCはかなり苦手意識があるので,KUPCで疲れたと言い訳をしてサボるつもりでいましたが,魔が差して出てしまいました.結果暖まったのでびっくりしています.

11/3 バチャ

祝日なのでICPC前最後のチーム練習をやります.2007年の国内予選のバチャをやります.ABCDFの5完.良い.相変わらず僕がCでカスみたいなバグを埋め込んで苦しみましたが,代わりにFを通せたのでうれしいです.Eは僕はほとんど見ていませんが,通せるわけがない問題(AOJ-ICPCで1000点)なので解けなくても問題ありません.仮の人君はその後通したらしいです.まじですごいよ.

精進

去年のICPC前はyosupo judge埋めにハマって全然やるべき精進をしていなかったという反省があります.今年はしっかりと問題を解いて力を磨いていくことにしました.

AOJ-ICPCの黄色(450~550点)で星が付いているものを中心に埋めていきました.模擬国内が悔しかったのでその後猛烈なペースで問題を埋め始め,本番までの2週間で多分100問以上ときました.星の数はだいたい実装と考察の比率だと僕は思っているのですが,やっぱり星が多いやつは賢い鮮やかな解法があってとても面白かったです.また実装力を鍛えるために星が少ないやつもいくつか解きましたが,最悪な気分になる問題が多くて楽しかったです.

AOJにはAtCoderにはなかなか出てこない幾何,構文解析,フローが頻出なので,最初の方はなかなか慣れず苦しんでいましたが,段々とコツを掴んできました.特にフローの問題は,フローに対する感覚を養っておかないと嘘DPや嘘貪欲に走りがちですが,AOJ埋めを通してその感覚を磨いていくことができたと思います.実際,最近のABCで珍しく最小カットが出ましたが,僕はAOJのおかげで解けました.最高.

あとは,いくつかのアルゴリズムやデータ構造を履修しました.あまりマニアックなものではなく,そこそこ知られているものを意識的に選んでやりました.Aho-Corasick,中国剰余定理,重心分解,DP高速化テク(CHT,分割統治,Monge)を学びました.まあ多分出ないんですけどね.

ICPC前日は,AOJで易しめの問題をいくつか解いた他(あまりむずかしいのに手を出して解けないまま当日を迎えると気持ち悪いので),蟻本を通読して基本的なテクニックの復習をしました.蟻本は一通り読んだつもりだったのですが,特に上級編を雑に読んでいたことが発覚し,「なにこれ初めて知った」みたいなのが出てきて困ってしまいました.

当日

9時頃起床.いくつか問題を解きます.少し前から考えていたEndless BFSが解けました.あとNim without Zeroを午前中ずっと考えていましたが,昼飯を食っている間に解けました.実装して通ったので家を出て大学に向かいます.東北大学は今日から学祭をやっているらしいです.楽しんでいる人々を横目にAobayamaを登って参加場所に向かいます.

2時半頃に参加場所に集合.部屋がいくつかあって,7チームが3+3+1に分かれることができます.レート和が高いチームから場所を選んでいくことができるんですが,最初のAobayama_dropoutが3の部屋を選んだので,2番目の僕らは個室を占領します.いくら騒いでも大丈夫.部屋が決まったらリハーサルをしました.全完.やったね.

コンテスト開始が4:30なのでそれまでかなり暇です.僕は少し前から考えていたBit Operation Gameに取り組みました.解けました!嬉しい.もう心残りはないです.残り30分位になると,緊張がすごくなります.心臓バクバク.もう新しい問題には取り組まず,タイピングゲームで指を温めたり,深呼吸をして落ち着こうとしたりします.

僕がやるべきことは,Cをさっさと解いて,中盤にしっかり時間をかけて解くこと.がんばるぞい.

開始

今年もコドフォりました!去年と全く同じ展開じゃないか?ただ今年は復旧が早かったので1分ほどで始まります.

Cを見ます.待ってました,重実装問.俺はこのためにAOJを埋めてきたんだ.大まかな方針はすぐ立ちます.構文解析して木DP.全方位木DPっぽくも見えるんですが,制約的に各頂点から木DPしても間に合うと判断し,実装開始.AOJ埋めのおかげで構文解析パートはスムーズに書けます.次に木DPのパートですが,結構大変です.子に順番があるので,ちゃんと気にしてあげなければいけません.ABがすぐ通って,Cにhelpはいらないかと聞かれましたが,見通しは立っているので大丈夫と宣言し,先に進んでもらいます.

僕がダラダラと実装をしている間に,チームメイトがDを考えていました.チームメイトの会話に片耳突っ込んで聞いていましたが,仮の人君が天才をしてくれたらしく通ります.すげー.

僕は1時間ほどでCが書き終わったのでテストケースを実行します.永遠に終わりません.困ったので脳死でメモ化したら爆速になって笑ってしまいました.提出する直前になって,詰めきれてなさそうなパートに思い当たったのですが,まあ落ちたら考えればいいだろと思って投げると通りました.わーい.

(追記: 僕は制約を勘違いして2乗が通ると思いこれを書いたのですが,普通に2乗では全然間に合いません.それに,普通はメモ化したところで速くなるはずはないです.ただしこの問題では木の頂点の次数が小さいことから,メモ化するとオーダが落ちます.勘違いから意図せず実装が楽な方針を引き当てたので,本当に運が良かっただけです)

この時点で4位だったらしいです.まあどうせ中盤で抜かれるので,序盤速くてよかったねという感じにしかこの時点では思いませんでした.

むげんさんがEをやっていて,方針が立っているらしいので仮の人君と一緒にFを見ます.二部マッチングですか.DM分解が頭をよぎります.少し前のABCでこれを貼るだけが出たので,その際に履修したのですが,ちょっと似ていたので使えるんじゃないかと思いました.ただ僕は使い方が思いつかなかったのでそれはすぐ捨てました.ちょっと考えると,残余グラフを強連結成分分解していい感じに辺を追加すればいいということに思い当たります.細かいところを詰めるのを仮の人君に任せ,僕はSCCするところまで実装します.その間にEをやっていたむげんさんが何度も叫んでは実装をやり直すというのを繰り返していましたが,突然「既出です!」と叫んで解き切ってくれます.むげんさんが床を転げ回って喜びます.

(追記: DM分解を捨てたと言っていますが,DM分解です.というか残余グラフをSCCするというのはDM分解そのものです.僕が知っていたのは coarse decomposition *1で,ここで出てきたのは fine decomposition の方でした)

5完!しかし模擬では5完で学内4位だったので,まだ油断はできません.ところで,まだ全体5位らしいです.嘘だろ.

Fを解きたい.みんなでFを考えています.SCCまではあってそうなんですが(僕がSCCだと主張しているだけでチームメイトはあまりピンときていなかったらしい),そこからどうやって辺を張るのかを詰めるのが難しい.嘘解法を見つけては反例を発見するの繰り返し,なかなか詰まりません.一度サンプルがあったのでダメ元で投げてみましたが当然WA.難しい....

しかし残り10分位になってまだFがほとんど通っていないことを知り,まだ全体6位であることを知ると,勝利を確信します.

振り返り

5完,全体6位,学内1位,アジア進出です.

僕らが一番ビックリしています.30位くらいになってアジアに行ければ大成功だと思っていました.始まる前,「学内2位に入れなくても全体20位に入ればアジアいけますよ〜」「いや無理だろ」と冗談を言っていたんですが,それすら超えてくるとは....夢のようです.

かなり運が良かったと思います.Eまでが僕らが解ける程度の難易度で,Fに大きな崖があったこと.黄色3人で勝てるのはこの傾斜しかないです.そしてチームメンバーそれぞれが得意な問題を担当できたこと.後で聞いたんですが,Dは天才だったらしいです.天才は仮の人君にしか解けない.Eは,似た問題を解いたことがあるむげんさんにしか解けない.そして実装を鍛えてきた僕はだれも解きたくないCを1人でちゃんと通すことができました.まさに適材適所.

誰か1人がキャリーしたという感じではなく,3人全員で勝つことができたと思います.誰か1人欠けていても勝てなかったと思います.感謝の気持でいっぱいです.

コンテスト後,チームメイトとラーメンを食べました.うまかったぜ.

アジアもがんばってきます.

AGC050C - Block Game

最高に面白い

問題

atcoder.jp

解法

重要な観察として,すぬけ君が動ける範囲はBが来るたびに半分以下になるというのがある.最終的にすぬけ君が勝つ,すなわちすぬけ君の動ける範囲が1マス以上あるためには,Bを置く回数は\(O(\log n)\)回程度でなければならない.

すぬけ君が勝つような文字列を数え上げて,\(2^Q\)からそれを引くという方針が立つ.\(i\)回目の操作のあとすぬけ君の動ける範囲が\(k\)マスの時,\(i-k\)回目の操作後すぬけ君が\(2k\)マス動ける状態にあり,また\(i-k\)回目以降にBが来ない.知りたいのは\(n\)回目の操作のあとすぬけ君が1マス以上動ける文字列の数だから,
\[
dp[i][j] := \text{\(i\)回目の操作後すぬけ君が\(2^j\)マス以上動ける状態にある文字列の数}
\]
としてDPをすればよい.

遷移は,
\[
dp[i][j] \leftarrow dp[i][j] + \begin{cases}
dp[i - 2^j][j + 1] & \text{if \(s[i] = B\) and \(s[i - 2^j...i-1]\) doesn't contain B} \\
dp[i-1][j] & \text{if } s[i] = S
\end{cases}
\]
みたいな感じになる.?ならどっちもやる.Bが区間に含まれるかどうかの判定は累積和でやる.

提出

atcoder.jp

ついでにshortestもとっておいた.

感想

2冪だけ覚えておくだけでいいのすごく非直感的.めちゃくちゃ面白い問題だと思った.コンテスト中にこういうのが通せるようになりたいなあ

ABC129F - Takahashi's Basics in Education and Learning

問題

atcoder.jp

解法

公式解説は行列累乗で解いているが,僕はダブリングで解いた(まあやっていることが実質同じといえば,それはそう).

まず\(d\)桁の項が\(s_l\)から\(s_r\)であるとする.\(l, r\)は二分探索で求められる.

この区間での和は
\[
10^{d(r-l)}(A+Bl) + 10^{d(r-l-1)}(A+B(l+1)) + \dots + 10^d (A+B(r-1)) + (A+Br) \\
= A(1 + 10^d + \dots 10^{d(r-l)}) + B(r + (r-1)10^d + \dots + l10^{d(r-l)})
\]
これに,\(s_{r+1}\)以降の桁数だけ後ろに0をつける.これは\(d\)を逆順に見て適当に持っておけばいい.

\(A, B\)の係数を計算する.等比数列なので和の公式を使いたくなるが,今回modが素数とは限らないので,分母の逆元が計算できないことがある.だから公式を使わずに和を直接計算する.

1項ずつ見ていったら当然間に合わないのでダブリングをする.まず\(A\)の係数を考える.これは,
\[
tableA[i] := 1 + 10^d + \dots + 10^{d (2^i - 1)}
\]
として,
\[
tableA[i] = tableA[i - 1] (1 + 10^{d 2^{i-1}})
\]
というように更新していけばよい.

\(B\)の方は,頭をバグらせないようにちょっと変形する.
\[
r + (r-1)10^d + \dots + l10^{d(r-l)} = r (1 + 10^d + \dots 10^{d(r-l)}) - (0 + 10^d + \dots + (r - l)10^{d(r-l)})
\]
とする.右辺第1項は\(A\)の係数を使いまわす.第2項を計算する. これはまた式変形をすると分かるのだが,
\[
tableB[i] := 0 + 10^d + \dots + (2^i - 1)10^{d (2^i - 1)}
\]
として,
\[
tableB[i] = tableB[i - 1] + (tableB[i - 1] + 2^{i-1}tableA[i - 1])10^{d2^{i-1}}
\]
とすればよい.

これらの前計算した値を使って答えが求まる.

atcoder.jp

感想

考察はかなり自然に思いつけるから,実装テクニックを問う問題だと思うんだけど,ダブリングも行列累乗も有名なテクだからこんなにdiffが高くつくのはなんでだろうな...という気持ち.

行列累乗の方がきれいなので,思いつけるようになりたいなあ.

ARC102-E Stop. Otherwise...

形式べき級数楽しい

問題

atcoder.jp

解法

サンプルを見ると明らかに解には対称性があるので,\(i \leq K + 1\)の時のみ考える.

解を多項式の係数に帰着する.\(x^n\)の係数が,数字を\(n\)個選んだ時の組み合わせの数という風に.

まず,\(a + b = i\)なる\(a, b\)に関しては,片方選べばもう一方は選べないが,その選んだ方は何回選んでもいいので,\(1 + 2x + 2x^2 + \dots\)を掛けることに対応する.このようなペアは,\(i = 2k\)のとき\(k - 1\)個,\(i = 2k + 1\)のとき\(k\)個ある.

\(2a = i\)なる\(a\)に関しては,それは一度選んだらもう選べないので,\(1 + x\)を掛けることに対応する.このような\(a\)は,\(i = 2k\)の時に1個ある.

どのように選んでも問題ない数に関しては,\(1 + x + x^2 + \dots\)を掛けることに対応する.このような数は\(i = 2k\)のとき\(K - 2k + 1\)個,\(i = 2k + 1\)のとき\(K - 2k\)個ある.

以上より,解は
\[
\begin{cases}
(1 + 2x + 2x^2 + \dots)^{k-1} (1+x) (1 + x + x^2 + \dots)^{K-2k + 1} & \text{if } i = 2k \\
(1 + 2x + 2x^2 + \dots)^{k} (1 + x + x^2 + \dots)^{K-2k} & \text{if } i = 2k + 1 \\
\end{cases}
\]
の\(x^N\)の係数となる.

ここで,\(1 + 2x + 2x^2 + \dots = \frac{1+x}{1-x}\),\(1 + x + x^2 + \dots = \frac{1}{1 - x}\)となるから,上の式は\(i\)の偶奇にかかわらず
\[
\frac{(1 + x)^k}{(1 - x)^{K - k}}
\]
と書き直せる.

分子は二項定理,分母は負の二項定理を用いてそれぞれ展開し,\(x^N\)の係数を求めればいい.計算量は\(O(KN)\)

終わり.下の提出はNTTを使っているんですが,これはなんか勘違いしていてNTTを使わないといけないと思っていた.求めるのは1項だけなのでNTTはいりません.

atcoder.jp


形式的べき級数は頭を使わなくても機械的な計算で答えが出るのでうれしい.

yosupo judge データ構造全部解いた

yosupo judge のデータ構造を全部解いた.解法と感想をメモする.

judge.yosupo.jp

Associative Array

実はこれ解いていない(タイトル詐欺).AC取るだけならmap使うだけなので,さすがにいいかなという感じだった.見逃してください.†X-fast trie†を書いたらこれでverifyする.書くとは言ってない.

Unionfind

union find っていろいろな実装法があって奥が深い.僕のは普通にunion by size,経路圧縮.

Static Range Sum

累積和,Fenwick tree,セグ木,disjoint sparse table, sqrt treeを比較した.Fenwick tree < 累積和 < セグ木 < sqrt tree < DST だった.Fenwick treeが累積和より速いのは本当に謎.

Static RMQ

セグ木,スパーステーブル,DST,Sqrt treeを比較した.セグ木 < sqrt tree < スパーステーブル < DST だった.線形のやつもそのうちやりたい.

Point Add Range Sum

セグ木.僕はモノイドの構造体を渡す設計にしている.僕は結構意味論にこだわっていて,セグ木はモノイドを扱うデータ構造なので,モノイドを渡す方がしっくりくる.また,モノイドを定義する性質が一つのところでまとめられるので可読性が高いという利点がある.

Point Set Range Composite

セグ木.可換性を要求しない実装にする.

Range Affine Range Sum

遅延セグ木.遅延セグ木の設計も結構こだわりがある.テンプレートパラメータでモノイド,作用素モノイド,作用の3つを渡すようにしている.まずACLみたいに全部バラバラに渡すのは可読性が著しく損なわれると感じていて,関連する情報はまとめてモノイド構造体にした方がすっきりする.しかし,これらの3つ全部をまとめたものになにか名前がついているわけでもないので(ある?),これ以上まとめる必要はないと考え,この設計を採用している.また,作用に区間の長さを渡すかどうかも流派が分かれるポイントだが,僕は渡しません.代わりにモノイドを(値,区間の長さ)のペアにする.こうすることで何がうれしいのかというと,作用の分配則が成り立つということ.別に分配則を成り立たせることに意味はないが,こっちの方が美しい気がする.

また,僕のセグ木は非再帰だが遅延セグ木は再帰で書いている.遅延伝搬のパートを非再帰でやろうとすると頭がバグるので,わかりやすい再帰で書いている.可読性と効率のトレードオフは悩みポイントで,割とその時の気分.

一応動的構築版も作ったが,これいる?

Range Chmin Chmax Add Range Sum

Segment tree beats.実装はやばいが持ってると結構殴れるのでオススメ.

Range Kth Smallest

Wavelet matrix.このデータ構造何でもできるので本当にすごい.しかしこの問題を解く以外で使ったことはない.他にも Mo + set,Mo + Fenwick tree,Mo + binary trie,平方分割,領域木などでできそう.

Vertex Add Path Sum

HL分解.これも結構使える.

Vertex Set Path Composite

HL分解.これ頑張った.HL分解を非可換でやる方法に関する文献が見つからなかったので自己流.HL分解は気を付けて実装すると上から下に順番に計算するようにできる.これを使って,lcaからu,lcaからvまでそれぞれ計算して,lcaからuの値をひっくり返す(これをするには,2方向から計算した時の値をモノイドに覚えさせておけばいい).ここでlcaを2回数えているので,lcaの逆元を掛ければuからvまでを非可換で計算できる.もっといいやり方がありそう.

(2022/5/11 追記)

別にこんな面倒なことをしなくても,HLDの中身をちゃんと見ればできる.木を登るときに,パスの両端からの積を別々に持っておいて適当にswapしながら登る.最後にflipをする.詳細は
Heavy-Light Decomposition | sotanishy’s competitive programming library

Vertex Add Subtree Sum

オイラーツアー + セグ木.オイラーツアー,すごく実装が軽くて好き.HL分解の配列を使いまわして部分木クエリができるのもおもしろい.

Dynamic Sequence Range Affine Range Sum

遅延Treap.平衡二分探索木,いろいろ種類があって面白い.たくさんの種類コレクションして鑑賞したいな.

Dynamic Tree Vertex Add Path Sum

link/cut木.とても難しい.toptreeとかいうのでもできるらしいがさらにヤバいらしいので履修してない.

Dynamic Tree Vertex Set Path Composite

link/cut木.非可換であることに注意.

Dynamic Tree Vertex Add Subtree Sum

Euler tour木.解説記事が非常に少なく実装にてこずった.結局hotmanさんのをパクって参考にして書いた.個人的にはET木の方がLC木よりも圧倒的に理解しやすく(オイラーツアーを,動的にしただけ),パスクエリはLC木,部分木クエリはET木と使い分けるのがよさそう.あとdynamic connectivityでもできる.

Dynamic Tree Subtree Add Subtree Sum

ET木で遅延伝搬する.ET木と遅延ET木を分けているがどうせどちらも永遠に使わないので遅延の方だけ残すか.

Dynamic Graph Vertex Add Component Sum

offline dynamic connectivity.個人的にこれ結構気に入っている.時間方向にセグ木のような構造をもって,undoable UF を持ちながらDFSするとできるというやつ.このテクは一般化できて,要素の追加とundoができるが削除ができないデータ構造で削除をしたいときにこのテクが使える.ET木を書いたので,online dynamic connectivityもそのうちやるかも.

Set Xor-Min

binary trie.これ好き.

Line Add Get Min

Li Chao tree.個人的にCHTよりもLi Chao treeの方がわかりやすかった.

Segment Add Get Min

Li Chao tree.僕のは静的なのでクエリを先読みして端点を知る必要があるが,動的にするとオンラインで行ける.そこまで遅くならないらしいので,動的もそのうち整備しようかな.

Queue Operate All Composite

sliding window aggregation.これもかなり好き.実装も軽くてうれしい.

キューじゃなくてデックも作れるらしいが,どうやるかはわからない.キューと同じようにやると,回転しまくるケースを作ることができるので,クエリあたり最悪O(N)になってしまう.じゃあ半分だけ回せばよさそうで,実際半分だけ回転させる実装を見たことがあるが,本当にそれで償却O(1)になるのか?よくわからない.

あと,なんか,スタックを6個使うと最悪O(1)になると聞いたことがある.やばそう.

Static Range Inversions Query

これは何をverifyするためにあるんですか?座圧 + Mo + Fenwick tree でO(N√N logN).クエリを並べ替えるだけで計算量が落ちるってすごい.

Rectangle Sum

永続セグ木.もっといいやり方がたくさんあるが永続セグ木をverifyしたかった.座圧してから,すべてのx座標に対して一本ずつy方向のセグ木を持つことを考える.xでのセグ木をst[x]としたとき,st[x] = (st[x-1] にこのx座標での点を追加)という風にする.こうすることで,st[x].fold(l, r)を呼ぶと,x座標は 0 から x ,y座標は l から r までの領域内の点の総和が求められる.当然,普通のセグ木をN本も持ったらだめだが,st[x]はst[x-1]との差分だけ持っておけばいいので,永続セグ木が使える.ACは取れるが(少なくとも僕の実装は)めちゃくちゃ遅い.

Point Add Rectangle Sum

四分木.2Dセグ木のメモリ効率がいい作り方がよくわからなかったので調べたら四分木というのが出てきたので実装してみた.これはセグ木にセグ木を乗せるのではなく,セグ木をそのまま2次元に拡張したものみたいなイメージ.1次元のセグ木は右と左に子ノードがある二分木だが,四分木は4つの象限に子ノードを持つ.動的に構築すればメモリは大丈夫.これめっちゃ使いやすいと思ったが,よく考えると長方形領域取得の最悪計算量がめちゃくちゃ悪い(O(N log N)?).まあよほど意地悪な入力じゃない限りそこそこ高速っぽい.どうせ四分木を落とすためのケースなんて作ってないだろ.たださすがにコンテストで使うのは怖いのでそのうち2Dセグ木もちゃんとやる.

Persistent Queue

永続キュー.永続配列を使ったクエリO(log N)の実装で通した.

まず最初はSWAGのように永続スタック2本(永続スタックは非常に簡単に作れる)でO(1)の永続キューを作ろうとしたが,two_stacks_killerというテストケースに殺される.SWAGはO(N)の回転操作をたまにやることで償却O(1)になっているが,永続キューの場合は過去の回転する前のバージョンにアクセスされると回転をしまくることになってまずいんですかね(適当).

ダブリングでもO(log N)でできるらしい.永続配列使うよりこっちの方が速そう.他にもO(1)でやる方法がたくさんあるっぽい.

ところで,永続配列は何分木にするのが一番いいのだろう.僕はn分探索の中で2分探索が最強みたいなノリで適当に2分木をデフォルトにしている.log(n)分木にするのが一番いいと聞いたことがあるがなんでそれがいいのかよくわかっていない.getとsetで計算量が違くて,getは子の数が多い方がよくてsetは子の数が小さい方がよいっていう認識であってるのかな.

Persistent Unionfind

永続union find.永続配列をつかって UF をやるだけ.

永続データ構造は配列,セグ木,スタック,キュー,UFを持っている.永続データ構造のよく見る実装は,ノードのポインタを外から渡していろいろ処理するという感じだが,個人的にはポインタは隠して普通の配列やセグ木を扱うようなインターフェースにしたいと考えている.更新した後の返り値はそのクラスそのものにしている(たとえば,永続配列で要素を変更する関数は PersistentArray set(int i, T x); みたいな感じ).まあ実質ポインタを管理しているのと大して変わりないのだが.

参考にしたサイト

hotmanさんの記事

qiita.com

大体これを見て知らないデータ構造とかのヒントを得ました.ありがとうございます.

その他,多くの方々のライブラリや記事を参考にさせていただきました.名前を挙げるときりがないので挙げませんが,ありがとうございます.