noshi91のメモ

データ構造のある風景

クエリが整数の Convex Hull Trick の 凸 判定

概要

Convex Hull Trick で最小値取得クエリの  x が整数の時、アルゴリズム内で直線が不要かどうかチェックする部分を簡単に記述することが出来る。 特に、傾きと切片を掛けることでオーバーフローする問題を発生させない。

本文

Convex Hull Trick で、最小値取得クエリの  x の値の全てが整数とします。 簡単な下処理により、直線の傾きは distinct と仮定してよいです。  3 本の直線  \lbrace \text{line}\ i = a _ i x + b _ i \ | \  i = 0, 1, 2 \rbrace があり、 a _ 0 > a _ 1 > a _ 2 とします。 ここで、 a _ 1 x + b _ 1 が不要かどうか判定することを考えます。  3 本の直線が下図のような位置関係にあるとしましょう。 座標軸の目盛りは整数の座標を表しています。

f:id:noshi91:20210323190912j:plain

Convex Hull の字面に従うならば、 \text{line} \ 1 は最小値を取り得る、すなわち必要な直線です。 しかし、クエリが整数しかないという仮定の下、 \text{line} 1 は最小値を取り得ません。 この条件を正確に記述します。

\begin{align} f(ax+b, cx+d) := \max \lbrace k \ | \ ak+b \leq ck+d \rbrace \end{align}

と定義すれば、

\begin{align} \text{line}\ 1 \text{が最小値を取り得る} \Leftrightarrow f(\text{line}\ 0, \text{line}\ 1) < f(\text{line}\ 1, \text{line}\ 2) \end{align}

です。 ただし、 2 直線が交わる  x 座標では傾きの大きい方が最小値を取るとしています。

ところで  a > c の仮定の下

\begin{align} f(ax+b, cx+d) :&= \max \lbrace k \ | \ ak+b \leq ck+d \rbrace \\ &= \left \lfloor \frac {d - b} {a - c} \right \rfloor \end{align}

であるため、 f は整数除算で計算することが出来ます。 これにより、分母を払って判定するアルゴリズムにおけるオーバーフローの問題が発生しなくなります。 また、この条件はクエリが実数の場合の条件より強い条件なので、残された直線群の 凸 性が保たれ、最小値取得クエリのアルゴリズムは普通の物を使うことが出来ます。

{x | b^x = a (mod m)} の構造

概要

 a, b, m が与えられて  \lbrace x \mid b^x \equiv a \pmod m \rbrace の構造を観察します。

結論

 S := \lbrace x \mid b^x \equiv a \pmod m \rbrace について、次のうちいずれかが成り立ちます。

  •  S = \emptyset
  •  \exists\ 0 \le c \lt \log_2(m),\ S = \lbrace c \rbrace
  •  \exists\ n \mid \lambda(m),\ \exists\ t,\ S = \lbrace x \mid x \equiv t \pmod n  \rbrace
  •  \exists\ 0 \le c \lt \log_2(m),\ \exists\ n \mid \lambda(m),\ \exists\ t,\ S = \lbrace x \mid x \gt c \land x \equiv t \pmod n  \rbrace

証明

 m = \prod_{i=1}^{k} p_i ^ {e_i} とすると、中国剰余定理から  \mathbb{Z} / m\mathbb{Z} \simeq \prod_{i=1}^{k} \mathbb{Z} / p_i ^ {e_i} \mathbb{Z} i を 1 つ固定して考えると、

  •  p \mid b, a \neq 0 \pmod {p^e} のとき、 x は存在しないか、一意に定まる。
  •  p \mid b, a = 0 \pmod {p^e} のとき、 x \gt c の形の条件が付く。
  •  p \nmid b のとき、 x は存在しないか、 x = t \pmod {n} の形の条件が付く (ただし  n \mid \lambda(p^e))。

これらすべての条件を合わせた条件が、挙げた 4 つのいずれかに該当することは確認できる。

 \lambda : Carmichael function - Wikipedia

割れない時の行列式

概要

逆元が必ずしもない状況で行列式を計算する方法を紹介します

一般の可換環上の行列式

n91lib_rs/division_free_determinant.rs at master · noshi91/n91lib_rs · GitHub

除算を使わない  \Theta(n^4) のシンプルなアルゴリズムが知られています。 除算無しでも  o(n^3) になるらしいのですが、詳しく勉強してないので詳細は分かりません。

 {\rm mod} 合成数

掃き出し法のアルゴリズムを少し変えるだけで高速に計算できます。 掃き出し法は「ある行を使って別の行の先頭要素を 0 にする」を繰り返すアルゴリズムなので、この部分を何とかできればいいです。 そこで、0 にしたい列の要素を一旦 ( {\rm mod}\ m ではなく) ただの数だと思って、ユークリッド互除法を行います。 ユークリッド互除法は定数倍と加減算なので、 {\rm mod}\  m でも問題なく計算できます。 これを繰り返せば  O(n^3 \log(m))アルゴリズムが得られます。  {\rm gcd} を繰り返しとる操作は計算量が落ちるので、実際には  O(n^3 + n^2 \log(m)) になっています。

さらに計算量を落とすことも出来ます。 そもそもユークリッド互除法の部分は 2 要素間だけで計算できるので、最終的にどういう線形変換になったかだけを追って、最後に行ベクトル同士を足し引きすればよいです。 これにより  O(n^3 + n \log(m)) となります。

ユークリッド互除法が行えるなら何でもいいので、多項式環などでも使うことが出来ます。




----------- 以下、yukicoder のネタバレを含みます ----------



















yukicoder No.1303 Inconvenient Kingdom - ひとなので

この解法で掃き出し法が機能しているのは、元のグラフが連結なので行列式が 0 でない、などの性質を利用しています。 先ほども述べた通り多項式ユークリッド互除法が使えるので、それを用いれば ad-hoc な証明をしなくても AC を取ることが出来ます。

https://yukicoder.me/submissions/586286

余談ですが、私は本番でこれに気付かなかったので、一般の可換環用の  \Theta(n^4)アルゴリズムの方を使って AC しました。

ICPC2020 国内予選 参加記

Poyashi で参加して 6 完 4 位でした

開始前

リハーサル参加のためかなり早めに来たが、模擬に参加してるのでほとんど練習の必要がないことに気付く。 ライブラリを書きながら開始を待つ。

コンテスト開始後

記憶がかなり曖昧なので間違ってるかも。

A moyashi B pulmn C noshi で読む。 AB はすぐに解けるが C が普通に解けないので D を読んだら D も解けなくて困った。

moyashi pulmn に C を伝えるが解法が出ないため、ぶん回し覚悟で重いコードを書く。 サンプルですら異常な時間が掛かっており不安になるが、とりあえず裏で回しておくことにした。(後で通った)

D が分かったので実装する。この間に pulmn が F を解いていた(問題見てなかった)。

pulmn「E は 230

noshi「100 時間くらいかかる計算」

pulmn「fibonacci になった」

noshi「215 になった」

moyashi「そうじゃん、割り算間違えてた」

3 人で E やってもしょうがないだろということで moyashi と pulmn が同時に E をコーディングして、どちらかで通ればいいですねをした。

その間に noshi が G を考えていた。

E が通ったが moyashi pulmn は G が苦手らしいので H をやってもらう。

こんなのは最小費用流の双対なんですねと言いながら双対を取ったら、解が一致しないので苦戦する。 実は weight と capacity を swap してしまっていたことが分かり、双対は正しく取れていた(は?)。 結局方針があってるのかもよく分からなかった。

終了

感想

2019 予選, 2020 模擬で大敗していたので、良い感じの結果が出て満足した。 序盤の問題が通り終わったら自分は全ての問題に目を通すべきだったというのと、後半時間が余りまくったのでもう少し考察を速めたい。

私用メモ: 畳み込めるものまとめ

長さ  n 同士の列を畳み込むのが  \mathrm{O}(n ^ {2-\epsilon}) でできるやつ

 \displaystyle c _ k = \bigoplus _ {i \circ j = k} a _ i \otimes b _ j

 \circ  \oplus  \otimes 時間計算量 備考
 +  +  \times  \mathrm{O}(n \log (n))  \mathbb{C}, 高速フーリエ変換
 +  +  \times  \mathrm{O}(n \log(n) \log(\log(n))) 環, Schönhage–Strassen
 \cap  +  \times  \mathrm{O}(n \log (n)) 環, 高速ゼータ変換, メビウス変換
 \cup  +  \times  \mathrm{O}(n \log (n)) 環, 高速ゼータ変換, メビウス変換
 \Delta  +  \times  \mathrm{O}(n \log (n)) 環, 2 で割れる, 高速アダマール変換
 \amalg  +  \times  \mathrm{O}(n (\log (n)) ^ 2) 環, subset convolution
 \mathrm{gcd}  +  \times  \mathrm{O}(n \log (\log (n))) 環, 高速ゼータ変換, メビウス変換
 \mathrm{lcm}  +  \times  \mathrm{O}(n \log (\log (n))) 環, 高速ゼータ変換, メビウス変換
 \times \pmod{n}  +  \times  \mathrm{O}(n \log (n))  n素数,  \mathbb{C}, 高速フーリエ変換
 +  \max  \min  \mathrm{O}(n \sqrt{n \log (n)})
 +  \min  +  \mathrm{O}(n)  b が下に凸
 +  \min  +  \mathrm{O}(n \alpha(n))  b が上に凸
 +  \times  \mathrm{O}(n (\log(n)) ^ 2)  \mathbb{C}, ※多変数畳み込み
 +  \times  \mathrm{O}(n (\log(n)) ^ 2 \log(\log(n))) 環, ※多変数畳み込み

 \circ が集合の演算になっているものは、添え字を集合の整数表現と解釈して演算することを意味する。


実装など

n91lib_rs/subset_convolution.rs at master · noshi91/n91lib_rs · GitHub

n91lib_rs/max_min_convolution.rs at master · noshi91/n91lib_rs · GitHub

Concave Max Plus Convolution | Library

多変数畳み込み(切り捨て)のアルゴリズム – 37zigenのHP

Range Mode Query 空間 Θ(n) 構築 Θ(n√n) クエリ Θ(√n)

概要

Range Mode Query - data-structures

これをもう少しだけ詳しめに解説します。 面白いので好きなアルゴリズムです。

アルゴリズム

まずは平方分割して、ブロックの境界を端点とする区間 ( \Theta(n) 個ある) の全てについて最頻値とその頻度を計算します。 他にも前計算するものがありますが、説明が難しいのでクエリの処理方法を追いながらやっていきます。

区間  \lbrack l, r) の最頻値を計算したいとします。 まずは前計算したものを使って、暫定的な最頻値と頻度 (  =: freq ) を取得します。 この暫定解が更新されるとすれば、それはブロックからはみ出た左右の端数の部分に含まれるいずれかの値によるものになります。 従って、端数の部分の全要素について、その要素が  \lbrack l, r) に出現する回数が  freq を超えるか判定して、適切に更新する、というアルゴリズムがまず思い浮かびます。

f:id:noshi91:20201026134939j:plain

クエリを  {\rm O}(\sqrt n) にしたいので、それぞれの要素について  {\rm O}(1) で頻度を取得したくなりますが、 {\rm O}(n) の空間計算量でこれは絶望的です。

そこで見る方向を変えて、「左側の端数部分に含まれる各要素  x について、その  freq 個後ろの  x の位置が  r より手前にあるか」という判定を考えます。 この判定自体は、全ての数についてその数だけ抜き出した index の列を持っておけば空間計算量  \Theta(n) 、時間計算量  \Theta(1) で可能です。 これを用いた「 r より手前にあったら  freq が更新されるということなので、 freq を 1 増やして繰り返す」というアルゴリズムは、 \Theta(\sqrt n) の時間計算量になります。 なぜなら、ステップを 1 つ進めるたびに「 x では更新されなくなったので、次の数に移る」「 freq が 1 増える」のいずれかが発生し、 freq は最初の暫定解から  \Theta(\sqrt n) しか増え得ないからです。 右側の端数部分に含まれる数については逆に、 freq 個手前の  x の index が  l より後ろにあるかを判定すれば全く同様です。