Sparse TableにO(log N)で非冪等性演算range foldを載せてみた

## はじめに

ネタ記事です。該当の機能がほしいならDisjoint Sparse Tableを使ったほうが良いことを前提に読んでください。

はじめましての人ははじめまして。[alumi](https://x.com/alumi_02)です。
最近Sparse Tableについて学習しているのですが、その際にSparse Tableにも非冪等性の演算に関するrange foldが載せられそうだな～と思いつき実際に載せられたので、Sparse Tableの実装とともに紹介します。全国のSparse Table未学習者の参考になればうれしいです。

## Sparse Tableってなに？

静的配列に対する区間foldを、空間$O(N\log{N})$、前計算$O(N\log{N})$、区間foldクエリ$O(1)$で答えるデータ構造です。ただし静的配列の要素はfold演算$\oplus$について以下の性質を満たす必要があります。

- [結合性](https://en.wikipedia.org/wiki/Associative_property)：$(A \oplus B) \oplus C = A \oplus (B \oplus C)$
- [冪等性](https://en.wikipedia.org/wiki/Idempotence)：$A \oplus A = A$
    - ex. min, max, gcd, lcm, bitwise-and, bitwise-or, $\cap$, $\cup$

結合性と冪等性はfold演算を$O(1)$でこなすために必要な性質です。例えばある2つの区間fold値$A_0 \oplus A_1 \oplus A_2 \oplus A_3$と、$A_2 \oplus A_3 \oplus A_4 \oplus A_5$とを演算すると、

```math
\begin{eqnarray*}
&& (A_0 \oplus A_1 \oplus A_2 \oplus A_3) \oplus (A_2 \oplus A_3 \oplus A_4 \oplus A_5) \\
&=& A_0 \oplus A_1 \oplus (A_2 \oplus A_3) \oplus (A_2 \oplus A_3) \oplus A_4 \oplus A_5 \\
&=& A_0 \oplus A_1 \oplus A_2 \oplus A_3 \oplus A_4 \oplus A_5
\end{eqnarray*}
```

というふうに同じ区間に対する結果を圧縮でき、別の区間fold値を計算することができます。

つまり「いい感じの長さの区間fold値を前計算しておけば、ほしい区間fold値を前計算した2つの区間fold値の演算によって$O(1)$で求めることができる」というのがSparse Tableの基本コンセプトです。

### 前計算

前計算で作る「いい感じの長さの区間fold値」は下の図のように、任意の始点$j (\in [0,N))$から、$j$からの長さが操作対象の静的配列の要素数$n$を超えない2冪の長さの範囲のfold値となります。図の$L$は2冪となる範囲の長さの指数（$level$）であり、範囲の長さが4なら$level$は2となります。このようにすると、任意の始点$j$の長さ$2^{level}$の区間fold値は`table[level][j]`で参照できます。

![fig1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3665424/7e319536-1341-5080-2e9f-2671a6762c74.png)

tableの構成は各始点ごとに行うよりも、各$level$ごとに行ったほうが効率的です。つまり、

- $level0$（長さ1）のfold値は操作対象の静的配列そのものなのでそのままコピー
- $level1$（長さ2）の始点$j$のfold値は、$level0$の隣り合う2つ$j, j+1$を演算すれば得られる
- $level2$（長さ4）の始点$j$のfold値は、$level1$の距離$2$挟んで隣り合う2つ$j, j+2$を演算すれば得られる
- $level3$（長さ8）の始点$j$のfold値は、$level2$の距離$4$挟んで隣り合う2つ$j, j+4$を演算すれば得られる
- $leveli$（長さ$2^i$）の始点$j$のfold値は、$leveli-1$の距離$2^{i-1}$挟んで隣り合う2つ$j, j+2^{i-1}$を演算すれば得られる

という流れで構成できます（下図参照）。

![fig2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3665424/ee1209dd-57ab-815f-4406-7cced1cf2db9.png)

$j$の最大値は、0を始点にしたとき$n$を超えない2冪の指数となるので、`n.bit_length() - 1`となります（`floor(log2(n))`と同値）。0から数えると$level$の階層の高さ$h$は`n.bit_length()`となります。
tableの大きさはほぼ$n \times h \simeq n\log{n}$となります。厳密には$level0$が$n$、$level1$が$n-1$、$level2$が$n-3$、...$leveli$が$n-2^i+1$となるので$$\sum_{i=0}^{h-1}(n-2^i+1)=(n+1)h-2^h+1$$となります。$(n+1)h$の項に対して$2^h$は十分小さいので、空間計算量は$O(N\log{N})$となります。時間計算量も同様に$O(N\log{N})$となります。ただ$level$が上がるに連れて`table[level]`の要素数が2冪で減少していくのでスカスカになります（多分これが"[Sparse](https://eow.alc.co.jp/search?q=sparse)" Tableの名前の由来）。

### 任意区間fold値の計算

上記のように2冪の長さの範囲のfold値を前計算すると、任意区間$[l,r)$のfold値は次の2つの区間fold値によって計算できます。

1. 左端$l$から始まる$r-l$を超えない最大の長さの区間fold値
    - ある数$x$を超えない2冪の指数は`x.bit_length() - 1`で計算できます（ex. 7(111)→2、1729(11011000001)→10）。よって`level = (r - l).bit_length() - 1`と計算でき、1.を満たす区間fold値は`table[lebel][l]`で参照できます
2. 1. と同じ長さで、右端が$r-1$で終わるような区間fold値
    - 1. と同じ長さなので`level`は同じ。右端が$r-1$で終わり、かつ区間の長さが$2^{level}$の区間の左端は算数で$r - 1 - 2^{level} + 1 = r - 2^{level}$と計算でき、よって2.を満たす区間fold値は`table[lebel][r - 2**level]`で参照できます

ちなみに区間$[l,l+2^{level})$, $[r-2^{level},r)$によって区間$[l,r)$は隙間なく被覆されます。被覆されないとすると$l+2^{level} < r-2^{level} \Rightarrow r - l > 2 \cdot 2^{level}$となりますが、$level$は定義から$r - l$を超えない最大の2冪の指数なので矛盾してしまいます。
よって区間$[l,r)$のfold値は、`level = (r - l).bit_length() - 1`として、
```math
fold\_range(l,r) = table[level][l] \oplus table[level][r - 2^{level}]
```
と計算できます。

![fig3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3665424/721b5aed-aca3-c375-8dfc-d783a43b42db.png)

以下はここまでの話をまとめた実装です。

```python
class SparseTable:
    __slots__ = ("fold_op", "table")

    def __init__(self, init_arr: list, fold_op) -> None:
        self.fold_op = fold_op
        n = len(init_arr)
        h = n.bit_length()
        self.table = [init_arr[:]]
        self.table.extend([0] * n for _ in range(h - 1))

        for i in range(1, h):
            ct, pt = self.table[i], self.table[i - 1]
            for j in range(n - (1 << i) + 1):
                ct[j] = fold_op(pt[j], pt[j + (1 << (i - 1))])

    def fold_range(self, l: int, r: int):
        """Fold to the range [l, r) for an idempotent operation."""
        level = (r - l).bit_length() - 1
        return self.fold_op(self.table[level][l], self.table[level][r - (1 << level)])

spt = SparseTable(A, max)
spt = SparseTable(A, min)
```

tableの構成で`table = [[0] * n for _ in range(h)]`, `table[0] = A[:]`ではなく、
`table2= [A[:]]`, `table.extend([0] * n for _ in range(h - 1))`としているのは、後者のほうが若干早いからです。前者のほうがリーダブルではあるのでお好みで変えてください。

また`self.table.extend([0] * n for _ in range(h - 1))`としているのですべての$level$でn個分のスペースを確保していますが、以下のようにすると必要な分しか各`table[level]`に格納しないので、前述のSparse性から使用メモリ数を若干減らせます。

```diff_python
class SparseTable:
    __slots__ = ("fold_op", "table")

    def __init__(self, init_arr: list, fold_op) -> None:
        self.fold_op = fold_op
        n = len(init_arr)
        h = n.bit_length()
        self.table = [init_arr[:]]
-       self.table.extend([0] * n for _ in range(h - 1))
+       self.table.extend([] * n for _ in range(h - 1))

        for i in range(1, h):
            ct, pt = self.table[i], self.table[i - 1]
            for j in range(n - (1 << i) + 1):
-               ct[j] = fold_op(pt[j], pt[j + (1 << (i - 1))])
+               ct.append(fold_op(pt[j], pt[j + (1 << (i - 1))]))
```

## 非冪等性演算range foldを載せる方法

やっと本題。
冪等性を満たさない演算を以下非冪等性演算とします。
非冪等性演算は冪等性を満たさない（小泉構文）ので、区間foldする際には**拾ってくる区間がそれぞれ重ならない**ように拾ってくる必要があります。
幸い任意の始点における2冪の長さのfold値は前計算されているので、**セグメント木のfoldのようにlevelを遡りながら、範囲を超えない最大の長さの区間を結合**しまくれば行けそうです。つまり区間$[l,r)$のfold値がほしい場合、

1. `while l < r`
    1. $r-l$を超えない最大の長さの`level`を計算する。`level = (r - l).bit_length() - 1`とし、`res = fold_op(res, table[level][l])`で`res`を更新
    2. 見る区間の左端を更新する。$2^{level}$分は見たので`l += 2**level`として$2^{level}$分シフト
2. `return res`

として区間$[l,r)$のfold値が取得できます（下図参照）。

![fig4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3665424/4045248e-bfe8-f11a-8ad0-0ee7e525c5f4.png)

### 実装

```python
    def fold_range(self, l: int, r: int):
        """Fold to the range [l, r) for any associative operation."""
        res = None
        while l < r:
            level = (r - l).bit_length() - 1
            if res is None:
                res = self.table[level][l]
            else:
                res = self.fold_op(res, self.table[level][l])
            l += 1 << level
        return res
```

ただし、
- 各始点ごとに前計算している & 左端からのみfoldしているのでセグメント木よりも空間・時間計算量ともに劣る
- 加減やxorなら累積和・累積xorで区間foldに対応できる
- Sparse Tableとほぼ同じ空間・時間計算量で非冪等性演算foldできるDisjoint Sparse Tableがある
ので、実用上のメリットはないでしょう。悲しい。

## 結び

というわけでネタ記事に託つけたSparse Tableの解説でした。
あまり存在感のないSparse Table（wikipediaの記事すら無い）ですが、初学者にもちょうどいい難易度のデータ構造なので、もっと市民権を得て活用できる機会が増えてくれれば嬉しいなぁ...
よかったら左上のいいねボタン押してください🙇

## 参考

- [Sparse Table - data-structures](https://scrapbox.io/data-structures/Sparse_Table)
- [Range Minimum Query #アルゴリズム - Qiita](https://qiita.com/okateim/items/e2f4a734db4e5f90e410#sparse-table-algorithm)
- [Sparse Tableをそれとなく理解する #Python3 - Qiita](https://qiita.com/recuraki/items/0fcbc9e2abbc4fae5f62)
- [Sparse Table](https://aquarchitect.github.io/swift-algorithm-club/Sparse%20Table/)