えびちゃんの日記

えびちゃん(競プロ)の日記です。

Montgomery 乗算について

Montgomery(モンゴメリ) 乗算の話をします。

Barrett reduction について前に 紹介したことがある のですが、それと同じように除算命令を避けて剰余演算をする手法の一つです。

前提

以下では、法を $N$ とします。また、以下の条件を満たすような $R$ を用意します。

  • $R \gt N$ である。
  • $\gcd(R, N) = 1$ である。
  • $R$ による除算・剰余演算は簡単に行える。

ここでは、$N$ は奇数であるとし $R = 2^{32}$ としてみます。 $\floor{x / R} = x \gg 32$ とできるほか、$x$ を符号なし 64-bit 整数型 x で表せば、$\floor{x / R}$ は x の上位 32-bit、$x\bmod R$ は x の下位 32-bit に対応します。

このとき、法 $N$ での $R$ の乗法逆元 $R^{-1}$ が存在します。すなわち、ある $0\lt R^{-1} \lt N$ と $0\lt N' \lt R$ が存在して $RR^{-1}-NN' = 1$ が成り立ちます。 除算を介さずに $N'$ を求める方法は後述しますが、最初に一度求めれば使い回せるので、適当に互除法などで(除算を用いつつ)求めてもよいでしょう。

Algorithm REDC

まず、$0\le T\lt RN$ から $TR^{-1}\bmod N$ を求めるアルゴリズムを紹介します。ここで、$RN \gt N^2$ に注意しておくとよいかもしれません。

  • function $\text{\small REDC}(T)$
    • $m \gets ((T \bmod R)\cdot N') \bmod R$
      • note: $0\le m\lt R$。
      • note: $m \equiv TN' \pmod{R}$。
    • $t \gets (T + mN) / R$
      • note: $mN \equiv mNN' \equiv m\cdot(RR^{-1}-1) \equiv -m \pmod{R}$。
      • note: なので $T + mN \equiv 0 \pmod{R}$ であり $t$ は整数となる。
      • note: また $tR \equiv T \pmod{N}$ なので $t \equiv TR^{-1} \pmod{N}$。
      • note: $0\le T\lt RN$ と $0\le mN\lt RN$ より $0\le t\lt 2N$。
    • return if $t \lt N$ then $t$ else $t-N$

値の表現・各演算

ここで、$i$ ($0\le i\lt N$) を使って $i$ の属する剰余類を表すのではなく、$i$ を使って $i R^{-1} \bmod N$ の属する剰余類を表すことを考えます。 あるいは、$i$ の属する剰余類を表すために $i$ ではなく $iR\bmod N$ を用いると見ることもできます。

このとき、$(a\pm b)R \equiv aR\pm bR \pmod{N}$ なので、加減算については普段通りに行うことができます。 $0\le a\pm b\lt N$ の範囲を超えた場合は、REDC 同様に $N$ を足し引きして調整します。

また、乗算については、$aR\cdot bR \equiv abR^2 \pmod{N}$ であるため、$(ab)R \equiv \text{\small REDC}(aR\cdot bR) \pmod{N}$ として求められます。

さらに、$aR \equiv \text{\small REDC}(aR^2) \pmod{N}$ なので、予め $R^2 \pmod{N}$ を求めておくことで、以降は除算なしで初期化できます。

各値の前計算

$R^2$ の計算の際は、$(R \bmod N)^2 \bmod N$ としてもよいですが、$R^2 \equiv R^2-N \pmod{N}$ であり、符号なし 64-bit 整数を使うと $R^2-N$ を $-N$ と書けるので、それを利用するのもよいでしょう。

$N'$ については、$N' \equiv -N^{-1} \pmod{R}$ であるので、$N^{-1} \pmod{R}$ を求めることを考えます。 これは、Newton 法を用いると除算を避けることができます。

$n' \gets N$ で初期化し、$n' \gets (n' \cdot (2 - n\cdot n') ) \bmod R$ で更新するのを 4 回程度繰り返すと $n' = N'$ に収束します。 $3$ 以上 $2^{30}$ 未満の各奇数について、これでうまくいくのを 1 秒程度で確認できました。高速ですごいですね。

法が $R$ じゃなくて一般のもののときは Newton 法だとうまくいかないように見えたのですが、うまくいく条件はなんなのでしょう。ねむくて考えられません。

並列化

たとえば、メモリ上で連続した位置にある 8 個の mod int をまとめて乗算したいような状況はそれなりにあります。 こうしたとき、SIMD の命令をうまく使うことで、REDC などを並列して行うことができ、高速化が期待できます。 FFT を高速化したいときなどに有効そうです。

下位ワードをまるっと捨ててよいのとかは、複雑なことを気にしなくてよいがちなのでうれしいですね。

正規化のタイミング

REDC 内で $t$ の範囲によって $t-N$ にしたりしなかったりする処理がありますが、これをしなくても $2N$ 未満にはなってくれることは前述の通りです。

そこで、乗算のための REDC で正規化をしないことにすると、加減算でも $0\le i \lt N$ に収めなくてよくなり、$0\le i\lt 2N$ に収まるように調整してもよい感じにできます。

ただし、値の比較をする際に注意しないと $x \not\equiv x+N$ と判定してしまいがちなので、気をつけましょう。

注意

各剰余計算ごとに REDC で変換するとオーバーヘッドが大きそうです。 計算の最初に $i → iR \bmod N$ に変換し、計算はずっとその表現方法のまま行い、計算が終わって出力する段階で $i → iR^{-1}\bmod N$ に戻すのがよいでしょう。

参考文献

  • Montgomery, Peter L. “Modular multiplication without trial division.” Mathematics of computation 44, no. 170 (1985): 519-521. [PDF]

おわり

朝になりました。なんで?

log をつけずに解くのが好きな人向けの記事

\(\Theta(n \polylog(n))\) 時間で解くのは簡単だけど、実際には \(\Theta(n)\) 時間で解けるよーという問題はちょくちょくあります。 競プロ界隈では log を削ることに意欲的な人は一部だけな気がします*1が、話として知っておいても損はあまりない気がするので、挙げてみます。

ここではざっくり紹介するに留め、詳細は各自調べてもらうようなスタンスになっています。

紹介

RMQ

長さ \(n\) の配列 \(a = (a_0, a_1, \dots, a_{n-1})\) に対して、区間最小値 \(\min_{l\le i\lt r} a_i\) をクエリ形式で答えよ。ただし、更新クエリは来ない*2

セグ木を使うと \(\angled{\Theta(n), O(\log(r-l))}\) で解けるものの、クエリ時間も \(O(1)\) にしたいです。

これは、sparse table と MoFR (method of Four Russians) で \(\angled{O(n), O(1)}\) で処理でき、線形 RMQ*3とか \(\angled{O(n), O(1)}\)-RMQ などの言い方で言及されがちです。

CS166 スライド など。

中央値

長さ \(n\) の配列 \(a = (a_0, a_1, \dots, a_{n-1})\) に対して、\(a\) を昇順にソートした列 \(b\) を考える。\(a\) の中央値 \(b_{\floor{n/2}}\) を答えよ*4

ソートして真ん中を答えればよいですが、ソートのために最悪 \(\Omega(n\log(n))\) 時間かかってしまいます。

これは、median of medians と呼ばれるアルゴリズムで \(O(n)\) 時間で求めることができます。 AtCoderMedian of Medians と呼ばれる別の設定の問題が出題されたため、検索しにくい感じがあります。

Introduction to Algorithms などに載っています。

篩による素数判定

正整数 \(n\) が与えられる。長さ \(n\) の整数列 \(p = (p_1, p_2, \dots, p_n)\) であって、\(i\) が素数のとき \(p_i = 1\)、そうでないとき \(p_i = 0\) となるものを答えよ。

Eratosthenes の篩と呼ばれるアルゴリズムで \(O(n\log(\log(n)))\) 時間ですが、線形篩と呼ばれる \(O(n)\) のアルゴリズムが存在します。 配る DP っぽく遷移しながら篩っていくのですが、素因数分解の一意性によって遷移の個数が \(O(n)\) 個に抑えられるような遷移をします。

実際には Atkin の篩などの劣線形*5アルゴリズムもありますが、競プロ界隈ではあまり見ないです。他の篩アルゴリズムもあった気がします。

cp-algorithms.com

過半数

長さ \(n\) の配列 \(a = (a_0, a_1, \dots, a_{n-1})\) が与えられる。\(a\) のうち過半数を占める要素があればそれを出力し、なければその旨を報告せよ。 ただし、\(x\) が過半数を占めるとは、\(|\{i\mid a_i = x\}| \gt |\{i\mid a_i \ne x\}|\) を満たすことを言う。

ソートなり連想配列なりをしてカウントすれば求められはします。ソートしたり木ベースの連想配列をしたりすると log がついてしまいます。 あるいは、ハッシュベースの連想配列では最悪計算量が保証できなくなってしまいます。

Boyer–Moore majority vote algorithm というのがあります。 ある要素 \(x\) に対して、それと異なる要素 \(y\) をペアにして消していくのを繰り返すと、消えずに残った(高々一種類の)要素が過半数の候補となります。 そこで、その候補の個数を数えればよいです。

ペアにして消すのは先頭から順に見ていけばよいので \(O(n)\) 時間でできます。大小判定 (\(\lt, =, \gt\)) ではなく同値判定 (\(=, \ne\)) しか必要ないのも汎用的ですごいです。

Boyer, Robert S., and J. Strother Moore. “MJRTY—a fast majority vote algorithm.” In Automated Reasoning, pp. 105–117. Springer, Dordrecht, 1991.

wikipedia にも載っています。

行最小値

\(n\) 行 \(m\) 列の行列 \(a = (a_{i, j})\) が与えられる。ただし陽には与えられず、\(a_{i, j} = f(i, j)\) なる \(f\) を呼び出すことで値が得られるとする。 このとき、各 \(i\) に対して \(i\) 行目の最小値を持つ列 \(\argmin_j a_{i, j}\) を求めよ。

ただし、\(a_{i, j}\) は totally monotone であるとする。すなわち、\(i\lt i'\), \(j\lt j'\) に対して \(f(i, j) \gt f(i, j') \implies f(i', j) \gt f(i', j')\) とする。 簡単のため、\(a\) は distinct とする。

totally monotone より弱い条件である monotone の性質として、\(i \lt i' \implies \argmin_j a_{i, j} \le \argmin_j a_{i', j}\) というのがあります。 これを利用した簡単な分割統治のアルゴリズムによって \(O(n + m\log(n))\) 時間で解けます。

分割統治のアルゴリズム

まず、\(n/2\) 行目の最小値を持つ列 \(j' = \argmin_j a_{n/2, j}\) を \(O(m)\) 時間で求めます。 それより上の列については \(j'\) 列目以下のみ、下の列については \(j'\) 列目以上のみを探索すればよいです。

SMAWK algorithm というのがあります。

まず偶数行目についてのみ(うまいこと再帰して)求め、奇数行目については偶数行目での答えと total monotonicity を利用して線形時間で求めます。 全体で \(O(n+m)\) 時間です。

行列のアクセス順に関して、\(i\) 行目の最小値を求めるまでは \(i\) 列目 \(f(\bullet, i)\) の値を得られないという制約を追加した問題設定もあります。 これは、ある種の DP の言い換えと見ることができる(状態 \(i\) への遷移コストの最小値が \(i\) 行目の最小値に、状態 \(j\) に最小コストで遷移してから状態 \(i\) に遷移したときのコストが \(f(i, j)\) に対応する)のですが、これも \(O(n)\) 時間で求めることができます。LARSCH algorithm と呼ばれるアルゴリズムです。

  • Aggarwal, Alok, Maria M. Klawe, Shlomo Moran, Peter Shor, and Robert Wilber. “Geometric applications of a matrix-searching algorithm.” Algorithmica 2, no. 1 (1987): 195-208.
  • Larmore, Lawrence L., and Baruch Schieber. “On-line dynamic programming with applications to the prediction of RNA secondary structure.” Journal of Algorithms 12, no. 3 (1991): 490–515.

接尾辞配列

長さ \(n\) の文字列 \(s = s_0 s_1 \dots s_{n-1}\) が与えられる。各 \(i \in\halfopen{0}{n}\) に対して、\(s\) の接尾辞 \(s_{i\dots}\) を \(s_i s_{i+1} \dots s_{n-1}\) と定義する。また、\(s_{n\dots}\) は空文字列とする。

このとき、\(\{0, 1, \dots, n\}\) の順列 \(p = (p_0, p_1, \dots, p_n)\) であって、\(s_{p_0\dots} \lt s_{p_1\dots} \lt \dots \lt s_{p_n\dots}\) となるものを求めよ。

ここでの \(p\) が接尾辞配列 (suffix array, SA) と呼ばれるものです。 これは、愚直に接尾辞を毎回作りながらソートすると \(\Omega(n^2\log(n))\) 時間かかりそうです。 蟻本などには \(\Theta(n\log(n)^2)\) の方法が載っていたり、\(\Theta(n\log(n))\) 時間のアルゴリズムが知られていたりします。

アルファベットを \(\Sigma = \{0, 1, \dots, \sigma-1\}\) (\(\sigma = n^{O(1)}\)) と仮定します。すなわち \(s_i\in\Sigma\) です。 このとき、SA-IS (suffix array, induced sorting) と呼ばれるアルゴリズムがあり、\(O(n)\) 時間でできます。

Nong, Ge, Sen Zhang, and Wai Hong Chan. “Two efficient algorithms for linear time suffix array construction.” IEEE transactions on computers 60, no. 10 (2010): 1471-1484. や CS166 のスライド など。

SA-IS の登場以前にもいくつか SA の線形時間構築のアルゴリズムは存在していたようです。

最深共通祖先

頂点 \(r\) を根とする \(n\) 頂点の木が与えられます。頂点対 \( (u, v)\) の最深共通祖先をクエリ形式で答えよ。ここで、\(u\) の祖先でもあり \(v\) の祖先でもある頂点のうち、\(r\) から最も遠い頂点を \( (u, v)\) の最深共通祖先 (lowest common ancestor, LCA) という。

Euler tour を \(O(n)\) 時間で求め、区間最小値に帰着する方法は競プロ界隈でも有名です。ここで線形 RMQ を使えばよく、\(\angled{O(n), O(1)}\) で解けます。

クエリの先読みができる状況では、別のアルゴリズムもあります。

rsk0315.hatenablog.com

Decremental neighbor query

集合に対する以下のクエリを処理せよ。

  • \(A \subseteq \{0, 1, \dots, n-1\}\) が与えられ、\(S \gets A\) で初期化。
  • 与えられた \(i\) に対し、\(A \gets A \setminus\{i\}\) で更新。
  • 与えられた \(i\) に対し、\(\min\,\{j\in A\mid j\ge i\}\) を出力。
  • 与えられた \(i\) に対し、\(\max\,\{j\in A\mid j\le i\}\) を出力。

std::setstd::collections::BTreeSet 系の、平衡木などを使えば \(\angled{O(n), O(\log(n))}\) で処理できます。 あるいは、van Emde Boas tree などを使うと \(\angled{O(n), O(\log(\log(n)))}\) などになります。

rsk0315.hatenablog.com

\(\angled{O(n), O(1)}\) になります。\(O(1)\) 時間で更新できるのってすごいですよね。

Level ancestor

頂点 \(r\) を根とする \(n\) 頂点の木が与えられます。頂点と非負整数の対 \( (u, k)\) に対して level ancestor をクエリ形式で答えよ。ここで、\(u\) の祖先のうち、\(r\) からの距離が \(k\) である頂点を答える問題を level ancestor (LA) という。

ダブリングで求める方法や HL 分解で求める方法などが有名そうですが、最悪 \(\angled{\Theta(n), \Theta(\log(n))}\) や \(\angled{\Theta(n\log(n)), \Theta(\log(n))}\) です。

ladder algorithm や jump-pointers algorithm と MoFR を組み合わせることで、\(\angled{O(n), O(1)}\) のアルゴリズムが得られます。

Bender, Michael A., and Martin Farach-Colton. “The level ancestor problem simplified.” In Latin American Symposium on Theoretical Informatics, pp. 508–515. Springer, Berlin, Heidelberg, 2002.

SWAG

モノイド \( (S, \circ, \mathrm{id}_{\circ})\) がある。これに関する以下のクエリを \(n\) 個処理せよ。

  • キューを空で初期化。\(q\gets ()\)。
  • 与えられた \(x\in S\) をキューの末尾に追加。\(q \gets q \concat (x)\)。
  • キューの先頭を削除。\(q = (q_0, q_1, \dots, q_{n-1})\) のとき \(q \gets (q_1, \dots, q_{n-1})\)。
  • キューの \(\circ\)-fold を出力。\(q = (q_0, q_1, \dots, q_{n-1})\) のとき \(\mathrm{id}_{\circ} \circ q_0 \circ q_1 \circ \dots \circ q_{n-1}\)。

これは、two-stack queue と呼ばれる、スタック二つを用いてキューを実装する手法で、\(\circ\)-fold の値も(累積和のように)管理することで、ならし \(O(1)\) 時間で可能です。代わりに two-stack deque を使うことで、両方向の出し入れが可能です。

蟻本にある、いわゆる「スライド最小値」より広い問題設定に対応できます*6

とはいえ、蟻本のスライド最小値(と、その直前にあるスタックによる最大長方形)で使うような、大小関係を利用する発想をできるようになっておくのも大事そうです。えびちゃんは苦手です。

SWAG (sliding window aggregation) という呼称は、競プロ界隈では上記の問題を解くデータ構造の意味合いで使われがちですが、下記の文献などを見た感じだと、もっと広い問題設定(あるいはデータ構造のインタフェース?)の意味合いで使われていそうです。

Hirzel, Martin, Scott Schneider, and Kanat Tangwongsan. “Sliding-window aggregation algorithms: Tutorial.” In Proceedings of the 11th ACM International Conference on Distributed and Event-based Systems, pp. 11–14. 2017.

Gray code の利用

\(n\) 要素の集合に対して、ビット全探索で \(2^n\) 通りの集合を列挙して各々に対して \(\Theta(n)\) 時間で処理をするケースは頻出です。 このとき、\(\Theta(n\cdot 2^n) = \Theta(2^n\log(2^n))\) 時間かかっています。

ところで、\(2^n\) 通りの集合を列挙する順序のうち、以下を満たすようなものが常に存在し、かつ簡単に得ることができます。

  • \(\emptyset\) から始まる。
  • 最後に出力した集合に、1 つの要素を追加または削除することで次の集合を得る。

たとえば、3 要素であれば 000001011010110111101100 といった具合です。

元々 \(\Theta(n)\) 時間かけていた処理が、差分のみの更新が \(O(1)\) 時間で可能なものであれば、上記の列挙順を利用することで \(\Theta(2^n)\) 時間になります。

±1 配列の転倒数

長さ \(n\) の配列 \(a = (a_0, a_1, \dots, a_{n-1})\) が与えられる。配列 \(a\) の転倒数 \(|\{(i, j) \mid i \lt j \wedge a_i \gt a_j\}|\) を答えよ。

ただし、\(a_i \in\N\) であり、\(a_0 = 0\) かつ \(|a_i - a_{i+1}| \le 1\) とする。

一般の配列 \(a\) に対する転倒数は、平面走査なり分割統治なりで \(O(n \log(n))\) 時間で求められ、有名です。

上記の制約がある場合は、maspy さんの記事 のような方法で \(O(n)\) 時間にできます。

0/1 配列に対する転倒数は、\(w\)-bit 整数に \(w\) 個の要素を詰めて与えられるとすると、\(O(n/w\cdot \log(w)^2)\) 時間でも解けそうです。もしかしたらもっと速くなる? わかりません。

rsk0315.hatenablog.com

周期検出

初期値 \(x_0 \in S\) (\(|S| \lt \infty\)) と関数 \(f: S\to S\) が与えられる。ただし、\(f\) は \(O(1)\) 時間で計算できる設定とする。

\(i \gt 0\) に対して \(x_i = f(x_{i-1})\) で定義するとき、ある \(\lambda\gt 0\) に対して \(x_{\mu} = x_{\mu + \lambda}\) と書けるような最小の \(\mu\ge 0\) と、それを満たす最小の \(\lambda\) を求めよ。

\(n = \mu+\lambda\) とおきます。

連想配列を用いると \(\Theta(n\log(n))\) 時間などで、空間も \(\Theta(n)\) になります。

tortoise and hare algorithm を用いることで、\(O(n)\) 時間、\(O(1)\) 空間で可能です。前述の majority voting のように、同値判定だけで済むのもおもしろいです。

周期検出 2

\(n\) 頂点の functional graph、すなわち \(n\) 本の有向辺 \(i, f(i)\) (\(0\le i\lt n\)) があるグラフが与えられる。 このとき、各 \(0\le i\lt n\) について、\(x_0 = i\) としたときの上記の問題の \( (\mu, \lambda)\) を求めよ。

がちゃがちゃやると \(O(n)\) 時間でできます。あんまりおもしろくはないかも。

その他メモ

同じ値で \(k\) 回の二分探索を行う際、愚直には \(\Theta(k \log(n))\) 時間になるが、fractional cascading を使うと \(O(k+\log(n))\) 時間に落とすことができる。

https://cs.brown.edu/courses/cs252/misc/resources/lectures/pdf/notes08.pdf

配る DP で区間加算が必要になるとき、双対セグ木を使わずに imos 法を on-the-fly で更新していくことで log をつけずにできる。えびちゃんはよくバグらせる。かわいそう。

union-find は十分な回数クエリするとならし \(O(1)\) 時間になるので、そういう場合は BFS するのと比べて損せずに済む。

文字列アルゴリズム界隈や簡潔データ構造界隈、線形時間前処理かつ定数時間クエリじゃないと気が済まないみたいな印象がある。

yosupo.hatenablog.com

CHT で仮定を置くと log が消えるやつ。こっちの方が有名そうな印象はある。

\(n\) 以下の逆元の列挙。最近あんまり話題にならないね。

Fibonacci heap を使うと Dijkstra 法や Prim 法の \(O(m\log(n))\) が \(O(m+n\log(n))\) になるやつ。

個数制限ありナップサック問題を log なしで解くやつ。最近 ABC でそんな感じでやっても log がつかない(?)やつがあったような?

なんか忘れてるかも。思い出したら書く。

おわり

いくつ実装したことがありましたか? いろいろな log なしアルゴリズムを実装していろいろな知見を得ましょう〜〜。

log ありで困ることがあるかと言われると微妙な感はありますが、log を消すアルゴリズムたちは(競プロであまり見慣れていないせいもあってか)天才に感じることが多く、面白いがちです。分割統治や MoFR などはそれ自体は典型的ですが、それらが使える形に変形するのは天才がちな気がします。

「これが log なしで解けるのか〜〜」という感情が生えるきっかけになれたらうれしいです。

ところで、\(n \log(n)\) が \(n\) になると log が消えた感があるのに、\(n\) が \(n/\log(n)\) になると log が消えた感はないですね。\(\log(n)\) で割れる系アルゴリズムもあまり有名でない感じがします。あるいは \(\log(\log(n))\) 系のものもあまり有名でない気がしますね。

*1:尺取り法で解けても二分探索が楽ならそうしてしまいがち。

*2:これがあるとソートが \(O(n)\) でできてしまい下界に反してしまういつものやつ。

*3:定数時間で答えられないものは RMQ ではないという思想を感じる命名

*4:\(\floor{n/2}\) としておくと、順序は定義されているが平均が定義されていない型でも扱いやすい。

*5:\(o(n)\) のこと。

*6:「スライド最小値」は特定の実装ではなく問題設定に関する名前だと思われますが、蟻本にある deque による実装を指して言う人もいそうです。

大まかな方針はわかっているはずなのに細かい実装でバグらせるのかわいそう

えびちゃんかわいそう。

\(O(1)\) 階差分を取ると imos 法で書ける形になるような区間加算クエリ*1を一生バグらせました。 特に、DP をしながら on-the-fly で imos 法の遅延を解消していくやつを毎回バグらせている気がします。

こういうのはコンテストで出たときにその場でどうにかしようとすると頭が間に合わなくなるので、コンテストのない間に整理しておくのがよさそうです。

null-mn.hatenablog.com

この手のセグ木を貼れば済むというのに後から冷静になって気づいたのですが、セグ木を貼れてもえびちゃんの頭がこわれていることの解決にはなりません。あと線形時間で解きたいです。

todo: 整理できるか、別のものでバグらせてかわいそうになったら更新します。

cf. 昔の記事 rsk0315.hatenablog.com

*1:区間に等差数列を足す、のようなものです。

ワードサイズの転倒数を O(log(w)²) 時間で求めるよ

これいる?

問題設定

\(w\)-bit 整数を 0/1 配列として見たときの転倒数を考えます。

すなわち、\(x = b_0 \cdot 2^0 + b_1 \cdot 2^1 + \dots + b_{w-1} \cdot 2^{w-1} = b_{w-1}\ldots b_1 {b_0}_{(2)}\) (\(b_i \in \{0, 1\}\)) として、\(i\lt j\) かつ \((b_i, b_j) = (1, 0)\) なる \((i, j)\) の個数を求めたいです。

たとえば、\(w = 16\) で \(x = 0010{,}0111{,}0110{,}0101_{(2)}\) とすると、求める転倒数は \(39\) です(右側が下位ビットであることに注意)。

方針

ビット並列手法の常套手段を使います。

rsk0315.hatenablog.com

転倒数は、平面走査ベースではなく分割統治で求める方法を使います。

atcoder.jp

一般の場合の 解説スライド

ここでは値が 0/1 だけなのでもっと簡単。実装例

すなわち、マージするときに、“左の子の 1 の個数” と “右の子の 0 の個数” の積を求めます。 なので、統治しながら 0/1 の個数を管理しつつ、それらの積も計算したいです。

解法

上記のビット並列の記事中の count ones 同様に、0/1 の個数を計算します。

x = abcd efgh ijkl mnop, !x = ABCD EFGH IJKL MNOP として、

abcd efgh ijkl mnop
ABCD EFGH IJKL MNOP

下位桁の 1 と上位桁の 0 の個数の積を求めるので、

    0b0d 0f0h 0j0l 0n0p
.*) 0A0C 0E0G 0I0K 0M0O
-----------------------
    0?0? 0?0? 0?0? 0?0?

のようになります。ここで、.* は上下に対応するブロックごとの積とします。たとえば、左の ? から順に b * A, d * C, ..., p * O の値が入ります。 この答え 0?0? 0?0? 0?0? 0?0? も count ones 同様に隣同士で足し合わせておきます。

続いて、count ones 同様にサイズ 2 のブロック内での 0/1 の個数がわかり、それらの積を求めたいです。

    00cd 00gh 00kl 00op
.*) 00AB 00EF 00IJ 00MN
-----------------------
    0??? 0??? 0??? 0???

同様のことを繰り返します。

    0000 efgh 0000 mnop
.*) 0000 ABCD 0000 IJKL
-----------------------
    0??? ???? 0??? ????
    0000 0000 ijkl mnop
.*) 0000 0000 ABCD EFGH
-----------------------
    0??? ???? ???? ????

これら各段での ???..? の総和が答えになります。

問題は、block-wise product .* です。ここが(SIMD 演算などで(あるいは、そういう基本命令を仮定することで))、\(O(1)\) 時間でできるのであれば、全体で \(O(\log(w))\) 時間になります。

たとえば pmullw という命令を使うと(特定のビット幅においては)できそうです。

この問題においては、0/1 の個数の最大値が \(w\) であることから、\(O(\log(w))\)-bit 程度で表せることを使うと、(乗算の筆算アルゴリズムを並列に行うことで)\(O(\log(w))\) 時間で計算でき、全体で \(O(\log(w)^2)\) 時間にできます。

例を示します。

    0000 0abc 0000 0def
.*) 0000 0101 0000 0110
-----------------------
    0000 0abc 0000 0000  # (0abc * 1, 0def * 0) << 0
    0000 0000 0000 def0  # (0abc * 0, 0def * 1) << 1
    000a bc00 000d ef00  # (0abc * 1, 0def * 1) << 2

0/1 に応じて 0abc/0000 なり 0def/0000 なりを並列に選んでくるのは、前述のビット並列の記事にあるような手法で \(O(1)\) 時間ででき、それを \(O(\log(w))\) 回足すので \(O(\log(w))\) 時間ということですね。

もしかしたら賢くやれば普通に \(O(1)\) 時間で .* を求められるかもしれません?

結局、コードとしては次のような感じになります。

fn main() {
    let x = 0x_6A6A6A12_BC4441D8_AA0EA523_D52ED8DC_u128;
    assert_eq!(inversion_u128(x), 2187);
}

const MASK_1_U128: u128 = 0x_5555_5555_5555_5555_5555_5555_5555_5555;
const MASK_2_U128: u128 = 0x_3333_3333_3333_3333_3333_3333_3333_3333;
const MASK_4_U128: u128 = 0x_0F0F_0F0F_0F0F_0F0F_0F0F_0F0F_0F0F_0F0F;
const MASK_8_U128: u128 = 0x_00FF_00FF_00FF_00FF_00FF_00FF_00FF_00FF;
const MASK_16_U128: u128 = 0x_0000_FFFF_0000_FFFF_0000_FFFF_0000_FFFF;
const MASK_32_U128: u128 = 0x_0000_0000_FFFF_FFFF_0000_0000_FFFF_FFFF;
const MASK_64_U128: u128 = 0x_0000_0000_0000_0000_FFFF_FFFF_FFFF_FFFF;

/// 128 要素の 0/1 配列の転倒数を求める。
fn inversion_u128(x: u128) -> u32 {
    // 各ブロックの転倒数が入る。
    let mut res = 0;
    // 各ブロックの 0 の個数が入る。
    let mut zero = !x;
    // 各ブロックの 1 の個数が入る。
    let mut one = x;

    let mask = [
        MASK_1_U128,
        MASK_2_U128,
        MASK_4_U128,
        MASK_8_U128,
        MASK_16_U128,
        MASK_32_U128,
        MASK_64_U128,
    ];
    for i in 0..7 {
        let m = mask[i];
        let w = 1 << i;
        if i > 0 {
            // 転倒数を隣同士でくっつけ、ブロックサイズを大きくしておく。
            res = (res >> w & m) + (res & m);
        }
        let bit = (2 * i + 1) as u32;
        // 上位の 0 と下位の 1 の個数の積 (.*) を足す。
        res += product_u128(zero >> w & m, one & m, w << 1, bit);
        // 0/1 の個数のブロックサイズを大きくする。
        one = (one >> w & m) + (one & m);
        zero = (zero >> w & m) + (zero & m);
    }
    res as u32
}

/// `lhs .* rhs` を求める。ただし、`lhs`, `rhs` は `bit`-bit 整数で、
/// サイズ `block` のブロックごとに分けられているとする。
/// `2 * bit - 1 <= block` を前提とする。
/// ここでは、block = [2, 4, 8, 16, ...] のとき bit = [1, 3, 5, 7, ...] である。
fn product_u128(lhs: u128, rhs: u128, block: u32, bit: u32) -> u128 {
    let mask = match block {
        2 => return lhs & rhs,  // 1 桁の積は bit-wise and でよい。
        4 => 0x_1111_1111_1111_1111_1111_1111_1111_1111,
        8 => 0x_0101_0101_0101_0101_0101_0101_0101_0101,
        16 => 0x_0001_0001_0001_0001_0001_0001_0001_0001,
        32 => 0x_0000_0001_0000_0001_0000_0001_0000_0001,
        64 => 0x_0000_0000_0000_0001_0000_0000_0000_0001,
        128 => return lhs * rhs,  // ブロックが一つなので、普通の積でよい。
        _ => unreachable!(),  // 上記以外は呼ばないとする。
    };

    let mut res = 0;
    for i in 0..bit {
        // `rhs & (mask << i)` で、乗数の 0/1 を取得し、`(それ >> i)` に対して
        // `!(!0 << block)`、すなわち `block` 個の連続する `1` を掛けることで、
        // 必要なマスクを得る。そのマスクを `lhs` にかけ、`res` に足していけばよい。
        let cur = ((rhs & (mask << i)) >> i) * !(!0 << block);
        res += (cur & lhs) << i;
    }
    res
}

Playground

感想

思いついたので書いてみたら一発でテストが通ってびっくりしました。愚直が間違ってなかったらいいなって思います。

その他

サイズ \(n\) の 0/1 配列を \(n/w\) 個のワードとして扱って転倒数を求めるとします。 各ワードの転倒数が \(O(n/w\cdot\log(w)^2)\) 時間で求まります。 また、今見ているワードより左にある 1 の個数の総和と、今見ているワードにある 0 の個数の総和の積を足していきますが、これは(0/1 の個数を求めるパートを含めて)\(O(n/w\cdot\log(w))\) 時間で求まります。

よって、全体の転倒数が \(O(n/w\cdot\log(w)^2) = o(n)\) 時間で求められます。

以下にある ±1 配列の転倒数と同じ方法でも \(O(n)\) 時間で求められますが、0/1 配列だともっと高速になるということですね。

maspypy.com

.* パートが \(O(1)\) になるとか別の考察で高速になるとかがあると、\(\log(w)\) ひとつぶんくらいはさらに速くなるかもしれません。

おわり

ねむいです。

ビット並列手法の基礎と基礎以外など

思い立ったので書いてみます。

導入

たとえば、次のような問題を考えます。

64-bit 整数 n が与えられます。n を 2 進法で表した際の 1 の個数を数えてください。

愚直には、ループを回して 1 ビットずつ確認すればよいです。

fn count_ones(n: u64) -> u32 {
    let mut res = 0;
    for i in 0..64 {
        if n >> i & 1 == 1 {
            res += 1;
        }
    }
    res
}
// Rust においては n.count_ones() というメソッドが用意されているが、今は触れない
// GCC における __builtin_popcountll() などに相当する処理

しかし、各ループにおいて(64-bit 整数は 64 個のビットをまとめて演算できるにもかかわらず)1 ビットずつしか見ていないのはもったいないです。 そこで、複数ビットをまとめて扱うことで高速化を図る手法があり、ビット並列 (bit-level parallelism) などと呼ばれています。

前提

\(w\)-bit 整数における算術演算 (+, -, *, /, %) とビット演算 (&, |, ^, <<, >>, !) と比較演算 (==, !=, <, <=, >, >=) を \(O(1)\) 時間で行えるものとします(/ は切り捨て除算、! はビット反転、>> は論理シフト。オーバーフローは \(2^w\) を法として計算)。 また、\(w\)-bit 整数を用いてメモリアクセス(配列の添字として使う)も \(O(1)\) 時間でできる(ランダムアクセス可能、という言い方をする)とします。 さらに、コンパイル時に決まっている定数を読み出すのも \(O(1)\) 時間でできるとします。

上記の \(w\) のことを word size と呼びます。 word size に収まる整数 \(x\in\{0, 1, \dots, 2^w-1\}\) のことを word と呼んだりします。 word に対するこれらの演算を \(O(1)\) 時間でできるとする計算モデルを word RAM model と呼んだりします。

たとえば 64-bit 整数の文脈においては \(w = 64\) となりますが、\(w = O(1)\) という意味ではないことに注意が必要です。

word size が定数でないという話に関して

入力サイズを \(n\) とします。word が \(n\) 個ある配列が与えられるのに相当しますが、これに対して \(w\)-bit 整数でランダムアクセスをするためには、(\(w\)-bit 整数で表せる個数が \(2^w\) 個なので)\(2^w \ge n\) である必要があり、\(w \ge \log_2(n)\) となります。

対象の問題の入力サイズに応じて word size が変わるのは不自然に感じるかもしれませんが、メモリが 128 KB しかないようなマシンで \(10^9\) 要素ある問題を解くのは無理があるというようなイメージをするといいかもしれません?

というわけで、入力サイズに応じて word size も大きくなるため、\(w = O(1)\) ではありません。オーダーとしては \(w = \Omega(\log(n))\) となりますね。

もう少し詳しくは以下の記事で触れています。

rsk0315.hatenablog.com

なので、この記事の話は定数倍高速化という意味合いではなく、\(\Omega(\log(n))\) 倍高速化などと見なせると思っています。 \(O(n^2/w) = O(n^2/\log(n))\) とかです。

以下、演算子の優先順位などは Rust の気持ちで書いているので、x >> i & 1 == 0 などと書いていますが、C++ に移植する際は (x >> i & 1) == 0 などとする必要があります。

基礎演算

note: word は \(\{0, 1, \dots, 2^{w-1}\}\) の符号なし整数で、オーバーフローの際は \(2^w\) を法として合同な値になるようにします。下 \(w\) bits のみに切り捨てると思ってもよいですね。

単項マイナス

x に対して単項演算 -x を考えます。-x == 0 - xx + (-x) == 0 が成り立ってほしいです。

これは、!x + 1 と等しくなります。 x + !x == 111...1 となり、これに 1 を足すと 0 になるためですね*1

!x + 1 が欲しくなったときに -x に置き換えてよい(逆も然り)ことを覚えておくと吉です。

最も右の諸々に関する演算

足し算において、答えのうちのある桁が(繰り上がりによって)右の桁の影響を受けることはありますが、左の桁の影響を受けることはありません。 こうした事情から、ある演算を足し算などの組み合わせで再現したいとき、左の桁の影響を受けるような演算を再現することはできません。

そこで、右側に関する演算を考えます。 たとえば、「最も右にある 1 のみからなる整数 (0100110000000100) を得る」のような演算は、「自分より右に 1 があるかどうか」で決まり、左の桁の影響は受けないので計算できそう、といった具合です。「最も左にある 1 のみからなる整数 (0100110001000000) を得る」は、「自分より左に 1 があるかどうか」に影響されるので無理そうです。

やや天下り的ですが、x に対して、演算 & | ^ と引数 x - 1 x + 1 -x の演算表を考えてみます。

具体例

x = 0101_1100

x - 1 x + 1 -x
0101_1011 0101_1101 1010_0100
x & 0101_1000 0101_1100 0000_0100
x | 0101_1111 0101_1101 1111_1100
x ^ 0000_0111 0000_0001 1111_1000

x = 1010_0011

x - 1 x + 1 -x
1010_0010 1010_0100 0101_1101
x & 1010_0010 1010_0000 0000_0001
x | 1010_0011 1010_0111 1111_1111
x ^ 0000_0001 0000_0111 1111_1110

x = 0000_0000

x - 1 x + 1 -x
1111_1111 0000_0001 0000_0000
x & 0000_0000 0000_0000 0000_0000
x | 1111_1111 0000_0001 0000_0000
x ^ 1111_1111 0000_0001 0000_0000

x = 1111_1111

x - 1 x + 1 -x
1111_1110 0000_0000 0000_0001
x & 1111_1110 0000_0000 0000_0001
x | 1111_1111 1111_1111 1111_1111
x ^ 0000_0001 1111_1111 1111_1110

x - 1 x + 1 -x
x & 最も右にある 1 を 0 にする 最右桁から続く 1 を 0 にする 最も右にある 1 以外を 0 にする
x | 最右桁から続く 0 を 1 にする 最も右にある 0 を 1 にする 最も右にある 1 より左全てを 1 にする
x ^ 最も右にある 1 より左全てを 0、右全てを 1 にする 最も右にある 0 を 1 に、それより左を 0 にする 最も右にある 1 を 0 に、それより左を 1 にする
  • exercise: x & (!x - 1) x | (!x - 1) x ^ (!x - 1) について考えてみよ。
  • exercise: x が 2 べき(ある \(k\) が存在して \(2^k\) と表せる)かを判定する式を考えてみよ。x == 0 は 2 べきでないことに注意せよ。

その他簡単な演算

x において連続する 1 が存在するか?というのは、x & (x >> 1) != 0 などで判定することができます。遷移が特殊な bit DP とかで役に立つかもしれません。

典型的な発想の紹介

以下、例のために 64 個のビットを書くのは大変なので、\(w = 16\) くらいとします。また、4 桁ごとに _ で区切って書くとし、2 進法であるとします。

にこにこで合わせるやつ(分割統治)

例を挙げながら説明します。

count ones

冒頭の例です。1 の個数を数えます。

例として、入力を x = 0010_1011_1100_0111とします。答えは 9 ということになります。

まず、これを 2 ビットずつに区切ったときの 1 の個数を数えます。区切ると [00, 10, 10, 11, 11, 00, 01, 11] なので、[0, 1, 1, 2, 2, 0, 1, 2] が欲しいということです。 2 ビットごとの左側と右側に分けた word を用意します。

0010_1011_1100_0111  # x
-------------------
0010_1010_1000_0010  # x & 1010_1010_1010_1010
0000_0001_0100_0101  # x & 0101_0101_0101_0101

これの桁を揃えて足します。

   0001_0101_0100_0001  # (x & 1010_1010_1010_1010) >> 1
+) 0000_0001_0100_0101  # x & 0101_0101_0101_0101
----------------------
   0001_0110_1000_0110  # [0, 1, 1, 2, 2, 0, 1, 2]

これにより、欲しかった値が得られました。この値を改めて x と置きます。

同様の処理を繰り返し、今度は 4 ビットずつに区切ったときの 1 の個数を数えます。 [0, 1, 1, 2, 2, 0, 1, 2] の隣同士を合わせればよく、[1, 3, 2, 3] が欲しい値です。

0001_0110_1000_0110  # x
-------------------
0000_0100_1000_0100  # x & 1100_1100_1100_1100
0001_0010_0000_0010  # x & 0011_0011_0011_0011

また揃えて足します。

   0000_0001_0010_0001  # (x & 1100_1100_1100_1100) >> 2
+) 0001_0010_0000_0010  # x & 0011_0011_0011_0011
----------------------
   0001_0011_0010_0011  # [1, 3, 2, 3]

同様に繰り返します。

   0000_0001_0000_0010  # (x & 1111_0000_1111_0000) >> 4
+) 0000_0011_0000_0011  # x & 0000_1111_0000_1111
----------------------
   0000_0100_0000_0101  # [4, 5]
   0000_0000_0000_0100  # (x & 1111_1111_0000_0000) >> 8
+) 0000_0000_0000_0101  # x & 0000_0000_1111_1111
----------------------
   0000_0000_0000_1001  # [9]

上記により、目標通り 9 が得られました。

この枠組みで言えば、元々の入力は「1 ビットごとに区切ったときの 1 の個数」を表していると見ることもできますね。

fn popcount(mut x: u16) -> u32 {
    x = ((x & 0xAAAA) >> 1) + (x & 0x5555);
    x = ((x & 0xCCCC) >> 2) + (x & 0x3333);
    x = ((x & 0xF0F0) >> 4) + (x & 0x0F0F);
    x = ((x & 0xFF00) >> 8) + (x & 0x00FF);
    x as u32
}

計算量は、\(O(\log(w))\) 時間です。

64-bit の場合・高速化

以下のようなコードが考えられます。

fn popcount(mut x: u64) -> u32 {
    x -= (x >> 1) & 0x_5555_5555_5555_5555;
    x = (x & 0x_3333_3333_3333_3333) + ((x >> 2) & 0x_3333_3333_3333_3333);
    x = (x + (x >> 4)) & 0x_0F0F_0F0F_0F0F_0F0F;
    x += x >> 8;
    x += x >> 16;
    x += x >> 32;
    (x & 0x7F) as u32
}
  • x -= (x >> 1) & 0x_5555_5555_5555_5555;
    • 元々 x == (x & 0x_AAAA_AAAA_AAAA_AAAA) + (x & 0x_5555_5555_5555_5555)
    • すなわち x == (((x >> 1) & 0x_5555_5555_5555_5555) << 1) + (x & 0x_5555_5555_5555_5555)
      • (x & 0x_AAAA_AAAA_AAAA_AAAA) >> 1 == (x >> 1) & 0x_5555_5555_5555_5555 なので
    • 欲しいのは ((x >> 1) & 0x_5555_5555_5555_5555) + (x & 0x_5555_5555_5555_5555)
    • なので差分の (x >> 1) & 0x_5555_5555_5555_5555 を引けばよい
      • a << 1a * 2 a + a と等しいので
  • x = (x & 0x_3333_3333_3333_3333) + ((x >> 2) & 0x_3333_3333_3333_3333)
    • これは & より先に >> をするように変形しただけ
  • x = (x + (x >> 4)) & 0x_0F0F_0F0F_0F0F_0F0F
    • (x & 0x_0F0F_0F0F_0F0F_0F0F) + ((x >> 4) & 0x_0F0F_0F0F_0F0F_0F0F)
    • この時点で 4 つブロックの個数を数えているので、最大でも x == 0x_4444_4444_4444_4444
    • 4 + 4 == 8 < 16 なので、繰り上がっても左の桁に影響を与えないため、マスクを取るのは計算後 1 回のみでよい
  • x += x >> 8 x += x >> 16 x += x >> 32
    • 最終的な答えが高々 64(かつ各ブロックはそれ未満)なので、7-bit ぶん以上離れたところに影響せず、マスクを取らなくても平気
  • x & 0x7F
    • 7-bit だけ得ればよい

reverse

上位ビットと下位ビットを逆順に並び換えます。

例として x = 0100_1101_0110_0001 とすると、1000_0110_1011_0010 が欲しい値です。

隣り合う二つを入れ換えます。

0100_1101_0110_0001  # x
-------------------
0000_1000_0010_0000  # x & 1010_1010_1010_1010
0100_0101_0100_0001  # x & 0101_0101_0101_0101
   0000_0100_0001_0000  # (x & 1010_1010_1010_1010) >> 1
|) 1000_1010_1000_0010  # (x & 0101_0101_0101_0101) << 1
----------------------
   1000_1110_1001_0010

or (|) の代わりに + を使ってもよさそうです。同様に、隣り合う 2 ビットのグループを入れ換えます。

1000_1100_1000_0000  # x & 1100_1100_1100_1100
0000_0010_0001_0010  # x & 0011_0011_0011_0011
   0010_0011_0010_0000  # (x & 1100_1100_1100_1100) >> 2
|) 0000_1000_0100_1000  # (x & 0011_0011_0011_0011) << 2
----------------------
   0010_1011_0110_1000

同様に 4 ビットで続けます。

0010_0000_0110_0000  # x & 1111_0000_1111_0000
0000_1011_0000_1000  # x & 0000_1111_0000_1111
   0000_0010_0000_0110 # (x & 1111_0000_1111_0000) >> 4
|) 1011_0000_1000_0000  # (x & 0000_1111_0000_1111) << 4
----------------------
   1011_0010_1000_0110

8 ビットで同様に。

1011_0010_0000_0000  # x & 1111_1111_0000_0000
0000_0000_1000_0110  # x & 0000_0000_1111_1111
   0000_0000_1011_0010  # (x & 1111_1111_0000_0000) >> 8
|) 1000_0110_0000_0000  # (x & 0000_0000_1111_1111) << 8
----------------------
   1000_0110_1011_0010

1000_0110_1011_0010 が得られました。めでたい。

suffix parity

パリティというのは、あるビット列のうち 1 の個数の偶奇を表す値です*2。 たとえば、0011_0101 は偶数個の 1 があるのでパリティ01011_1100パリティ1 です。

さて、各ビットに対して、自分より右(自分含む)のパリティを求めることを考えます。 累積和の xor 版を各ビットについて行うと言った方がわかりやすいかもしれません。

以下のようなビット列を考えます(a, b, ... などは 01 を表すとします)。

abcd efgh ijkl mnop  # x

これに対し、x << 1 との xor を考えます。

abcd efgh ijkl mnop  # x
bcde fghi jklm nop0  # (x << 1)

これを改めて x とおき、x ^ (x << 2) を考えます。

abcd efgh ijkl mnop  # x
bcde fghi jklm nop0  # ...x
cdef ghij klmn op00  # (x << 2)
defg hijk lmno p000  # ...(x << 2)

縦に見ると abcd... となっていることがわかります。同様に繰り返すことで全体の suffix parity が得られます。

fn suffix_parity(mut x: u16) -> u16 {
    x ^= x << 1;
    x ^= x << 2;
    x ^= x << 4;
    x ^= x << 8;
    x
}

シフトの方向を逆にすることで、prefix parity も計算できます。

二分探索

most-significant set bit (msb)

一番左の 1 の位置を求めます。

例として、入力を x = 0000_0010_1101_0110 とします。x >> 9 == 1 となることから、9 が答えです。

コーナーケースについて

入力が 0 の場合には、状況に応じて適当に処理することにします。 「どこまで左に行っても 1 が現れない」ということで \(\infty\) や \(w\) を返したり、「x >> (i + 1) == 0 を満たす最小の非負整数 i」と定義することで 0 を返したり、あるいは例外を投げたり、未定義ということにしたり。

0 でない場合の返り値は \(0\) 以上 \(w\) 未満であることも考慮しておくとよいかもしれません。

答えを 0 で初期化しておきます。 まず、左半分の 8 ビットに 1 があるかを調べます。

0000_0000_0000_0010  # x >> 8

x >> 8 != 0 なので、最上位の 1 より右に 8 つ以上のビットがあることになります。 なので、答えに 8 を足して、次は x >> 8 の msb を求めます。

この値を改めて x とし、該当の 8 ビットのうちの左半分の 4 ビットに 1 があるかを調べます。

0000_0000_0000_0000  # x >> 4

x >> 4 == 0 なので、最上位の 1 より右には 4 つ未満のビットしかないことになります。 同様に繰り返します。

0000_0000_0000_0000  # x >> 2

x >> 2 == 0 なので 2 つ未満しかなし。

0000_0000_0000_0001  # x >> 1

x >> 1 != 0 なので 1 つ以上あり。答えに 1 を足す。

よって、答えは 8 + 1 == 9 とわかりました。

least-significant set bit (lsb)

最も右の 1 の位置を求めます。

x のうち、最も右のビットのみからなる整数は x & -x で得られます。なので、これに対して msb を求めればよいです。

de Bruijn sequence の利用

de Bruijn sequence というのがあります*3。 種類数 \(k\) の文字集合 \(\Sigma\) における order \(n\) の de Bruijn sequence とは、長さ \(n\) の \(\Sigma\) 上の文字列が全て(循環してもよい)部分文字列として現れる文字列のことです。長さ \(k^n\) のものが常に存在することが 知られています

たとえば、\(\Sigma = \{0, 1\}\) における order \(4\) の de Bruijn sequence として、0000100110101111 があります*4。 左から順に 0000, 0001, 0010, 0100, 1001, 0011, 0110, 1101, 1010, 0101, 1011, 0111, 1111, 1110, 1100, 1000 です。

exact log

\(x = 2^k\) が与えられたときに \(k\) を返す(\(x\) が 2 べきでない入力は未定義)関数を考えます。

de Bruijn sequence を \(x\) 倍したときの上位 \(\log(w)\)-bit を見ることを考えます。 たとえば、x == 128 のとき (0b_0000_1001_1010_1111 * 128) >> (16 - 4) == 0b_1101 == 13 です。 * x<< k と同じであることや de Bruijn sequence の性質に注意すると、(\(w\) 種類ある)x の値によって (0b_0000_1001_1010_1111 * x) >> 12 の値は相異なることがわかります。

そこで、長さ \(w\) の配列を用意して上記の対応を覚えておくことで、\(O(1)\) 回の演算で求めることができます。 上記の例では、\(128 = 2^7\) なので a[13] = 7 となるような配列 a ということです。

ここで利用する de Bruijn sequence ですが、左シフトしたときには右側が 0 で埋められる関係で、上位桁が十分な個数の 0 から始まっている必要があります。たとえば 0b_1111_0101_1001_0000 を使うと、01110011 などが出現せず 0000 が複数回出現してしまいます。 その点にさえ注意すればどれを使ってもよいですが、使う列によって用いる配列が異なるので注意しましょう。

値が小さい場合の計算

\(w\)-bit の整数は、「\(\sqrt{w}\)-bit の整数が \(\sqrt{w}\) 個ある」と見なすこともできます。ここで、\(w\) は平方数だと仮定したりします*5。 あるいは、\(i\)-bit 整数が \(\floor{w/i}\) 個あるとも見なせますね。

一つの整数を小さい整数の配列と見ることでまとめて演算したり、一つの小さい整数を複数個にコピー(一つの整数に格納)してまとめて演算したりする手法があります。

左シフトと加算を複数回行う代わりに、一回の乗算でまとめて行う手法がよく使われます。ここの章にあるものは \(O(1)\) time のものばかりです。

distribute

4-bit 整数 abcd を入力とします(a, b, c, d0, 1 のいずれか)。

0000_0000_0000_abcd  # x

これを、以下の形に変換したいです。

abcd_abcd_abcd_abcd

これは、以下のように計算できます。

   0000_0000_0000_abcd  # x
   0000_0000_abcd_0000  # x << 4
   0000_abcd_0000_0000  # x << 8
+) abcd_0000_0000_0000  # x << 12
----------------------
   abcd_abcd_abcd_abcd

数式で書くと \(x + x \cdot 2^4 + x \cdot 2^8 + x \cdot 2^{12}\) であり、\(x \cdot (2^0 + 2^4 + 2^8 + 2^{12})\) となります。 よって、以下の乗算一回の形で計算できます。

   0000_0000_0000_abcd  # x
*) 0001_0001_0001_0001  # (1 << 0) | (1 << 4) | (1 << 8) | (1 << 12)
----------------------
   abcd_abcd_abcd_abcd

これにより、\(\sqrt{w}\)-bit の整数を \(\sqrt{w}\) 個に分配することなどができます。これが基本的な演算・考え方となります。

count ones

distribute された状態で入力が与えられるとします。

abcd_abcd_abcd_abcd

これに対して、各ブロックごとに相異なる一つだけが取り出されるようにマスクをかけます。

   abcd_abcd_abcd_abcd
&) 0001_0010_0100_1000
----------------------
   000d_00c0_0b00_a000

さらに、これらの桁を揃えて足します。

   000d_00c0_0b00_a000  # x
   000c_00b0_0a00_0000  # x << 3
   000b_00a0_0000_0000  # x << 6
+) 000a_0000_0000_0000  # x << 9

これは、先ほど同様に乗算に帰着できて、x * 0b_0000_0010_0100_1001 とできます。あとは、>> 12 をして桁を揃えれば終わりです。 \(O(\sqrt{w})\)-bit の整数であれば、popcount (count ones) が \(O(1)\) 時間でできるということです。

is-zero

\(\sqrt{w}\) 個 の \(O(\sqrt{w})\)-bit 整数 abcd efgh ijkl mnop が与えられます。

abcd_efgh_ijkl_mnop  # x

これの各々に対し、0 であれば 0000、それ以外であれば 0001 を返す関数を考えます。たとえば、0010_1011_0000_1000 であれば 0001_0001_0000_0001 を返します。

各々の最上位ビットとそれ以外に分けておきます。

a000_e000_i000_m000  # x_hi
0bcd_0fgh_0jkl_0nop  # x_lo

x_hi についてはほぼできあがっているので、x_lo について考えます。 0111_0111_0111_0111 との足し算をし、繰り上がり(キャリー)を見ればよいです。

先の例であれば

   0010_0011_0000_0000  # x_lo
+) 0111_0111_0111_0111
----------------------
   1001_1010_0111_0111

これと 1000_1000_1000_1000 の and を取り、x_hi との or を取ればよいですね。

   1000_1000_0000_0000  # x_lo & 1000_1000_1000_1000
|) 0000_1000_0000_1000  # x_hi
----------------------
   1000_1000_0000_1000

あとは >> 3 とかをして桁を揃えましょう。

from bits

\(\sqrt{w}\) 個のビット a b c d が与えられます。

000a_000b_000c_000d

たとえば上記の is-zero の返り値などがこの形式です。ここから一つの \(\sqrt{w}\)-bit 整数 abcd を作りたいです。

distribute と同様です。

   000a_000b_000c_000d
   0000_b000_c000_d000
   0000_0c00_0d00_0000
+) 0000_00d0_0000_0000
-------------------
   000a_bcdb_cd0c_d00d

適当に乗算に帰着して、その後にシフトすればよいです。

msb

\(\sqrt{w}\)-bit 整数の msb(最も左にある 1 の位置)を求めます。 返り値は、1 が存在すれば \(0, 1, \dots, \sqrt{w}-1\) のいずれか、存在しなければ \(\sqrt{w}\) とします。

以下の自明な前処理を行うことで、以降は返り値が \(0, 1, \dots, \sqrt{w}-2\) のケースを仮定します。

  • \(0\) との比較によって、答えが \(\sqrt{w}\) かの判定は \(O(1)\) time でできる。
  • 最上位ビットを調べることで、答えが \(\sqrt{w}-1\) かの判定は \(O(1)\) でできる。

最上位ビットは 0 と仮定してよいので、入力は 0bcd とし、distribute しておきます。

0bcd_0bcd_0bcd_0bcd

ここで、「bcd001 以上のとき(b の左の桁に)繰り上がりが発生する値」「bcd010 以上のとき繰り上がりが発生する値」「bcd100 略」を考えます。それぞれ 0111 0110 0100 です。これとの足し算を行い、キャリービットとのマスクを取ればよいです。

   0bcd_0bcd_0bcd_0bcd
+) 0000_0111_0110_0100
---------------------
&) 0000_1000_1000_1000

あとは、count ones を行えばよいです(010 以上であれば 001 以上であるなどの関係から従います)。

reverse

上位 \( (\sqrt{w}-1)\)-bit と下位 1-bit に分けます。

入力を abcd とします。

shift

0000_0000_0000_0abc

distribute(シフト幅は異なる)

0000_abca_bcab_cabc

mask

0000_a000_b000_c000

distribute

00c0_0000_0000_0000
000b_000c_0000_0000
0000_a000_b000_c000
-------------------
00cb_a00c_b000_c000

shift

0000_0000_0000_0cba

or

0000_0000_0000_dcba

parity

\( (\sqrt{2w}-1)\)-bit のパリティを付与した \(\sqrt{2w}\)-bit 整数を返します。 たとえば、_101_10111101_1011_110_01010110_0101 です。

\(w = 32\) とし、入力を _abc_defg とします。

0000_0000_0000_0000_0000_0000_0abc_defg

distribute and mask

0000_0000_0000_0000_0000_0000_0abc_defg
0000_0000_0000_0000_00ab_cdef_g000_0000
0000_0000_000a_bcde_fg00_0000_0000_0000
0000_abcd_efg0_0000_0000_0000_0000_0000
defg_0000_0000_0000_0000_0000_0000_0000
---------------------------------------
d000_a000_e000_b000_f000_c000_gabc_defg

ここで新しいことをします。やりたいことは、a b ... g の xor を 7-bit 目に入れることです。

各桁は、\(d \cdot 2^{31}\), \(a \cdot 2^{27}\), \(e \cdot 2^{23}\), ..., \(c \cdot 2^{11}\), \(g \cdot 2^7\) です。 これらを \(15\cdot 2^7\) で割ったあまりを考えます。\(2^{31} = 2^4\cdot 2^7\cdot 2^{20} = (15+1)\cdot 2^7\cdot 2^{20}\) なので、\(d\cdot 2^{31} \bmod (15\cdot 2^7) = d\cdot 2^7\) です。他の桁についても同様に \(a\cdot 2^7\), \(e\cdot 2^7\) などになります。

よって、上記の値を \(15\cdot 2^7\) で割ったあまり(の下位 8-bit)が所望の値となります。

\(\sqrt{w}\)-bit の特定の値が含まれるかを探します。その値を distribute したものと xor することで帰着できるので、\(0\) を探します。 C 言語における strlen のような関数をイメージするといい場合があるかもしれません。

検索パターンは \(\sqrt{w}\)-bit ではなくてもよいです。

abcd_efgh_ijkl_mnop

is-zero をします。

000r_000s_000t_000u

1111 を掛けることで複製します。distribute のシフト幅が違う版です。

rrrr_ssss_tttt_uuuu  # x

このとき、たとえば 0000_1111_1111_1111 とか 1111_0000_1111_1111 とかになっています。

このうち、一番右の 0 の位置にのみ 1 があり、他が 0 である整数(0001_0000_0000_00000000_0001_0000_0000)を計算し、それに対して exact log をすればよいです。

上記の計算は、!x & (x + 1) などでできます。

実装例へのリンクがついたツイートがありました。

応用

O(1)-time msb

\(\sqrt{w}\)-bit の msb が \(O(1)\) time でできることを利用し、\(w\)-bit の msb も \(O(1)\) time にできます。

abcd_efgh_ijkl_mnop

is-zero と from bits をしておきます。

0000_0000_0000_rstu

これに対して msb を求めると、元の整数における「一番左の 1 があるブロックの位置 \(i\)」がわかります。 続いて、そのブロックの \(\sqrt{w}\)-bit 整数について msb \(j\) を求めることで、全体の答えが \(i\cdot\sqrt{w}+j\) とわかります。

pext

parallel extract の略です。二つの引数 src mask を受け取り、src のうち mask で立っている箇所のビットを集めてくる演算です。例は以下のような感じです。

      abcd_efgh_ijkl_mnop  # src
pext) 1010_0001_0111_0010  # mask
-------------------------
      0000_0000_0ach_jklo

mask で立っていない src のビットには関心がないので、(適宜 and を取るなどして)そこは 0 であるとします。

      a0c0_000h_0jkl_00o0  # src
pext) 1010_0001_0111_0010  # mask
-------------------------
      0000_0000_0ach_jklo

各ビットについて、右にいくつかずつ動かしており、いくつなのかを調べます。

  • o: 1
  • l: 3
  • k: 3
  • j: 3
  • h: 4
  • c: 8
  • a: 9

これは、そのビットの位置より右にある(mask での)0 の個数であることがわかります。 なので、流れとしては以下のようなことをしたいです。

  • src のビットのうち、移動幅が奇数(2 で割って 1 あまる)のものを 1 動かす
  • src のビットのうち、移動幅が 4 で割って 2 あまるものを 2 動かす
  • src のビットのうち、移動幅が 8 で割って 4 あまるものを 4 動かす
  • ...

移動幅は、「mask において自分より真に右にある 0 の個数」に対応しますが、扱いやすくするため、「自分より右(自分含む)にある 1 の個数」(の差分)の形に変換しておきます。これは、mk = !mask << 1 で得られます。 移動幅が奇数のものは mk における suffix parity を計算することで得られ、それと mask との and を計算することで、ずらすべきビットがわかります。

1010_0001_0111_0010  # mask
1011_1101_0001_1010  # mk = !mask << 1
1001_0100_1111_0110  # mp = suffix_parity(mk)
1000_0000_0111_0010  # mv = mp & mask

これに従って src mask を動かします。

a0c0_000h_0jkl_00o0  # src
1010_0001_0111_0010  # mask
1000_0000_0111_0010  # mv
-------------------
0ac0_0000_00jk_l0o0  # (src ^ (src & mv)) | ((src & mv) >> 1)
0110_0001_0011_1001  # (mask ^ mv) | (mv >> 1)

残り動かすべき幅を確認しておきます。

  • o: 1 → 0
  • l: 3 → 2
  • k: 3 → 2
  • j: 3 → 2
  • h: 4
  • c: 8
  • a: 9 → 8

さて、mk のうち suffix parity が奇数のものは動かしたので、mk & !mp で更新しておきます。 この状態で suffix_parity(mk) を求めることで、移動幅が 4 で割って 2 あまるものを求めることができます。

あとはこの流れを繰り返すだけですが、やや難しいかもなので、考えながら手を動かしてみるのがよいかもしれません。

fn pext(src: u64, mask: u64) -> u64 {
    let mut x = src & mask;
    let mut mk = !mask << 1;
    for i in 0..6 {
        let mp = suffix_parity(mk);
        let mv = mp & mask;
        mask = (mask ^ mv) | (mv >> (1 << i));
        let t = x & mv;
        x = (x ^ t) | (t >> (1 << i));
        mk &= !mp;
    }
    x
}

計算量は \(O(\log(w)^2)\) 時間ですが、mask が共通であれば mv の値は同じになるので、そういう状況では前処理 \(O(\log(w)^2)\) 時間のクエリ \(O(\log(w))\) 時間です。

pdep

parallel deposit で、pext の逆です。下位ビット(mask1 の個数ぶん)を、mask の位置に配置します。

      abcd_efgh_ijkl_mnop # src
pdep) 1010_0001_0111_0010  # mask
-------------------------
      j0k0_000l_0mno_00p0

pext のループを逆順に辿ります。mv で使う値を予め求めて配列に入れて、後ろから操作する感じになると思います。計算量はこれも \(\angled{O(\log(w)^2), O(\log(w))}\) です。

たとえば、「x のうち右から \(i\) 番目の 1 はどこか?」という際に、pdep(src = 1 << i, mask = x) として(必要に応じて exact log も使いつつ)求めることもできます。

permute

予め与えられた順序にビットを並べ換えます。

たとえば、x = abcd_efghp = [2, 4, 1, 5, 3, 6, 0, 7] であれば aceg_dhfb です(h 側が 0 とし、i 番目のビットを p[i] ビット目に移動させる。たとえば p[0] = 2 なので、hf のあった位置に移動する)。

これは結構(説明するのが)難しいです。

まず、p を表すのに \(w\)-bit 整数 \(w\) 個の配列を使うのは無駄が大きいです。各要素は \(w\) 未満なので \(\log(w)\)-bit で足りてしまうためです。 そこで、p[i] の \(j\) ビット目が \(j\) 要素目の \(i\) ビット目にくるような \(\log(w)\) 要素の \(w\)-bit 整数の配列 p' を用います。

[2, 4, 1, 5, 3, 6, 0, 7]  # p
[0, 0, 1, 1, 1, 0, 0, 1]  # p'[0]
[1, 0, 0, 0, 1, 1, 0, 1]  # p'[1]
[0, 1, 0, 1, 0, 1, 0, 1]  # p'[2]

実際には p' の要素たちは \(w\)-bit 整数なので、(配列では 0 要素目が左、整数では 0 ビット目が右なのに注意し)[1001_1100, 1011_0001, 1010_1010] の形であるとします。

この各要素を mask と見て pext など*6を繰り返すことで、基数ソートのように左右に振り分けることで達成可能です。

p' の形式で与えられるとすると、計算量は \(\angled{O(\log(w)^3), O(\log(w)^2)}\) です。

combination

蟻本に載っている、\(k\) 個の 1 がある \(w\)-bit 整数を昇順に列挙するやつです。

rsk0315.github.io

別で書いたのでリンクだけ載せます。

おきもち

「pdep やら popcnt やら CPU 命令で用意されているんだから、自前で実装する必要はないのでは」みたいなのは、まぁそうなんですが、魔法としてそれらを使って生きているのもつまらないかなと思いました。基本的な演算による実装を知ることで面白くなったり、欲しくなった演算が CPU 命令にないときにも応用が効くかなという気持ちがあります。

それから、この手のややガチャガチャした手法に対して「早すぎる最適化は諸悪の根源」みたいなことを言ってくる人がいそうな気もしますが、うるさいな以外の感情がありません。手法のお勉強をしているだけなのでそっとしておいてほしいです。

実際、競プロにおいて(これらの演算自体や、これらを高速化することが)あまり役立つことはないというのはわかった上でやっています。

別の O(1)

\(0\) 以上 \(n\) 未満の各整数に対して count ones をしたいなら、

let mut dp = vec![0; n];
for i in 1..n {
    dp[i] = dp[i >> 1] + (i & 1);
}

などの DP で、一つあたり \(O(1)\) time で得られるほか、ならし計算量でよくある解析(繰り上がりの回数)を用いる方法もありそうです。 DP はいろいろ応用できそうです。

  • exercise: 上記の DP を使って各整数の lsb を求めてみよ。
  • exercise: 同様に msb, reverse, suffix parity を求めてみよ。
    • 書いたことないけどたぶんできるはず

参考文献

おわり

思い立って書き始めた結果、連休が溶けてしまいました。助けてほしいです。

*1:cf. 1 の補数 (ones’-complement)、2 の補数 (two’s-complement)。

*2:より一般に、何らかの偶奇を指す言葉としても使われるかもしれません。

*3:「de Bruijn はオランダ人数学者で,ド・ブランと発音するが,ド・ブルインと発音する人が多い.」らしいです。

*4:一意ではないです。

*5:あるいは、ワードふたつで扱うことにしてもよいかもしれません。

*6:右に寄せる版と左に寄せる版を別途用意するなりします。

リアクティブ問題のデバッグに関して

ABC 269 EABC 244 C など、リアクティブ・インタラクティブ形式の問題は ABC でもたまに出題されます。典型 90 :: 53 にもありますね。

普通の形式の問題であれば、oj などを使ってデバッグしている人でも、リアクティブ形式だと手作業で実行する人は多いのではないでしょうか。

実際には(これくらいの難易度帯の)リアクティブのジャッジを書くのは、(それくらいの難易度帯に挑戦するくらいの)競プロ er にとっては造作もないことかと思います。 とはいえ、デバッグ用の情報をいちいち出力したりするのは面倒な上、全ての問題で共通するパートでしょうから、そういうのは予め準備しておきたいです。

つくった

% sniff_reactive ./e -- ./e-judge test/e/sample-1.in
 +--------------------------------------+-------------------------------------+
 | contestant                           | judge                               |
 +--------------------------------------+-------------------------------------+
 |                                      | 3                                   |
 | ? 1 2 1 3                            |                                     |
 |                                      | 2                                   |
 | ? 1 3 1 3                            |                                     |
 |                                      | 2                                   |
 | ? 1 3 1 2                            |                                     |
 |                                      | 2                                   |
 | ? 1 3 1 3                            |                                     |
 |                                      | 2                                   |
 | ! 3 3                                |                                     |
 |                                      | AC                                  |
 +--------------------------------------+-------------------------------------+

e が提出用のプログラム、e-judge がジャッジのプログラムです。test/e/sample-1.in がテストケースで、e-judge はこれを読んでケースの情報を取得し、それに基づいて e とやり取りをする想定です。

% sniff_reactive solution_program -- judge_program

の形式で実行します。sniff は(通信を)盗聴するとかの意味合いがあるようなので、そういう名前にしてみました。

solution_program は複数コマンドでもよく、judge_program は(上記と違って)自分でケースを生成しながらやり取りするタイプでもよいです。

% # examples
% sniff_reactive python3 e.py -- cargo run --bin e-judge
% sniff_reactive ./e -- ./e-judge 'test/e/sample-1.in'

test/e/sample-1.in の内容は、(自分で用意する)ジャッジが読みやすい形式ならなんでもよいです。 ABC 269 E の例では、たとえば以下のような感じでよいでしょう。

3
0 1 0
1 0 0
0 0 0

ケースに関して

別で e-gen などを用意して、テストケースの生成をするのがよいでしょう。oj が便利です。

% oj g/i ./e-gen

ジャッジが自分でケースをランダム生成するタイプだと、前に試して WA になったケースなどを再現させるのに手間がかかるので、分離させておく方がえびちゃん的には好みです。 ランダム生成をするのに手間がかかると、こうした作業をするまでのハードルが高くなってしまうので、楽にできるような環境を整えておくのがおすすめです。 えびちゃんは、rand_gen! というマクロを用意しています。

rand_gen! {
    rng: ChaCha20Rng;    // 乱数生成器を指定

    n in 10_usize..=20;  // n in [10, 20] を生成
    a in [1_i32..10; n]; // a_i in [1, 10) を n 個生成
}

これは、リアクティブ問題に限らず便利です。

つかう

Gist へのリンクを貼っておきます。

  • sniff_reactive.py
    • Pythonスクリプトだけ置いています。適切な場所に置いて python3 sniff_reactive.py solution -- judge のように使ってもよいでしょう。よしなに。
    • インストールしている Python のバージョンが古いと動かないかもしれません。バージョンを上げるという解決策があります。
  • e-gen.rs
    • ABC 269 E のジェネレータです。Rust 製なので、Rust に詳しくない人は自分の好きな言語で書いてもよいでしょう。
    • えびちゃんのライブラリに依存しているので、Cargo.toml も添えておきます。
      • Rust に詳しくない人は読み飛ばしていいです。
  • e-judge.rs
    • ABC 269 E のジャッジです(AtCoder のサーバで動いているという意味ではないです)。
    • これも Rust なので、好みの言語で作り直してもよいでしょう。

gist.github.com

試していないのでわかりませんが、Windows だと動かないかもしれません。Windows 以外の PC を買うという解決策があります。

意図的なバグは入れていないはずですが、「このツールを使って WA/RE になったから提出しなかったのに、コンテスト後に投げてみたら AC だった」のような苦情は受け付けません。

おわり

ねこちゃんです。

えびちゃんの AtCoder ユーザ解説まとめ

振り返りたくなったのでついでに書きます。

新しくなる順。

  • J - 長い長い文字列/Long Long String
    • ほぼ愚直にシミュレートするのが、解析すると妥当な計算量とわかって面白かったので書いた
    • \(x \bmod y\) が \(x\) の半分になるので高々 \(\log(x)\) 回で \(1\) になる
  • Prefix K-th Max
    • 区間 \(k\)-最小値が並列二分探索でできるのに気づいて、面白い気がしたので書いた
    • ついでに WM の宣伝
  • Divisible Substring
    • 公式解説では後ろから見る方針が書かれていたが、前からで解けたので書いた
    • 内容自体は、変数分離のよくあるやつ
  • Sequence Query
    • WM の宣伝
  • Yamanote Line Game
    • Rust の proconioインタラクティブで使って変にしている人を複数見たので書いた
    • コード例も載せたので、自分もここからコピペして使うようになった
  • Inverse of Permutation
    • おふざけではあるが、何らかの気づきを得る人がいるかもしれないと思って書いた
    • SA-IS の base case が、入力が distinct(逆順列を返す)のときであるやつすき
  • LCS
    • LCS がよくある DP 以外で解けるのが面白かったので書いた
    • 古い論文で \(O(n\log(n)+\delta^2)\) 時間が述べられていたが、SA-IS + LCPA の線形時間構築 + 線形 RMQ で \(O(n+\delta^2)\) 時間になった
  • Cutting Woods
    • decremental neighbor query を使えてうれしかったので書いた
    • \(n = 10^9\) bits くらいなら \(n/w\) words くらい持っても平気なやつすき
  • Together Square
    • 線形篩を使えてうれしかったから書いた
    • ぽかいこちゃんが \(\Theta(\sqrt{N})\) 時間の 解説 を書いてくれた
  • Circumferences
    • union-find に対して、全頂点ペアぶんだけ \(\Theta(n^2)\) 回クエリすれば \(\alpha(n)\) がつかなくなる(解析側の問題)のすきなので書いた
    • BFS でやるのと比べて最悪時のオーダーで損していないということ
  • Intersection
    • ビット演算で思考停止できて実装も軽めだったので書いた
  • Long Bricks(★5)
    • 区間代入クエリをしたいが、遅延セグ木もない tree set もないみたいな状況で実装を比較的軽めにやりたいみたいな話があった
    • でも遅延セグ木を勉強すればいいという話もある
  • Slimes
    • 初心者の頃に名前は知っていて挫折していた Hu–Tucker algorithm を書けたので紹介を書いた
    • 中身を書けていなくてすみません
  • Frog 3
    • 黄色になった後くらいに話には聞いて挫折していた LARSCH algorithm を書けたので紹介を書いた
    • 中身はそのうち書きたいです
  • Moving Piece
    • 公式解説にある bonus に加えて、順列の制約も外せてうれしかったので書いた
    • functional graph にモノイドを載せるところまではすぐ思いついたのに、一日くらいダブリングが頭から抜けていて失格だった
      • 問題固有の考察に引っ張られていたりしたのもよくなかった
  • Split?
    • ビット演算で処理できてうれしいと解説を書きがちになってしまう
    • A 問題にこういう解説を書くとノイズになってしまうかもしれません?
  • Index × A(Continuous ver.)
    • 公式解説と雰囲気が違ったのと、CHT 入門のお気に入り問題の式変形だったので書いた
    • 添字をちゃんと考えましょうみたいな記事 の直後だったのにフォロヮが添字を適当にやっていたので書きたくなってしまった
  • Submask
    • pdep とか pext みたいな CPU 命令が有名でない気がしたので書いてみてしまった
    • 競プロにおいては、こういうのをやるより学ぶべきことがある感はあるが、そう言っていると一生学ばない系な気もする
    • こういうのを知っていると 何らかの文脈 で高速化が図れることもあるかも

何らかのことを思いついてうれしくなったり、積んでいたのを書けてうれしくなったり、バグらせている人を見てかわいそうになったりすると書きたくなってしまうみたいです。 A 問題など低難易度帯の解説は、ある程度成長した人の早解きテク用として有用な場合があるかもしれませんね。

自分用のリンク

atcoder.jp

他の人も見れる用のページが無いのは若干不便な気もするものの、あってもうれしい局面が少ない気もしますね。