javaで記述した　TimSortよりも 高速 かつ 安定・安全 な qsort (Sqsort)

最近のソートは、クイックソートに代わってマージソートが主流のようです。
最大の理由は「安定」性にあると思われます。しかし、マージソートは
　（１）要素の移動回数が多い
　（２）同じキー値が複数現れる場合でも速度向上が十分には望めない
などの欠点があります。
#####このことから「安定」な qsort (Sqsort) (比較ソート) を作りました。
#####最大の特徴は、キー値に同じ値が現れる場合に極めて高速なことです。
Sqsortはjavaで記述してあり、同じくjavaで記述したTimSortより高速です。
大雑把にいうと、Sqsortは補助配列(bas2)と対象配列(base)で
以下のような移動を繰り返します。

```
   bas2:...............  →  bas2:5555......77787  →  bas2:..........77787
   base:357358257257237  Ａ base:332223.........  Ｂ  base:3322235555.....
```

プログラム(約150行 classファイル 3164byte) は最後に記載します。
TimSort(約900行 classファイル 8801byte) よりかなりコンパクトです。
処理時間を測定したものを以下に記載します。

```
    ソート名   キーの種類  　要素数　 　繰返し回数　処理時間(秒)
     TimSort    d=-3      e=10000     R=2000     C=8.21 
     Sqsort     d=-3      e=10000     R=2000     C=5.65 
                                                            
     TimSort    d=100     e=10000     R=2000     C=6.69 
     Sqsort     d=100     e=10000     R=2000     C=2.81 
                                                            
     TimSort    d=2       e=10000     R=6000     C=7.41 
     Sqsort     d=2       e=10000     R=6000     C=2.09 
 
      d=-3:キーの重複なし　d=100:キーが100種類　d=2:キーが男女など         
      繰返し回数：処理時間の測定誤差を減らすためにソートを繰り返した回数   
```

すべての場合で、Sqsort は TimSort より高速です。
より詳細な測定結果・ベンチマークプログラム・TimSort/Sqsortのソースコードなどは
　　　　http://ww51.tiki.ne.jp/~srr-cake/qsort/sqs18j/ 　　にあります。

「安定」とは、キー値が同じ時、ソート前の順序をソート後も保持することです。
「安全」とは、「悪意のある配列により、異常に長い処理時間を発生させることができない」
　　　　　　　こととします。

Sqsort では「安全」を実現するために、ピボットの候補を
乱数(System.nanoTime())により決定しています。
これにより、最悪計算時間O(n^2)は変化しませんが、
元々その発生確率は非常に低いものです。
これまで長い間、C言語のnewlibライブラリのqsort()の実装において
最悪計算時間への配慮がなかったことが、その低さを示していると思います。
今回、最悪計算時間の発生確率は変わりませんが、
「悪意のある配列」を作ることは不可能になっていると思います。


以上３種類のソートを比較すると次の表のようになります。

| 　　　　 | TimSort | Sqsort　 | qsort 　 |
|:--------|:-------:|:--------:|:--------:|
| 速さ　　 |   △    |    ○     |    ◎    |
| 作業領域 |   △    |    △    |    ◎    |
| 安定性　 |   ○     |    ○     |    ×     |
| 安全性　 |   ○     |    ○     |    ×     |

また、TimSortでは配列に内在するソート済み部分列を有効利用しますが、
Sqsortにその機能はありません。
代わりに、Sqsortは同じキー値が複数回現れる場合に高速です。
TimSortにこの機能はありません。
すなわち、県名別にソートしたり、男女を分けるためにソートする場合に
極めて高速です。
　

Sqsortでは、最悪計算時間や安全性の問題を、
乱数による確率の問題として扱おうとしています。
これについては、以下の記事で解説します。
　　「乱数でピボットを選択した qsort の最悪計算時間 について」
　　　 https://qiita.com/t-kawa/items/3dfba10c407d1aede95d

　

以下が Sqsort のソースコードです。

```
/***************************************************/
/*         Sqsort   (stable qsort)                 */
/*                                                 */
/* by 河村 知行 (kawamura tomoyuki)    2021.2.27   */
/*             t-kawa@crux.ocn.ne.jp               */
/***************************************************/

// Sqsort は安定な比較ソート。TimSortより高速で,作業領域は2n(TimSortはn),安全性はほぼ同じ
/*
対象配列(base)と補助配列(bas2)(baseと同じ大きさ)を用いて安定なソートを行う。
以下のように Ａ,Ｂ 2段階で要素を移動する。ここで「5」はピボット(分割要素)を表している。

   bas2:...............  →  bas2:5555......77787  →  bas2:..........77787
   base:357358257257237  Ａ  base:332223.........  Ｂ  base:3322235555.....

移動Ａで、ピボットと同じ要素をbas2の左端に集め、5より小さい要素をbaseの左端に集め、
5より大きい要素をbas2の右端に集める。このとき77787は初期順序(元順)の逆順になっている。
移動Ｂの完了時点で5555はすでにソート完成状態の位置に収まっている。
Ｂ直後、5555以外の要素は、base上では元順に、bas2上では逆順に並んでいる。
次に332223の{所属フラグ(base/bas2)・先頭位置・末尾位置}をスタックに保存する。
これ以降、77787に対しbaseとbas2の意味を入れ替えて同じ処理を繰り返す。
部分配列が長さ15以下になれば挿入ソートを実施して部分配列のソートを完了し、次に
スタックから{フラグ(base/bas2)・先頭位置・末尾位置}を取り出して処理を続ける。
スタックが空になればソートを終了する。
*/
#define C(x,y)  {           x=(y);       }  /* ソート対象配列の要素のコピーに用いる */
#define S(x,y)  {T tmp=(x); x=(y); y=tmp;}  /* ソート対象配列の要素のスワップに用いる */

#define CMP(x,y)       cmp.compare((x),(y))
#define med3p(x,y,z)  (CMP(x,y)<=0 ? (CMP(y,z)<=0 ? y : (CMP(x,z)<=0 ? z : x)) : \
                                     (CMP(y,z)>=0 ? y : (CMP(x,z)<=0 ? x : z)) )

#define I(x) {x++;} /* 配列base,bas2の添字を増やす */
#define D(x) {x--;} /* 配列base,bas2の添字を減らす */

import java.util.Comparator;

class Sqsort {

static int L,R;        // L,Rは、baseまたはbas2上のソート対象となる領域[L～R]の左端と右端
static boolean flagB;  // [L～R]領域がbase上にあるときtrue、bas2上にあるときfalse
// popメソッド内で L,R,flagB へ代入できるように、sortメソッドのローカル変数からここへ移動した

//#define  base0  base  /*これを有効にすると、配列base0を複写しない。 32bitJVMでは有効の方が速い*/

@SuppressWarnings({"unchecked"})
public static <T> void sort( T[] base0, Comparator<? super T> cmp ) {//base0:ソート対象である配列
                                                                     //  cmp:要素の比較をする関数
 int l,r,m; //ptrB[L～R]間の要素をptrB[L～l],ptrR[r～R],ptrR[L～m]へ振り分ける。l領域 r領域 m領域
 T[] ptrB, ptrR; //flagBのとき ptrB==base ptrR==bas2、!flagBのとき ptrB==bas2 ptrR==base。
 int n,t,p;      // nは、現在ソート中の[L～R]領域の要素数。tは作業用。pはL～R間を移動する添字
 final long rndm = (System.nanoTime() & 0xFFFFFFFFFFFFL); //nanoTime()を乱数発生器として使用

 class StkElt {             //スタックをローカルクラスで実装
   boolean sB; int sL,sR;
   void push(boolean b, int l, int r) {      sB=b;     sL=l;     sR=r;}
   void pop ()                        {flagB=sB;     L=sL;     R=sR;  }
 }
 final int STK_LEN=32;      // base0.lengthの最大値は2^31-1なので３２で十分
 StkElt stk[] = new StkElt[STK_LEN];
 for (int i=0; i<STK_LEN; i++) stk[i] = new StkElt();
 int top = 0;               // スタックポインタの初期値
   
 int Z = base0.length;      // このsortメソッドに渡されたソート対象配列の要素数
#ifndef base0               /* 64bit JVM では、コピーした配列をソートした方が速いので次行を使用 */
 T[] base = (T[]) new Object[Z];   System.arraycopy(base0, 0, base, 0, Z);
#endif
 T[] bas2 = (T[]) new Object[Z];  //baseと対になって使用。ただし、m領域のコピー先は常にbase

 flagB=true; L=0; R=Z-1;

 //ここまでは、sortの１回の呼び出して１度だけ実行される

 for(;;) {
   for(;;) {//LOOP:
     if (R<L)  {                                 break;/*goto nxt;*/} //要素なし
     if (L==R) {{if (!flagB) C(base[L],bas2[L])} break;/*goto nxt;*/} //要素数１
     
     n = (R - L + 1);      /*要素数 n>=2 */

     if (n <= 15) { //挿入ソート   結果はbase上に作る   n==2でも動く  「15」は実験で得られた値
       if (flagB) { //ソート領域base[L～R]は 元順 になっている
         if (CMP(base[L],base[L+1])>0) S(base[L],base[L+1])
         for (r = L+2; r <= R; r++) {
           if (CMP(base[r-1],base[r])<=0) continue;
           T tmp=base[r]; C(base[r],base[r-1])
           for (l=r-2; L<=l && CMP(base[l],tmp)>0; l--) C(base[l+1],base[l])
           base[l+1]=tmp;
         }
       }else{       //ソート領域bas2[L～R]は 逆順 になっている
         if (CMP(bas2[R-1],bas2[R])<0) {C(base[L],bas2[R-1]) C(base[L+1],bas2[R])  }
         else                          {C(base[L],bas2[R])   C(base[L+1],bas2[R-1])}
         for (r = R-2; L <= r; r--) {
           t=L+(R-r);
           if (CMP(base[t-1],bas2[r])<=0) {C(base[t],bas2[r]) continue;}
           T tmp=bas2[r]; C(base[t],base[t-1])
           for (l=t-2; L<=l && CMP(base[l],tmp)>0; l--) C(base[l+1],base[l])
           base[l+1]=tmp;
         }
       }
       break; /*goto nxt;*/
     }

     if (flagB) {ptrB=base; ptrR=bas2;} else {ptrB=bas2; ptrR=base;}

     {//３点処理（３要素からピボットを選択　高速化のため９点処理や２７点処理の版もある）
       int rnd = (int)(rndm % ((n >> 1) - 2)); //nanoTime()を乱数発生器として使用
       T ll=ptrB[L+rnd], rr=ptrB[R-rnd],       //悪意ある攻撃を無効化するため、要素を乱数で選択
         mm=ptrB[L + (n >> 1)]; //中央の要素   //３点の内、定点が２点では悪意ある配列が作れてしまう

       T mmm=med3p(ll,mm,rr);                  //３点を比較して、mmmにピボットをセットする

     // l:ピボットより小さい要素(小要素)を次に置く所  r:大要素を次に置く所  m:同要素を次に置く所
     // ptrR[L～m]:左端からm領域を作る  ptrR[r～R]:の右端からr領域を作る(bas2上のr領域は逆順となる)
     // ptrB[L～l]:左端からl領域を作る  ptrB[L～R]:この全要素を左端からl領域,r領域,m領域へ振分ける
       l=m=p=L; r=R; 
       do {T pp=ptrB[p]; if ((t=CMP(pp,mmm)) < 0) {C(ptrB[l],pp) I(l)}
                         else if          (t > 0) {C(ptrR[r],pp) D(r)}
                         else                     {C(ptrR[m],pp) I(m)} I(p)} while(p<=R);
     }

     assert L<m;                    //m領域には少なくとも要素１個が配置される
     assert (l-L)+(m-L)+(R-r)==n;   //l領域,r領域,m領域の要素数の合計は ptrB[L～R]の要素数に等しい

     if (m-L == 1) {        //m領域をbase上のあるべき位置に移動する
       C(base[l],ptrR[L])   //要素数が１のときは、単純に代入（これは高速化のため。省略可能）
     }else{
       if (flagB) {
         System.arraycopy(bas2, L, base, l, m-L);                 //bas2上のm領域は常に元順
       }else{
         if (m<=l) {
           p=l; do {D(m) C(base[p],base[m]) I(p)} while (L<m);    //base上のm領域は常に逆順
         }else{
           for (t=m-1, p=l; p<t;) {S(base[p],base[t]) I(p) D(t)}  //m領域とその移動先が重なるとき
           for (p=m, m=l; L<m;) {D(m) C(base[p],base[m]) I(p)}
         }
       }
     }

     if (l-L < R-r) {stk[top++].push(!flagB,r+1,R); /*flagB=flagB;*/ R=l-1;} /*左を先にソート*/
     else           {stk[top++].push( flagB,L,l-1); flagB=(!flagB);  L=r+1;} /*右を先にソート*/
   }//goto LOOP;

   //nxt:                   /*スタックが空になったとき終了する*/
   if (top == 0) {
#ifndef base0               /* 64bit JVM では、コピーした配列をソートした方が速いので次行を使用 */
     System.arraycopy(base, 0, base0, 0, Z);
#endif
     return;
   }
   stk[--top].pop();   /*スタックに積んである値を取り出して、処理を繰り返す*/
 }
} // sort メソッドの終了
} // Sqsort クラスの終了
```