yukicoder No.1922 Separate and Attach

問題

https://yukicoder.me/problems/no/1922

解法

 Q = (1, 2, \ldots, N) として一般性を失いません。

 k 回の操作を行ったとして、値  i j 回目の操作で  S に属したなら  b _ {i, j}=0 T に属したなら  b _ {i,j}=1 とします。

操作前の値  i の位置を  x _ i k 回の操作後の値  i の位置を  y _ i とします。

また、各  i に対して、 b _ {i, k} b _ {i, k-1} \cdots b _ {i, 1} を 2 進数として解釈した値を  v _ i とします。

このとき、任意の  i, j\ (i\neq j) に対して次が成り立ちます。

  •  v _ i \neq v_j ならば、 v _i, v_j の大小関係と  y _ i, y _ j の大小関係は一致する
  •  v _ i = v_j ならば、 x _i, x _ j の大小関係と  y _ i, y _ j の大小関係は一致する (i.e. 元の順番を保つ)

従って、 i \lt j \Rightarrow y _ i \lt y _ j を満たす  y を与える  v であって、 v _ n が最小であるようなものを構築すれば良いです。

これは次のような貪欲法により達成されます。

  •  v _ 1 \leftarrow 0 とする。
  •  i = 2, \ldots, n に対して、
    •  x _ {i - 1} \lt x _ i ならば、 v _ i \leftarrow v _ {i - 1} とする。
    •  x _ {i - 1} \gt x _ i ならば、 v _ i \leftarrow v _ {i - 1} + 1 とする。

実装

https://yukicoder.me/submissions/831860

スマホコーディングなのでインデントがめちゃくちゃです。覚えてたら後で直します。

問題の一般化について

2 個の部分列を取り出して連結する操作ではなく  m 個の部分列を取り出して連結する操作である場合も、操作を  m 進数で解釈することで同様に解くことができます。

Kth term of Linearly Recurrent Sequence

問題

judge.yosupo.jp

動機

Bostan-Mori のアルゴリズムで $\displaystyle \lbrack x ^ N \rbrack \dfrac{P}{Q}$ を高速に計算できると知っていても Kth term of Linearly Recurrent Sequence を解くにはややギャップがあると思ったので、書きました。新規性は何もありません。

解法

当然ですが、$\displaystyle P(x) := \sum _ {i = 0} ^ \infty a _ i \cdot x ^ i$ とすれば $\lbrack x ^ N\rbrack P$ が答えです。入力からは $a _ 0, \ldots, a _ {d - 1}$ しか分からないので、$\displaystyle Q(x) := 1 - \sum _ {i = 1} ^ d c _ i \cdot x ^ i$ と定義して式 $(1)$ を用いて計算します。

$$\lbrack x ^ N\rbrack P = \lbrack x ^ N \rbrack\dfrac{P\cdot Q}{Q} \tag{1}$$

ここで、$P\cdot Q$ に対して次の補題が成り立ちます。

補題
$n \geq d$ なる任意の整数 $n$ に対して、$\lbrack x ^ n\rbrack P\cdot Q = 0$ が成り立つ。

補題の証明
$n \geq d$ なる整数 $n$ を任意に固定する。$\lbrack x ^ n\rbrack P \cdot Q$ を定義通りに計算することで次を得る。 $$\lbrack x ^ n\rbrack P \cdot Q = a _ n - \sum _ {i = 1} ^ d c _ i a _ {n - i}.$$ $a$ の定義より $i \geq d$ なる整数 $i$ に対して $\displaystyle a _ i = \sum _ {j = 1} ^ d c _ j a _ {n - j}$ が成り立つので、これを $i = n\geq d$ として適用することで $\lbrack x ^ n\rbrack P \cdot Q = 0$ が従う。$\blacksquare$

補題より $P\cdot Q = (P\cdot Q) \bmod x ^ d = ((P\bmod x ^ d)\cdot Q)\bmod x ^ d$ として計算できます。$P\bmod x ^ d$ および $Q$ は入力から直ちに求まり、$P\cdot Q$ を $\Theta(d \log d)$ 時間で得ることができます。あとは Bostan-Mori のアルゴリズムを適用することで、全体 $\Theta(d\log d \log N)$ 時間で問題を解くことができます。

ICPC 2022 Yokohama Regional 参加記 (Bu-Bu-Du-Ke / suisen)

昨年組んでいただいたゆきのんさんが引退されたので、kenchoくんに入ってもらってsuisen, otera, kenchoの 3 人でチーム「Bu-Bu-Du-Ke」として ICPC 2022 Yokohama Regional に参加しました。昨年はアジア予選に進出することは叶わず、今年に関しても国内予選や模擬国内/地区予選ではずっと2桁順位だったのですが、アジア予選で 7 位という一番良い結果を残せました。特にsuisen, oteraは今年でICPC引退で、最後のチャンスで力を十分に発揮できてよかったです。株式会社オプトさんから7位の企業賞まで頂けて嬉しかったです。

終結

チームの役割分担は

  • suisen: 重実装の人柱・典型・ライブラリ
  • otera: 実装・天啓。必要十分条件エスパー
  • kencho: Java使いなので考察専門。特にパズルと幾何とad-hoc

みたいな感じです。得意分野がばらけていてバランスはいいと思ってます

本番前々日

ライブラリが完成したので印刷をした。Y.Y.さん対策のDM分解が入っていたり、絶対に使わないであろう特性多項式のライブラリが入っていたりする。ライブラリ無制限の文言を見てライブラリを全部持ってくるくらいしようかと思ってたけど、写経用に短く書き直すのと書き直しで壊れてないかのテストを書くのが大変すぎて泣く泣く断念...

ライブラリ一覧

3時間だけ本番に近い状態でチーム練をしようという話になっていた。VS Codium に拡張機能が何も入っていないという噂を聞いて急遽 CLion をインストールして練習したが、よくわからないお節介機能 (慣れると便利なんだろうけど) でめちゃくちゃになっていた。

18時開始予定だったけどなんだかんだ開始が遅れたり感想戦をしたりしていたら終電になっていて危なかった。

生活習慣の調整が下手だったので眠気が来なくて、27時くらいにやっと寝れた気がする。

本番前日

チームメイトと中華街で合流して、昼食を摂って会場に向かった。朝が得意ではないので12時くらいに到着するように新幹線を取っていたけど、2人は早起きして横浜を満喫しててちょっと後悔。

会場に早く着いたので、チーム紹介ドキュメントを眺めていた。うちのドキュメントは実行可能なBrainf*ckのプログラムになっていてコピペして実行してもらおうという魂胆だったけど、よく考えると印刷物はコピペできないので...

↓にコピペ用のドキュメントを貼っておくので、良かったら https://copy.sh/brainfuck/ などで実行してみてください。

>>>>>>>>>++++++[-<+++++<+++++++++++++<+++++<+++++++++++++<+++++<++++++<++++
+++++++++<+++++<+++++wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww++++++++>>>>>>>>>]<+
+<-<++<-<++<+<<++<->-* Bu-Bu-Du-Ke (Kyoto University) *------------------.>
>>.........<+++++++++* > kencho                       *<...>...............
.....................* > otera                        *<.....<...<<.>>>..<<
+++++++++++++++++++++* > suisen                       *<+.<<..........>>+++
++++++.>--.>.........* > suibaka (Coach)              *<...................
...<++.>>+.-.........wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww+.<<--.<<.>>>..>-...
.........+.<...............................<----.<<<............>>+++++++++
++++++.<.>>>.<++++++.<<<...-------.>>>>....>++++++++++++++++++++++.>>...<++
+.<<..<...>....<....>....wwwwwwwwwwwwwwwwwwwwwwwwwww..............>.>>..+.>
....<<<<<.>>>>-....<+wwwwMw!+=.=<<~..~::~~<..>.. ??MMwww...<....<<.>>>.<.>>
>>...<.<<..<..>>>>.MM**~ ...:~``7H=<<====<=<``~~:::<==<?MM<.>>>>..>>>---.++
+..<--.<<...>>>--Mv*-.+++.~  .;d\<.XXXXAzzwHm   <  ~:*;*==?M<..............
....>>>...<.<<..M*.<=:~```  .WWmgkJZXzkwXXXI_*dY^     `~**>>M.>>...<.<<..<.
>>>>...<++.<<..M*<=*;*   .=>!>>(=sjXkZXkzXUr?TTXQaJ(>   .?*=.M.>+.<<<<<<<.>
>>.<--.>>>>...+M*(;*~.   `    dHH-.77TUWWWX.(~(Hmmg ~    .*;=M...<--.<-----
-.<...>>>...<<<Myz<<~<~    < ?H _7YH= .gHY"Wa>_J9\  `   _>~(zM>>...>>>...<<
<<<<<.>>>>.....Mblli=.~.. .!`  TH=.d"%  .dm==(Y` `     ..(JllM.....>>>...+.
.-......<-.<<...Mslllz(=<+.>_    ?=   ~*`        _ _>>(JllllM>...<<++.<..<+
.>>>>...<<<<<<.>M5zzOz>>.1>>>---.___   _   _ __+++((zzztlltlM..........+.<<
<<--.>..<++.>>>>-.Mzl..zIzz<?<?--1-=(((((((=-zzzzOOOllllOzOM.<..<.>>>>...<<
+.<...<--.>>>>.....MzOzl.?.=zzIzzzzz.==.==<+?=++=lllzOZOllM.<<...<++.>>>>..
.........<<+++++++++MelvOwzz?++==???=?++=+=++.zzzZZOltlllM<<.>..>---------.
>>..<<<<<<<.>>>>...>++Mz+=.zOVrwyzzz>>zzzzvOz..llllltlllM<--.<<....>>>.....
....<<<<<<<------.>.>>>.Mx>?--?=-.>>z=z1z=...==lzzzttOM<<----.<...<--.>>>>.
...<<<<<<<-----.>>>>.<++.>Mae>zz>>OOOOtOOOOOOOOzllOgMM..+.<<.<..<.>>>>-...<
<----.<.................<.>>>Mzwzzzzz>...-zzzzudM------.>-....<+++++++....<
<++++++.<..>>>...<<-----.<...>MwwwwwwwwwwwwwwwwM>>...<<<<<<<+++++.>.>>>.>>>
....>...<<<<<.>>+++++++.>>....<<--------.<.>>>---.+++..<.<<...>>>---.+++...
<<<<<-.>>.................<-.>>>>...<<+.<...<+.>>>>....<.<<.<.>>>>...<+++++
+.<<..<.>>>>...-------.<<<<<<.>>>.>>>+++++++....+..-......<<<<<<<.>>>>..>>>
---.+++...+.<<++++++++++.>>-.....<<<<<+++.>>.................>>>---.+++...+
..-.......<<---------.<.<.>>>>...+.<<+++++++++.>>-.....-------.<<<<<<.>>.>>
>>+++++++..........>+++.-.--...<<<<<.>>>>.....<<<<<<<++.>>>>>>>....<<------
-.<.................>++++.>>..........<<<<+++++++++++.>>----.<...>>>+++++++
.-------.....<<<<-----.>>>>...-------.<<<<<<..>>>.<------.>>>-.<+++++++.>>+
++++++......+.<<++++++.<<.>................................>---------.>>-..
.>+.-....<<<<<.>>>>.+.<++++.<<<<<.>>>.<.>>>>-...........<<---------.<......
.........................>>>...+.<<+++++.<..<.>>++++.>>-....<+++++++.<<<<<.
>>>.<--.>>>>...<+++++++.<<...<.>>>>....<<<<++.>............................
..>>>...>+++.---.<<<<<.---.>>>>....<<<<<<<++.>>>>.........<+++..>.<<<.>>>.<
------.>>>>...<<----.<....>>>---.+++...<<---.<.<++++++.>>>>...<<-.<..<.>>>>
...<<<<<.>>................<.>>>>...+.>++++++++.<-....<<<+++.++++.-------..
....<.---.>>>>.....+.<<<<+.<<.>>>.>>>----.+++..+.<<<.....<+.>>>>-...<<.<.<.
>>>>...<<++++.<..<.>>>>..+.<<<.................<+.>>>>-........<<<<<-------
---.>>.......<.>>>>..<<<+++.>--.<---.<<++.>>>>>..+.<<<<.<<.>>>.>++++++.>>-.
.<<<+++.---.....>.>>...<<<+.-.>.>>..-------.<<<...>.>>+++++++..<<<+++.---..
..<--.>>>>.........>---.<<<<..<++.>>>>........<<<<<------.>>......<.>>>>...
<<<+.-..<.>>>>...<<<<<+++++++++++++.<.>>>.>>>...<++++++.<<....<.>>>>...<+++
+++.<<..>>>...<<-------.<...>>>...<.<<.................<.>>>>...<<.<<--.>>>
>....+.<<<<++.>....>>>-...........<<+++.<<<<.>>>.>>>...<.<<..<.>>>++++.>...
<<<<<<<----.>>>>..<.>>>>...<<<<<--.>>..<.>>>>...<----.<<.................<-
.>>>>...<<---.<..>---------.>>...+.<<<<-.>...>>>-...<.<<<<<<-.>>>>+.<<+++.>
>>++++++++++++...<<<<.>>>-.>>>...+..-.....<.<<...<++.>>>>...<<+++++++++++++
.<<.>>>++++.>....<----.<<.................<------.>>>>...++++++++++++++++.<
<<...<.>>>+++++++....>+++++.<<<..>>...>.<<<<++++++.>..<.>>>.+.<<<--.<<.>>++
.>>>-........<<<<<<+.>>>>+.-.....>>-----.+++++.........>-----.<<<..........
.......<------.>>>...>.<<<....<++++++.>>>..<<<<<<.>>>-------.>..<+++++.>>>.
......-----.>+.<<<<<<..

リハーサルでは、ジャッジシステムは模擬地区予選と変わらなさそうだったので最低限の確認だけやったあとはタイピング練習をした。

ノートPCしかまともに使ったことが無かったのでストロークの深いキーボードを打つのは本当に辛くて、数十行のテンプレートを写経するのに10-15分くらい掛かっていて不安が残った。4年半JIS配列でずっとやってきたのを突然US配列に変更したというのもあって、本当にボロボロだった。

US配列の練習としては、本番まで時間が無かったので練習初日にABCのA,B,C問題でまだ解いていないものを全部実装して無理やり矯正するなどしていた。

US配列の練習

以降もJIS配列を一切封印して練習して今までの7-8割くらいの速度が出るようになって安心していたんだけど、いざリハーサルでキーボードを触ってみると物理的に形と大きさが違うキー (Enter とか \ とか) とかストロークの深さとかの諸々に全く対応できなかった。

コンテスト

コンテスト開始 ~ AB 2 完 (49:00)

コンテスト開始時の役割は予め話し合っていて、

  • suisen: 環境のセットアップ・テンプレート写経
  • otera: 前から順番に読んで、テンプレ写経が終わり次第A,Bを書けそうなら書く
  • kencho: 後ろから順番に読んで考察

という分担で開始。CLionのプロジェクトをC++14で作ってしまったり案の定テンプレート写経がとんでもなく遅れてしまったりで結局oteraにPCを渡したのが多分15分時点ぐらいだった。その後Aが通って、Bの考察も終わっていて書けそうとのことだったのでそのままBを実装してもらった。1ペナしつつも49分で2完。

この間suisenはCDEを読んで

  • Cはフローの見た目をしてるけど何もわからない
  • Dは「右端または左端を揃える」&「上端または下端を揃える」で探索すれば筋力で解けそう?みたいなことを言っていた。動かす駒が右端から左端 (または上端から下端) に動いている場合に壊れることを指摘してもらって、よく分からなくなったので放置
  • Eはシンプルに何もわからない

という感じだった。kenchoは後ろから問題を読んで概要をまとめつつ考察も軽くやっている感じだったと思う。

~ ABG 3 完 (80:00)

その後suisenがGを読んだ (Fはsuisenがテンプレート写経でもたついている間にoteraが読んでいた)。簡単な木クエリの問題ですんなり解けた (つもりになった) ので、PC をもらって実装を始めた。辺uvを張った時にできるサイクルに含まれる頂点であって入り口から一番近い頂点を求めるパートを lca の多数決を取るやつでできると勘違いしていてサンプルが合わなかった。よく考えると入口を根とした木の lca(u, v) が欲しい頂点と分かって、修正してGをAC (80分時点)。lcaと距離計算のためにHLDを写経していて、これが無駄になるとそこそこの痛手になるところだったのでヒヤッとした。書き始めた時点では0ACだったのでFAを取りたかったが、迷走している間に3チームくらい通っていて結局FAから9分遅れだった。

~ ABDFG 5 完 (167:00)

その後はACが出ているDEFを考えていた。oteraにDの概要を話すと辞書順で高々2個ずつ取ってそれぞれ対応させる解法を生やしてくれた。実装が大変めなのでsuisenが書くことに。元の配置を回転する方針を選んでしまったせいで最後の移動の復元で頭を壊してしまい、炎上。4通りの回転それぞれに対して丁寧に場合分けを手打ちしていたのがどう考えても悪手だった。後にhenoさんから目標を回転すると簡単ということを聞いて、大反省...

自分がDをバグらせている間、oteraはEを、kenchoはFの考察をしていた。Dのデバッグを終えて提出しようとすると、submitボタンが反応しなくなるトラブルが起きた。スタッフの方を呼んで対処してもらい、submitに成功してDが通った (155分)。スタッフの方が見守る中の提出で、通ったときに拍手をもらえたのはかなりレアな経験そう?このトラブルの影響でBu-Bu-Du-Keのコンテスト時間が1分伸びた。

Dが通った後にFの進捗を聞くと微妙な反応で、考察を共有してもらってみると何と解けていたので実装してAC (167分時点)。本人曰く正当性に自信が無かったらしい?3完から一気に5完になって順位が跳ねた。

~ ABDEFG 6 完 (237:00)

Fが通ったあとは順位表の雰囲気的にEは通すべき問題っぽかったのでsuisenはEを考えていた。(Iの個数)=(Cの個数)<(Pの個数)の場合を考えてみると、

  • (Iの個数)=(Cの個数) の条件は、Iを+1に、Cを-1に、Pを0に置き換えた数列の累積和  S について  S _ l = S _ r が成り立つことと言い換えられる
  • (Cの個数)<(Pの個数) の条件は、Iを0に、Cを+1に、Pを-1に置き換えた数列の累積和  T について  T _ l \gt T _ r が成り立つことと言い換えられる

で、区間和クエリに帰着されて Fenwick Tree で log 1 つで数えられそうとなった。(Cの個数)=(Pの個数)<(Iの個数)の場合と(Pの個数)=(Iの個数)<(Cの個数)の場合も同様で、結局置き換えのパターン 3 種類に対して累積和を作って、累積和の値ごとに座圧した Fenwick Tree を持てば元の問題も log 1 つで数えられることが分かった。

大量の座圧をガチャガチャやっていたら当然バグらせて、悩んでいる間にJの考察が終わったらしくoteraに実装を渡してEの紙デバッグを始めた。結局バグの原因は座圧ガチャガチャパートではなくFenwick Treeにあった。Fenwick Treeくらいなら写経するよりそらで書いた方が速いので練習でも毎回そらで書いていて、今までそれでバグらせたことは無かったので気づくのが遅れてしまった。

修正したコードをまだバグってるだろうな~と思いつつ提出すると一発で合ってくれて、6完 (237分時点)。

~ コンテスト終了

Jの実装をoteraに任せてsuisen, kenchoでKの考察をした。想定解に近いところまでは考察できていたが、計算量を落としきれなかった。その後Jの実装をsuisenが引き継ぐも、線分の交差判定のやり方がまずかったのかサンプル2が合わないままコンテスト終了。

おわりに

オンサイト開催になって本当に良かった。準備していただいた方や交流していただいた方、チームを組んでくれたoteraくんとkenchoくん、ありがとうございました。

AOJ 1430 Even Division

$$ \gdef\V{\mathcal{V}} \gdef\set#1{\lbrace #1 \rbrace} \gdef\card#1{\vert #1 \vert} \gdef\sqbrack#1{\lbrack #1 \rbrack} \gdef\induced#1#2{#1\sqbrack{#2}} $$

これは 帰ってきた AOJ-ICPC Advent Calendar 2022 - Adventar 15日目の記事です。

2021年の横浜Regionalで使われた問題のネタバレを大いに含みます。練習に使う予定がある方は注意してください。

問題

https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1430

以下要約:

$N$ 頂点 $M $ 辺の連結な単純無向グラフ $G$ が与えられる。$V = \set{0, \ldots, N - 1}$ の部分集合族 $\V = \set{V _ 1, \ldots, V _ k} \subseteq 2 ^ V$ が以下の 4 つの条件を全て満たすとき、またそのときに限って「$\V$ は $G$ の even division である」という。

  1. $\V$ は $V$ の 分割 である。
  2. 任意の $V _ i \in \V$ について $\card{V _ i}$ は 正の偶数 である。
  3. 任意の $V _ i \in \V$ について、$V _ i$ による 誘導部分グラフ $\induced{G}{V _ i}$ は連結 である。
  4. $\V$ の 細分 $\V'$ であって $\V' \neq \V$ かつ上記 1,2,3 を満たすものは存在しない。

$G$ の even division を 1 つ求めよ。

制約

  • $2 \leq N \leq 10 ^ 5$
  • $N$ は偶数
  • $N - 1 \leq M \leq 10 ^ 5$
  • 与えられるグラフは連結であり、多重辺や自己ループを含まない。

解法

方針: $G$ の DFS 木 $T$ を取り、ボトムアップに成分を取り除くことを繰り返して even division を構築する。

まず、$T$ の部分木であって頂点数が偶数のものがあればその部分木は切り離してもよい (i.e. 条件1,2,3 を満たすことは保たれる) ので、切り離す。操作のイメージは図 1 を参照。以降の図においては断りなく木を構成する辺を黒で、後退辺をで描く。

図1. 偶数サイズの部分木の切り離し

この操作を可能な限り繰り返すことで $T$ はいくつかの DFS 木 $T _ 1, \ldots, T _ l$ に分割される。各 DFS 木 $T _ i$ は次の性質 $(\star)$ を満たす。

性質 $(\star)$
頂点数が偶数の連結グラフであり、自身を除いた全ての部分木のサイズは奇数である。

$T _ i$ 毎に独立に分割が求まればよいので、以降は $T$ が $(\star)$ を満たすと仮定する。

$T$ の後退辺 $\set{x, y}$ に注目する。後退辺は必ず先祖/子孫関係にある2頂点を結ぶので、$x$ は $y$ の先祖であるとして一般性を失わない。図2のような、$xy$ パスに沿った分割の操作を考える。

図2. 後退辺の除去

図2において青で囲った成分の1つ $C$ に注目する。$C$ において部分木のサイズが変化するのは根のみであり、奇数サイズの子部分木を取り除いているので $\card{C}$ は偶数である。$C$ 自身以外の部分木のサイズは奇数であったので、$C$ は $(\star)$ の性質を満たしている。

続いて、橙で囲った成分 $C'$ に注目する。後退辺 $\set{x, y}$ を木を構成する辺として用いることで、$C'$ は再び DFS 木となり、$(\star)$ の性質を満たす。

従って、この分割に関して次が成り立つ:

  • 後退辺が 1 つ減る。
  • 得られた全ての連結成分は $(\star)$ を満たす。

これを可能な限り繰り返すことで得られた成分たちには後退辺が存在しない。即ち、最終的には森が得られる。

各成分が $(\star)$ を満たすように分割を繰り返したので、森の各連結成分もまた $(\star)$ を満たす。ここで、以下の事実が重要である。$V(T)$ は $T$ の頂点全体の集合を表す。

重要な事実
木 $T$ が $(\star)$ を満たすならば、$\set{V(T)}$ は even division である。

証明

$V(T)$ の任意の 2 分割 $\set{A,V(T) \setminus A}$ が条件 1,2,3 のいずれかに違反することを示せば十分。$\set{A,V(T) \setminus A}$ が条件 1, 2, 3 を満たすと仮定する。

条件 2 より $A\neq \emptyset$ かつ $V(T) \setminus A\neq \emptyset$ が必要である。従って、条件 3 を満たすような分割は $T$ の辺を 1 本削除することに対応する。しかし、辺を 1 本削除することで得られる $2$ つの成分のサイズは $(\star)$ の性質より必ず奇数となり、条件 2 に違反する。$\blacksquare$

以上で even division が得られた。この手続きは全体 $O(N + M)$ 時間で可能であり、十分高速である。

実装

実装は結構大変で、私はかなりバグらせました。

以下の実装は2重DFSで書かれていて、DFS木を構築する外側のDFSにおいて偶数サイズの部分木を発見したら後退辺を除去する内側のDFSを呼んでいます。

外側はともかく内側のDFSを正確に処理するのは結構大変で、全域森を管理を上手くやる必要があります。ここでは親へのリンク $\mathrm{par}$ を用意して、常に $\set{\set{v, \mathrm{par}\sqbrack{v}} \mid \mathrm{par}\sqbrack{v} \neq -1}$ が全域森の辺集合と一致するように $\mathrm{par}$ を更新しています。$\mathrm{par}$ を使えばある辺が全域森に属するかどうかを簡単に判定できます (実装例の is_tree_edge がそれ)。

実装する上で特に気を付けるべき点は以下の3つくらいでしょうか。

  • 先祖方向の後退辺は適切に skip する。実装例では深さの配列 dep を用意してその大小で先祖方向 or 子孫方向の判定をしています。
  • 内側のDFSで後退辺が指す先が既に切り離された頂点である (i.e. 異なる成分に属している) 可能性があることに注意する。
  • $\mathrm{par}$ の更新を忘れない。

提出リンク

#include <iostream>
#include <utility>
#include <vector>

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int n, m;
    std::cin >> n >> m;

    std::vector<std::vector<int>> g(n);
    for (int i = 0; i < m; ++i) {
        int u, v;
        std::cin >> u >> v;
        g[u].push_back(v);
        g[v].push_back(u);
    }

    std::vector<int> dep(n);      // 深さ。後退辺の方向 (先祖方向 or 子孫方向) の判定に用いる
    std::vector<int> par(n, -1);  // 全域森の管理。{v, par[v]} が全域森の辺

    auto is_tree_edge = [&par](int u, int v) { return par[u] == v or par[v] == u; };

    std::vector<int8_t> visit(n);    // 既に訪問した頂点かどうか
    std::vector<int8_t> removed(n);  // 既に切り離された頂点かどうか

    std::vector<std::vector<int>> ans;
    auto collect = [&](int r) {  // 全域森で r が属する成分の頂点を集めて答えに追加
        auto &cmp = ans.emplace_back();
        auto collect_rec = [&](auto collect_rec, int cur, int prv) -> void {
            cmp.push_back(cur);
            removed[cur] = true;
            for (int nxt : g[cur]) if (is_tree_edge(cur, nxt) and nxt != prv) {
                collect_rec(collect_rec, nxt, cur);
            }
        };
        collect_rec(collect_rec, r, -1);
    };

    auto remove_even_subtrees = [&](auto remove_even_subtrees, int root) -> int {
        visit[root] = true;
        int sub_size = 1;
        for (int v : g[root]) if (not visit[v]) {
            par[v] = root;
            dep[v] = dep[root] + 1;
            sub_size += remove_even_subtrees(remove_even_subtrees, v);
        }

        if (sub_size % 2 == 0) {
            par[root] = -1; // 全域森から辺を削除 (root を根とする部分木を切り離す)
            auto remove_backward_edges = [&](auto remove_backward_edges, int u, int par_u) -> void {
                for (int v : g[u]) if (v != par_u) {
                    if (is_tree_edge(u, v)) {
                        remove_backward_edges(remove_backward_edges, v, u);
                    } else if (not removed[v] and dep[v] > dep[u]) { // 子孫方向の削除されていない後退辺
                        // 全域森から辺 {v, par[v]} を削除して辺 {v, u} を追加 / uとvの間にある頂点 r を列挙
                        for (int r = std::exchange(par[v], u); r != u;) {
                            // 辺 {r, par[r]} を削除
                            int nxt_r = std::exchange(par[r], -1);
                            collect(r);
                            r = nxt_r;
                        }
                    }
                }
            };
            remove_backward_edges(remove_backward_edges, root, -1);
            collect(root);
        }
        return sub_size;
    };
    remove_even_subtrees(remove_even_subtrees, 0);

    std::cout << ans.size() << '\n';
    for (const auto& cmp : ans) {
        std::cout << cmp.size();
        for (int v : cmp) std::cout << ' ' << v;
        std::cout << '\n';
    }
}

Static Range Inversions Query

前置き

judge.yosupo.jp

を解いたのでメモ。前計算 $O(N\sqrt{N})$、クエリ $O(\sqrt N)$ で解きます。

問題

長さ $N$ の整数列 $A=(A _ 0, \ldots, A _ {N-1})$ が与えられる。以下のクエリを $Q$ 個処理せよ:

  • l r : $A$ の連続部分列 $A _ l A _ {l+1} \cdots A _ {r-1}$ の転倒数を計算する。

記号の定義

  • 列 $A$ の転倒数を $I(A)$ と書く
  • 列 $A,B$ に対して、$A _ i \gt B _ j$ なる $(i, j)$ の個数を $I(A,B)$ と書く
  • 列 $A,B$ をこの順番に連結してできる列を $AB$ と書く
  • 列 $A$ の連続部分列 $A _ l A _ {l + 1} \cdots A _ {r - 1}$ を $A _ {\lbrack l : r)}$ と書く

解法

区間クエリなので Segment Tree などの利用が思い浮かびますが、モノイドではあるものの演算を高速に計算することが難しいので、平方分割による解法を考えます。

列 $A$ を長さ $L = \Theta(\sqrt{N})$ のブロックで分割します。簡単のため $L \mid N$ を仮定して整数 $K = N / L= \Theta(\sqrt{N})$ を定めます。実装においても $\infty$ を末尾に $L$ 個未満だけ加えることで $L \mid N$ とできます。

さて、クエリ区間とブロックの関係は次の $2$ パターンが考えられます。

タイプ1: クエリ区間が $2$ つ以上のブロックに跨るケース

タイプ2: クエリ区間が $1$ つのブロックに含まれるケース

$I(XYZ)$ に関して、次が成り立ちます。

$$I(XYZ) = I(XY) + I(YZ) - I(Y) + I(X, Z). \tag{1}$$

タイプ 1 は式 $(1)$ を用いて計算し、タイプ 2 に関しては式 $(1)$ を変形した次の式 $(2)$ を用いて計算します。

$$I(Y) = I(XY) + I(YZ) - I(XYZ) + I(X,Z). \tag{2}$$

式 $(1),(2)$ の右辺に現れる $I(X,Z)$ 以外の項は $A$ のある連続部分列の転倒数となっています。ここで、これらの連続部分列の左端または右端の少なくとも一方がブロックの境界と一致していることが重要です。

なぜなら、左端がブロックの境界と一致するような連続部分列の個数は $\displaystyle \sum _ {i = 0} ^ {K - 1} (K - i) L = \Theta(K ^ 2 L) = \Theta(N\sqrt{N})$ 個しかないためです。同様にして、右端がブロックの境界と一致するような連続部分列の個数も $\Theta(N \sqrt{N})$ 個です。

従って、左端または右端がブロックの境界と一致するような合計 $\Theta(N \sqrt{N})$ 個の連続部分列に対する転倒数を全て前計算してしまえば、式 $(1), (2)$ の右辺のうち $I(X,Z)$ 以外の項は全て簡単に得ることができます。

ぱっと思いつく前計算の方針として、Binary Index Tree で各値の個数を管理しながら列を伸ばしていくという方法があります。しかし、この方法では前計算に $\Theta(N \sqrt{N} \log N)$ 時間掛かってしまうので、工夫する必要があります。

以下、この前計算を $O(N\sqrt{N})$ 時間で行う方法を説明します。

左端がブロックの境界と一致するような連続部分列の転倒数を前計算

まず、長さ $L$ 以下のものに関しては愚直に計算してしまっても $\Theta(L ^ 2 K) = \Theta(N\sqrt{N})$ 時間で済みます。当然、Binary Index Tree などを用いた高速化を行っても構いません。その場合は $\Theta(KL \log L) = \Theta(N \log N)$ 時間です。

長さが $L$ より大きく、左端がブロックの境界と一致するような連続部分列 $A _ {\lbrack i\cdot L : (i + 1)\cdot L + k)}\; (k\gt 0)$ の転倒数の計算を考えます。

$A _ {\lbrack i\cdot L : (i + 1)\cdot L + k)}$ を次のように $3$ つに分解します。

$$A _ {\lbrack i\cdot L : (i + 1)\cdot L + k)} = A _ {\lbrack i\cdot L : (i + 1)\cdot L)} A _ {\lbrack (i + 1)\cdot L : (i + 1)\cdot L + k - 1)} A _ {(i + 1)\cdot L + k - 1}.$$

$(1)$ は一般に成り立つのでこれを適用すると次のようになります。ただし、表記の簡単のため $\mathrm{invL}(i, l) := I(A _ {\lbrack i \cdot L: i \cdot L + l)})$ と定めました。

$$\begin{aligned} \mathrm{invL}(i, L + k) &= \mathrm{invL}(i, L + k - 1) + \mathrm{invL}(i+1, k) - \mathrm{invL}(i + 1, k - 1) \\ &+I(A _ {\lbrack i\cdot L : (i + 1)\cdot L)}, A _ {(i + 1)\cdot L + k - 1}). \end{aligned}$$

ここで、$c _ {i, k} := I(A _ {\lbrack i\cdot L : (i + 1)\cdot L)}, A _ {(i + 1)\cdot L + k - 1})$ は $A _ {i\cdot L + j} \gt A _ {(i + 1)\cdot L + k - 1}$ なる $0 \leq j \lt L$ の個数と一致します。従って、$A _ {\lbrack i\cdot L : (i + 1)\cdot L)}$ の要素を予めソートしておき、$A _ {(i + 1)\cdot L + k - 1}$ の値が大きい順に $k$ を走査して尺取り法を行うことで、全ての $k$ に対する $c _ {i, k}$ を計算できます。

この方法でネックとなるのは「$A _ {(i + 1)\cdot L + k - 1}$ の値が大きい順に $k$ を走査して」という部分で、ここで毎回ソートしてしまうと全体で $\Theta(N \sqrt{N} \log N)$ 時間掛かってしまいます。そこで、$i$ の降順に計算することにすれば、$i+1$ のときの $A _ {\lbrack (i + 1) \cdot L : KL)}$ のソート結果と $A _ {\lbrack i \cdot L : (i + 1)\cdot L)}$ のソート結果をマージソートの要領で併せることで、$A _ {\lbrack i \cdot L : KL)}$ のソート結果を $O(N)$ 時間で得られます。

以上より、すべての $c _ {i, k}$ の計算を全体 $O(N \sqrt{N})$ 時間でできたので、$\mathrm{invL}$ の計算も全体 $O(N \sqrt{N})$ 時間でできました。

右端がブロックの境界と一致するようなものに関しても同様にできるので、省略します。


さて、最後に残ったのは式 $(1), (2)$ における $I(X,Z)$ の計算です。ここで重要なのは、いずれの場合においても $X$ および $Z$ は $1$ つのブロックに包含されているということです。

即ち、各ブロックをソートした列を前計算しておき、$X$ を包含するブロックと $Z$ を包含するブロックの間で尺取り法を行うことで $I(X,Z)$ を $O(L) = O(\sqrt{N})$ 時間で計算できます。

なお、ソート列の前計算に関して、尺取り法において要素が $X$ や $Z$ に含まれるかの判定を行う必要があるので、添字の情報を付加した組 $(A _ i, i)$ をソートした列を覚えておく必要があることに注意します。

実装

$L \mid N$ となるように番兵を置くことで、左端/右端がブロック境界と一致するような連続部分列の転倒数の計算を共通化できて実装が楽になります。(右端の時はreverseして比較関数を逆転すればOK)

judge.yosupo.jp

第11回 アルゴリズム実技検定

114:52 で全完しました。

A

うさぎは $\dfrac{X}{A} + C$ 秒、かめは $\dfrac{X}{B}$ 秒掛かるので分母を払って整数で比較。

B

連想配列に入れてカウント。

C

$10 ^ 9$ より大きい値は $10 ^ 9 + 1$ に潰してもよいので、$\min(N ^ k, 10 ^ 9 + 1)$ を前から求めていくと overflow しない。

D

Union Find で同一視されるグループを求めておく。すべての $i$ に対して $S_i$ と $T_i$ が同じグループに属しているなら Yes で、それ以外は No

E

$(v,\ldots,2,1,2,\ldots,v)$ は $2v - 1$ 項なので、操作 $v$ が終わった時点での数列 $A$ の長さは $\sum _ {i = 1} ^ v (2 i-1) = v ^ 2$ と分かる。

従って、$N$ 項目は操作 $\lfloor \sqrt{N-1} \rfloor + 1$ で追加されることが分かるので、あとは添え字を合わせれば OK。

F

愚直にシミュレートして $O (H (W+N) )$ で間に合う。

G

Union Find で連結成分数を求めて、これがちょうど $1$ なら Yes で、それ以外は No

H

$\mathrm{dp}[i][w_a][w_b]$ を 「$1,\ldots, i$ 番目までのお菓子を見て、$1$ つめのナップサックに入れたお菓子の重さの合計が $w_a$ 以下、$2$ つめのナップサックに入れたお菓子の重さの合計が $w_b$ 以下となるように詰める場合の価値の最大値」として定めれば状態数 $O(NAB)$、遷移 $O(1)$ なので十分高速。

I

すぬけ君の位置が $\mathrm{pos} _ s$、荷物の位置が $\mathrm{pos} _ a$ の状態 $(\mathrm{pos} _ s, \mathrm{pos} _ a)$ を頂点にしたグラフで BFS をすればよい。

J

Pythondatetime モジュールを使うとよい。

K

$\# \{S \text{の部分列} \} + \# \{T \text{の部分列} \} - \# \{S \text{と} T \text{の共通部分列} \}$ として計算する。

$\# \{S \text{の部分列} \}$ は https://judge.yosupo.jp/problem/number_of_subsequences をやればよく、共通部分列の場合はこれを拡張して状態を添え字のペアにすれば OK。

L

最長共通部分列の DP に文字を変更した回数を追加で状態として持てばよい。

M

$2$ 人のスコアが両方とも $i\; (0 \lt i \lt N)$ となった次のゲームで $\dfrac{1}{2}$ の確率で逆転が起こる。従って、求める答えは次で表される。

$$ \sum _ {i = 1} ^ {N - 1} \binom{2i}{i} \cdot \dfrac{1}{2 ^ {2i + 1}} $$

N

最小カットなので、頂点を in と out に分けて最大流。

O

点 $p_i = (x_i ,y_i)$ の寄与を考える。$\arg p _ j \in [\arg p _ i, \arg p _ i + \pi)$ なる $j$ の集合を $S$、それ以外の $j$ の集合を $T$ とすれば、点 $p _ i$ の寄与は次のように書ける。

$$\begin{aligned} \sum _ {j \in S \sqcup T} \vert x _ i y _ j - y _ i x _j \vert &= x_i \left(\sum _ {j\in S} y _ j \right) - y _ i \left(\sum _ {j\in S} x _ j \right) \\ & - x_i \left(\sum _ {j\in T} y _ j \right) + y _ i \left(\sum _ {j\in T} x _ j \right) \end{aligned}$$

偏角ソートしておけば $S, T$ はそれぞれ定数個の区間の和となるので、予め偏角を座圧して BIT に $x _ j$ や $y _ j$ の和を乗せることで答えを差分更新できる。

任意次元Fenwick Treeとその実装

$$ \gdef\defarrow{\overset{\text{def}}{\Longleftrightarrow}} $$

注: 計算量解析が全体的に怪しいです。疑いながら読んでください(?)

問題

$D$ 次元空間の格子点 $\bm{x}\in \Z ^ D$ の重み $w(\bm{x})$ が全て $0$ で初期化されている。以下の種類のクエリが合計で $Q$ 個飛んでくるので、順に処理してください。

  • $\mathrm{add}(\bm{x}, v)$ : 格子点 $\bm{x}$ の重みに $v$ を加算する。即ち、$w(\bm{x})\leftarrow w(\bm{x}) + v$ とする。
  • $\mathrm{sum}(\bm{x}, \bm{y})$ : $\displaystyle \sum_{\bm{x} \leq \bm{z} \lt \bm{y}} w(\bm{z})$ を出力する。

ここで、$\bm{x}=(x_1,\ldots,x_D)$ および $\bm{y}=(y_1,\ldots,y_D)$ に対して、

$$\bm{x} \lt \bm{y} \defarrow \forall i \in \{1,\ldots,D\}.\; x_i\lt y_i$$

と定義します。その他の不等号 ($\leq,\geq,\gt$) に関しても同様に要素ごとの比較で定義します。


$D$ を定数とみなせば、この問題はクエリ先読みと $D$ 次元 Fenwick Tree により $O(Q (\log Q) ^ D)$ 時間で解くことができますが、一般の $\mathrm{sum}$ クエリに関しては定数倍として少なくとも $2 ^ D$ が掛かります (後述)。

なお、$\bm{x}=-\bm{\infty}\coloneqq (-\infty,\ldots,-\infty)$ であるような場合の $\mathrm{sum}$ クエリはこの $2 ^ D$ を落とすことができます。そのため、このケースを特別視して以降は $\mathrm{sum}(-\bm{\infty},\bm{y})$ を $\mathrm{prefix\_sum}(\bm{y})$ と書きます。

解法

$D=1$ の場合

この場合は、$\mathrm{add}$ クエリで現れる座標を先読みして座標圧縮を行うことで、通常の Fenwick Tree により全体 $O(Q\log Q)$ 時間で解くことができます。

$D\gt 1$ の場合

概要

$D=1$ の場合と同様に $\mathrm{add}$ クエリで現れる座標の $1$ 次元目を座標圧縮しておきます。そして、各ノードに $D-1$ 次元 Fenwick Tree を乗せることで $D$ 次元 Fenwick Tree を構築します。これにより、次で述べるような方法で各種クエリに答えることができます。

$3$ 次元 Fenwick Tree

上図は $D$ 次元 Fenwick Tree のノード $N_i$ が $D-1$ 次元 Fenwick Tree であり、即ち再帰的な構造を持っていることを示しています。

クエリ処理

一旦 $D$ 次元 Fenwick Tree の構築が完了したと思うことにして、$\mathrm{add}, \mathrm{sum}$ クエリの実装イメージを示します。通常のFenwick Treeを知っていれば難しくないと思います。

# 1-indexed
# n : ノード数
# N[i] := Fenwick Tree の i 番目のノード

def add((x1, ..., xD), v):
  index = compress(x1) # x1 を座標圧縮して得られる添字
  while index <= n:
    N[index].add((x2, ..., xD), v)
    index += index & -index

def sum((x1, ..., xD), (y1, ..., yD)):
  if x1 >= y1:
    return 0

  result = 0
  
  i = max_index_less_than(y1) # y1未満の最大の座標に対応する添字
  while i > 0:
    result += N[i].sum((x2, ..., xD), (y2, ..., yD))
    i -= i & -i

  i = max_index_less_than(x1) # x1未満の最大の座標に対応する添字
  while i > 0:
    result -= N[i].sum((x2, ..., xD), (y2, ..., yD))
    i -= i & -i

  return result

以降の計算量解析では $D$ は定数とし、また $\log$ の底は $2$ とします。

$\mathrm{add}$ クエリの時間計算量 $T _ {\mathrm{add}}(D)$ は

$$\begin{aligned} T _ {\mathrm{add}}(1) & = O(\log Q), \\ T _ {\mathrm{add}}(D) & \leq (\log Q + 1) \cdot (T _ {\mathrm{add}}(D-1) + O(1)) \quad (D \gt 1) \end{aligned}$$

より $T_{\mathrm{add}}(D) = O( (\log Q) ^ D)$ です *1

$\mathrm{prefix\_sum}$ クエリの時間計算量 $T _ {\mathrm{psum}}(D)$ は、sum において $2$ つ目の while を skip した場合の計算量なので

$$\begin{aligned} T _ {\mathrm{psum}}(1) &= O(\log Q), \\ T _ {\mathrm{psum}}(D) &\leq (\log Q + 1) \cdot (T _ {\mathrm{psum}}(D-1) + O(1)) \quad (D \gt 1) \end{aligned}$$

より $T _ {\mathrm{psum}}(D) = O( (\log Q) ^ D)$ です。

$\mathrm{sum}$ クエリの時間計算量 $T _ {\mathrm{sum}}(D)$ は

$$\begin{aligned} T _ {\mathrm{sum}}(1) &= O(\log Q), \\ T _ {\mathrm{sum}}(D) &\leq (2\log Q + 1) \cdot (T _ {\mathrm{sum}}(D-1) + O(1)) \quad (D \gt 1) \end{aligned}$$

より $T _ {\mathrm{sum}}(D) = O( (\log Q) ^ D)$ です。ただし、上式から分かるように、定数倍として少なくとも $2 ^ D$ が掛かります *2。実用する際は意識するとよいかもしれません。

構築

$\mathrm{add}$ クエリで与えられる格子点の列 $\bm{x} ^ 1,\ldots,\bm{x} ^ k$ をクエリ先読みにより得られているとします。

格子点 $\bm{x} ^ i = (x ^ i _ 1, \ldots, x ^ i _ D)$ の $1$ 次元目の値 $x ^ i _ 1$ を集めてできる (多重でない) 集合 (形式的には $\{ p \mid \exists i \in \{1, \ldots, k\} \; \text{s.t.}\; x ^ i _ 1 = p \}$) の要素を昇順に並べてできる列を $P=(p_1, \ldots, p_n)$ とすると、構築する $D$ 次元 Fenwick Tree のイメージ図は次のようになります。

格子点 $\bm{x} ^ i$ の情報を持つべきノード $N_j$ は $\mathrm{add}$ クエリを処理する場合と同様にして列挙できるので、以下のようにして構築パートを実装できます。

# 1-indexed

points = empty_list()

# クエリ先読みでFenwick Treeに点を登録する
def register_point((x1, ..., xD)):
  points.append((x1, ..., xD))

# p[i] は上で p_i として説明した変数
p = None

# 座標圧縮
def compress(x):
  l = 0
  r = p.size() + 1
  while r - l > 1:
    m = floor_div(l + r, 2)
    if p[m] < x:
      l = m
    else:
      r = m
  return r

def build():
  # pointsの1次元目の値のリスト
  coordinates = empty_list()
  for (x1, ..., xD) in points:
    coordinates.append(x1)
  # coordinatesの要素を昇順に、重複なしで並べたリスト
  p = sorted_unique_list(coordinates)
  n = p.size()

  # 再帰終端
  if D == 1:
    return

  for (x1, ..., xD) in points:
    # p[index]=x1 なるindex
    i = compress(x1)
    # addクエリと同様にして i の親を辿っていく
    while i <= n:
      # N[i]は(D-1)次元Fenwick Tree
      N[i].register_point((x2, ..., xD))
      i += i & -i

  for i in {1, ..., n}:
    N[i].build()

構築にかかる時間計算量を考えます。各点に対して、最も内側の while が回る回数は $O( (\log Q) ^ {D - 1} )$ 回です。従って、$p$ の計算にかかる時間計算量はソートがボトルネックとなり

$$ O( (Q (\log Q) ^ {D - 1} ) \log (Q (\log Q) ^ {D - 1} )) = O(Q (\log Q) ^ {D - 1} (\log Q + D \log \log Q)) $$

と評価できます。$D$ は定数なので、これは $O(Q (\log Q) ^ D)$ です。また、座標圧縮を行う部分の時間計算量は $O(Q \cdot ( (\log Q) ^ {D - 2} \log Q)) = O(Q (\log Q) ^ {D - 1})$ です。

以上より、構築にかかる時間計算量は $O(Q (\log Q) ^ D)$ です。while が回る回数から、空間計算量は $O(Q (\log Q) ^ {D - 1})$ と評価されます。

実装

ここまで python 風の擬似コードを書いておいてなんですが、C++ による実装例を以下に示します。

折り畳み

  • 動作は保証しません
  • 説明とは異なり全て0-indexedです
#include <algorithm>
#include <array>
#include <cassert>
#include <limits>
#include <vector>

namespace compressed_fenwick_tree {
    namespace detail {
        // 座標圧縮
        template <typename CoordinateType>
        struct Compressor {
            using coordinate_type = CoordinateType;

            Compressor() = default;

            void register_coordinate(const coordinate_type& x) {
                xs.push_back(x);
            }

            int build() {
                std::sort(xs.begin(), xs.end());
                xs.erase(std::unique(xs.begin(), xs.end()), xs.end());
                return xs.size();
            }

            // x以上の最小の値の添字(存在しない場合はxs.size())
            int operator()(const coordinate_type& x) const {
                return std::lower_bound(xs.begin(), xs.end(), x) - xs.begin();
            }
        private:
            std::vector<coordinate_type> xs;
        };
    }

    /**
     * @brief D>1の場合のD次元Fenwick Tree
     *
     * @param T 要素型
     * @param op 二項演算の結果を返す関数
     * @param e 単位元を返す関数
     * @param inv 逆元を返す関数
     * @param D 次元
     * @param Index 座標の型
     */
    template <typename T, T(*op)(T, T), T(*e)(), T(*inv)(T), std::size_t D = 1, typename Index = int>
    struct CompressedFenwickTree {
        using value_type = T;
        using index_type = Index;
        // 点を表す型
        using point_type = std::array<index_type, D>;
        // D次元の直方体領域を表す型
        using cube_type = std::array<std::pair<index_type, index_type>, D>;
        // ノードの型
        using node_type = CompressedFenwickTree<value_type, op, e, inv, D - 1, index_type>;

        CompressedFenwickTree() = default;

        /**
         * @brief addクエリで渡される点を予め与えておく
         *
         * @param p 登録する点
         */
        void register_point(const point_type& p) {
            compressor.register_coordinate(p[0]);
            points.push_back(p);
        }

        /**
         * @brief Fenwick Treeを構築する (クエリを投げる前に必ず呼ぶ)
         */
        void build() {
            nodes.assign(n = compressor.build(), node_type{});
            for (const auto& p : points) for (int k = compressor(p[0]) + 1; k <= n; k += k & -k) {
                nodes[k - 1].register_point(tail(p));
            }
            points.clear();
            points.shrink_to_fit();
            for (auto& t : nodes) t.build();
        }

        /**
         * @brief D次元直方体領域に含まれる点の重みの総和を計算する
         *
         * @param cube D次元直方体領域
         */
        value_type sum(const cube_type& cube) const {
            return op(sum(cube[0].second, tail(cube)), inv(sum(cube[0].first, tail(cube))));
        }

        /**
         * @brief 点の重みに加算する
         *
         * @param p 重みを加算する点 (必ずregister_pointで予め渡した点を与える。それ以外の点に対する動作は未定義)
         * @param val 加算する値
         */
        void add(const point_type& p, const value_type& val) {
            for (int r = compressor(p[0]) + 1; r <= n; r += r & -r) {
                nodes[r - 1].add(tail(p), val);
            }
        }
    private:
        // ノードの数
        int n;
        // 1次元目の値を座標圧縮するための構造体
        detail::Compressor<index_type> compressor;
        // addクエリが来る点のリスト
        std::vector<point_type> points;
        // ノード
        std::vector<node_type> nodes;

        /**
         * @brief D次元直方体領域に含まれる点の重みの総和を計算する。(1次元目に下限が無い場合)
         *
         * @param head_r 1次元目の上限
         * @param p 2次元目以降を表すD-1次元直方体領域
         */
        value_type sum(const index_type& head_r, const typename node_type::cube_type& p) const {
            if (head_r == std::numeric_limits<index_type>::min()) return e();
            value_type res = e();
            for (int r = compressor(head_r); r; r -= r & -r) {
                res = op(res, nodes[r - 1].sum(p));
            }
            return res;
        }

        /**
         * @brief array<point_type>またはarray<cube_type>の第二要素以降を取り出す
         *
         * @param p 点またはD次元直方体領域
         */
        template <typename X>
        static constexpr auto tail(const X& p) {
            return tail_impl(p, std::make_index_sequence<D - 1>{});
        }
        template <typename X, std::size_t... Seq>
        static constexpr auto tail_impl(const X& p, std::index_sequence<Seq...>) {
            using return_type = std::conditional_t<
                std::is_same_v<X, point_type>,
                typename node_type::point_type, // Xがpoint_typeの場合
                typename node_type::cube_type   // Xがcube_typeの場合
            >;
            return return_type{ p[1 + Seq]... };
        }
    };

    /**
     * @brief 1次元Fenwick Tree
     *
     * @param T 要素型
     * @param op 二項演算の結果を返す関数
     * @param e 単位元を返す関数
     * @param inv 逆元を返す関数
     * @param Index 座標の型
     */
    template <typename T, T(*op)(T, T), T(*e)(), T(*inv)(T), typename Index>
    struct CompressedFenwickTree<T, op, e, inv, std::size_t(1), Index> {
        using value_type = T;
        using index_type = Index;
        using point_type = index_type;
        using cube_type = std::pair<index_type, index_type>;
        using node_type = value_type;

        CompressedFenwickTree() = default;

        /**
         * @brief addクエリで渡される点を予め与えておく
         *
         * @param p 登録する点
         */
        void register_point(const point_type& p) {
            compressor.register_coordinate(p);
        }

        /**
         * @brief Fenwick Treeを構築する (クエリを投げる前に必ず呼ぶ)
         */
        void build() {
            nodes.assign(n = compressor.build(), e());
        }

        /**
         * @brief 半開区間内の重みの総和を計算する
         *
         * @param l 区間の左端 (閉)
         * @param r 区間の右端 (開)
         */
        value_type sum(const index_type& l, const index_type& r) const {
            return op(sum(r), inv(sum(l)));
        }

        /**
         * @brief 半開区間内の重みの総和を計算する
         *
         * @param range 半開区間
         */
        value_type sum(const cube_type& range) const {
            return sum(range.first, range.second);
        }

        /**
         * @brief 点の重みに加算する
         *
         * @param p 重みを加算する点 (必ずregister_pointで予め渡した点を与える。それ以外の点に対する動作は未定義)
         * @param val 加算する値
         */
        void add(const point_type& p, const value_type& val) {
            for (int r = compressor(p) + 1; r <= n; r += r & -r) {
                nodes[r - 1] = op(nodes[r - 1], val);
            }
        }
    private:
        int n;
        detail::Compressor<index_type> compressor;
        std::vector<node_type> nodes;

        /**
         * @brief prefix query
         *
         * @param p 区間の右端 (開)
         */
        value_type sum(const point_type& p) const {
            if (p == std::numeric_limits<index_type>::min()) return e();
            value_type res = e();
            for (int r = compressor(p); r; r -= r & -r) {
                res = op(res, nodes[r - 1]);
            }
            return res;
        }
    };

    namespace detail {
        template <typename T, T(*e)()>
        T inv_dummy(T) { return e(); }
    }

    /**
     * @brief 任意次元Fenwick Treeで、sumクエリがprefix sumのみの場合
     *
     * @param T 要素型
     * @param op 二項演算の結果を返す関数
     * @param e 単位元を返す関数
     * @param D 次元
     * @param Index 座標の型
     */
    template <typename T, T(*op)(T, T), T(*e)(), std::size_t D = 1, typename Index = int>
    struct CompressedFenwickTreePrefixQuery : private CompressedFenwickTree<T, op, e, detail::inv_dummy<T, e>, D, Index> {
        using base_type = CompressedFenwickTree<T, op, e, detail::inv_dummy<T, e>, D, Index>;

        using value_type = typename base_type::value_type;
        using index_type = typename base_type::index_type;
        using point_type = typename base_type::point_type;

        using base_type::base_type;
        using base_type::register_point;
        using base_type::build;
        using base_type::add;

        /**
         * @brief prefix sumを計算する
         *
         * @param p 領域の右端
         */
        value_type sum(const point_type& p) {
            return base_type::sum(make_cube(p));
        }

    private:
        using cube_type = typename base_type::cube_type;

        cube_type make_cube(const point_type& p) {
            static constexpr index_type neg_inf = std::numeric_limits<index_type>::min();
            if constexpr (D == 1) {
                return cube_type{ neg_inf, p };
            } else {
                cube_type cube;
                for (std::size_t i = 0; i < D; ++i) cube[i] = { neg_inf, p[i] };
                return cube;
            }
        }
    };
} // namespace compressed_fenwick_tree

使用例

結び

多分どこか間違えてるので指摘してください。

*1: $T _ {\mathrm{add}}(D) \leq \sum _ {i = 0} ^ {D} c ^ D _ i (\log Q) ^ i$ なる係数列 $c ^ D$ を $Q$ の値に依存しないように $D$ が小さい順に決めることができて、これに対して十分大きな係数 $C _ D$ を取って $\forall x \gt 0.\; \sum _ {i = 0} ^ {D} c ^ D _ i x ^ i \leq C _ D x ^ D$ を満たすようにできます。これは当然 $D$ に (のみ) 依存しますが、今回の解析では $D$ は定数なので $C_D$ も定数です。従って、$D$ を定数とすれば $T _ {\mathrm{add} }(D) = O( (\log Q) ^ D)$ と結論付けられます。そして、これは $T _ {\mathrm{sum} }(D)$ や $T _ {\mathrm{prefix\_ sum} }(D)$ に関しても同様です。

*2:この $2 ^ D$ は $\mathrm{prefix\_sum}$ を用いて包除原理で超直方体領域内の和を計算していることで掛かる係数と解釈できます。