####◎旧来の qsort の不便さ
C言語の著名なライブラリに、newlib と glibc があります。
その中の qsort 関数は２０年以上使用されていると思うのですが、
いくつか不便な点があります。
１．newlib
newlib の qsort は、直接ソートのみで間接ソートへの自動切り替えがありません。
そのために、元配列の要素が大きいとき、低速となります。
（間接ソート：各要素へのポインタ配列を作り、これをソートする。最後に元配列の要素を移動する。）
２．glibc
glibc の qsortは、実体はマージソートです。
そのために、キーが（社員番号のような）重複がないときは高速ですが、
キーに（男・女のような）重複があるときは、低速です。
また、マージソートなので安定ソートです。
しかし、作業領域の確保(malloc)に失敗したときは、
普通のクイックソートを行います。そのため、「安定」を保証できません。
今回公開する ss18 は、高速な安定ソートです。
「不安定で良いから、とにかく高速なソート」を望むなら、
ss15 (http://ww51.tiki.ne.jp/~srr-cake/qsort/qs15/) をお薦めします。
newlib・glibc 意外はあまり調べていません。より高速な比較ソートがあれば教えて下さい。
####◎作業領域について
ss18 は作業領域(オーダー$O(n)$)(=補助配列)を必要とします。
glibc の qsort は安定性が保証されないのに同じ作業領域$O(n)$を使用します。
筆者としては、主メモリが安価な時代ではこれでよいと考えています。
####◎計算時間について
qsortでは、最悪計算時間が$O(n^2)$となります。
しかし、ss18 はピボットの選定を丁寧に行っているので、
（偶然または運悪く発生する）相性の悪いソート対象配列はほぼないと
考えています。
ただし、悪意を持って作られたソート対象配列では計算時間が$O(n^2)$
となり得ます。この対策としては、「ピボットの選定をより丁寧に行う」
ことがあります。（ただ、現在のバージョンでは実装していません）
具体的には、「プロセス内での最初の ss18 の呼び出し時に、time()関数などで
再現性のない乱数を発生させて、この値でピボットを変動させる」などがあります。
この方法は安定ソートである ssort で特に有効です。ピボットが変動しても、
ソート後に対象配列は常に同じになるからです。qsort では異なる可能性があります。
#クイックソートを改造して安定な比較ソート stable sort (ssort)(ss18) を作成しました
以前に公開した ss14 とは全く別のアルゴリズムです。
よりコンパクトで、より高速です。大雑把にいうと、
ss18は補助配列(ptr2)と対象配列(ptr)で以下のような移動を繰り返します。
ptr2:............... → ptr2:5555......77787 → ptr2:..........77787
ptr :357358257257237 Ａ ptr:332223......... Ｂ ptr:3322235555.....
アルゴリズム解説のための簡易版(要素サイズを８バイトに固定)(ss18c4d)の
処理時間を測定したものを以下に記載します。
簡易版のプログラムは最後に記載します。
キーの種類 要素数 サイズ 繰返し回数 処理時間(秒)
qs_glibc d=-3 e=10000 s=8 R2463 T=5.74
ss18c4d d=-3 e=10000 s=8 R2463 T=4.73
qs_glibc d=100 e=10000 s=8 R7035 T=15.47
ss18c4d d=100 e=10000 s=8 R7035 T=6.25
qs_glibc d=2 e=10000 s=8 R24630 T=35.19
ss18c4d d=2 e=10000 s=8 R24630 T=2.35
d=-3:キーの重複なし d=100:キーが100種類 d=2:キーが男女など
繰返し回数：処理時間の測定誤差を減らすためにソートを繰り返した回数
qs_glibc は、Cのライブラリ "glibc" のqsort関数を性能測定用に改造したものです。
qs_glibc は、はmalloc()に失敗しない限りマージソートとして動作します。
マージソートは安定なソート(https://ja.wikipedia.org/wiki/ソート) です。
すべての場合で、ss18c4d は qs_glibc より高速です。
正規版(ss18e8b)は、 http://ww51.tiki.ne.jp/~srr-cake/qsort/ss18/ にあります。
readme2.txt に正規版のベンチマークテストの方法を載せました。
正規版では、qsnewlibの処理時間も測定しています。
( http://ww51.tiki.ne.jp/~srr-cake/qsort/ss18/bench-sample.txt を参照 )
qsnewlib は、Cのライブラリ "newlib" のqsort関数です。安定でないソートです。
なので、qsnewlib が最速になって当然なのですが、間接ソートへの
自動切り替えがないため要素サイズが大きい場合に qsnewlib は遅くなります。
qsnewlib(非安定) qs_glibc(安定) ss18e8b(安定) を比べると
64bitシステムでは、概ね ss18e8b が最速です。
/********** ssort (stable sort) ss18 *******/
/*
ss18 の簡易版 (ss18は安定な比較ソート。マージソートやss14より高速。)
対象配列(ptr)と補助配列(ptr2)(ptrと同じ大きさ)を用いて安定なソートを行う。
以下のようにＡＢ2段階で要素を移動する。ここで「5」はピボット(分割要素)を表している。
ptr2:............... → ptr2:5555......77787 → ptr2:..........77787
ptr :357358257257237 Ａ ptr:332223......... Ｂ ptr:3322235555.....
移動Ａで、ピボットをptr2の左端に集め、5より小さい要素をptrの左端に集め、
5より大きい要素をptr2の右端に集める。このとき77787は初期順序の逆順になっている。
移動Ｂのとき、もしも5555が逆順ならば、Ｂ直後の5555が正順になるように移動する。
この時点で5555はすでにソート完成状態の位置に収まっている。
次に332223の{先頭位置・末尾位置・正順逆順}をスタックに保存する。
これ以降、77787に対して同じ処理を繰り返す。部分配列が長さ2以下になれば部分ソートを終了して
スタックから{先頭位置・末尾位置・正順逆順}を取り出して処理を続ける。
スタックが空になれば処理を終了する。
（以上は簡易版なので、ピボットの収納位置が正規版とは異なっている）
*/
#include <stdio.h>
//typedef long unsigned int size_t;
void *malloc(size_t); void free(void*); void exit(int);
static void assert(int x, char *s) //必要な場合は、下の行のコメントを取り去る
{ /* if (x) return; fprintf(stderr, "++++ %s ++++\n", s); exit(1); */ }
#define F(x) (*((long long int *)(x)))
#define COPY(x,y) { F(x)=F(y); } /*配列の要素のコピー処理*/
#define SWAP(x,y) {long long int v=F(x); F(x)=F(y); F(y)=v;} /*配列の要素のスワップ処理*/
//上の3行は要素の移動に関するもの 別なものに書き換え可能です
//ここから下は、ss18によるソートプログラムです
typedef struct { char *LLss, *RRss; int Rv; } stack_node; /*L,Rを積むスタックの構造体*/
#define PUSH(llss,rrss,rv) {top->LLss = (llss); top->RRss = (rrss); top->Rv = (rv); ++top;} /*積む*/
#define POP(llss,rrss) {--top; llss = top->LLss; rrss = top->RRss; rev = top->Rv;} /*戻す*/
int ssort( void *base, size_t nel, size_t size, int (*cmp)(void *a, void *b) )
/* base : ソートしようとする配列へのポインタ */
/* nel : 配列baseの要素数 */
/* size : 配列baseの要素の大きさ（バイト単位） */
/* cmp : 要素の大小比較をする関数へのポインタ */
{
char *l,*l0,*r,*r0,*L,*R; // l0,r0,L,Rは配列の左右の起点。 L～R間の要素をl,r,mへ振り分ける
char *ptr,*ptrE,*ptr2; // ptr,ptr2 は L～R,l0～r0 配列の先頭 送受は毎回入れ替わる
char *m,*m0,*p; // pはL～R間を移動 *mにピボットを貯める
int t, rev; // tは作業用 rev==0のときL～R間は正順
long long int Z; // Zは配列間の距離(正負ありうる)
stack_node stack[32], *top = stack; /*現在のところは３２で十分*/
assert(sizeof( int) == 4, "sizeof( int) == 4");
assert(sizeof(long long int) == 8, "sizeof(long long int) == 8");
assert(size == 8 , "簡易版のなので、sizeは8に固定とする" );
//簡易版なので、baseのアライメントはsizeof(char*)になっているとする
ptr2 = malloc( nel * size ); //ptr2はソートのための補助配列
if (ptr2 == NULL) {assert(0,"mallocに失敗した"); return 1;}
ptr = base;
Z = ptr2 - ptr; //Zは２個の配列・補助配列間の距離、正負両方あり
ptrE = ptr + nel * size; // ptrE はL～Rが対象配列・補助配列いずれにあるか、判定に使用する
L=ptr; R=&ptr[size*(nel-1)]; rev=0;
goto LOOP;
//ここまでは、ssortの１回の呼び出して１度だけ実行される
nxt:
if (stack == top) {free(ptr2); return 0;} /*スタックが空になったとき終了する*/
POP(L,R); /*スタックに積んである値を取り出す*/
LOOP:
if (R<L) { /*要素なし*/
goto nxt;
}
if (L==R) { /*要素数１*/
if (ptr<=L && L<ptrE) { } // L<ptr2 の判定が普通だが、ptr2<base のシステムも存在する
else {COPY(L-Z,L)}
goto nxt;
}
if (L+size==R) { /*要素数２*/
if ((t=cmp(L,R))<0) if (ptr<=L && L<ptrE) { }
else {COPY(L-Z,L) COPY(R-Z,R)}
else if (t>0) if (ptr<=L && L<ptrE) {SWAP(L,R) }
else {COPY(L-Z,R) COPY(R-Z,L)}
else if (ptr<=L && L<ptrE) if (rev==0) { }
else {SWAP(L,R) }
else if (rev==0) {COPY(L-Z,L) COPY(R-Z,R)}
else {COPY(L-Z,R) COPY(R-Z,L)}
goto nxt;
}
if (ptr<=L && L<ptrE) {l=l0=L ; r=r0=R+Z; m=m0=L+Z;} // l : ピボットより小さい要素を次に置く場所
else {l=l0=L-Z; r=r0=R-Z; m=m0=L; } // r : ピボットより大きい要素を次に置く場所
// m : ピボットと同じ 要素を次に置く場所
COPY(m,L) m+=size; //簡易版のなので、単純に先頭の要素をピボットとする
for (p=L+size; p<=R; p+=size) {
if ((t=cmp(p,m0))<0) {COPY(l,p) l+=size;}
else if (t>0) {COPY(r,p) r-=size;}
else {COPY(m,p) m+=size;}
}
assert((l-l0)+(m-m0)+(r0-r)==(R-L+size),"(l-l0)+(m-m0)+(r0-r)がおかしい");
if (rev==0) {
p=l; do {COPY(p,m0) m0+=size; p+=size;} while (m0<m); //m0～mには少なくとも１個ある
if (l-l0 < r0-r) {PUSH(r+size,r0,1) L=l0; R=l-size; rev=0;} /*左から先にソートする*/
else {PUSH(l0,l-size,0) L=r+size; R=r0; rev=1;} /*右から先にソートする*/
}else{
p=l; do {m-=size; COPY(p,m) p+=size;} while (m0<m); //m0～mには少なくとも１個ある
if (l-l0 < r0-r) {PUSH(r+size,r0,0) L=l0; R=l-size; rev=1;} /*左から先にソートする*/
else {PUSH(l0,l-size,1) L=r+size; R=r0; rev=0;} /*右から先にソートする*/
}
goto LOOP;
}