noshi91のメモ

データ構造のある風景

ABC339-F Product Equality の確率を保証する

 \lbrack 10 ^ 9 , 2 \times 10 ^ 9 \rbrack から一様ランダムに素数  P を選ぶ *1 A _ i A _ j \neq A _ k となる  1 組の  ( i , j , k ) を固定し、誤判定する、すなわち  A _ i A _ j \equiv A _ k \pmod P となる確率を求める。 これは  A _ i A _ j - A _ k P で割り切れることと同値である。

 \lbrack 10 ^ 9 , 2 \times 10 ^ 9 \rbrack の範囲に素数 47374753 個存在する *2 。 一方、 0 \lt \lvert A _ i A _ j - A _ k \rvert \leq 10 ^ { 2000 } であるから、 A _ i A _ j - A _ k は範囲内の素因数を  222 個以下しか持たない。 よって一様ランダムに選んだとき、誤判定の確率は  222 / 47374753 \leq 10 ^ { - 5 } 以下である。

 4 つの素数を一様ランダムかつ独立に選ぶと、全ての素数で同時に誤判定する確率は  \left ( 10 ^ { - 5 } \right ) ^ 4 = 10 ^ { - 20 } 以下である。 独立に選ぶ部分が重要で、独立な事象が同時に発生する確率はその積となるという性質を用いて確率を抑えている。

 1 つのケースで最大  N ^ 3 \leq 10 ^ 9 組の  ( i , j , k ) に対して判定を行うので、 1 回以上誤判定する確率は  10 ^ { - 20 } \times 10 ^ 9 = 10 ^ { - 11 } 以下である。 この部分は union bound ( Boole's inequality - Wikipedia ) と呼ばれる直感的にはほとんど明らかな不等式を用いた。 この  N ^ 3 組では同じ  P を使っているため、これらの事象は独立ではない。 よって  1 - \left ( 1 - 10 ^ { - 20 } \right ) ^ {10 ^ 9} と計算するのは誤りである

テストケース間では  P を取り直すのでそれぞれは独立となるが、独立性を気にせず union bound を使った方が簡単になる。 全体で  1000 ケースある場合、 1 ケース以上 WA となる確率は  10 ^ { - 11} \times 1000 = 10 ^ {- 8} 以下である。  1000 ケースという想定が甘いかどうかと、 10 ^ {- 8 } という確率が許容可能かどうかは各自で判断するしかない。 人生が終わるまでにジャッジできる量、宇宙線アルゴリズムが破損するより低い確率、などで無理やり上限を付けることはできるかもしれない。

(おまけ) 乱択に工夫をするときに気を付けること

本問で独立に素数を選んだ部分で、既に選んだ素数を選ばないようにするという工夫が考えられる。 同じ素数を選ぶことは得にならないのでこれは改善になっているが、十分多い候補から素数を選ばなくてはならない都合上、同じ素数を選ぶ確率自体がとてつもなく小さいため効果が無視できるほどしかない。

また、乱択繋がりで rolling hash の基数に原始根を選ぶという工夫がある。 ランダムケースに対しては改善になっているかもしれないものの、もし作問者が参加者のそのような考えを見抜いて原始根を集中的に落とすケースを用意した場合は改悪となる *3

工夫によりアルゴリズムの挙動が複雑になると、確率の証明を付けることが難しくなりやすい。 そもそも証明により安全性が保障されたシンプルなアルゴリズムがあるのだから、それをそのまま使えばよい。

*1: 範囲内から一様ランダムに整数を選び、素数判定することを繰り返す

*2: WolframAlpha に PrimePi[2*109]-PrimePi[109-1] と入力する

*3: 乱数の出目は作問者には想定のしようがないことに注意

FFT の回数を削減するテクニック集

載ってないテクニックを募集しています。

前提

 \omega _ n := \exp \left ( 2 \pi i / n \right ) と置く。 以降  i と書いたらそれは虚数単位ではなく添え字である。

数列  \left ( a _ i \right ) _ { i = 0 } ^ { n - 1 } に対し、その DFT  \left ( \hat a _ j \right ) _ { j = 0 } ^ { n - 1 }  \hat a _ j = \sum \limits _ { i = 0 } ^ { n - 1 } \omega _ n ^ { i j } a _ i で定義する。

基本事項

以下の事実は式変形で確認できる。

関係式  \left ( 1 \right )

 \left ( a _ i \right ) _ { i = 0 } ^ { 2 n - 1 } に対し、 \left ( e _ i \right ) _ { i = 0 } ^ { n - 1 } , \left ( o _ i \right ) _ { i = 0 } ^ { n - 1 }  e _ i = a _ i + a _ { n + i } , o _ i = \omega _ { 2 n } ^ i \left ( a _ i - a _ { n + i } \right ) で定義する。 このとき

\begin{align} \hat a _ { 2 j } &= \hat e _ j \cr \hat a _ { 2 j + 1 } &= \hat o _ j \end{align}

が成り立つ。 すなわち、長さ  2 n の DFT の偶数部分と奇数部分は、それぞれが長さ  n のある列の DFT として与えられる。

参考: Cooley-Tukey Algorithm

関係式  \left ( 2 \right )

 \left ( a _ i \right ) _ { i = 0 } ^ { 2 n - 1 } に対し、 \left ( e _ i \right ) _ { i = 0 } ^ { n - 1 } , \left ( o _ i \right ) _ { i = 0 } ^ { n - 1 }  e _ i = a _ { 2 i } , o _ i = a _ { 2 i + 1 } で定義する。 このとき

\begin{align} \hat a _ j &= \hat e _ { j \bmod n } + \omega _ { 2 n } ^ j \hat o _ { j \bmod n } \end{align}

が成り立つ。 すなわち、長さ  2 n の DFT は、元の列の偶数部分と奇数部分それぞれの DFT を基に計算できる。

参考: https://hcpc-hokudai.github.io/archive/math_fft_001.pdf

DFT の使いまわし

もっとも基本的なテクニック。 同じ列の DFT を複数回使うときは  1 度計算すれば十分である (それはそう)。 やや気付きにくいケースとして、IDFT を行った列を後々もう一度 DFT し直している場合がある。 多項式 2 乗を求めるときや、有理式を通分して足し合わせるときなど、あらゆる場面で使う。

middle product

 \left ( a _ i \right ) _ { i = 0 } ^ { n + m - 1 }  \left ( b _ i \right ) _ { i = 0 } ^ { n } を畳み込んだ列を  c として、 \left ( c _ i \right ) _ { i = n } ^ { n + m - 1 } だけが必要となることがある。 DFT を用いた畳み込みは巡回畳み込みであるため普通は巡回しないよう十分な長さを取るが、この場合はいくらか巡回しても問題なく、長さが  n + m 以上の巡回畳み込みで足りる。 特に  n = m = 2 ^ k の場合が頻出で、計算量が半分になる。 ワイルドカードマッチング、転置原理を用いて導出されたアルゴリズム、オンライン畳み込みなどで用いる。

 2 倍の長さの DFT の計算

 \left ( a _ i \right ) _ { i = 0 } ^ { n - 1 } について  \hat a が求まっているとして、 a を長さ  2 n に伸ばして  0 埋めした列  \displaystyle \left ( b _ i \right ) _ { i = 0 } ^ { 2 n - 1 } , b _ i = \begin{cases} a _ i \quad \left ( i \lt n \right ) \cr 0 \quad \left ( i \geq n \right ) \end{cases} に対する  \hat b を計算したい場面がある。 関係式  \left ( 1 \right ) より  \hat b _ { 2 j } = \hat a _ j であり、 \left ( c _ i \right ) _ { i = 0 } ^ {  n - 1 } , c _ i = \omega _ { 2 n } ^ i a _ i に対して  \hat b _ { 2 j + 1 } = \hat c _ j である。 よって  \hat c を計算するだけでよい。 subproduct tree を作るときなどに用いる。

 f \left ( - x \right ) の DFT の計算

 \left ( a _ i \right ) _ { i = 0 } ^ { 2 n - 1 } について  \hat a が求まっているとして、 \displaystyle \left ( b _ i \right ) _ { i = 0 } ^ { 2 n - 1 } , b _ i = \left ( - 1 \right ) ^ i a _ i に対する  \hat b を計算したい場面がある。  \hat b _ j = \sum \limits _ { i = 0 } ^ { 2 n - 1 } \omega _ { 2 n } ^ { i j } \left ( - 1 \right ) ^ i a _ i = \sum \limits _ { i = 0 } ^ { 2 n - 1 } \omega _ { 2 n } ^ { i \left ( j + n \right ) } a _ i = \hat a _ { \left ( j + n \right ) \bmod 2 n } であるから、 \hat a から  \Theta \left ( n \right ) で計算できる。 Bostan-Mori 法で用いる。

 f \left ( x ^ 2 \right ) の DFT の計算

 \left ( a _ i \right ) _ { i = 0 } ^ { n - 1 } について  \hat a が求まっているとして、 \displaystyle \left ( b _ i \right ) _ { i = 0 } ^ { 2 n - 1 } , b _ i = \begin{cases} a _ { i / 2 } \quad \left ( i \equiv 0 \pmod 2 \right ) \cr 0 \quad \left ( i \equiv 1 \pmod 2 \right ) \end{cases} に対する  \hat b を計算したい場面がある。 関係式  \left ( 2 \right ) より  \hat b _ j = \hat a _ { j \bmod n } であるから、 \hat a から  \Theta \left ( n \right ) で計算できる。

前節の内容と組み合わせると以下の事実が分かる。 次数  2 n - 1 多項式  g について  g \left ( x \right ) g \left ( - x \right ) = f \left ( x ^ 2 \right ) が成り立っているとき、 g の DFT の前半と後半を掛け合わせると  f の DFT が得られる。

 x f \left ( x \right ) の DFT の計算

 \left ( a _ i \right ) _ { i = 0 } ^ { n - 1 } について  \hat a が求まっているとして、 \displaystyle \left ( b _ i \right ) _ { i = 0 } ^ { n } , b _ i = a _ { \left ( i - 1 \right ) \bmod n } に対する  \hat b を計算したい場面がある。  \hat b _ j = \sum \limits _ { i = 0 } ^ { n } \omega _ { n } ^ { i j } a _ { \left ( i - 1 \right ) \bmod n } = \sum \limits _ { i = 0 } ^ { n } \omega _ { n } ^ { \left ( i + 1 \right ) j } a _ i = \omega _ { n } ^ j \hat a _ j であるから、 \hat a から  \Theta \left ( n \right ) で計算できる。

偶数部分 (奇数部分) だけ IDFT

未知の  \left ( a _ i \right ) _ { i = 0 } ^ { 2 n - 1 } について  \hat a が分かっているとして、 \left ( e _ i \right ) _ { i = 0 } ^ { n - 1 } , e _ i = a _ { 2 i }  \hat e が求めたい場面がある。 関係式  \left ( 2 \right ) と同様に、 a の偶数部分と奇数部分をそれぞれ  e , o と定義する。 関係式  \left ( 2 \right ) より、 0 \leq j \lt n について  \hat a _ j = \hat e _ j + \omega _ { 2 n } ^ j \hat o _ j , \hat a _ { n + j } = \hat e _ j + \omega _ { 2 n } ^ { n + j } \hat o _ j = \hat e _ j - \omega _ { 2 n } ^ j \hat o _ j であるから、 \hat e _ j = \left ( \hat a _ j + \hat a _ { n + j } \right ) / 2 である。 よって  \hat e  \Theta \left ( n \right ) で求まり、それを IDFT することで  e も求まる。  o に関するものも同様に計算できる。

低次部分 (高次部分) の DFT の取り出し

未知の  \left ( a _ i \right ) _ { i = 0 } ^ { 2 n - 1 } について  \hat a が分かっているとして、 \left ( \mathrm { lo } _ i \right ) _ { i = 0 } ^ { n - 1 } , \mathrm { lo } _ i = a _ i に対する  \hat { \mathrm { lo } } が求めたい場面がある。  \left ( \mathrm { hi } _ i \right ) _ { i = 0 } ^ { n - 1 } , \mathrm { hi } _ i = a _ { n + i } と定義する。 関係式  \left ( 1 \right ) より、 \hat a の偶数部分は  \mathrm { lo } + \mathrm { hi } の DFT である。 また、 \hat a の奇数部分を IDFT することで  \omega _ { 2 n } ^ i \left ( \mathrm { lo } _ i - \mathrm { hi } _ i \right ) が求まる。 ここから  \mathrm { lo } - \mathrm { hi } が求まるため、これを DFT して  \hat a の偶数部分と足し合わせることで  2 \mathrm { lo } の DFT が求まる。  \hat { \mathrm { hi } } も同様に計算できる。

長さ  2 n の IDFT ののち長さ  n の DFT を行うことと比べて  2 / 3 倍の時間計算量になる。

 f \left ( x ^ { - 1 } \right ) や反転した列の DFT

 \left ( a _ i \right ) _ { i = 0 } ^ { n - 1 } について  \hat a が求まっているとして、 \displaystyle \left ( b _ i \right ) _ { i = 0 } ^ { n } , b _ i = a _ { \left ( - i \right ) \bmod n } に対する  \hat b を計算したい場面がある。  \hat b _ j = \sum \limits _ { i = 0 } ^ { n } \omega _ { n } ^ { i j } a _ { \left ( - i \right ) \bmod n } = \sum \limits _ { i = 0 } ^ { n } \omega _ { n } ^ { - i j } a _ i = \hat a _ { \left ( - j \right ) \bmod n } であるから、 \hat a から  \Theta \left ( n \right ) で計算できる。

IDFT の計算において回転子を  \omega _ n から  \omega _ n ^ { - 1 } に取り換えた DFT を計算するがこれも同じ式になるので、FFT の実装において IDFT 用に回転子の切り替えを可能にしておく必要はないことが分かる。 もっとも、FFT に凝っている人は DFT の結果を bit 反転された状態で管理すると思われ、その場合 IDFT は別で実装することになるので回転子を一致させる旨味は少ないだろう。

 x f \left ( x \right ) の DFT の計算と組み合わせると  a を反転した列の DFT も  \Theta \left ( n \right ) で求まる。

push_back / suffix の xor 基底クエリ Θ(log(A))

概要

最初空の数列  A があり、以下のクエリを処理するアルゴリズムを説明する。

  •  x が与えられ、 A の末尾に  x を追加する。
  •  i が与えられ、 \left \lbrace A _ i , A _ { i + 1 } , \dots \right \rbrace の xor に関する基底  \left ( b _ 0 , b _ 1 , \dots \right ) を出力する。
    • ただし  x の highest set bit を   \operatorname { bsr } \left ( x \right ) := \max \left \lbrace k \mathrel { } \middle \vert \mathrel { } x _ k = 1 \right \rbrace と定義したとき  \operatorname { bsr } \left ( b _ j \right )  j について単調減少である。

時間計算量はいずれも  \Theta \left ( \log \left ( \max A \right ) \right ) である。

アルゴリズム

 0 \leq x \lt 2 ^ w のクエリだけが与えられるとして、各クエリを  \Theta \left ( w \right ) で処理するアルゴリズムを説明する。 多少工夫すれば  w を各クエリ時点での  \log \left ( \max A \right ) に合わせることもできる。

基本となる操作は  l \lt r に対する  A _ l \leftarrow A _ l \oplus A _ r という操作である。  A のいかなる suffix に着目してもこの操作は基本変形になっているから、この操作だけを行っている限りは  \operatorname { span } \left \lbrace A _ i , A _ { i + 1 } , \dots \right \rbrace は変化しない。

 A の添え字の列  \left ( i _ 0 , i _ 1 , \dots , i _ { w - 1 } \right ) を管理する。 これらは以下の不変条件を維持する。

  •  \operatorname { bsr } \left ( A _ { i _ k } \right ) = k
  •  i _ k として現れない  i について  A _ i = 0

 A が最初からいくつかの要素を持っていてもクエリには影響しないため、最初  A = \left ( 2 ^ 0 , 2 ^ 1 , \dots , 2 ^ { w - 1 } \right ) として  i _ k = k とすれば条件を満たす。  A = \left ( A _ 0 , A _ 1 , \dots , A _ { n - 1 } \right ) の末尾に要素を追加した場合、以下の処理を実行し不変条件を回復する。

\begin{align} & j := n \cr & \mathtt { for } ~ k := w - 1 ~ \mathtt { to } ~ 0 \cr & \quad \quad \mathtt { if } ~ \left ( A _ j \right ) _ k = 0 \cr & \quad \quad \quad \quad \mathtt { continue } \cr & \quad \quad \mathtt { if } ~ j \gt i _ k \cr & \quad \quad \quad \quad \left ( j , i _ k \right ) \leftarrow \left ( i _ k , j \right ) \cr & \quad \quad A _ j \leftarrow A _ j \oplus A _ { i _ k } \end{align}

アルゴリズムの正当性は以下の事実に着目することで確認できる。

  • ループの始まりでは  \operatorname { bsr } \left ( A _ j \right ) \leq k であり、ループの終わりでは  \operatorname { bsr } \left ( A _ j \right ) \lt k である。
  • 不変条件が、 A _ j \neq 0 となり得ることを除き常に満たされている。
  •  A を変更する操作は、 l \lt r に対する  A _ l \leftarrow A _ l \oplus A _ r の形でしか行われていない。

処理の終わりでは  A _ j = 0 となっているから、不変条件が完全に満たされる。  \left \lbrace A _ i , A _ { i + 1 } , \dots \right \rbrace の基底は  \left \lbrace A _ { i _ k } \mathrel { } \middle \vert \mathrel { } i \leq i _ k \right \rbrace であるから、ただ出力すればよい。

実際には  A を陽には管理せず、 \left ( i _ k , A _ { i _ k } \right ) を管理すると空間計算量が抑えられる。

応用

区間 xor 基底クエリ

静的な区間 xor 基底クエリが  \Theta \left ( \left ( N + Q \right ) \log \left ( \max A \right ) \right ) で処理できる。  A _ 0 , A _ 1 , \dots , A _ i まで処理したときの  \left ( i _ k , A _ { i _ k } \right ) を全て保持しておけば、オンラインクエリにも対応可能。

実装例

木上のパス xor 基底クエリ

静的な木上のパス xor 基底クエリが  \Theta \left ( N \log \left ( \max A \right ) + Q \left ( \log \left ( \max A \right ) \right ) ^ 2 \right ) で処理できる。 自由に根を決め、根から各頂点までのパスについての  \left ( i _ k , A _ { i _ k } \right ) を計算しておく。 これは各  v について、根から  v の親へのパスに対する  \left ( i _ k , A _ { i _ k } \right )  a _ v を push すればいいので  \Theta \left ( N \log \left ( \max A \right ) \right ) で計算できる。  uv パスに対するクエリにおいては、 \operatorname { LCA } \left ( u , v \right ) から  u , v それぞれへのパスの xor 基底をマージすればよい。 それぞれは根から  u , v へのパスの suffix であるから本稿のアルゴリズムで求めることができる。

実装例 (LCA の計算に  \Theta \left ( \log \left ( N \right ) \right ) アルゴリズムを使っているので計算量は少し増えている)

拡張

アルゴリズムを観察すると、要素の追加が末尾である必要はないことが分かる。 すなわち、最初  S = \emptyset として、以下のクエリを処理できる。

  •  j , x が与えられ、 S \leftarrow S \cup \left \lbrace \left ( i , x \right ) \right \rbrace
  •  i が与えられ、 \left \lbrace x \mathrel { } \middle \vert \mathrel { } \exists \left ( j , x \right ) \in S . ~ j \geq i \right \rbrace の xor に関する基底  \left ( b _ 0 , b _ 1 , \dots \right ) を出力する。
    • ただし  \operatorname { bsr } \left ( b _ k \right )  k について単調減少である。

木上の Disjoint Sparse Table

概要

 n 頂点の木の各頂点にモノイド  M の元  a _ v が乗っているとして、指定されたパス上の値の総積を計算するクエリが、前計算  \Theta \left ( n \log \left ( n \right ) \right ) 、クエリ  \Theta \left ( 1 \right ) で処理できる。

アルゴリズム

重心分解を行う。

木の重心を  c とし、各  v \in V に対して、 c v パスから  c を除いたものに対するクエリ結果を計算しておく *1 u v パスが  c を含む場合、 u v パスに対するクエリはこの前計算を用いて  \Theta \left ( 1 \right ) で計算できる。  c を含まない場合に対処するため、木から  c を取り除くことで得られる部分木たちに対し再帰的に同様の前計算を行う。

クエリとして  u , v が与えられた際に適切な  c を発見するためには、重心分解によって得られた木構造において LCA を計算すればよい。 (Disjoint) Sparse Table を用いた、前計算  \Theta \left ( n \log \left ( n \right ) \right ) 、クエリ  \Theta \left ( 1 \right ) LCA があるためそれを用いる。

良い性質

一般の  x , y \in M に対する  x y の計算はクエリの際に  1 度行われるだけで、前計算はすべて  x \in M, v \in V に対する  x a _ v , a _ v x の計算を用いている。 よって、この種の計算が通常より高速に行える状況では計算量が改善する。 例: 各頂点に非負整数が定められており、パス上の非負整数に対する xor 基底をクエリで求めたい場合。

これは Disjoint Sparse Table も持つ性質である。 そもそも Disjoint Sparse Table はこのデータ構造をパスグラフに対して適用したものと一致する。 したがってこのデータ構造は Disjoint Sparse Table の一般化であると見ることもできる。

LCA を計算する部分の対案

本稿のデータ構造でクエリを行う際に  c を計算する部分は、Disjoint Sparse Table においては整数の xor と bsr で行っていた。 1/3 重心分解などを用いて LCA を計算したい木を二分木にし、根から  v までのパスの左右を 01 でエンコードすれば、同様に xor と bsr で LCA を計算できる。 正確には LCA ではなく LCA の深さが求まるが、そもそも前計算のデータを (頂点, 深さ) に対してデータを対応させるような形式で管理しておけばよい。

*1: モノイドが非可換なら双方向の向きについて計算しておく

ICPC Yokohama 2023 参加記

リハーサル

MLE を見ようとして大量にメモリを取ったら PC が落ちたりして面白かった。

本番

ある程度関わった問題を何となく時系列っぽく並べる

問題文と解説: Regional | ICPC 2023 Asia Yokohama Regional

B

親の顔より見た平均の式変形。 c の小ささを利用しようかと考えたが数秒で思いつかなかったので却下した。

D

ほとんどやるだけ。連続部分列の高速な一致判定が欲しくなったので  \Theta(n  ^ 2) で LCP のテーブルを作っておいた。

H

potato と一緒に考察する。 しばらく悩んでいたが、potato が「最小カットじゃないか?」と発言したのでそこから数秒で解けた。 最小カットで非自明な問題を久々に見た気がするが、すぐに解けるべきだった。

K

線分の端点が正方形の内側にあるようなクエリのことが脳から抜け落ちていて、 x 軸と  y 軸のそれぞれに対して座標を特定しようとするとクエリがギリギリ溢れるななどと考えていた。 その結果アルゴリズムも複雑になり、実装に手間取ったうえにバグらせてしまった。

J

最初の方で potato から「min plus の畳み込みは高速になるか?」と聞かれ、いかにも担当っぽい問題だったので貰う。 木上の最小費用流になることはすぐに分かり既知のアルゴリズムを適用しようとしたのだが、頂点の湧き出しに対するコストが  f m ^ 2 なので辺の本数が  2 乗になってしまう。 こういうときはアルゴリズムの動作を観察して適切に高速化すると (例えば等差数列加算とか) 上手く行くことが多いのだが、今回は等差数列の insert だったので解けそうで解けない。 得意分野のはずという気持ちでかなり時間を消費してしまったのだが、結局最後まで残ってしまう。 最小費用流なので多分何かしらの貪欲だろうということで potato Hemimor に任せているとそれっぽい貪欲が提案されたので実装してもらい、通った。

結局、等差数列の inset を愚直にやっても上手く計算量が抑えられていたらしい。 この分野に関連するアルゴリズムはかなり触っており自信があったので、ショックだった。

C

とりあえずポリアはする。

回転を同一視しない数え上げを考えると、カタラン数を少し弄ったようなものが出てくる。 ラグランジュ反転の思い出しに自身が無かったのと、最後の方まで畳み込みなしで押せる自信が無かったので FPS で書いておく。

周期が決まっている場合の数え上げだが、当初スタックに色を積んでいく方針で考えていたので  1 周期で全部消えるケースしかないんじゃないかと予想した。 この方針で実装するとサンプル  1 が合わないので見てみると、消えないケースがある。 サンプルすら合わない解法を実装し始めてしまうのは悪い癖なのだが直りそうもない。 しばらく考えると、隣り合う同じ色を消す方針を思いつき、奇数長回文が残る場合が必要十分と分かった。 これも FPS で立式できるので書いて出してみると TLE。

ここからの高速化の作業は悪夢だった。 まず FPS の立式を偶奇で分けることで長さを半分にする (TLE)。 FPS の inv を定数倍が良好な実装に変える (TLE)。 同じ列の DFT を繰り返し計算している部分を使いまわすようにする (TLE)。 その他あり得そうな部分をひたすら高速化して回る (TLE)。

結局この方針は十分に高速化された FFT でギリギリ通せるラインのものだったようだ。 SpeedStar にはどのみち負けていたのだが、食らいつけなくて残念。

総評

不調気味だった気がするが  3 位には落ちなかったし、 1 位は好調でも難しかっただろうと思うので、全完できなかったことだけが心残りになった。 海外に行くのが嫌いなので playoff 確定でげんなりしているが頑張りたい。

懇親会

例の大きい名札を付けてうろうろし、色々な人に話しかけてもらった。 初めて会う人もたくさんいたので嬉しい。

懇親会でアルゴリズムの話をする難しさをこの手のイベントの度に感じていて、「最近良い感じのアルゴリズムありましたか?」から切り出せる相手があまり多くない。 次の機会があれば、hotman さんに倣って懇親会で話せるアルゴリズムのお品書きを作っておくべきかもしれない。 fenwick tree の最悪計算量を半分にする方法や、競プロで使えそうにない最近のフローアルゴリズム群などについて話した。

終わりに

準備に携わったスタッフの皆様、ありがとうございました。

Aliens DP における二分探索の色々

とにかく使いたい人向け

スコアとペナルティがどちらも整数とする。 *1

ペナルティを  p にしたときの最適値を  g(p) と置く。

事実  1

ペナルティが  p の時の最適解がペナルティを受ける回数としてあり得るものは、区間  \lbrack l _ p , r _ p \rbrack になっている。

事実  2

最小化 (ペナルティの度にスコアに  + p ) なら、 r _ {p + 1} = l _ p = g(p + 1) - g(p) である。

最大化 (ペナルティの度にスコアに  - p ) なら、 r _ {p + 1} = l _ p = g(p) - g(p + 1) である。

使い方

 g(p) の計算の際にペナルティ回数の復元は不要で、 g(p + 1) - g(p) (最大化なら逆) を返せばいい。

詳しめの解説

概要

Aliens DP などで二部探索をするときに用いる複数の手法を紹介する。

手法 利点 欠点
ペナルティ回数の管理 分かりやすい。ペナルティ回数を求める際に計算時間がほとんど増えない問題では最速。 ペナルティ回数の管理が面倒。tie-breaking に気を使わないと壊れる。
傾き二分探索 回数の管理が不要な中では最も分かりやすく、三分探索より速い。 最速ではない。
三分探索 凸じゃないケースで誤魔化せる可能性が傾き二分探索より高いかもしれない。 傾き二分探索のほぼ下位互換。
黄金分割探索 回数の管理で計算時間が嵩む場合最速。 ほぼライブラリ化前提、その場で書くのは苦しい。

問題設定

下に凸な関数  f \colon \mathbb{Z} \to \mathbb{Z} \cup \lbrace + \infty \rbrace, ~ f \neq + \infty があり、 g \colon \mathbb{Z} \to \mathbb{Z} \cup \lbrace - \infty \rbrace, ~ p \mapsto \inf \limits _ x \left ( f(x) - px \right ) と定義する *2 g(p) は、 y = f (x) のグラフに傾き  p の接線を引いた際の  y 切片の値と解釈することもできる。 我々の目標は、任意の  p \in \mathbb{Z} に対して  g \left ( p \right ) が高速に計算できるとき、与えられた  \tilde x \in \mathbb{Z} に対する  f \left ( \tilde x \right ) を計算することである。

ペナルティ回数の管理

 \inf \limits _ x \left ( f(x) - px \right ) を達成する最小の  x x ^ {\ast} _ p とする。  g \left ( p \right ) の評価の際に  x ^ {\ast} _ p を同時に計算できる場合を考える。

 x ^ {\ast} _ p p に対して単調増加であるから、二分探索を用いて  x ^ {\ast} _ {\tilde p} \leq \tilde x \lt x ^ {\ast} _ {\tilde p + 1} を満たす  \tilde p を発見できる。 すると  \tilde x \in \mathop{ \mathrm{arg\,min} } \limits _ x \left ( f(x) - \tilde p x \right ) であるから、 g \left ( \tilde p \right ) = f \left ( \tilde x \right ) - \tilde p \tilde x から  f \left ( \tilde x \right ) の値が分かる。

 x ^ {\ast} _ p の定義を  \inf を達成する最大の  x とした場合も、 x ^ {\ast} _ {\tilde p - 1} \lt \tilde x \leq x ^ {\ast} _ {\tilde p} を満たす  \tilde p を発見すれば同様である。 一方で、 \min を達成する  x を何か  1 つ計算できるとした場合、この方法では  f \left ( \tilde x \right ) を特定できない *3

傾き二分探索

 \mathop{ \mathrm{arg\,min} } \limits _ x \left ( f(x) - p x \right )  \mathop{ \mathrm{arg\,min} } \limits _ x \left ( f(x) - \left ( p + 1 \right ) x \right ) の共通部分は  1 元集合であり、 \mathop{ \mathrm{arg\,min} } \limits _ x \left ( f(x) - \left ( p + 0.5 \right ) x \right ) と一致する。 証明は容易なので省略する。

上記の事実から、 \left \lbrace x ^ {\ast} _ {p + 0.5} \right \rbrace = \mathop{ \mathrm{arg\,min} } \limits _ x \left ( f(x) - \left ( p + 0.5 \right ) x \right ) と定義すると以下の式が成り立つ。

\begin{align} g ( p ) &= f \left ( x ^ {\ast} _ {p + 0.5} \right ) - p x ^ {\ast} _ {p + 0.5} \cr g ( p + 1 ) &= f \left ( x ^ {\ast} _ {p + 0.5} \right ) - \left ( p + 1 \right ) x ^ {\ast} _ {p + 0.5} \end{align}

ここから  x ^ {\ast} _ {p + 0.5} = g ( p ) - g ( p + 1 ) が分かるため、 g の評価ができるだけで  x ^ {\ast} _ {p + 0.5} を得ることができる。 後は二分探索を用いて  x ^ {\ast} _ {\tilde p - 0.5} \leq \tilde x \leq x ^ {\ast} _ {\tilde p + 0.5} を満たす  \tilde p を発見すれば、ペナルティ回数の管理の節と同様である。

本節の内容は  \mathop{ \mathrm{arg\,min} } \limits _ x \left ( f(x) - p x \right ) = \left \lbrack x ^ {\ast} _ {p - 0.5} , x ^ {\ast} _ {p + 0.5} \right \rbrack を意識すると理解しやすいかも知れない。

また、ペナルティ回数の管理を用いて  x ^ {\ast} _ {p + 0.5} を直接求めることで本節の内容を実行しても良い。 その場合、 \mathop{ \mathrm{arg\,min} } \limits _ x \left ( f(x) - \left ( p + 0.5 \right ) x \right ) = \mathop{ \mathrm{arg\,min} } \limits _ x \left ( 2 f(x) - \left ( 2 p + 1 \right ) x \right ) を用いると整数の範囲で計算できる。

三分探索

 g の定義より、任意の  p について  g(p) \leq f \left ( \tilde x \right ) - p \tilde x である。 これを移行して  f \left ( \tilde x \right ) \geq g(p) + \tilde x p としたのち両辺を  p についての  \sup をとると、 f \left ( \tilde x \right ) \geq \sup \limits _ p \left ( g(p) + \tilde x p \right ) を得る。  f \left ( \tilde x \right ) \lt +\infty として、 f の凸性から  g(p) = f \left ( \tilde x \right ) - p \tilde x を満たす  p が存在するため前述の式は等号で成立し、

\begin{align} f \left ( \tilde x \right ) = \sup \limits _ p \left ( g(p) + \tilde x p \right ) \end{align}

である。  g は線形関数の  \inf で定義されるから上に凸であり、従って  p \mapsto g (p) + \tilde x p も上に凸であるから、三分探索を用いて最大化すればよい。

傾き二分探索の節の内容は、 p \mapsto g (p) + \tilde x p の最大化を傾きの二分探索で行っているのと本質的に同じである。 よって、傾き二分探索の節の内容を押し進めて三分探索に行き着くことも可能である。

黄金分割探索

三分探索の節の内容をそのまま黄金分割探索にする。

問題設定の変種

 f \colon \mathbb{R} \to \mathbb{R} \cup \lbrace + \infty \rbrace が下に凸かつ狭義単調増加で、 g (p) = \inf \limits _ x \left ( f(x) - px \right ) が高速に計算できるとする。 与えられた  k に対し、 f \left ( \tilde x \right ) = k を満たす  \tilde x を求めよ。

 f が離散ではなく連続になっていることに注意。 これは  f \left ( \tilde x \right ) = k となる  \tilde x が存在した方が議論が簡潔になるためであり、離散の設定なら線形に連続に拡張してしまえば良い。

ペナルティ回数の管理, 傾き二分探索

 x ^ {\ast} _ {p} が分かれば  g (p) = f \left ( x ^ {\ast} _ p \right ) - p x ^ {\ast} _ p から  f \left ( x ^ {\ast} _ p \right ) が分かるため、二分探索を行うことができる。

三分探索, 黄金分割探索

まず、 f を平行移動することで  k = 0 とする。 すると、 f(x) のグラフと  x 軸との交点の  x 座標を求める問題となる。

 f x 軸との交点における傾き  \tilde p を求めることを考える。  h \colon \mathbb{R} _ {\gt 0} \to \mathbb{R} \cup \lbrace + \infty \rbrace h(p) = - \dfrac{g (p) }{p} で定める。  h(p) は、 f (x) に引いた傾き  p の接線の  x 切片である。  f が単調かつ凸であるから、 h は単峰で  \tilde p において最小値を取り、三分探索を適用することで求まる。 このとき  \tilde x = h \left ( \tilde p \right ) である。

正確には、 f x 軸との交点において微分可能とは限らないため、 \tilde p は任意の劣勾配ということになる。

注意

  • 二分探索や三分探索の初期値は、問題毎に考えて行うしかない。
  •  g(p) = - \infty の場合  x ^ {\ast} _ p などが定まらないが、これも問題毎に処理する。ただし  f の定義域が有界区間であることが多く、その場合  g は常に有限になる。また、 g(p) = - \infty となるのは  p が大きすぎるか小さすぎるかのときなので、 x ^ {\ast} _ p を適切に  \pm \infty にすることでも回避できる。

余談

本稿で定義した  g の符号を反転した関数は  f の凸共役と呼ばれ、 f ^ {\bullet} で表される。 凸共役に関する各種の定理を認めると、本稿の内容は大幅に簡略化される。

関連する記事

Aliens DP で辺の本数が区間で指定される場合 - noshi91のメモ

*1: スコアが実数ならペナルティも実数にしないと駄目で、tie-breaking は自由にやってよい

*2: p がペナルティというよりはボーナスになっているが、argmin が単調減少ではなく単調増加の方が見た目が分かりやすいためである。

*3: 正確には、次節の内容と組み合わせると特定できる

省メモリな区間 add 区間 min

概要

長さ  n の数列  a に対する区間 add と区間 min を時間計算量  \Theta ( \log ( n ) ) で処理する、空間計算量  cn のデータ構造を説明する。 ただし、一連の操作における  \lvert a _ i \rvert の最大値を  X としたとき、 \lbrack -2X , 2X \rbrack の範囲の数を管理するために必要な空間計算量を  c としている。

例: 長さ  10 ^ 5 の数列に  \pm 10 ^ 9 の加算クエリが  10 ^ 6 回ほど行われる状況では、 64 bit 整数  10 ^ 5 個のメモリで処理できる。

構造

簡略化のため、 n 2 の冪の場合をまず考える。 Segment Tree と同様に、配列を用いて二分木の構造を作る。 まず、各ノードに対して対応する区間内の最小値を考える (この値は実際には保持しない)。 続いて、各ノードに対して左子  - 右子を計算し、これを保持する。 葉には子が無いので、データは保持しない。 すると下図のような構造となる。

この右側のデータに加えて全体の最小値も持つことで、Segment Tree に対するものと似たようにクエリの処理が可能となる。 例えば区間 add では、データを適切に更新しつつ、着目しているノードの区間の最小値の増分を返すような関数を設計すると良い。 区間  \lbrack l , r )  x を加算するクエリは以下のようになる

  1. 着目しているノードの区間 \lbrack l , r ) に完全に含まれているなら、 x を返す。
  2. 着目しているノードの区間 \lbrack l , r ) と交差しないなら、 0 を返す。
  3. 左子と右子に対して再帰的に加算クエリを行う。それぞれの返り値から着目しているノードの値を更新し、返り値を計算する。

 n 2 の冪でない場合も、空間計算量が  2n の Segment Tree と同様に完全二分木でない二分木構造を作ればよい。 ノードの個数は  n - 1 になるため、全体の最小値の管理と併せて  n 個の整数を管理することになる。