####◎旧来の 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;
}