乱数でピボットを選択した qsort の最悪計算時間 について

####開発に至る経緯

javaのsortとしてTimSortが採用されたと聞いた。TimSortはマージソートが元に
なっているので速度的に問題があると思った。そこで、クイックソートを元にした
　 Sqsort 　 https://qiita.com/t-kawa/items/4be381efb74d44fc991c を開発した。
初代のSqsortは、TimSortより速く、安定だが、最悪計算時間で負けていた。
そこで「乱数を用いてピボットを選択」を導入して、最悪計算時間の改善を図った。
これにより、Sqsort の最悪計算時間が O(n^2) になる確率は
実用上ほぼ０であると考える。これを解説するのが本稿の目的である。

本稿は証明ではなく、「実用上問題ない」ことの解説である。

####理想的な解決方法（解析的方法）
「乱数でピボットを選択したqsort」は
最悪計算時間や安全性の問題を、確率の問題として扱おうとしている。
ユニーク乱数による長さ n の配列をクイックソート(qsort)するとき、
要素の比較回数の平均は C × n × log2(n) である。
ソート１回における比較回数が C × n × log2(n) の a 倍 以上 になる確率 P ( n, a ) は
どのような数式になるのだろうか？
この問題に解析的な答えがあるらな、P( 1000000, 2.0 ) などを手(て)計算して、
その値が十分に小さいこと示せば、「実用上問題ない」ことを示せたことになる。
しかし、解析的な答えは作成不可能のように思われた。

　　（「ユニーク乱数による配列」とは、ソートキーに同じものが１つもなくて、
　　　　キーの並びを乱数で決定した配列のことである。）

####シミュレーションによる解決方法
qsortの最悪計算時間には、次の3つの場合がある。
　１．n! 通りある要素の並びから選ばれた配列が、本当に本当に運悪く最悪である場合。
　２．"並び" と ピボット選択方法 の相性が悪い場合。（並びが既昇順＋ピボット先頭固定など）
　３．特定のqsortアルゴリズムを攻撃するために作成された「悪意のある配列」である場合。

この内、「１．」は qsort の本質的な特徴であり、改善することはできない。
一方、「ピボットの位置を乱数を用いて決定」することにより
「２．」と「３．」の出現確率は、「１．」と同じ出現確率になると考えられる。

以下では「乱数でピボットを選択」を導入した Sqsort について解説する。
これにより、「１．２．３．」全ての場合で最悪計算時間を実用上 O(n×log(n)) にできると考える。

ソートキーがユニーク乱数で、要素数が n=1000、ソートの繰返し回数が R=1万回 のとき、
1万回中の i 回目のソートにおける"要素の比較関数"の呼出し回数(以後,比較回数)を ci とする。
また、計算時間は ci に比例し、ci の平均は n×log(n)（＝g）に比例するものとする。

以上を前提とするとき、ci/g (i=1,2,…,1万) は i や n が変化してもほぼ一定の値を取る（1.0に近い値）。
しかし、かなりの変動がある。R が大きいほど変動幅も大きくなる。
c1/g, c2/g, …, c1万/g を計算して ci/g の分布を取ると次のようになった。(n=1000, R=10000の場合)

|ci/g区間 | [0.9,1.0) | [1.0,1.1) | [1.1,1.2) | [1.2,1.3) | [1.3,1.4) | [1.4,1.5) | [1.5,1.6) |
|:-------:|:---------:|:---------:|:---------:|:---------:|:---------:|:---------:|:---------:|
|出現頻度 |   434     |   4907    |   3727    |    825    |    92     |    14     |   　  1　 |

これを n=1000、R=1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10 について計算したのが表１である。
ここで 1e4 は 10^4 を表している。 @1.2 は区間 [1.2,1.3) を表している。

|ci/g区間| @0.8 | @0.9 | @1.0 | @1.1 | @1.2 | @1.3 | @1.4 | @1.5 | @1.6 | @1.7 | @1.8 | @1.9 | @2.0 |
|:-------|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|
| R=1e3  |      |   47 |  454 |  397 |   89 |   12 |    1 |      |      |      |      |      |      |
| R=1e4  |      |  434 | 49e2 | 37e2 |  825 |   92 |   14 |    1 |      |      |      |      |      |
| R=1e5  |      | 45e2 | 48e3 | 38e3 | 82e2 | 10e2 |   98 |    3 |    1 |      |      |      |      |
| R=1e6  |    2 | 46e3 | 49e4 | 38e4 | 82e3 | 98e2 |  827 |   44 |    4 |    # |      |      |      |
| R=1e7  |   20 | 45e4 | 49e5 | 37e5 | 82e4 | 99e3 | 86e2 |  555 |   32 |    3 |    # |      |      |
| R=1e8  |  250 | 45e5 | 49e6 | 37e6 | 82e5 | 99e4 | 85e3 | 59e2 |  299 |   13 |    1 |      |      |
| R=1e9  | 23e2 | 45e6 | 49e7 | 37e7 | 82e6 | 99e5 | 85e4 | 58e3 | 33e2 |  158 |   10 |    1 |      |
| R=1e10 | 24e3 | 45e7 | 49e8 | 37e8 | 82e7 | 99e6 | 85e5 | 58e4 | 33e3 | 1513 |   59 |    5 |    # |

　　　　　　表１　n=1000 での比較回数のバラツキ(ci/g)の分布 (空白欄は"0") (行の横合計=R)

表１の各行の右端の "1" に着目する。（ "#" は "1" の出現が予測される位置である。）
ソートをR回繰返したとき、"1" は比較回数(バラツキci/g)が最大となる区間を表している。
R=1e4=1万の場合は、@1.5 の下に "1" があるので
ソートを1万回(R)繰返した場合、比較回数 @1.5 × n × log(n) 回を必要とするソートが
1万回(R)の中に平均で約 "1" 回あることを示している。
逆に言うと、「 @1.5 の変動を目撃したければ約1万回の繰返し(R)を覚悟しろ」と言うことである。


表１の "1"（"#"を含む）に着目し、8個ある "1" を出現させるための繰返し回数(R)
を１行で表すと次のようになる。(n=1000の場合)

|ci/g区間    |   @1.4    |   @1.5    |   @1.6    |   @1.7    |   @1.8    |   @1.9    |   @2.0    |
|:----------:|:---------:|:---------:|:---------:|:---------:|:---------:|:---------:|:---------:|
|繰返し回数R |   1e3     |   1e4     |   1e5     |   1e6     | 1e7～1e8  |   1e9     |   1e10    |

さらに、これを n=1000=1e3 だけでなく n=1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8 
について計算したのが表２である。

|ci/g区間| @1.4 | @1.5 | @1.6 | @1.7 | @1.8 | @1.9 | @2.0 | @2.1 | @2.2 |  …  | @2.6 |
|:-------|:-----|:-----|:-----|:-----|:-----|:-----|:-----|:-----|:-----|:-----|:-----|
| n=1e2  | 1e3  | 1e4  | 1e5  | 1e6  | 1e7  | 1e7  | 1e89 | 1e10 | 1e11 |  …  | 予想 |
| n=1e3  | 1e3  | 1e4  | 1e5  | 1e6  | 1e78 | 1e9  | 1e10 |      |      |  …  |(1e16)|
| n=1e4  | 1e4  | 1e4  | 1e56 | 1e78 | 1e9  |      |      |      |      |  …  |(1e17)|
| n=1e5  | 1e3  | 1e35 | 1e67 | 1e8  |      |      |      |      |      |  …  |(1e18)|
| n=1e6  | 1e3  | 1e4  | 1e57 |      |      |      |  ※  |      |      |  …  |(1e18)|
| n=1e7  | 1e3  | 1e45 | 1e67 |      |      |      |      |      |      |  …  |(1e18)|
| n=1e8  | 1e34 | 1e5  | 1e6  |      |      |      |      |      |      |  …  |(1e18)|

　　　表２　各n,各区間に対応する "1" が出現するはずの繰返し回数(R)　(@2.6の列は予想値)
　　　　　　1e35 は 10^3～10^5 を表している。ただし、1e1* は 10^1* を表している。

この表の値を qsort で計算しようとすると時間がかかり過ぎて計算できない。
「ソートせずに比較回数(ci)だけ」を計算するプログラムsss(約100倍高速)
http://ww51.tiki.ne.jp/~srr-cake/qsort/sqs18j/main_sss.c を作成・使用している。
この表の空白欄は "0" ではない。sssでも時間がかかり過ぎて計算できない部分である。

####計算できない部分（結論）
計算できない部分も予想することはできる。表２の@2.6の列は予想値である。
ただし、表２の "※" などの位置に小さな値(1e5など)が現れると予想はできなくなってしまう。
@2.6列はあくまで予想であって根拠はない。

n=1e3 行 で @2.6 列 の欄の予想回数(R)は、1e16 以上と予想される。
予想回数R(1e16回) のどの 比較回数(ci) も概ね n × log(n) 以上である。
よって、この場合の "総" 比較回数は
　　予想回数R × n × log(n)　≒　1e16 × 1000 × 10　＝　1e20
以上となる。
比較1回あたりの計算時間を10ナノ秒とすると、
　　予想回数実行に必要な計算時間　＝　1e20比較 × 10ナノ秒／比較　≒　3万年
以上となる。

一般に、qsortの平均比較回数(＝C×n×log(n))の C は 1.2～1.3 である。
@2.6 は C のほぼ２倍である。そのため、偶然に平均の２倍の計算時間に
遭遇するには、n=1000 の場合で "３万年" 待つことが "予想" される。
計算機の故障頻度などを考慮すれば、３万年に１度の現象はほぼ出現しない現象といえる。
そして、仮に出現したとしても、平均の "たった" ２倍の計算時間で正常にソートは終了する。
このことは、「１．」の　n! がいかに 大きな数 で、
計算時間 O(n^2) を必要とするケースがいかに 稀 であるかを示している。
以上により、「乱数でピボットを選択したqsort」の最悪計算時間は
O(n^2) でなく、O(n×log(n))として実用上問題ないと思われる。

####悪意のある配列について
「悪意のある配列」に対しても「乱数でピボットを選択したqsort」は有効である。
与えられた配列がたとえ「悪意のある配列」であったとしても、
乱数でピボットを選択することにより、「悪意のある並び順」が無効になるからである。

####乱数を用いた Sqsort と qsort のソート結果の再現性
Sqsort はその処理の先頭で nanoTime() を１度だけ呼び出して乱数を作っている。
そのため、Sqsortの実行ごとに異なる乱数でピボットを選択する。
その結果、２つの同じ配列をソートしたとき、２つの"処理過程"は同じにはならない、
しかし、Sqsortは "安定" なソートなので "処理結果"(２つの配列) は同じになる。（再現性あり）

一方、従来のqsortに「乱数でピボットを選択」を導入した場合、
２つの同じ配列をソートしたとき、"処理結果" は異なる。（再現性なし）
qsortが "非安定" なため、ソートキーが同じ要素の最終位置が変動するからである。
「最悪計算時間」や「悪意のある配列」に関して、このqsortは良好な動作をする。
しかし、この再現性のなさは悪い動作といえる。