noshi91のメモ

データ構造のある風景

競プロにおける約数の個数の見積もり

概要

 n \leq N を満たす  n の正整数の約数の個数の最大値は  N ^ {1/ 3} 程度と思っておけば競プロの計算量の見積もりで便利そう。

実験

 d(n) n の正の約数の個数として、 \displaystyle r = \frac{\max _ {1 \leq n \leq N} d(n)}{N ^ {1 / 3}} と定義し、 N r の関係を見る。

つまり、 N ^ {1 / 3} と見積もると実際 *1 にはその  0.1 倍から  3.5 倍程度に収まることになる。

例えば、 n \leq 10 ^ {12} なら  n の約数の個数の  2 乗が  10 ^ {8} より小さいくらいになって  1 sec に間に合いそうと分かる。

*1:競プロだと  N \leq 10 ^ {18} がほとんどだろうという仮定の下

アフィン部分空間の共通部分の計算

2022/05/02 記事中の記号の使い方やその他の軽微な修正をしました。

概要

 \mathbb{F} _ p ^ d のアフィン部分空間が  \mathrm{span}\lbrace a _ 1, a _ 2, \dots, a _ n \rbrace + b の形で与えられるとき、 2 つのアフィン部分空間の共通部分を時間計算量  \Theta( (n + d) d ^ 2 + d \log(p) ) で計算する。

アフィン部分空間

ベクトル空間  V の部分集合  S は、部分空間  U \subseteq V b \in V によって  S = \lbrace u + b \mid u \in U \rbrace と表現できるとき、アフィン部分空間であると言う。 アフィン部分空間は部分空間の一般化になっており、後述するが競プロでも出現する概念である。  S に対して  U は一意だが  b は一般には一意でない。

アフィン部分空間の  2 通りの表現

 \mathbb{K} を体として、 S \mathbb{K} ^ d のアフィン部分空間とする。 定義から明らかに、 A \in \mathbb{K} ^ {d \times n}, b \in \mathbb{K} ^ d によって  S = \lbrace Ax + b \mid x \in \mathbb{K} ^ n \rbrace と表現できる。 ここで  n S = \lbrace u + b \mid u \in U \rbrace としたときの  U \subseteq \mathbb{K} ^ d の次元である。

(ドット積についての) 直交補空間を考えると、 \lbrace Ax \mid x \in \mathbb{K} ^ n \rbrace = \lbrace y \mid y \in \mathbb{K} ^ d , Cy = 0 \rbrace となる  C \in \mathbb{K} ^ {(d - n) \times d} が存在する。 すると、 S = \lbrace y + b \mid y \in \mathbb{K} ^ d, Cy = 0 \rbrace = \lbrace y + b \mid (y + b) \in \mathbb{K} ^ d, C(y + b) = Cb \rbrace であり、諸々の記号を置きなおすと結局  S = \lbrace y \in \mathbb{K} ^ d \mid Cy = d \rbrace という表現を得る*1。 逆に、 \lbrace y \in \mathbb{K} ^ d \mid Cy = d \rbrace と表現される集合は、空であるかアフィン部分空間であることもすぐに確かめられる。

アフィン部分空間の共通部分の計算

 S _ 1, S _ 2 をアフィン部分空間として、 S _ 1 = \lbrace y \in \mathbb{K} ^ d \mid C _ 1 y = d _ 1 \rbrace, S _ 2 = \lbrace y \in \mathbb{K} ^ d \mid C _ 2 y = d _ 2 \rbrace と表現できるとする。 すると  \displaystyle S _ 1 \cap S _ 2 = \left \lbrace y \in \mathbb{K} ^ d ~ \middle | \pmatrix{C _ 1 \cr C _ 2} y = \pmatrix{d _ 1 \cr d _ 2} \right \rbrace であり、容易に共通部分の表現が得られる。 この式から分かるように、これらのアフィン部分空間の共通部分は空であるかアフィン部分空間となる。  \lbrace Ax + b \mid x \in \mathbb{K} ^ n \rbrace \lbrace y \in \mathbb{K} ^ d \mid Cy = d \rbrace の表現の相互の変換は、直交補空間の計算と同様に行うことができる (参考: Parity-check matrix - Wikipedia)。

実用例

注意: cp-unspoiler この問題のネタバレを含みます。





















G - Xor Cards

 C _ i を、 A _ i, B _ i \mathbb {F} _ 2 ^ {30} の元として見てこの順に連結したベクトルと定義する。 問題は以下のように整理される。

長さ  N の非負整数列  (C _ i) と、非負整数  K が与えられる。  C から  1 つ以上選んで xor を計算し、前半  30 bit が  K 以下になるようにしつつ、後半  30 bit を最大化せよ。

 K 以下の非負整数全体は、「下位  t bit は自由で、残りは決まっている」という形の集合の高々  \log _ 2(K) + 1 個の直和に分解される。 この分解のそれぞれについて独立に計算し、最後に最大値を取ればよい。

 1 つ以上選ぶという条件を一旦  0 個以上に緩和すると、「 C から  0 個以上選んで xor を取ることで得られる数全体」と「下位  t bit は自由で、残りは決まっている数全体」はいずれも ( \mathbb{F} _ 2 ^ {60} の) アフィン部分空間になっている。 したがって、これらの共通部分  S を計算した後、 S の元で後半  30 bit が最大となるものを計算すればよい。 実際には、 C _ i の定義を  B _ i, A _ i の順に連結することにして、単に  \max S を求めるのが簡単だろう。  \max S = 0 かつ  (C _ i) が線形独立である場合に限り最大値を  -1 に修正することで、 1 つ以上選ぶ条件をカバーする。

 \max \lbrace Ax + b \mid x \in V \rbrace は、 A の列ベクトルを掃き出して簡約化し、 b に対して基底を貪欲に足していけばよい。

実装例

簡略化のため、 K \leftarrow K + 1 として以下を未満に言い換えている。

Submission #31329771 - Monoxer Programming Contest 2022(AtCoder Beginner Contest 249)

*1: Cy + d = 0 の形式とどちらが綺麗なのかは分からない。

Xorshift を用いた Zobrist hashing を衝突させる方法

概要

Xorshift によって生成された乱数列について、いくつかの場所の xor を取ると seed に依存せずに  0 になる。

Xorshift

Xorshift - Wikipedia

説明は Wikipedia に丸投げする。 本記事では簡略化のため、実装例の一番上にある周期  2 ^ {32} - 1 の実装について議論する。

Zobrist hashing

 N 個の要素があるとして、要素  i に乱数  r _ i \in \mathbb{Z} _ {\geq 0} を割り当てる。 このとき、集合  S の hash  h(S) \displaystyle \bigoplus _ {i \in S} r _ i で定める。 ただし  \oplus は bitwise xor とする。 この hashing の方法を Zobrist hashing と呼ぶ。

hash の衝突

異なる集合  T _ 1, T _ 2 について  h(T _ 1) = h(T _ 2) となることを衝突と呼ぶ。 Zobrist hashing の場合  h(T _ 1) = h(T _ 2) \Leftrightarrow h(T _ 1 \mathbin{\triangle} T _ 2) = 0 であるから、 S \neq \emptyset かつ  h(S) = 0 となるような  S が分かれば衝突を引き起こすことができる。

結論から言うと、 r _ i を Xorshift が  i 番目に出力した値として設定した場合、seed に依存せず  h(S) = 0 となるような  S \subseteq \lbrace 0, 1, \dots, 32 \rbrace が存在する。 Xorshift の乱数列は途中から切り出しても何らかの seed を与えた Xorshift の乱数列であるから、 r _ i の設定が行われるまでに何度か乱数生成されていても衝突することになる。 思いつく回避方法を挙げる

  •  r _ i の設定を行う際に 1 つ飛ばしや逆順にするなどの *1 変化を加える。ソースコードを読まれた場合対応される可能性がある。
  • Xorshift における shift の幅などを変える。ソースコードを読まれると対応される上、パラメータを間違えると乱数の質が落ちる。
  • Xorshift 以外の乱数生成アルゴリズムを用いる。

 S の存在の証明、及び構成方法

32 bit 整数をベクトル空間  V := \mathbb{F} _ 2 ^ {32} の元として解釈する。 整数の xor が  V における  + に対応することに注意せよ。 このとき、Xorshift における演算は全て  V 上の線形変換となっている。 従って、線形写像  A: V \to V によって Xorshift は以下のように表現できる (酷く厳密性に欠ける表現ではあるが)。

V Xorshift() {
    static V y = seed;
    y = A(y);
    return y;
}

つまり、Xorshift の乱数列は初項を  t として  t, A(t), A ^ 2(t), \dots と表現される。  A の最小多項式 p(X) \in \mathbb{F} _ 2 \lbrack X \rbrack とすると、定義より  p(A)(t) = 0 である。  S := \lbrace i \mid \lbrack X ^ i \rbrack p(X) = 1 \rbrace とすると、 S は求める集合になっている。

Cayley-Hamilton の定理より、 p として最小多項式ではなく固有多項式を用いても良い。 周期  2 ^ {128} - 1 の実装例を用いた場合でも、全く同様の議論で衝突を起こすことができる。

Wikipedia の 32 bit の実装を衝突させる例

固有多項式を用いて、 S = \lbrace 0, 6, 9, 14, 15, 17, 18, 19, 20, 21, 32 \rbrace を得た。 seed を自由に変えて実行してみると、毎回 hash が衝突していることが分かる。

[C++] gcc 11.1.0 - Wandbox

追記

 \mathbb{F} _ 2 \lbrack X \rbrack における  2 乗は  X の指数を  2 倍したものに等しいため、Xorshift を  1 つ飛ばしにしてもそのまま衝突する。 また、 \mathrm{rev}(p) p は reverse について不変なので、Xorshift を逆から読んでも衝突する。

*1:この変化は役に立たない。追記を参照せよ。

重心分解で 1 点更新区間取得

概要

木上の等高線集約クエリ - suisen のブログ この記事で未解決になっていた変種  1, 1' を解決する。

 M を可換モノイドとし、 N 頂点の木  T の各頂点には  M の元が書き込まれているとする。

  •  v \in  V(T) に書かれている値を  x \in M に書き換える。
  •  v \in V(T) から距離  a 以上  b 未満の頂点に書かれた値の総和を出力する。

以上のクエリを前計算  \Theta(N \log(N)) 1 クエリ  \Theta( (\log(N) ) ^ 2) の時間計算量で処理できる。 ただし  M の演算は  \mathrm{O}(1) で可能とする。

また、同様の手法によって区間作用  1 点取得も処理できる。 JOI 2021/2022 春合宿 day3-2 がこの問題だが、時間計算量的に通るかどうかはよく分からない。

距離が辺の重み付きで定義される場合も前計算が  \Theta(N (\log(N) ) ^ 2) かかるものの同様に処理できる。

アルゴリズム

概要で述べた記事の内容を概ね理解していることを前提とする。

 T から重心を取り除いたときに分かれる部分木たちを  L とする。 従来の手法で問題となっていたのは、重心の次数が  \lvert L \rvert になってしまい、ここの処理に時間がかかるという点である。 そこで、 L を一度に重心に付ける  \lvert L \rvert 分木ではなく、 V(T) と対応しないノードを間に挟み二分木を作る。

f:id:noshi91:20220327041409j:plain

上図において実線の三角形は  L の元であり、点線の長方形は  V(T) に対応しない新たに追加されたノードである。 書き込まれた数は後述するアルゴリズムにおける重みである。

新たに追加されたノードについても、その子孫の頂点たちを深さ順に並べて Segment Tree で管理する。 また、 L の元については従来の重心分解同様、再帰的に重心を取り除き木構造を作成し、得られた根付き木を  U とする。 この状態で更新クエリを処理することを考えると、 U において  v \in V(T) から根に至るのパスの長さが  \mathrm{O}(\log(N) ) になっているかどうかが問題となる。  L の元の重みをその部分木の頂点数として  L を葉とする Huffman coding の木を作ることで、それが達成される。 具体的には以下のアルゴリズムを実行する。

  1.  Q := L とする。
  2.  \lvert Q \rvert \geq 2 の間、以下を繰り返す。
    •  Q から重みが最小及び  2 番目に小さい元を取り出し、それらを子とするノードを作成する。 作成されたノードの重みを子の重みの和と定義し、 Q に追加する。

 U において  v \in V(T) から根に至るパス  P の長さを考える。  x \in V(U) の重み  w(x) を、 U において  x の子孫である  T の頂点の個数と定義する。  P において  v の次に現れる  T の頂点を  x とすると、 v- x パスの長さは  c + c \log(w(x) / w(v) ) \cdots \bigstar 以下である。 これは  v の親を  2 つ辿ると  w 2 倍以上になるためであり、Huffman coding の木を作るアルゴリズムから証明できる。  P において  T の頂点が現れる回数は重心分解自体の性質により  \mathrm{O}(\log(N) ) であり、 \bigstar の第  1 項の和は  \mathrm{O}(\log(N) ) となる。  \bigstar の第  2 項の和は telescoping sum の形であるから  \mathrm{O}(\log(N) ) となる。

雑記

厳密な定式化はできていないが、領域木でできるクエリは全部できるのではないかと思う (形が同じなので)。 裏を返すと、一般の区間更新区間取得は難しそうに見える。

ICPC2021 Yokohama Regional 感想

チームメイトに投げる場面が多く自分の実力の低下を感じた。 K を検索したのが一番の貢献だったはずだが通せなかったため無に 来年度どうするかあまり考えていない (参加はしたい)

今年は特に問題が面白いと感じたので upsolve 予定

周期関数の極小点を求めるアルゴリズム

定義

 f: \mathbb{Z} \to \mathbb{R}, ~ p \in \mathbb{Z} _ {\gt 0} \forall x. ~ f(n) = f(n + p) を満たすとする。  n を与えると  f(n) \mathrm{O} (1) の時間計算量で取得できるとして、 f(n - 1) \geq f(n), ~ f(n) \leq f(n + 1) を満たす  n \Theta(\log ( p ) ) の時間計算量で求めるアルゴリズムが存在する。

アルゴリズム

一般の数列を三分探索すると極小点が求まるとは限らない - noshi91のメモ

上記の記事の後半で紹介しているアルゴリズムを、 (l, m, r) := (0, p, 2p) で初期化して、 f に対して適用すればよい。

実用例

 2 次元平面上に点群があり、ベクトル  v の方向に最も遠い点を求めることを考える。 点群の凸包を  P、その頂点を順に  a _ 1, a _ 2, \dots として  f(n) := v \cdot a _ n と定義する。  f P の頂点数を周期とする関数であり、 f の極大点が求める点である。

一般の数列を三分探索すると極小点が求まるとは限らない

概要

三分探索を行うと極小点が  1 つも求まらないような数列を構成した。 また、一般の数列に対して  \Theta(\log (N)) の時間計算量で極小点を  1 つ求めるアルゴリズムを説明する。

定義

 A = (A _ 1, A _ 2, \dots, A _ N) ~ (N \geq 1) を実数列とする。 整数  i ~ (1 \leq i \leq N) A の極小点であることを、 A _ {i - 1} \geq A _ {i}, ~ A _ {i} \leq A _ {i + 1} を満たすことと定義する。 ただし仮想的に  A _ 0 = A _ {N + 1} = + \infty とする。

三分探索

 A が丁度  1 つの極小点を持つとき、三分探索を用いると  A \Theta(\log (N)) 回アクセスすることでその極小点を求めることができる。 三分探索は以下のようなアルゴリズムである。

最初、 (l, r) := (0, N + 1) とする。 「 A の極小点  i l \lt i \lt r を満たす」という条件を維持したまま  r - l を狭めていく。

  •  \textrm{(i)}  r - l = 2 の場合

     (l + r) / 2 が極小点であり、アルゴリズムを終了する。

  •  \textrm{(ii)}  r - l \geq 3 の場合

     (l', r') := (\lfloor (2l + r) / 3 \rfloor, \lfloor (l + 2r) / 3 \rfloor ) とする。  (l, (2l + r) / 3, (l + 2r) / 3, r) は隣接項の差が  1 以上であるから、それぞれに floor 関数を適用した  (l, l', r', r) は狭義単調増加列になることに注意。

    •  \textrm{(ii - i)}  A _ {l'} \lt A _ {r'} の場合

       (l, r) \leftarrow (l, r') と更新して繰り返す。

    •  \textrm{(ii - ii)}  A _ {l'} \gt A _ {r'} の場合

       (l, r) \leftarrow (l', r) と更新して繰り返す。

 \textrm{(ii)} が起こるたびに  r - l \lceil 2(r - l) / 3 \rceil 以下になるから、時間計算量は  \mathrm{O}(\log (N)) である。

三分探索で極小点が求まらないような一般の数列

 N = 26 の例を図示した。 数列ではなく関数のグラフになっているが、適当に離散化して考えて欲しい。

f:id:noshi91:20220313203519j:plain

 3 回のイテレーションの後  (l, r) = (9, 17) となり、極小点が存在しないことが分かる。

一般の数列に対して、極小点を  1 つ求めるアルゴリズム

黄金分割探索は、一般の数列に対しても  1 つ極小点を求めることができる (詳細略)。 本記事では、黄金分割探索とは異なるアルゴリズムを説明する。

最初  (l, m, r) := (0, \lceil N / 2 \rceil, N + 1) とする *1。 「 l \lt m \lt r かつ  A _ l \geq A _ m, ~ A _ m \leq A _ r」を満たすという条件を保ちつつ、 r - l を狭めていく。

  •  \textrm{(i)}  r - l = 2 の場合

     l + 1 = m = r - 1 であり、不変条件から  m が極小点である。アルゴリズムを終了する。

  •  \textrm{(ii)}  r - l \geq 3 の場合

     (l ^ {\prime}, r ^ {\prime}) := (\lfloor (l + m) / 2 \rfloor, \lceil (m + r) / 2 \rceil) とする。  l \leq l' \lt m \lt r' \leq r である。

    •  \textrm{(ii - i)}  A _ {l'} \lt A _ m の場合

      不変条件と併せて  A _ l \geq A _ m \gt A _ {l'} である。 特に  A _ l \neq A _ {l'} より  l \lt l' も満たされるため、 (l, m, r) \leftarrow (l, l', m) と更新して繰り返す。 f:id:noshi91:20220313211811j:plain

    •  \textrm{(ii - ii)}  A _ {m} \gt A _ {r'} の場合

       \textrm{(ii - i)} の場合と対称である。  (l, m, r) \leftarrow (m, r', r) と更新して繰り返す。

    •  \textrm{(ii - iii)}  A _ {l'} \geq A _ m, ~ A _ m \leq A _ {r'} の場合

       (l, m, r) \leftarrow (l', m, r') と更新して繰り返す。 f:id:noshi91:20220313212122j:plain

 w := \max (m - l, r - m) と定義すると、 \textrm{(ii)} が起こるたびに  w \lceil w / 2 \rceil 以下になるから、時間計算量は  \mathrm{O}(\log (N)) である。

*1:実は  m = 1 としてもほとんど問題ないが、イメージのしやすさのため  m は中央付近に取った