radix heapの計算量解析



## はじめに

どうも。aqualengthです。[^1]

この記事は、[「木 Advent Calendar 2023」](https://adventar.org/calendars/9821)の22日目です。
昨日の記事は[blood](https://atcoder.jp/users/blood)の[C言語で木構造を実装をするときのテクニック](https://qiita.com/aqualength/items/8f857d24eeb74d0ed981)でした。
明日の記事は[31536000さん](https://adventar.org/users/43927)の(ここに記事のタイトルを入力)です。

償却計算量の解析は最悪計算量や平均計算量の解析と比較して難しいと思います。[^2]
そこで、今日は練習としてradix heapの計算量を解析します。

なお、ちゃんとした償却計算量の解析は初めてなので、間違ってたらすみません。(わりとお気持ち寄りの解説なので、数学的には厳密じゃないと思います)


## 目次

- [はじめに](#はじめに)
- [目次](#目次)
- [前提条件](#前提条件)
- [radix-heapの提供する操作](#radix-heapの提供する操作)
- [償却計算量の解析について](#償却計算量の解析について)
- [ポテンシャル関数の定義](#ポテンシャル関数の定義)
- [pushの解析](#pushの解析)
- [decrease-keyの解析](#decrease-keyの解析)
- [popの解析](#popの解析)
- [考察](#考察)
- [ところで](#ところで)

## 前提条件

- この記事で扱う計算量は時間計算量のことを指します。
- radix heapの実装は[3種類ほど](https://dl.acm.org/doi/abs/10.1145/77600.77615)ありますが、今回は[よすぽさんの実装](https://yosupo.hatenablog.com/entry/2015/04/03/224649)を解析します。
- データ数をN,データの最大ビット長をkとします。
- 四則演算およびビット演算は定数時間で完了するものとします。
- bsr()は定数時間で完了するものとします。そういうことにしておいてください[^3]

## radix-heapの提供する操作

以下の3つの操作を提供するものとします。

- push(i,n) … 管理番号をiとする値nをヒープに追加します。
- pop(*i,*n) … 最小値を削除し、その管理番号をi,nの参照先に返します。
- decrease-key(i,n) … 管理番号iの値をnまで減少させます。

具体的な実装は省略します。
[よすぽさんの実装](https://yosupo.hatenablog.com/entry/2015/04/03/224649)を基に解析するので特に問題はありません。

<!---
## 実装



C言語による32ビット実装は以下の通りです。
なお、動作検証はしていません。

``` c
#define MAX 100100 
#define bsr(n,ans) do{ans=0;if(n<=0){ans=-1;n=0;}if(n>>16){ans+=16;n>>=16;}if(n>>8){ans+=8;n>>=8;}if(n>>4){ans+=4;n>>=4;}if(n>>2){ans+=2;n>>=2;}if(n>>1){ans+=1;n>>=1;}}while(0)

#define radixheap_init(heap) do{heap.node[0].value=34;}while(0)
struct heap_node{
	unsigned value;
	unsigned next; //同バケット内の次の要素へのインデックス
	unsigned back; //同バケット内の前の要素へのインデックス
	unsigned idx;
};
typedef struct radix_heap{
	struct heap_node node[34+MAX];
}radixheap;

void init(radixheap *heap){
	heap->node[0].value=34; //ヒープ内のデータ数
	heap->node[0].next=0; //最後に取り出された値
}

void push(radixheap *heap,unsigned i,unsigned x){
	unsigned last=heap->node[0].value;
	
	heap->node[0].value++;
}

```

なお、解析は[よすぽさんの実装](https://yosupo.hatenablog.com/entry/2015/04/03/224649)を基に行うため、C言語での実装は一旦忘れてください。

--->

## 償却計算量の解析について

直感で計算量を理解する以外の方法として、ポテンシャル関数を使った解析がよく使われます。

ポテンシャル関数を使った解析は[別の記事](https://trap.jp/post/883/)を見てもらう方が早いですが、直感で理解する方法を数学的に説明したものと思って大丈夫です。

ポテンシャル関数はおそらく[天から下りてくるもの](https://dictionary.goo.ne.jp/word/%E5%A4%A9%E4%B8%8B%E3%82%8A%E7%9A%84/)であるため、償却計算量の解析は難しいです。

## radix-heapの性質

[よすぽさんの記事](https://yosupo.hatenablog.com/entry/2015/04/03/224649)によるとradix heapにはいくつかの性質が成り立ちます。

解析に必要な部分を抜き出すと

- v[i+1]の値は必ずv[i]の値より大きい
- 再振り分けされたv[i]の値は必ずv[j](j<i)へ行く

です。


加えて

- bsr(x^last)はbsr(x)以下になる

という性質も重要です。

## ポテンシャル関数の定義

バケットvについて、各バケットv[i]に入ってる要素数をn[i]とします。
ポテンシャル関数を
$$
\sum_{i=1}^{33} i \times n[i]
$$
と定義することにします。

## pushの解析

push一回で実際にかかる計算量はbsrを計算してリストに追加するだけです。
bsrの計算もリストへの追加も$O(1)$[^3]なので、合計で$O(1)です。

ポテンシャル関数は最大でbsr(x)==kだけ増加するので、償却計算量は
$O(k)$です。

## decrease-keyの解析

decrease-key一回で実際にかかる計算量はpushと同じく$O(1)$です。

「v[i+1]の値は必ずv[i]の値より大きい」という性質よりポテンシャル関数は増えることはなく、最高でも0しか増加しません。

よって、償却計算量は$O(1)$です。

## popの解析

### v[0]が空でない場合

v[0]内の値は全て最小値なので、v[0]の値を一つ削除すれば良く、実際の計算量は$O(1)$

ポテンシャル関数は増加しないため、償却計算量は$O(1)$です。

### v[0]が空の場合

ここで、v[i]が空でない最小のiを考えます。

v[i]が空でない最小のiを検索するのに$O(i)$
また、v[i]の最小値を求めるのに$O(n[i])$
v[i]の値を全て再配置するのに$O(n[i])$の計算量がかかるため、
実際の計算量は$O(i+n[i])$です。

ポテンシャル関数はv[i]の最小値を削除することでiだけ減少し、
「再振り分けされたv[i]の値は必ずv[j]\(j<i)へ行く」という性質から再配置によって最低でもn[i]以上減少します。
合計でi+n[i]以上のポテンシャルが減少するため、
償却計算量は$O(1)$です。


## 考察
解析の結果、radix heapのpushにかかる計算量は$O(k)$であることが分かりました。
二分ヒープやフィボナッチヒープのpushが$O(1)$[^4]なことを考えるとpushの計算量が意外と大きいことが分かると思います。
蟻本などで紹介されているジェネリックダイクストラ法[^5]だとpushの回数がpopと同程度に大きくなるため、ジェネリックダイクストラ法ではpushに時間のかかるradix heapを活かすことが出来ないことが分かります。[^6]

反面、decrease-keyの計算量は$O(1)$であるため、[decrease-keyを使った(オリジナルに近い)ダイクストラ法](https://rsk0315.hatenablog.com/entry/2019/10/29/151823)では効果を発揮することが期待されます。

popの計算量は$O(1)$と破格の計算量ですが、popの回数はpushの回数以下でなければならないため、全体としてはあまり影響はありません。

## ところで

[~~非負整数Xを10進法で表したときの各桁の数字の和をf(X)とします。~~](https://atcoder.jp/contests/iroha2019-day1/tasks/iroha2019_day1_h)

radix heapってヒープなんですかね。
ヒープって親の値が子の値より小さい木/森って印象なんですけど、これが木かと言われると…いや、確かに閉路のない無向グラフではあるんですけど。
百歩譲ってこれが木だとして、親の値が子の値より小さいという条件を満たしているわけでもない(というか親が定義されてない)のでヒープと言い張るには不適切じゃないかなと思います。


## 追記

ポテンシャル関数を
$$
\sum_{i=1}^{33} i \log{i} \times n[i]
$$
と定義すればbsr()が対数時間でも解析できそうな気がします。

そのときの計算量は
- push … $O(k \log{k})$
- decrease-key … $O(\log{k})$
- pop … $O(1)$
になると思います。[^7]


[^1]: bloodという名前で活動する場合もあります
[^2]: 極限に馴染みのない方であれば平均計算量の方が難しいかもしませんし、素数判定など最悪計算量解析が難しいアルゴリズムもあります。
[^3]: 実際には対数時間かかりますが、対数時間で解析したらうまく計算できませんでした。32bit範囲だと5回の加算で済むため、logは定数です。(たすけて)
[^4]: 二分ヒープの場合は償却計算量ではなく平均計算量で保証されることに注意が必要です。
[^5]: 勝手に命名しました。decrease-keyをpop+pushで代用したようなダイクストラ法のことです。(あるいは、BFSのキューをヒープに変えたようなダイクスラ法とも)
[^6]: [twīt](http://zpdic.ziphil.com/dictionary/fez-ja?enableSuggestions=false&ignoreCase=false&mode=both&orderDirection=ascending&orderMode=unicode&page=0&text=tw%C4%ABt&type=prefix)を検索してみると、ジェネリックダイクストラ法で使用してかえって速度を低下させたという声がちらほらあります。
[^7]: 誰か解析してください。俺は力尽きました