miscalc のブログ

主に競プロの話をします

floor に関する等式の証明

最近気づいたので書きます（有名かも？）

有名な等式として \displaystyle \left\lfloor \frac{N}{ab} \right\rfloor = \left\lfloor \frac{\left\lfloor \frac{N}{a} \right\rfloor}{b} \right\rfloor（ N, a, b は正整数）があります。これを証明します。

整数 n と実数 r に対して

\begin{align} n \leq r \iff n \leq \lfloor r \rfloor \tag{ * } \end{align}

であることを利用します。 n を整数として

\begin{align} &\phantom{{}\iff{}} abn \leq N \tag{1} \\ &\iff bn \leq \frac{N}{a} \tag{2} \\ &\iff n \leq \frac{N}{ab} \tag{3} \end{align}

ですが、 (3) の後でだけ ( * ) を適用すると \displaystyle n \leq \left\lfloor \frac{N}{ab} \right\rfloor が得られ、 (2), (3) の両方の後で ( * ) を適用すると \displaystyle n \leq \left\lfloor \frac{\left\lfloor \frac{N}{a} \right\rfloor}{b} \right\rfloor が得られます。このふたつが同値で、右辺が整数であることから、右辺は等しいです。 \square

雑にいうと、操作の途中で floor をとったりとらなかったりしてよい、という感じです。 \displaystyle \left\lfloor \sqrt{\frac{N}{i}} \right\rfloor = \left\lfloor \sqrt{\left\lfloor \frac{N}{i} \right\rfloor} \right\rfloor なども同様に示せますね。

XOR の基底　事実・できることまとめ

間違っている部分やもっと簡単にできる部分などがあれば教えていただけると助かります。

この記事での記法など

  • XOR 演算子は \oplus で書く。

  • 非負整数の（有限）集合 S に対して全体の XOR \displaystyle \bigoplus _ {s \in S} s を \mathrm{XOR}(S) とも書く。

  • 非負整数の（有限）集合 S に対して X(S) = \{ \mathrm{XOR}(S) \mid T \subseteq S \}（何個か選んで作れる XOR 全体の集合）と書く。

  • S を \mathbb{F} _ 2 行列とみたときの基底（XOR 基底）の 1 つを B(S) と書く。

  • S の要素はワードサイズに収まることを仮定する。ワードサイズを w と書く。

本題

[1] 作れる XOR の種類数

|X(S)| = 2 ^ {|B(S)|}

[2] XOR を作る方法の数

c(x) := \# \{ T \subseteq S \mid \mathrm{XOR}(T) = x \} とする。

  • x \notin X(S) のとき、 c(x) = 0
  • x \in X(S) のとき、 c(x) = 2 ^ {|S| - |B(S)|}

つまり、作れる XOR はどれも、作る方法の数が同じである。*1

F - Limited Xor Subset

No.2672 Subset Xor Sum - yukicoder

[3] 作れる XOR のうち、ある bit が立っているものの個数

  • その bit が立っている要素が S にない場合、明らかに 0

  • そうでない場合、その bit が立っているものと立っていないものは半々である。つまり、 \# \{ x \in X(S) \mid \mathrm{test}(x, j) = 0 \} = \# \{ x \in X(S) \mid \mathrm{test}(x, j) = 1 \} である（ x の j bit 目を \mathrm{test}(x, j) と表した）。

    • 理由： B(S) から選ばれる要素のうち、 j bit 目が立っているものの個数の偶奇で分かれる。 j bit 目が立っている要素を 1 つ固定し、それを選ぶ・選ばないを flip することで全単射が作れる。

[4] chmin で XOR 基底を求めるアルゴリズム（いわゆる noshi 基底）

このアルゴリズムを実行すると、msb*2 が distinct な基底が得られる。行列の言葉でいうと、（行を降順ソートすれば）階段行列になる。 S の要素 1 つを処理するのに O(w) 時間。

[5] 作れる XOR のうち k 番目に小さいもの

[4] のアルゴリズムを少し変形すると、msb が distinct なだけでなく

  • msb は他の要素で立っていない（すなわち、 b_i, b_j \in B(S), b_i \neq b_j ならば \mathrm{test}(b_j, \mathrm{msb}(b_i)) = 0）

を満たす基底が得られる（具体的なアルゴリズムはこの節の一番下にある問題の解説を参照）。行列の言葉でいうと、（行を降順ソートすれば）簡約行列になる。基底を計算するパートの計算量は [4] と同じ。

この基底をソートする（昇順に b_0, b_1, \ldots, b_{m-1} とする）。すると、各 b_i を選ぶか選ばないかを 01 で表したものが k の 2 進表記と対応する。すなわち、 X(S) の要素を昇順に x_0, x_1, \ldots, x _ {2 ^ m - 1} とすると、 \displaystyle k = \sum _ {i = 0} ^ {m-1} c_i 2 ^ i \: (c _ i \in \{ 0, 1 \}) のとき \displaystyle x_k = \bigoplus _ {i = 0} ^ {m-1} c_i b_i となる。クエリ O(w) 時間。

verify に使える・より詳細な解説が書いてある：G - Partial Xor Enumeration

[6] 作れる XOR のうち v 以下で最大のもの（lower_bound や upper_bound のようなもの）

単純な二分探索で自明に O(w ^ 2) 時間になるが、もっと高速になる。

[5] の基底をとる。基底を降順に見ていき、「選んでも v を超えないなら選ぶ」ことを繰り返す。クエリ O(w) 時間。

[7] 奇数個・偶数個選んで作れる XOR

  1. 今あるどの bit よりも上の余っている bit を 1 つ選び、全部の値で立てておく。その bit が 0 なら偶数個、 1 なら奇数個選んだといえる。 X(S) を昇順に並べたとき、 [0, 2 ^ {|B(S)|-1}) 番目はこの bit が 0 に、 [2 ^ {|B(S)|-1}, 2 ^ {|B(S)|}) 番目はこの bit が 1 になる（[3] の事実も参照）。

  2. X _ {\mathrm{ev}}(S) := \{ \mathrm{XOR}(T) \mid T \subseteq S, |T| \equiv 0 \pmod 2\}（偶数個選んで作れる XOR の集合）と書く。 S = \{ s _ i \} _ {0 \leq i \lt n} に対して S' := \{ s _ 0 \oplus s _ i \} _ {1 \leq i \lt n} とすると、 X(S') = X _ {\mathrm{ev}}(S) となる。

No.803 Very Limited Xor Subset - yukicoder

ARC のネタバレ注意

E - Rearrange and Adjacent XOR

[8] 作れる XOR のうち、 x との XOR が k 番目に小さいもの

  1. 今あるどの bit よりも上の余っている bit を 1 つ選び、 x のその bit を立てたものを x' とする。 S' := S \cup \{ x'\} を考えると、 X(S') のうち先程新たに立てた bit が立っているもの、すなわち昇順で [2 ^ {|B(S')|-1}, 2 ^ {|B(S')|}) 番目のものが対応物（[3] の事実も参照）。

  2. 最大だけなら、後で扱う問題の解説 に書いてある方法で求まる：[4] または [5] の基底を降順に見て、 \mathrm{chmax}(x, x \oplus b _ i) を繰り返す。
    任意の k 番目で同様のことができるかは不明。

[9] 作れる XOR のうち、 x との XOR が v 以下で最大のもの

[8] 1. と [6] を組み合わせるとできる。

例題

G - Xor Cards

非負整数列 (s _ i) _ {0 \leq i \lt n} と (t _ i) _ {0 \leq i \lt n} が与えられる。 I \subseteq \{ 0, 1, \ldots, n-1 \} を I \neq \emptyset, \displaystyle \left( \bigoplus _ {i \in I} s _ i \right) \leq v のもとで選べるとき、 \displaystyle \bigoplus _ {i \in I} t _ i の最大値を求めよ（あるいは、そのような選び方はないと判定せよ）。 n \leq 1000, s _ i, t _ i, v \lt 2 ^ {30}

（記号をこの記事用に変えています）

いったん I = \emptyset を許して考える。

u _ i := s _ i \cdot 2 ^ {30} + t _ i 、 U := \{ u _ i\} _ {0 \leq i \lt n} とすると、求めるものは

  • X(U) の (v + 1) \cdot 2 ^ {30} 未満の要素のうち下 30 bit が最大のもの

となる。 (v + 1) \cdot 2 ^ {30} 未満となる境界を [6] の方法で求める。 r _ 0 番目が (v + 1) \cdot 2 ^ {30} 未満で最大であるとすると、条件を満たす選び方に対応する 2 進数が、区間 [0, r _ 0 ] になっている。つまり、境界を求めるのに使った [5] の基底の下 30 bit を取り直して ( t' _ i ) _ {0 \leq i \lt |B(U)|} とすると、求めるものは

  • I' \subseteq \{ 0, 1, \ldots, |B(U)|-1 \} を、対応する 2 進数が区間 [0, r _ 0 ] にあるように選ぶときの \displaystyle \bigoplus _ {i \in I'} t' _ i の最大値

となる。（ここで |B(U)| \leq 60 であることが計算量的に効いてくる。）（ここで先程許した I = \emptyset を考える。 I = \emptyset 以外で条件を満たせないのは、 r _ 0 = 0 で、 I' = \emptyset に対応する I が I = \emptyset の自明な 1 通りしかない場合。[2] より、 n = |B(U)| かどうかで判定できる。）

f(j, r, x) := ( t' _ i ) _ {0 \leq i \leq j} から、対応する 2 進数が区間 [0, r ] にある選び方をしたとき、選んだ整数と x の総 XOR の最大値 とする。答は f(|B(U)|-1, r _ 0, 0) である。 j 番目を選ぶかどうかで場合分けして再帰する。

  • \mathrm{test}(r, j) = 0 のとき

    • j 番目を選んでしまうと r より大きくなるので選べない。 f(j-1, r, x) へ。

  • \mathrm{test}(r, j) = 1 のとき

    • j 番目を選ばないとき

      • 区間 [0, 2 ^ j - 1 ] に対応。つまり ( t' _ i ) _ {0 \leq i \leq j-1} から好きに選んだときの選んだ整数と x の総 XOR の最大値 となる。これは [8] なので解ける。

    • j 番目を選ぶとき

      • 区間 [2 ^ j, r ] に対応。つまり t _ j を必ず選ぶとした上で、 ( t _ i ) _ {0 \leq i \leq j-1} から区間 [0, r - 2 ^ j ] にある選び方をしたときの選んだ整数と x の総 XOR の最大値 となる。 f(j-1, r - 2 ^ j, x \oplus t _ j) へ。

やっていることは公式解説と多分同じ（再帰をループで書いているだけ）。計算量も同じく O(nw + w ^ 3) になる。

*1:教えていただきありがとうございます https://twitter.com/Sophia_maki/status/1766876738024300696

*2:most significant bit: 立っている最上位ビット

Range Chmin Point Get

書いた動機

range chmin でググるとなぜか beats の話しかヒットしないので

区間 chmin 一点取得 だけならもっと簡単にできます。（当たり前ですが chmin を chmax に替えても可）

ACL の遅延セグ木を使う

ABC179 解法 • knshnbのブログ の F のようにすればできます。

区間積などを呼び出すと壊れますが、呼び出さないことで解決

双対セグ木

双対セグメント木 - HackMD あたりを読みましょう（面白いです）

問題例

JMO 2017 予選 9 を母関数で解いてみた

競プロにもありそうな教育的な問題だと思いました。

問題

問題リンク

1, 2, \ldots, 2017 の並べ替え \sigma = (\sigma(1), \sigma(2), \ldots, \sigma(2017)) について、 \sigma(i) = i となる 1 \leq i \leq 2017 の個数を F(\sigma) とする。すべての並べ替え \sigma について F(\sigma) ^ 4 を足し合わせた値を求めよ。

解答

a _ {n, m} を、 F(\sigma) = m であるような 1, \ldots, n の順列 \sigma の個数とします。

m = 0 のときは撹乱順列の個数で、包除原理から \displaystyle a _ {n, 0} = \sum _ {k=0} ^ {n} (-1) ^ k \binom{n}{k} (n-k)! = n! \sum _ {k=0} ^ {n} \frac{(-1) ^ k}{k!} とわかります（母関数による導出もあります：指数型母関数入門 – 37zigenのHP）。

一般の m に対しては、一致する箇所 m 個を選んだ後、残り n-m 個の場所で撹乱順列を作ると考えて、 \displaystyle a _ {n,m} = \binom{n}{m} a _ {n-m, 0} = \frac{n!}{m!} \sum _ {k=0} ^ {n-m} \frac{(-1) ^ k}{k!} となります。

答えは \displaystyle \sum _ {m=0} ^ n m ^ 4 a _ {n, m} = n! \sum _ {m=0} ^ n \frac{m ^ 4}{m!} \sum _ {k=0} ^ {n-m} \frac{(-1) ^ k}{k!} （で n = 2017 としたもの）です。

母関数を使って考察します。まずは撹乱順列の個数 \displaystyle \sum _ {k=0} ^ {n-m} \frac{(-1) ^ k}{k!} の母関数を考えましょう。 \displaystyle \frac{(-1) ^ k}{k!} = \lbrack x ^ k \rbrack e ^ {-x} であり、撹乱順列の個数はその累積和なので \displaystyle \frac{1}{1 - x} をかけたものが母関数になります。つまり \displaystyle \sum _ {k=0} ^ {n-m} \frac{(-1) ^ k}{k!} = \lbrack x ^ {n-m} \rbrack \frac{e ^ {-x}}{1 - x} です。

次に \displaystyle \frac{m ^ 4}{m!} の母関数を考えましょう。このままでは考えにくいので、 m ^ 4 = m + 7m(m-1) + 6m(m-1)(m-2) + m(m-1)(m-2)(m-3) と変形します（このように冪乗を下降冪の級数で表したときの係数は 第 2 種スターリング数 です）。すると

\begin{align*} \frac{m ^ 4}{m!} &= \frac{1}{(m - 1)!} + \frac{7}{(m - 2)!} + \frac{6}{(m - 3)!} + \frac{1}{(m - 4)!} \\ &= \lbrack x ^ {m-1} \rbrack e ^ x + 7 \lbrack x ^ {m-2} \rbrack e ^ x + 6 \lbrack x ^ {m-3} \rbrack e ^ x + \lbrack x ^ {m-4} \rbrack e ^ x \\ &= \lbrack x ^ m \rbrack ( x e ^ x + 7 x ^ 2 e ^ x + 6 x ^ 3 e ^ x + x ^ 4 e ^ x) \end{align*}

となります（負の階乗の逆数や負冪の係数は 0 とみなしておきます）。よって答えは

\begin{align*} &\phantom{{}={}} n! \sum _ {m=0} ^ n \frac{m ^ 4}{m!} \sum _ {k=0} ^ {n-m} \frac{(-1) ^ k}{k!} \\ &= n! \sum _ {m=0} ^ n \lbrack x ^ m \rbrack (x + 7 x ^ 2 + 6 x ^ 3 + x ^ 4) e ^ x \cdot \lbrack x ^ {n-m} \rbrack \frac{e ^ {-x}}{1 - x} \\ &= n! \lbrack x ^ n \rbrack \frac{x + 7 x ^ 2 + 6 x ^ 3 + x ^ 4}{1 - x} \\ &= n! \lbrack x ^ n \rbrack (x + 8 x ^ 2 + 14 x ^ 3 + 15 x ^ 4 + 15 x ^ 5 + 15 x ^ 6 + \cdots) \\ &= 15 \cdot 2017! \end{align*}

となります。（補足：2 行目から 3 行目は FPS の積の定義、3 行目から 4 行目は \displaystyle \frac{1}{1 - x} 倍が累積和）

競プロ

2017 を n に、 4 を d に一般化したときの答えを \mathrm{ans}(n, d) とすると、 \displaystyle \mathrm{ans}(n, d) = n! \sum _ {k = 0} ^ {\min(n, d)} S(d, k) です（ S(\cdot, \cdot) は第 2 種スターリング数）。これを（mod 998 とかで）計算する競プロの問題だと思うことにします。[多項式・形式的べき級数] 高速に計算できるものたち | maspyのHP を参考にすると、次のことがわかります。

  • D 固定、クエリ n に対して \mathrm{ans}(n, D) を求める：「第 2 種 Stirling 数」の項の「前者」の方法で S(D, 0), \ldots, S(D, D) を O(D \log D) で前計算して累積和をとることで、各クエリについて（階乗の計算を除いて） O(1) で求まる。

  • N 固定、 d = 0, 1, \ldots, D \ (D \leq N) に対して \mathrm{ans}(N, d) を求める： S(d, k) がすべての k = 0, 1, \ldots, d について足されることから Bell 数になって、 O(D \log D) で求まる（「Bell 数」の項も参照）。

  • N 固定、 d = 0, 1, \ldots, D に対して \mathrm{ans}(N, d) を求める：「第 2 種 Stirling 数」の項の「後者」を参考にすると \displaystyle \sum _ {k=0} ^ {\min(N, D)} \frac{1}{k!} (e ^ x - 1) ^ k を D 次まで求めればよい。これは \displaystyle \sum _ {k=0} ^ {\min(N, D)} \frac{1}{k!} x ^ k を Polynomial Taylor Shift した後 e ^ x を代入すればよいので（「Polynomial Taylor Shift」および「特殊な合成」の項を参照）、（階乗の計算を除いて）全体 O(D \log ^ 2 D) で求まる。

Stern–Brocot Tree についてのメモ

この記事の目標

SBT についてざっくり理解して、次の機能を実装できるようになります。

おことわり

  • 新規性は特になく、既存の記事を自分用にまとめた程度の内容になっています。

  • 気持ち寄りの説明が多めで、議論がやや雑な部分があるかもしれません。

Stern–Brocot Tree について

定義

Stern–Brocot Tree (SBT) は、次のように定められる、（無限に続く）完全二分木です。

  • 各頂点は非負整数の 4 つ組 (p, q, r, s) を持ち、正の有理数 \dfrac{p+r}{q+s} が対応する。*1

  • 根は (0, 1, 1, 0) である（対応する有理数は \dfrac{1}{1}）。

  • (p, q, r, s) の左の子は (p, q, p+r, q+s)、右の子は (p+r, q+s, r, s) である。

性質

Stern–Brocot Tree には次の重要な性質があります。

  1. 二分探索木である。つまり、 (左の子孫の値) \lt (親の値) \lt (右の子孫の値) である。（cf. 加比の理）

  2. 根から右の子に a_0 回、左の子に a_1 回、 \ldots、( k-1 が奇数 ? 右の子 : 左の子) に a_{k-1} 回進んだ頂点が対応する有理数は、連分数 a _ 0 + \dfrac{1}{a _ 1 + \dfrac{1}{\ddots + \dfrac{1}{a _ {k-1} + 1}}} として表される。

  3. 正の有理数がすべて現れる。（2. からわかる。すべての正の有理数は連分数展開を持つので）

  4. 逆に、各頂点に対応する正の有理数は相異なる。（1. からわかる。あるいは、正の有理数の連分数展開は a_i \geq 1 \: (i \geq 1) に限定すれば一意なので、と思ってもよい）

  5. 根から始めて、有理数 \dfrac{p}{q} に対応する頂点に到達するまでのパスを考える。「右の子/左の子 に d 回進む」ことをひとまとまりの操作とみなしたときの操作回数（パスの連長圧縮の長さ、2. での k の値 *2）は O(\log(p+q)) である。（2. からわかる。連分数展開は実質的に互除法なので）

  6. 頂点 (p, q, r, s) の子孫の集合は、開区間 \left(\dfrac{p}{q}, \dfrac{r}{s} \right) になる（ \dfrac{1}{0} は \infty とみなす）。（気持ちとしては、左の子に d 回進むと \dfrac{p + (dp + r)}{q + (dq + s)}、右の子に d 回進むと \dfrac{(p + dr) + r}{(q + ds) + s} で、 d \to \infty とするとそれぞれ \dfrac{p}{q}、 \dfrac{r}{s}）

Library Checker - Stern–Brocot Tree

ENCODE_PATH

有理数から SBT 上のパス（の連長圧縮）を得る問題です。性質 2. を利用して連分数展開すればよいです。（性質 5. より、出力は多くなりません）

DECODE_PATH

SBT 上のパス（の連長圧縮）から有理数を得る問題ですが、素直に根から辿ると有理数だけでなく頂点の 4 整数 (p, q, r, s) を得ることができます（そう実装すると後々楽です）。（これも性質 5. より、入力は多くなりません）

LCA

2 つの有理数の LCA（に対応する有理数）を求める問題です。 2 つの有理数を ENCODE_PATH したのち、根から辿っていってどこまで一致するか見ればよいです。性質 5. より間に合います。

ANCESTOR

有理数 \dfrac{p}{q} の祖先であって深さ d の頂点（に対応する有理数）を求める問題です。これも \dfrac{p}{q} を ENCODE_PATH したのち、根から深さ d まで辿っていけば性質 5. より間に合います。

RANGE

有理数が与えられて、その子孫の下限と上限を求める問題です。性質 6. そのものです。ENCODE_PATH してから DECODE_PATH すると (p, q, r, s) が得られるので、それが答えです。

ABC333 G - Nearest Fraction

SBT 上を探索する話です。

より一般に、次の形で考えましょう。（この定式化は [2], [3] をかなり参考にしています）

関数 f \colon \mathbb{R} _ {\gt 0} \to \{ \mathrm{true}, \mathrm{false} \} は単調性がある。すなわち、ある \alpha \in \mathbb{R} _ {\gt 0} が存在して、 x \leq \alpha \implies f(x) = \mathrm{true}, x \gt \alpha \implies f(x) = \mathrm{false} といった形になっている（等号・不等号や true/false は入れ替わってもよい）。このとき、境界 \alpha を分母・分子が n 以下の有理数で近似したい。 \dfrac{p}{q} \leq \alpha \lt \dfrac{r}{s} *3 を満たす 0 \leq p,q,r,s \leq n のうち、 \alpha - \dfrac{p}{q} および \dfrac{r}{s} - \alpha が最小になるものを求めよ。

性質 1. を利用して、二分探索木を探索することを考えます。性質 6. より、探索の途中で頂点 (p, q, r, s) にいるとき、 (p, q, r, s) は \dfrac{p}{q} \leq \alpha \lt \dfrac{r}{s} を満たしています。よって、素直に二分探索木を探索して、 n を超えたら打ち切る、とすると O(n) で解けます。性質 5. を利用して、「その方向にどれだけ進むか」を二分探索 *4 すれば速くなります（log ふたつに見えて実は O(\log n) です。証明は [2]）。

実装で自分が苦労した細かい部分を解説します。まず、頂点 (p, q, r, s) に到達したとき、そこからはもう進まずに (p, q, r, s) を答えとする条件は p + r > n または q + s > n です。 (p + r, q + s) は、左および右に 1 個進んだときの (r, s) および (p, q) だからです。

同様に、進む深さ d を二分探索するときの d の上限 \mathrm{ng} は、（左に進むのであれば） dp + r > n または dq + r > n を満たす最小の d、というようにして決めることができます *5。

こういったことは分母・分子が n 以下という制限が強い場面でなければ（例：[2] のように、真の値が有理数で、かつ分母・分子がある値以下であることが示せる場合）気にしなくてもあまり問題はありませんが、今回の問題のような場面では出力が n を超えてしまうバグの原因になるので、注意して実装する必要があります。

実装例

クリックして展開

#include <iostream>
#include <vector>
#include <functional>
#include <cassert>
using namespace std;

// stern-brocot tree
namespace sbt
{
  // (p, q, r, s) から (isleft ? 左 : 右) に d 個進んだ頂点を返す
  template<class T>
  tuple<T, T, T, T> child(T d, T p, T q, T r, T s, bool isleft)
  {
    if (isleft)
      r += d * p, s += d * q;
    else
      p += d * r, q += d * s;
    return make_tuple(p, q, r, s);
  }

  // (p, q, r, s) の親を返す
  // (p, q, r, s) が根なら (0, 0, 0, 0) を返す
  template<class T>
  tuple<T, T, T, T> parent(T p, T q, T r, T s)
  {
    if (p == 0 && q == 1 && r == 1 && s == 0)
      return make_tuple(0, 0, 0, 0);
    if (p < r || q < s)
      r -= p, s -= q;
    else
      p -= r, q -= s;
    return make_tuple(p, q, r, s);
  }

  // 1/1 から p/q へのパスが右に a_1 回、左に a_2 回、… としたときの a を返す
  // (a_1, …, a_k + 1) は p/q の連分数展開になっている
  template<class T>
  vector<T> encode_path(T p, T q)
  {
    vector<T> a;
    if (p < q)
    {
      a.emplace_back(0);
      swap(p, q);
    }
    while (p != 1)
    {
      a.emplace_back(p / q);
      p %= q;
      swap(p, q);
    }
    if (!a.empty())
    {
      if (a.back() == 1)
        a.pop_back();
      else
        a.back()--;
    }
    return a;
  }

  // 1/1 から右に a_1 回、左に a_2 回、… 移動した頂点 (p, q, r, s) を返す
  // (a_1, …, a_k + 1) は (p+r)/(q+s) の連分数展開になっている
  template<class T>
  tuple<T, T, T, T> decode_path(const vector<T> &a)
  {
    T p = 0, q = 1, r = 1, s = 0;
    for (int i = 0; i < (int)a.size(); i++)
      tie(p, q, r, s) = child(a[i], p, q, r, s, i & 1);
    return make_tuple(p, q, r, s);
  }

  // p/q と r/s の LCA (f, g, h, k) を返す
  template<class T>
  tuple<T, T, T, T> lca(T p, T q, T r, T s)
  {
    vector<T> a = encode_path(p, q), b = encode_path(r, s);
    int n = min(a.size(), b.size());
    T f = 0, g = 1, h = 1, k = 0;
    for (int i = 0; i < n; i++)
    {
      T c = min(a[i], b[i]);
      tie(f, g, h, k) = child(c, f, g, h, k, i & 1);
      if (a[i] != b[i])
        break;
    }
    return make_tuple(f, g, h, k);
  }

  // p/q の祖先であって深さが d のノード (f, g, h, k) を返す
  // 存在しなければ (0, 0, 0, 0)
  template<class T>
  tuple<T, T, T, T> ancestor(T d, T p, T q)
  {
    vector<T> a = encode_path(p, q);
    T f = 0, g = 1, h = 1, k = 0;
    for (int i = 0; i < (int)a.size(); i++)
    {
      T c = min(d, a[i]);
      tie(f, g, h, k) = child(c, f, g, h, k, i & 1);
      d -= c;
      if (d == 0)
        break;
    }
    return d == 0 ? make_tuple(f, g, h, k) : make_tuple(0, 0, 0, 0);
  }

  // p/q の子孫の下限 f/g と上限 h/k を返す
  // p/q に対応する頂点 (f, g, h, k) でもある
  template<class T>
  tuple<T, T, T, T> range(T p, T q) { return decode_path(encode_path(p, q)); }

  // judge: [0, ∞] → {true, false}
  // judge は単調性を持つ、すなわち、ある実数 α を境に true と false が切り替わる
  // 分母・分子が max_value 以下の有理数のうち、α に最も近いもの (下側: p/q, 上側: r/s) を返す
  template<class T>
  tuple<T, T, T, T> search(const function<bool(T, T)> &judge, const T &max_value)
  {
    T p = 0, q = 1, r = 1, s = 0;
    const bool judge_01 = judge(0, 1), judge_10 = judge(1, 0);
    assert(judge_01 != judge_10);
    while (p + r <= max_value && q + s <= max_value)
    {
      const bool judge_now = judge(p + r, q + s);
      const bool isleft = judge_now ^ judge_01;
      T maxd;
      if (isleft)
        maxd = p == 0 ? (max_value - s) / q : min((max_value - r) / p, (max_value - s) / q);
      else
        maxd = s == 0 ? (max_value - p) / r : min((max_value - q) / s, (max_value - p) / r);
      T dl = 1, dr = 2;
      while (dr <= maxd)
      {
        auto [np, nq, nr, ns] = child(dr, p, q, r, s, isleft);
        if (judge(np + nr, nq + ns) != judge_now)
          break;
        dl = dr, dr = min(2 * dr, maxd + 1);
      }
      while (dr - dl > 1)
      {
        T dm = dl + (dr - dl) / 2;
        auto [np, nq, nr, ns] = child(dm, p, q, r, s, isleft);
        if (judge(np + nr, nq + ns) == judge_now)
          dl = dm;
        else
          dr = dm;
      }
      tie(p, q, r, s) = child(dl, p, q, r, s, isleft);
    }
    return make_tuple(p, q, r, s);
  }
}

// 使用例 (https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1208)
int main()
{
  while (true)
  {
    int p, n;
    cin >> p >> n;
    if (p == 0 && n == 0)
      break;
    const function<bool(int, int)> judge = [&](int x, int y) -> bool
    { return y * y * p < x * x; };
    auto [u, v, x, y] = sbt::search(judge, n);
    cout << x << "/" << y << " " << u << "/" << v << endl;
  }
}

参考文献

[1] Stern–Brocot tree - Wikipedia

[2] Editorial - AtCoder Beginner Contest 294

[3] 有理数上の二分探索（ざっくり導入編） - えびちゃんの日記

*1:有理数だけを対応させる定義もあります。私にとっては 4 つ組が対応すると思ったほうがわかりやすかったのでこの記事ではそちらを採用します

*2:ちょっと嘘で、 a_0 = 0 のとき 1 ずれます

*3:等号・不等号は f が満たす条件の等号・不等号で変わる

*4:指数探索すればこれより速いです。定数倍の違いだと思っています

*5:このあたりは「分母が n 以下、分子は気にしない」といった場面でも同じ考え方でできますね

「周期性」の問題の抽象化（ライブラリ）

概要

競プロでは「周期性」を利用して解く問題がたまに出ます。例：

  • D - Teleporter

    長さ N の数列 A = (A_1, A_2, \ldots, A_N) が与えられる（ 1 \leq A_i \leq N）。 x = 1 から始めて、 x \leftarrow A_x という置き換えを K 回行ったときの x の値を求めよ。 N \leq 2 \times 10 ^ 5, K \leq 10 ^ {18}

  • E - Sequence Sum

    漸化式 A_1 = X, A _ {n+1} = (A _ n ^ 2 \bmod M) で定まる数列 A に対し、 \sum_{i=1}^N A_i を求めよ。 N \leq 10 ^ {10}, 0 \leq X \lt M \leq 10 ^ 5

難しくはないのですが、毎回書くと意外と頭を壊しがちなので、抽象化してライブラリにしてみました。この記事ではそのライブラリについて書きます。なお、問題の解法自体は既知として進めます。

抽象化

これらの問題は次のように抽象化できます。

モノイド (\mathrm{Data}, \bullet) およびインデックスの集合 \mathrm{Index} があり、各インデックス x \in \mathrm{Index} にはデータ \mathrm{get \_ data}(x) \in \mathrm{Data} が対応づけられている。また、インデックス列 (x _ i) _ i が漸化式

\displaystyle x _ 0 = \mathrm{start}, x _ {i+1} = \mathrm{next \_ index}(x _ i)

で定まっている。これについて、先頭からの区間積

\mathrm{get \_ data}(x _ 0) \bullet \mathrm{get \_ data}(x _ 1) \bullet \cdots \bullet \mathrm{get \_ data}(x _ {K-1})

を求めるクエリを処理したい（ K がクエリで与えられる）。

これは、

  • x _ \lambda = x _ {\lambda + \mu} となる (\lambda, \mu) が存在する（特に最小のものをとる）。つまり、インデックスは ρ のような形で周期性を持ち、非循環部分の長さが \lambda、循環部分の長さが \mu である

  • a \in \mathrm{Data}, n \in \mathbb{N} に対して冪 a ^ n = \underbrace{a \bullet \cdots \bullet a} _ {n \ \mathrm{個}} を高速に求められる（簡単のため、この冪演算およびその他必要な演算は O(1) で行えるものとする）

という条件のもとで、前計算 O(\lambda + \mu)、クエリ O(1) で処理できる。

「冪が高速に計算できるモノイド」である必要があるのは、（ K が大きいとき）

\begin{align} &{\phantom{\bullet}} \left( \mathrm{get \_ data}(x _ 0) \bullet \cdots \bullet \mathrm{get \_ data}(x _ {\lambda - 1}) \right) \tag{非循環部分} \\ &\bullet \left( \mathrm{get \_ data}(x _ \lambda) \bullet \cdots \bullet \mathrm{get \_ data}(x _ {\lambda + \mu - 1}) \right) ^ n \tag{循環部分} \\ &\bullet \left( \mathrm{get \_ data}(x _ {\lambda + n\mu}) \bullet \cdots \bullet \mathrm{get \_ data}(x _ {K - 1}) \right) \tag{あまり} \end{align}

と変形して計算するためです*1。クエリ O(1) にできるのは、クエリで計算する区間積が 0 スタートのものと \lambda スタートのものしかなく前計算できるためです。

設計

コンストラクタ

Period<Index, Data, next_index, get_data, op, e, power> pe(start); *2

  • インデックスの型 Index（演算子 != が定義されている必要があります）
  • データの型 Data
  • 一つ後のインデックスを得る関数 Index next_index(Index x)
  • インデックスからデータを得る関数 Data get_data(Index x)
  • \mathrm{Data} 上の二項演算 Data op(Data a, Data b)
  • (\mathrm{Data}, \bullet) の単位元 Data e()
  • a ^ n = \underbrace{a \bullet \cdots \bullet a}_{n \ \mathrm{個}} を計算する関数 Data power(Data a, long long n)

を定義する必要があります。start は最初のインデックス x _ 0 です。

query

Data query(long long K)

\mathrm{get \_ data}(x _ 0) \bullet \mathrm{get \_ data}(x _ 1) \bullet \cdots \bullet \mathrm{get \_ data}(x _ {K-1}) を返します。

実装

  • 既に訪れたかどうかを配列で管理するのが最も思いつきやすいと思いますが、Index の型によっては連想配列を用いる必要が出てきてしまいます。フロイドの循環検出法 *3 を用いることでこれを回避しています。

  • 限界高速化とかはしていません。

#include <iostream>
#include <vector>
using namespace std;

template <typename Index, typename Data, Index (*next_index)(Index), Data (*get_data)(Index), Data (*op)(Data, Data), Data (*e)(), Data (*power)(Data, long long)>
struct Period
{
private:
  vector<Data> dat, prod_head, prod_tail;
  int lambda, mu;
  Data head, tail;

public:
  Period(Index start)
  {
    Index a = start, b = start;
    do
    {
      dat.emplace_back(get_data(a));
      a = next_index(a);
      b = next_index(next_index(b));
    } while (a != b);
    mu = dat.size();
    Index c = start;
    while (a != c)
    {
      dat.emplace_back(get_data(a));
      a = next_index(a);
      c = next_index(c);
    }
    lambda = (int)dat.size() - mu;

    prod_head.resize((int)dat.size() + 1);
    prod_tail.resize((int)dat.size() + 1);
    prod_head[0] = e(), prod_tail[lambda] = e();
    for (int i = 1; i < (int)prod_head.size(); i++)
      prod_head[i] = op(prod_head[i - 1], dat[i - 1]);
    for (int i = lambda + 1; i < (int)prod_tail.size(); i++)
      prod_tail[i] = op(prod_tail[i - 1], dat[i - 1]);

    head = prod_head[lambda];
    tail = prod_tail[lambda + mu];
  }

  Data query(long long k)
  {
    if (k <= (int)dat.size())
      return prod_head[k];

    Data mid = prod_tail[lambda + (k - lambda) % mu];
    return op(op(head, power(tail, (k - lambda) / mu)), mid);
  }
};

使用例

A_i をインデックス x_i にして、 \mathrm{get \_ data}(A_i) は A_i をそのまま返すようにします。二項演算 \bullet は整数の足し算 + にすればよいです。このとき power は掛け算（個数倍）になります。

using Index = long long;
using Data = long long;
int M;
Index next_index(Index x) { return x * x % M; }
Data get_data(Index x) { return x; }
Data op(Data a, Data b) { return a + b; }
Data e() { return 0; }
Data power(Data a, long long n) { return n * a; }

int main()
{
  long long N;
  int X;
  cin >> N >> X >> M;

  Period<Index, Data, next_index, get_data, op, e, power> pe(X);
  cout << pe.query(N) << endl;
}

  • D - Teleporter

    長さ N の数列 A = (A_1, A_2, \ldots, A_N) が与えられる（ 1 \leq A_i \leq N）。 x = 1 から始めて、 x \leftarrow A_x という置き換えを K 回行ったときの x の値を求めよ。 N \leq 2 \times 10 ^ 5, K \leq 10 ^ {18}

各操作の後の x の値をインデックス（兼データ）にします。二項演算 op がややトリッキーです。直感的には「最も右のもの」としたいですが、このままだと単位元がないので困ってしまいます。単位元 e を「意味のないデータ」にしておいて、op は「意味のあるデータのうち、最も右のもの（なければ『意味のないデータ』）」を返すようにするとよいです。

using Index = int;
using Data = int;
vector<int> A;
Index next_index(Index x) { return A[x]; }
Data get_data(Index x) { return x; }
Data e() { return -1; }
Data op(Data a, Data b) { return b == e() ? a : b; }
Data power(Data a, long long) { return a; }

int main()
{
  int N;
  long long K;
  cin >> N >> K;
  A.resize(N);
  for (int i = 0; i < N; i++)
  {
    cin >> A[i];
    A[i]--;
  }

  Period<Index, Data, next_index, get_data, op, e, power> pe(0);
  cout << pe.query(K + 1) + 1 << endl;
}

その他の問題例

*1:単位元が必要なのは、それぞれの部分の長さが 0 になることがあるためです

*2:英語は cycle の方が正しいような気もしますが、閉路と紛らわしいので period にしました

*3:高速素因数分解のポラード・ロー法に使われていることのほうが有名？

JAG 夏合宿 2023 参加記

JAG 夏合宿 2023 に参加した。

ICPC のチームは Magical Fish で、メンバーは自分と magsta, confeito の 3 人。今回 confeito 君に合宿に誘われたので行くことにした。（本当は 3 人で行きたかったけど magsta 君が来れなかったみたいで残念）

1 日目

コンテスト

ぎりぎりに家を出たら電車が遅延して 10 分遅刻（すみません……）。

チームのもう 1 人として kumakuma さんに加わっていただいた。戦略としては、kumakuma さんは英語がかなり得意らしいので長めの問題文担当 + US 配列使いらしいので基本的に実装はしないで自分か confeito 君がやる、という感じになった。

confeito 君が A を読む。自分は B を読もうとしたが長そうだったので、kumakuma さんと交代して C を読むことにした。C を読むと一瞬で解けた（というか ABC で既出なのを思い出した）ので書いてすぐ通る。A も簡単らしいので書いてもらい、すぐ通る。

B は大変らしいということで、シンプルそうな E を読む。一瞬むずそうに見えたけど最適戦略が貪欲っぽい？ ということで 2 人に相談、正しそうということになる。細かい解法がすぐに出てこなかったが kumakuma さんが二分探索でできるというのを出してくれて、たしかにとなる。実装して AC。後で考えたら（ソート後は）差分更新で線形でいけるけど、にぶたんのほうがバグらせにくそう。

その後は F ができそうという話になる。DAG の最大パスは EDPC にあって、一般の有向グラフだと最大パスは解けないけど今回の問題は「ウォークを取ってきてそこに含まれる distinct な頂点の数」という解釈になる（はず）なので強連結成分の中を好きなだけぐるぐる回れるから、SCC して頂点に重みがついた DAG とみなした後同様の DP で解ける。とすると出次数が 1 以下という制約は何に使うんだ？ と一瞬疑ったけどそのまま書き始めてしまった（これがよくなかった）。kumakuma さんの蟻本から SCC を写経（蟻本が昔の本すぎて範囲 for すら使ってなくてやや困った）、DP はメモ化再帰が楽そうだからそれをすることにして書き終えた。サンプルも試して提出すると TLE、よく見たらメモ化になっていない。提出、TLE。この辺で N \leq 10 ^ 6 と TL 1 sec に気づく。やめてくれ〜　入力を scanf にしたけどだめで、そもそも SCC で再帰 2 回するのが重くて出次数 1 以下を活かした定数倍がより高速な解法が必要なのでは？ という話になって一旦後回しに。

K が解けるらしい。問題概要と解法（DP）を kumakuma さんから聞く。えそれ DP じゃだめでダイクストラでは？ と一瞬なったけど話し合ったところ自分が勘違いしていて、DP であっている。サッと実装して AC。

confeito 君が G 解けたというので実装を任せて、kumakuma さんと 2 人で F を詰める。出次数 1 以下なことからサイクルから辺が出ていくことはない、よって普通に DP した後サイクルだったらサイクル長を足した値を考えればよさそうで、queue を使うサイクル検出を DP しながらやる、という解法が生えた。定数倍がまだ不安だが、これ以上のものは出なかった。G の実装がまだ終わらないらしいので順位表で解かれている J も見る。既視感があって、ABC に O(N ^ 2 \log N) 想定で既出 + 誰かが O(N \sqrt{N}) で解ける旨を言っていたような気がするところまで思い出したけど、肝心の O(N \sqrt{N}) 解法を知らず、ちょっと考えたけどよくわからず、諦めた。

G サンプルが合わないらしい。問題概要と解法を聞くとかなり茨の道なことをやっていて、さらに勘違いがあったらしく MLE しそう。一旦諦めて F をやることに。実装して出したけど TLE、グラフを持つのが vector<vector<int>> ではなく vector<int> でいい、とかもやったけどだめだった。結局 ACEK の 4 完で 17 位、激冷え。

コンテスト後

G は DP の値に First / Second / Draw を持たせることばかり考えていたけど、実は普通に得点の差を評価値にしてよくてこれで状態数が削れる。これは気づくべきだった……。

F は confeito 君によると JOI 典型らしく、じゃあ F と G の担当逆にすべきだったねという話になった。サイクルが高々 1 個なのを利用するの、なるほどね。（しかしそこまでの定数倍高速化を要求するのはちょっと不快だなあという気持ちもあり）

夕食をとったり解説を聞いたりした後は部屋で交流した。入浴を済ませた後談話室で ABC に出た。机が低かったり他の団体の声が響いたりしてあまり集中できなかった。SSRS さんのタイピングが速すぎてびびっていた。

部屋に戻って、作問の話とかをしつつ就寝。

2 日目

コンテスト

1 日目チームを組んでいただいた kumakuma さんは 2 日目からソロ参加（元々ソロ参加の予定だったが 1 日目 PC を忘れたらしい）ということで、2 日目は漁師さんとチームを組んでいただいた。

昨日の反省を活かして、ある程度難しい問題は実装を始める前に人に相談しようということにした。

A を見るといかにも FPS 使いそうな 998 で、自分がやることに。ざっと眺めると他にも 998 がいっぱいあって、これは modint 書いたほうがよさそうという話になる。漁師さんに modint を写経してもらっている間に A が解けて AC。

順位表で解かれている C を見る。簡単。confeito 君と 2 人で詰めて AC。

confeito 君が B を、漁師さんが D, E を読んでいて、両方からいろいろ聞いた。confeito 君によると B はいけそうらしいので D, E を考えるけどよくわからず。

順位表で解かれている H を読んで問題概要を伝えると confeito 君から天才考察が返ってくる。コンビネーションを H でも B でも使うっぽいので H の実装ついでに写経する。あとは和が s で積が p な 2 数のパートだけちょっと詰めて実装、AC。

confeito 君が B を実装している間に順位表で解かれている K を漁師さんと 2 人で考える。超頂点作ってダイクストラかなーとかいろいろ考えたけどよくわからず。B を通して戻ってきた confeito 君（FA、すごい！）に助けを求めると、値に (\mathrm{cost}, \mathrm{from}) の配列を持つ DP をすると遷移回数がいい感じに抑えられるらしい。たしかにじゃん。実装してもらって AC。この人天才か？

このあたりで順位がだいぶよかった。順位がいいということは逆に解く問題の見極めを順位表に頼れないということでもある。残っている問題で唯一 AC が出ている F を考える。括弧列まではすぐ言い換えられて、カタラン数の証明が思い出せればちょっと進展がありそうだと思ったけど思い出せなくて苦しんでいた。（実際はこれがわかっても解けなかった気がする）

J に AC が出る。読むとゲーム。通したのは上位チームではない。とすると実験エスパーで解けるやつじゃねということになって、漁師さんが実験を書いてくれた。自分は F を考え続けるけど進展なし。

実験が書けたらしいので結果を見にいく。長方形に A, B をプロットすると見やすいという話になって（ここ地味に偉い）、見ると大きいところでは規則性がありそう。 3 \times 3 と 3 \times 5 がコーナーなので H,W \leq 5 なら愚直、そうでなければ規則性、のコードを漁師さんが書いて提出。WA。実験結果を H,W \leq 7 とかじゃなくて H \times W \leq 50 とかで表示すると、 2 \times \mathrm{BIG} のケースが処理できていないことに気づく。ここを自分が詰めて書いてもらって AC。

F を考え続けるけど進展がない。

confeito 君が L 解けた気がするらしい。聞くと MCF でできそうとのこと。これを正しいと思ってしまった（実際は嘘だったけど、指摘できなかった。とてもよくない）。漁師さんに kactl の MCF を写経してもらうと、謎の機能を使っていてコンパイルが通らない。こうならないように持ち込むライブラリはちゃんと把握しておくべきですね。Gifted Infants の MCF を見ると変なことはしていなさそうでコンパイルが通った。いざ実装、というところで confeito 君が解法の誤りに気づく。聞くとたしかに嘘だ。残り 10 分くらい、ここから立て直すのは厳しい。結局 ABCHJK の 6 完で 5 位。だいぶいい順位。この日は confeito 君がかなり強かったという印象。

コンテスト後

解説を聞くと F も L も激ヤバでびっくり。チャンスがあったとしたら D, E なのかなあ。まあ今の実力だとこの 6 完が限度という気もする。

談話室でボドゲをした後、部屋で ARC に出た。B が通せなくて激冷え・青落ちで険しい気持ちになった。

その後談話室で Kite_kuma さんと将棋を指した。tokusakurai さんとも指したかったが時間がなくて残念。また機会があればよろしくお願いします。

部屋では 1 行ずつ作問をして遊んだ。滅茶苦茶な問題になっても面白いし、それに混じってたまに良さげな雰囲気の問題が生えても面白いので、レクとしていいかもしれない。nonon さんの参加記 に問題が載っているのでよければ考えてみてください。

3 日目

コンテスト

3 日目はチームが組めず、2 人でやることになった。

confeito 君が A を、自分が B を読む。A は簡単ですぐ通る。B も簡単なはずだが、実装にちょっと手間取ってしまった。

C の概要を聞く。「過半数が同じ」が条件ということの正当性を 2 人で確認。confeito 君曰くこういう数え上げは自信ないらしいので自分が考えることに。その間 I を考えてもらう。

C を立式すると解けた気になる。998 が他にもあるので modint とコンビネーションを写経して（結果的には modint 書くべきだったかは怪しい）、実装。この modint を普段使いしていないせいでなぜかコンパイルが通らなかったりして手間取る。実装すると合わない。よく確認すると立式が間違っていることに気づく。この式だと今の解法では解けないので一旦後回しに。

I を考える。いろいろ議論すると ( を prique にためて ) がきたら「まだペアになっていない ( と組ませる」「もうペアになった ) と交代」のどっちかをやるという感じになりそう。たぶん正しいとは思ったがちょっと自信がなかった。confeito 君が実装。

その間に順位表で解かれている J を見るとド典型。C とか I より先にこっちやればよかった。

I を提出、WA。実装に細かいミスがあるのか方針がそもそも嘘なのかよくわからないので、コードを印刷してリカバリしてもらっている間に J を実装。実装中に DP だとうまくいかなくて 01BFS にする必要があることに気づく。実装重くないか心配されたけど、B でも BFS をやったこともあり大丈夫だった。AC。この間に I のバグに気づいてもらって I もすぐ AC。

あとは C, E と心中するしかないかなという感じ。C はあれから進展がなく、E を 2 人で考える。valid な順番だと区間を広げていく感じになることに confeito 君が気づいて、それを逆から見て立式すると \mathrm{dp}(i, j) := 左から i 個、右から j 個取ったときのなんかの積の総和 という感じになりそう。あとはこれが線形になればいいんだけど、わからない。自分は C と E を行ったり来たりしたが、最後までわからなかった。結局 ABIJ の 4 完で 14 位、激冷え 2。

コンテスト後

2 人だと 3 人のときに比べて明らかに思考が凝り固まっているのを感じた。やっぱり 3 人いるべきだなあ。

C 分割統治はなるほどなあ。ただもしコンテスト中に通すなら二項係数こねこねをしないといけなかった気がする。

E はもっとうまい考え方をすると自然に線形になるっぽい。2 乗を高速化することに固執しすぎたのがよくなかったかも。

さいごに

コンテスト時間が 3 日間合計で 3 + 5 + 5 (+ 1.67 + 2) = 16.67 時間というとても濃い合宿で、競プロモチベがとても高まった。来年は絶対に国内予選を通過したい。

最後になりますが、このような素晴らしい会を運営してくださった JAG の皆様や関係者の方々、また交流していただいた皆さん、誠にありがとうございました。