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}$ を用いて包除原理で超直方体領域内の和を計算していることで掛かる係数と解釈できます。

AtCoder 橙/Codeforces IGMになりました

はじめに

両者とも人生最後の上方向の色変になる可能性があるので書きます。今まで色変記事は書いたことがなかったので、黄->橙というよりは今までにやったことを振り返る感じです。

うれしい!

うれしい!

精進・コンテスト参加

全体

Solved (全体)

AtCoder

レーティング遷移 (AtCoder)

2019年10月にAtCoderに初めて参加しました。青になるまではかなり順調だったんですが、そこから先でかなり苦労しました。今年1月に2371を踏んだ時はこのまま橙になれるかなと思っていましたが、このあとかなりグダグダして半年以上かかりました。

AtCoderに関しては、始めてから現在までABC/ARC/AGCはほとんど全て参加しています。

Solved (AtCoder)

Rated

AtCoderでの精進は↑の通りです。Ratedを中心に結構埋めたと思います (水・青は100%、黄・橙は1問欠け、赤↑は4割と少し)。Streakを繋ぐのは苦手なので意識してやったことはほとんどありません。試験管の問題に関しては、解いていて苦な問題もありますが面白い問題もたくさんあるという印象です。

Unratedコンの問題は全然埋められていないので、その辺りも今後手を出すかもしれません (CF Div.1の方が優先度は高いかも)。

参考までに、ABC/ARC/AGCの埋め具合を貼っておきます。画像がかなり縦に長いので注意してください。

ABCの埋め具合 画像の縦の長さに制限があるみたいなので分割しています。
8問時代 Gは埋めるとしてH/Exも埋めたいんですが、単純に難易度が高くて埋まりません。ABCなので解けないなら解説open推奨なのかな。
Solved (ABC 8問時代)
6問時代
Solved (ABC 6問時代)
4問時代 Dで唯一解けてないのは D - 25個の整数 です。また今度通したい
Solved (ABC 4問時代)

ARCの埋め具合 6問になってからのFはほとんど手が出てません。最近は全然復習できていないので、やっていきたいです
Solved (ARC)

ACCの埋め具合 後半の問題も古いものであれば比較的手が出やすい印象です
Solved (AGC)

Codeforces

レーティング遷移 (Codeforces)

私は競技を始めた当初はJavaを使っており、CodeforcesはTLが厳しいがちだという噂を聞いていたので中々参加できていませんでした *1。ただ、ICPCに参加するにあたってチームメイトが2人ともC++使いだったこともありC++への乗り換えを本格的に始めていたので、2021年4月にようやくCodeforcesに手を出しました。この頃はAtCoder 2000-2100の所謂 "黄溜まり" でさまよっていた頃で、刺激が欲しかったというのもあります *2

同年6月にGrandmaster (薄赤)になったまでは良かったんですが、それから5ヶ月ほどサボってしまい、また11月頃に復帰して今に至るという感じです。時間帯の問題もあってDiv.1すら毎回参加とはいかないんですが、ほどほどに参加しています。

Solved (Codeforces)

Codeforcesでの精進は↑の通りです。そこそこの問題数を解いているように見えますがこれは過去のDiv.3埋めが大きく寄与していて、実はDiv.1に関してはほとんど埋められていません。最近はECRを埋めたりしていますが、そろそろDiv.1から逃げるのをやめないとなぁと思っています。

CodeChef

今年7月に初めて参加して4回コンテストに参加しただけです。この4回に限った印象としては特に悪い印象はなく、問題も面白かったです。

ライブラリ整備について

suisen-cp.github.io

↑は私の競技用ライブラリです。結構量が多い方だと思いますが、あまり他の方のライブラリを見たことがないのでそうでもないかもしれません。

ライブラリ整備が好きなのでLibrary Checkerなんかもそこそこ埋めているんですが、レートに直結している感覚はありません。特にAtCoderに関しては、ARC以上 (特にAGC) では「考察 >> (ライブラリを含む) 実装」な問題を選んで出しているというような内容を含むこどふぉブログを読んだことがある気がします (ソースは忘れました)。Codeforcesに関しては、例えばHLDなんかはそこそこの頻度で貼っているので、基本的なものに関しては役に立っていそうです。

言語の変更について

Codeforcesのところでも言及したんですが、2021年4月~5月ごろを境にメインで使用する言語をJavaからC++へと変更しました。当初は

  • Javaで書いた大量のライブラリを捨てるのが勿体ない
  • 新しい言語で最低限必要なライブラリを準備しないといけない
  • そもそも言語を覚えないといけないので、一時的にパフォーマンスが落ちるかも

みたいなことを考えて足踏みをしていましたが、1つ目に関しては単なる気分の問題なので良いとして、2つ目に関してはAtCoder Libraryが導入されたこと、3つ目に関してはAtCoderが考察>>実装な問題を好むということで、思ったより障壁は小さかったという印象です。

おわりに

書きたいことがあまり思い浮かばなかったので淡白になってしまいました。何か気になることがあれば聞いてもらえれば追記するかもしれません。

*1:これは私が勝手に足踏みしていただけで、JavaPythonCodeforcesに参加している人も一定数います。参加を推奨しないということではありません。

*2:レーティングシステムの性質なのか、Codeforcesではレートが乱高下しやすいので刺激は十分に得られると思います。

ICPC 国内予選 2022 参加記(Bu-Bu-Du-Ke)

suisen・otera・kenchoの3人でチーム「Bu-Bu-Du-Ke」として参加し、全体12位・学内3位の成績で国内予選を通過しました。自分とoteraさんにとっては最終年だったので、何とか突破出来て嬉しいです。

予選前

Heno Worldの権利が無ければ突破は基本的に出来るんじゃないか、みたいなことを考えていました。が、まだ権利があったらしく、KUB1ppのことも考えると全体20位は必須かな〜という話をしていました。

予選本番

A:suisen、B:otera、C:kenchoの割り当てで開始。Aは簡単だったので丁寧に丁寧に出力ファイルを確認して提出してAC。Dに取り掛かる。

oteraさんのBの読解が大変そうで、Cを通したkenchoさんとoteraさんの2人で問題の理解を擦り合わせていた。

その間自分はDの考察をしていて、queueの先頭であることが確定している集合 S を持ってmin SとT[i] の大小関係によって決められそうなことが分かった。実装するもサンプルが合わず、デバッグ開始。バグを見つけて1ケース合うたびに他が壊れる、みたいなことを繰り返して辛くなっていたら、oteraさんから「Dまだあんまり通ってないですよ」みたいなことを聞いて安心。最終的に、min S=T[i]の場合、今までpopした値の最大値Mに対してM>T[i]ならT[i]から始まるqueueを作ってはいけないというケースが漏れていることに気付いて修正、S[i]とT[i]を書き間違えるバグで時間を溶かしつつも何とかAC。デバッグをしてる間にoteraさんがBを通して、kenchoさんとEを考えていた。4完の速度が結構速く、蓋を開けてみるとこの時点で予選通過が確定していた。

Eを2人に任せて、構文解析の練習をしていた自分がFに取り掛かっていた。「位置iから始まる部分で、左からl文字、右からr文字消す場合の辞書順最小」を求める再帰を書こうとする。構文解析パートは練習の成果もあり5分足らずで書けたが、l,rの管理でめちゃくちゃバグらせてサンプルすら合わない。

自分がバグと格闘していると、oteraさんからEは無限場合分けをすれば通せるということを聞く。kenchoさんは場合分けを網羅できるように色んなテストケースを作ったり、oteraさんのコードのデバッグを手伝ったりしていた。oteraさんは200行の実装を書ききり、2ペナしつつも何とかEが通る。kenchoさんの用意したテストケースがかなり役に立っていたっぽい。感謝…

Eが通るまでの間、oriKoU3も4完でEへの提出がありかなり怖かった。4完時点のペナルティはうちの方が小さく、oriKoU3もEで2つペナを出していたので、うちがEを通した時点でoriKoU3がEを通しても順位が逆転しないことが確定。これは通ったんじゃないかという話をしていたが、今考えるとF以降を通されると逆転する可能性は残っていたはず。

Eが通るまでの間、自分はというと、oteraさんが全部手打ちで場合分けを書いているのを応援しながら、ずっとデバッグしていて辛くなっていた。Fは出力の長さすら合っておらず、バグの原因が判明しないままコンテスト終了。

Fが通せなかったのは個人的に悔いが残る部分だが、予選通過の方が大事なのでok。

予選後

軽くコンテストの振り返りをした。最後僅かな時間でkenchoさんが考えていたGが正しかったらしい。結果論だけど、Fは放棄した方が良かったかなと思ったりした。

振り返りの後は近くの店でご飯を食べた。trineutronさんに遭遇し、話をした。

感想

とにかく予選突破できて嬉しい。横浜でも頑張ります。その頃には橙×3のチームになってたらいいな。

第10回 アルゴリズム実技検定 バチャ録

83:42+2ペナで完走しました.

A - 3枚のカード (5:03 (+5:03))

格好付けて minmax_element で出力しようとしたら色々おかしくなって原因を探ったりしてたら 5 分掛かった (?)

B - 花束 (6:00 (+0:57))

赤と青でそれぞれ花束を作ると誤読して $\lfloor X/A \rfloor + \lfloor Y/B \rfloor$ を出力するコードを書いたら合わない.よく読むと和じゃなくて $\min{\lfloor X/A \rfloor, \lfloor Y/B \rfloor}$ だった.

C - Go Further (8:47 (+2:47))

ちょっと面倒.$\lfloor \log _ {10} \alpha \rfloor \bmod 3$ で場合分け.

D - ハイスコア (11:06 (+2:19))

差分更新の問題だと思ってコードを書き上げてから愚直でいいことに気付いた.ただ,この問題は愚直より差分更新の方が楽そう.

E - 良い日付 (17:56 (+6:50))

こういう問題が出ると困ってしまう.Python の datetime にこういうのがあるので,検索して使い方を見ながら書いた.

答えの上限は 3000/03/03 なので,日付を $1$ つずつ試しても高々 $366\times 100 + 100$ 回の操作で終わる.余裕かと思って生 Python で出すと 601 ms で焦った.

F - 地図の塗り分け (20:50 (+2:54))

こういうのは $4$ 隣接のうち $[0,H)\times [0,W)$ に含まれるものだけを列挙する便利ライブラリを用意しておくと楽だと思います.

G - 方程式 (32:11 (+11:21))

まず二分法が思い浮かぶ.初めは $x$ 軸に接するケースをケアしないといけないと思ったけど,$f' = 0$ となるのは $x ^ 4 = -b/5a$ のときで,制約から $-b/5a \lt 0$ が分かるので実はケアしなくても良い.次にケアする必要があるのは $f(1) = 0$ や $f(2) = 0$ となるケース.これは,二分法の初期区間を $[1 + \varepsilon, 2 - \varepsilon]\ (\varepsilon = 10 ^ {-10})$ とすることで回避できる.

あとで気付いたけど,$f(1)=f(2)=0 \Rightarrow 31a+b=0$ の一方で $a,b\gt 0$ だから結局 $f(1)\neq 0$ または $f(2) \neq 0$ が成り立って,つまり符号を確定できるので $\varepsilon$ だけずらしたりしなくても良かったかも.

H - 連結成分 (35:12 (+3:01))

貼るだけ.

連結成分を取得できる Union Find | cp-library-cpp

BFSをすると,密な連結成分ばかり指定された場合に計算量が壊れるので注意.

I - 対称変換 (42:29 (+7:17, 1 WA))

$y$ 軸に平行な直線で対称移動すると,$S$ の $x$ 座標最小の点は $T$ の $x$ 座標最大の点に移るはずなので,対称軸が一意に定まる.$x$ 座標と $y$ 座標を入れ替えて $2$ 回試すことで $x$ 軸に平行な直線での対称移動を考慮できる.

ソート忘れで 1 WA.

J - 区間の期待値 (45:20 (+2:51))

$A$ をソートすると,寄与が簡単に書ける.総和は $\displaystyle \sum _ {i = 1} ^ {N} A _ i \cdot \binom{i - 1}{K - 1} - \sum _ {i = 1} ^ {N} A _ i \cdot \binom{N - i}{K - 1}$ になるので,これを $\displaystyle \binom{N}{K}$ で割ったものが答え.

K - 旅行計画 (50:39 (+5:19))

問題文で与えられている有向グラフ $G$ と,辺を逆向きに張った有向グラフ $H$ を構築する.$G$ における $1$ から $k$ までの最短距離と $H$ における $N$ から $k$ までの最短距離の和が答え.ただし,どちらか一方で到達不可能なら $-1$ を出力.

L - N mod M (57:12 (+6:33))

厄介なのは $\underbrace{1\cdots 1}_{d個} = \dfrac{10 ^ d - 1}{9}$ を $\mathrm{mod}\ M $ で計算したくなる点 ($\mathrm{mod}\ M $ において $9$ の逆元が存在するとは限らない).

ただ,これも繰り返し二乗法の要領で半分に割って再帰的に計算すれば $\Theta( (\log d) ^ 2)$ で計算出来る.$10 ^ d$ の計算も同時にやれば $\Theta(\log d)$ に落とせるが,この問題では $\Theta( (\log d) ^ 2)$ で十分だった.

M - ランキング (63:13 (+6:01) )

Binary Trie や平衡二分探索木などを使えばよい.Binary Trie で適当にやったら 2499 ms で冷や汗.

N - 400億マス計算 (65:28 (+2:15) )

畳み込んで非零の係数を数える.

O - 3-順列 (83:42 (+18:14, 1WA) )

$(A _ i, B _ i)$ に辺を張ったグラフを構築する.各連結成分に関して,それが二部グラフでなければ,その成分に属する頂点は全て $3$ で割った余りが $0$ でなければならない.一方で,二部グラフであれば,彩色時に色 $0$ で塗った頂点数を $X$,色 $1$ で塗った頂点数を $Y$ とすれば,次の $3$ つの可能性がある.

  1. $X+Y$ 個の頂点は全て $3$ で割った余りが $0$
  2. $X$ 個の頂点は $3$ で割った余りが $1$,$Y$ 個の頂点は $3$ で割った余りが $2$
  3. $X$ 個の頂点は $3$ で割った余りが $2$,$Y$ 個の頂点は $3$ で割った余りが $1$

あとは,部分和問題なので $\Theta(N ^ 3)$ の DP が生える.bitset で高速化出来る形なので,$w$ をワードサイズとして $\Theta(N ^ 3 / w)$ へと計算量を削減でき,これは十分高速.

彩色パートのBFSで実装をミスっていたのと,$1,\ldots, N$ を $3$ で割った余りを集計するパートで間違えて $0,\ldots, N-1$ を $3$ で割った余りを集計するミスがあり,1WA.同時に両方見つけられて良かった.

仮想関数を回避するテク [C++]

(追記 2022/06/22 22:35)

本記事の内容は Curiously recurring template pattern (CRTP) と呼ばれるものらしいです.より正確な情報や関連する情報を得たい場合は "Curiously recurring template pattern" や "CRTP C++" などと調べると良さそうです.

(追記終)

本文

競プロのライブラリを書いていると,継承されることを想定した基底クラスでデフォルトの動作を指定して,継承先で動作を微妙に変更したいことがあります (例えば平衡二分探索木において,添字アクセス・区間和取得・区間作用・反転などをサポートするかどうかで回転などにより木構造を変更した場合に更新すべき情報が変化します).

ベースのコードを書いた後はコピペして色んなバージョンを作ってもよいのですが,後でバグが発覚したり,インタフェースを変更したくなったときの修正箇所が多くて辛くなりがちです.

この問題への対処として仮想関数を用いる方法がありますが,私は使ったことがないのでここでは解説しません (できません).1 つ言えることとしては,仮想関数の呼び出しは通常の関数よりもオーバーヘッドが大きくなるということです.競プロライブラリにおいて速度は重要なので,これも避けたいです.

そこで,次のようなトリックを用います.

#include <iostream>

// Derived は base を継承することを想定
template <typename Derived>
struct Base {
    static void f() {
        // デフォルトの動作
        std::cout << "default f()" << std::endl;
    }
    static void g() {
        // デフォルトの動作
        std::cout << "default g()" << std::endl;
    }
    // 継承先の f() を呼ぶ関数
    static void call_f() {
        Derived::f();
    }
    // 継承先の g() を呼ぶ関数
    static void call_g() {
        Derived::g();
    }
};

// 全てデフォルトの動作をするクラス
struct Default : Base<Default> {};

// f の定義だけ上書きしたクラス
struct Derived : public Base<Derived> {
    using base_type = Base<Derived>;

    // f の定義を上書き
    static void f() {
        // デフォルトの動作を呼び出す (必要なければ呼ばなくても OK)
        base_type::f();
        // 差分の動作を記述
        std::cout << "extended" << std::endl;
    }
};

int main() {
    std::cout << "Default::call_f():" << std::endl;
    Default::call_f();
    std::cout << "Derived::call_f():" << std::endl;
    Derived::call_f();
    std::cout << "Default::call_g():" << std::endl;
    Default::call_g();
    std::cout << "Derived::call_g():" << std::endl;
    Derived::call_g();
    return 0;
}

このコードの実行結果は次のようになります ([C++] gcc 11.1.0 - Wandbox から確かめられます).

Default::call_f():
default f()
Derived::call_f():
default f()
extended
Default::call_g():
default g()
Derived::call_g():
default g()

何をやっているかですが,継承されることを想定した基底クラス Base を定義し,テンプレート引数に継承先のクラス Derived を与えるようにします.継承先で上書きされた関数を呼び出したい場合は Derived::f() のように書けばよいです.

継承先のクラス Derived では Base<Derived> を継承し,定義を上書きしたい関数のみを書きます.例えば上の例では,関数 f の定義を上書きすることで Derived::call_f() 内では Derived 内で定義した f が呼び出されてくれます.また,上書きする必要のない関数に関しては,何も書かなくてもよいです (Default::call_f()Derived::call_g() などの結果を見ると分かると思います).

使用例

私の赤黒木ライブラリはこのテクを使用して実装されています.継承先の各クラスは基底クラスに対して比較的差分を小さく抑えた実装になっていると思います.

おわりに

そこそこ面倒くさいので,もう少し簡単なやり方があったら教えて欲しい.C++に詳しい人いませんか.