バレンタインデーにquadsortがやってきた・・・けど

quadsortとは？

それはバレンタインデーにやってきた

「なんか、qsortよりも良い結果がでたっぽいんだけど！」という話がでてきて、みんな色々検証をしていた。

もとスレッドは「モダンってなんだよ！」とか「それならこっちのアルゴリズムとか見たほうが良いんじゃ？」とか、なかなかエキサイティングしているっぽい。

ソースコードとか技術的な説明はこちら。

性能検証結果

 quadsort: sorted 1000000 elements in 0.074616 seconds. (random order)
    qsort: sorted 1000000 elements in 0.101402 seconds. (random order)
 quadsort: sorted 1000000 elements in 0.000912 seconds. (forward order)
    qsort: sorted 1000000 elements in 0.027091 seconds. (forward order)
 quadsort: sorted 1000000 elements in 0.004904 seconds. (reverse order)
    qsort: sorted 1000000 elements in 0.025820 seconds. (reverse order)
 quadsort: sorted 1000000 elements in 0.018378 seconds. (random tail)
    qsort: sorted 1000000 elements in 0.043725 seconds. (random tail)

quadswap

従来の2要素のスワップでは、２入力に対し、一時変数１つを利用していた。

従来の2要素スワップ

 if ( val[0] > val[1] ) {
  tmp[0] = val[1]
  val[0] = val[0]
  val[1] = tmp[0]
 }

これに対して、quadswapでは２つのステージに分け、４入力に対して、一時変数４つで処理する。

第１ステージ

第１ステージでは、2要素ずつの最小、最大を求める。

tmp[0] = min( val[0], val[1] )、
tmp[1] = max( val[0], val[1] )、
tmp[2] = min( val[2], val[3] )、
tmp[3] = max( val[2], val[3] )、

第１ステージ

if (val[0] > val[1])
{
    tmp[0] = val[1];
    tmp[1] = val[0];
}
else
{
    tmp[0] = val[0];
    tmp[1] = val[1];
}

if (val[2] > val[3])
{
    tmp[2] = val[3];
    tmp[3] = val[2];
}
else
{
    tmp[2] = val[2];
    tmp[3] = val[3];
}

第２ステージ（前半）

第２ステージは更に２つのサブステージに分けて考えられる。

まず、順番が一意に確定できるパターンがある。

tmp[1](= max( val[0], val[1] ) が、tmp2 以下であるならば、必然的に４要素の順番は以下で確定する。

tmp[0] = min ( val[0],val[1] )
tmp[1] = max ( val[0],val[1] )
tmp[2] = min ( val[2],val[3] )
tmp[3] = max ( val[2],val[3] )

一方、tmp[0](= min( val[0], val[1] ) が、tmp3以上であるならば、必然的に４要素の順番は以下で確定する。

tmp[2] = min ( val[2],val[3] ) ※ tmp[3]より絶対小さい
tmp[3] = max ( val[2],val[3] )
tmp[0] = min ( val[0],val[1] ) ※ tmp[1]より絶対小さい
tmp[1] = max ( val[0],val[1] )

第２ステージ前半

if (tmp[1] <= tmp[2])
{
    val[0] = tmp[0];
    val[1] = tmp[1];
    val[2] = tmp[2];
    val[3] = tmp[3];
}
else if (tmp[0] > tmp[3])
{
    val[0] = tmp[2];
    val[1] = tmp[3];
    val[2] = tmp[0];
    val[3] = tmp[1];
}
else

第２ステージ（後半）

次に、順番に決めるパターンを解く。

tmp[0](= min( val[0], val[1] ) が、tmp2 以下であるならば、最初の２要素は以下で定まる。

val[0] = tmp[0] = min ( val[0], val[1] )
val[1] = tmp[2] = min ( val[2], val[2] )

そうでなければ、順番は入れ替わる。

val[0] = tmp[2] = min ( val[2], val[2] )
val[1] = tmp[0] = min ( val[0], val[1] )

同様に、tmp[1](= max( val[0], val[1] ) が、tmp3 以下であるならば、最初の２要素は以下で定まる。

val[2] = tmp[1] = max ( val[0], val[1] )
val[3] = tmp[3] = max ( val[2], val[2] )

そうでなければ、順番は入れ替わる

val[2] = tmp[3] = max ( val[2], val[2] )
val[3] = tmp[1] = max ( val[0], val[1] )

第２ステージ後半

{
   if (tmp[0] <= tmp[2])
   {
       val[0] = tmp[0];
       val[1] = tmp[2];
   }
   else
   {
       val[0] = tmp[2];
       val[1] = tmp[0];
   }

   if (tmp[1] <= tmp[3])
   {
       val[2] = tmp[1];
       val[3] = tmp[3];
   }
   else
   {
       val[2] = tmp[3];
       val[3] = tmp[1];
   }
}

quadsort

基本的な考え方は、マージソートと同じである。正直、この部分に目新しさがないかも。。。

例えば３パスで処理する話が記載されています。うーん、これは普通のマージソートなのでは……

In the first row quadswap has been used to create 16 blocks of 4 sorted elements each. In the second row quadsort has been used to create 4 blocks of 16 sorted elements each. In the last row we're left with 1 block of 64 sorted elements.
（邦訳）最初のrow quadswapでは、それぞれ整列された4要素を持つ、16のブロックが生成される。
２番目のrow quadsortでは、それぞれ整列された16要素を持つ、4つのブロックが生成される。
最後のrowでは、整列された64要素を持つ１つのブロックが完成する。

追記(例えばどうなる？）

A<B<C<D とした場合、組み合わせは、24通り。

A-B-C-D, A-B-D-C 
A-C-B-D, A-C-D-B
A-D-B-C, A-D-C-B

B-A-C-D, B-A-D-C
B-C-A-D, B-C-D-A
B-D-A-C, B-D-C-A

C-A-B-D, C-A-D-B
C-B-A-D, C-B-D-A
C-D-A-B, C-D-B-A

D-A-B-C, D-A-C-B
D-B-A-C, D-B-C-A
D-C-A-B, D-C-B-A

これらが、第１ステージを経ると、以下のようにパターンが重複するようになる。

A-B-C-D, A-B-C-D ,
A-C-B-D, A-C-B-D
A-D-B-C, A-D-B-C

A-B-C-D, A-B-C-D
B-C-A-D, B-C-A-D
B-D-A-C, B-D-A-C

A-C-B-D, A-C-B-D
B-C-A-D, B-C-A-D
C-D-A-B, C-D-A-B

A-D-B-C, A-D-B-C
B-D-A-C, B-D-A-C
C-D-A-B, C-D-A-B

重複を除外すると、実際のパターン数は6パターン。

A-B-C-D(4)
A-C-B-D(4)
A-D-B-C(4)
B-C-A-D(4)
B-D-A-C(4)
C-D-A-B(4)

第２ステージ前半部では、8パターンは即時に決定できる。

A-B-C-D(4)　⇒　(tmp[1] <= tmp[2])  なので、 {A,B,C,D}
C-D-A-B(4)　⇒　(tmp[0] > tmp[3])　なので、 {A,B,C,D}　※入れ替える

第２ステージ後半部では、のこり16パターンをmin部とmax部とで分けて比較する。

A-C-B-D(4)　⇒ {A,*,B,*} + {*,C,*,D}
A-D-B-C(4)　⇒ {A,*,B,*} + {*,D,*,C}
B-C-A-D(4)　⇒ {B,*,A,*} + {*,C,*,D}
B-D-A-C(4)　⇒ {B,*,A,*} + {*,D,*,C}

{A,*,B,*) ならば、(tmp[0] <= tmp[2])であり、 {A,B,-,-}
{B,*,A,*) ならば、!(tmp[0] <= tmp[2])であり、 {A,B,-,-} ※入れ替える

{*,C,*,D) ならば、(tmp[1] <= tmp[3])であり、 {-,-,C,D}
{*,D,*,C) ならば、!(tmp[1] <= tmp[3])であり、 {-,-,C,D} ※入れ替える

よって、

A-C-B-D(4)　⇒ {A,*,B,*} + {*,C,*,D} = {A,B,C,D}
A-D-B-C(4)　⇒ {A,*,B,*} + {*,D,*,C} = {A,B,C,D}
B-C-A-D(4)　⇒ {B,*,A,*} + {*,C,*,D} = {A,B,C,D}
B-D-A-C(4)　⇒ {B,*,A,*} + {*,D,*,C} = {A,B,C,D}

実はquadsortは並列化しやすいハードウェアに最適化できる可能性が高いのでは？

ここで、第１＋第２ステージに注目する。このうち、val[x]はメモリへのload/store、tmp[x]はレジスタ操作になる。

// 第１ステージ
tmp[0] <- val[0] or val[0]
tmp[1] <- val[0] or val[1]
tmp[2] <- val[2] or val[3]
tmp[3] <- val[2] or val[3]
// 第２ステージ
// 前半
val[0] <- tmp[0] or tmp[2]
val[1] <- tmp[2] ot tmp[0]
// 後半
val[2] <- tmp[1] or tmp[3]
val[3] <- tmp[3] or tmp[2]

つまり、４要素を対象にして、LoadとStoreをまとめてやれるようになるアルゴリズムと言えるわけです。

例１）FPGAやらreconfigurable processorあたりが処理しやすいといえば、処理しやすいんじゃないのかなーという気がしていますね。

例２）１レジスタにvalの配列全部ぶち込んで、レジスタ内でまとめて入れ替えられるのではないかという気がしてくるわけです……

むしろ、RISC-Vなんかで４要素ソート専用の命令を作ってしまえば、圧倒的に高速になる気がしますね……

あと、マージソートも同じですが、マルチコアでぶん回すときに、圧倒的に処理しやすいですね。もしかすると、現状のイントロソートの代わりに全部quadsort構成のほうが良かったりする可能性も否定しきれないですね（検証が必要）。

結論

• 検証方法などはまだまだ詰める必要はありそう。
• 「４要素まとめて処理する」というアイデアは、実は悪くないかもしれない。

---

余談） qsortってquick sortじゃないよね？

さて、こういう議論で出てくるqsort()はかなりめんどくさい立ち位置にいる。

C99( Draft N1256)

The qsort function sorts an array of nmemb objects, the initial element of which is pointed to by base. The size of each object is specified by size.
(邦訳) qsort関数は、baseで指示された初期要素に対して、、nmemb objectsの配列をソートします。それぞれの要素の大きさはsizeで定義されます。

C11( これは Draftですが N1570 で確認 )

The qsort function sorts an array of nmemb objects, the initial element of which is pointed to by base. The size of each object is specified by size.
(邦訳) qsort関数は、baseで指示された初期要素に対して、、nmemb objectsの配列をソートします。それぞれの要素の大きさはsizeで定義されます。

Linux Programmer's Manual

qsort() 関数は、 nmemb 個の大きさ size の要素をもつ配列を並べ変える。 base 引き数は配列の先頭へのポインターである。

glibc実装

うん、イントロソートでございますなあ…

https://ja.wikipedia.org/wiki/%E3%82%A4%E3%83%B3%E3%83%88%E3%83%AD%E3%82%BD%E3%83%BC%E3%83%88

しかもスタックとかを自前でやってるから、結構面白い構造になっておりますね！！

で、ここで終わってたら「qsortはquicksortじゃないんだよ、」って断言できたんですが……

(BSD) Library Functions Manual

The qsort() function is a modified partition-exchange sort, or quicksort.
（邦訳）qsort()関数は、改良された部分交換ソート、quicksortです。

oh....