[C]クイックソートを実装してみた

# はじめに

　競技プログラミングなどで、ソートを使う場面が多いですが、自分でクイックソートを実装してみたくなったので実装してみました。<br>　まず簡単にクイックソートを実装し、それを改良してみます。

# クイックソートの仕組み

クイックソートは以下の操作を再帰的に行うことでソートを行っています。（以下左から右に昇順になるようソートした場合）
- ソートする範囲の中から一つ値を選ぶ（選ばれた値はpivotと呼ばれます）
- pivotの左側はすべてpivot以下の値、右側はpivot以上の値になるようにデータを交換していく
- pivotより右の範囲と左の範囲でそれぞれソートをする

# 実装
以上の仕組みを実装したものが以下になります。
```c
#define BOOL char

/* ソートは昇順にする return a<b にすると降順にソートができる */
BOOL cmp(int a,int b)
{
    return a>b;
}

/* 値を交換する */
void swap(int* i,int* j)
{
    int tmp = *i;
    *i = *j;
    *j = tmp;
    return;
}

/* 半開区間 [left,right) をソートする */
void quick_sort(int* left, int* right,BOOL (*cmp)(int a,int b))
{
    /* ソート範囲に値が１つしかない場合はreturn */
    if(right-left <= 1)return;
    
    /* 配列の真ん中の値をpivotとする */
    int *pivot = (left+((right-left)/2)); 
    
    /* 配列の両端から線形探査し、値を交換していく */
    int *l = left , *r = right-1;
    while(l<r)
    {
        while(cmp(*pivot,*l) && l < pivot)l++;
        while(cmp(*r,*pivot) && pivot < r)r--;
        swap(l,r);
        if(l ==  pivot) pivot = r;
        else if(r == pivot) pivot = l;
        if(l != pivot)l++;
        if(r != pivot)r--;
    }

    /* pivotを除き、残りの範囲をそれぞれソート */
    quick_sort(left,pivot,cmp);
    quick_sort(pivot+1,right,cmp);
}


```

# 改良する

　こちらのソートは基本的には$O(n$$log$$n)$で動作しますが、pivotをソート範囲の中央の値にしているため、仮に選んだpivotの値が毎回そのソート範囲の最大値ないしは最小値になってしまうと計算量が$O(n^2)$まで膨らんでしまいます。<br>　理想としてはpivotが毎回ソート範囲の中央値であることが理想になるので、このpivotを選ぶ方法を改良していきます。

# 改良案

　今回はpivotを選ぶやり方として、範囲の値を三分木の葉の要素とみなし、それぞれ3つある子の要素のうちの中央値を親ノードの値としていき、その結果の根の値をpivotとするようにします。実装は以下です。
 ```c
 #include <stdlib.h>
 
/* [left,right)の中の[l,r)について再帰的に計算する */
int* _cal_pivot(int* left, int* right,int* l, int* r,int size)
{
    /* [l,r)が[left,right)の範囲外のときNULL(ポインタ)を返す */
    if(right <= l)return NULL;

    /* 範囲に要素が1個だけのときそのポインタを返す */
    if(size==1)return l;

    /* 三分木の子の要素の中央値を再帰的に求める */
    int *m1 = _cal_pivot(left,right, l,l+size/3,size/3);
    int *m2 = _cal_pivot(left,right,l+size/3,l+((size/3)*2),size/3);
    int *m3 = _cal_pivot(left,right,l+((size/3)*2),r,size/3);

    /* ３つあるこの要素の内、一つでもNULLポインタが返ってきた場合はm1を返す */
    if(m3 == NULL)return m1;

    /* m1 >= m2 >= m3 になるように並び替える */
    if(*m1 < *m2)swap(m1,m2);
    if(*m1 < *m3)swap(m1,m3);
    if(*m2 < *m3)swap(m2,m3);
    return m2;
}

/* 半開区間 [left,rigth)  */
int* cal_pivot(int* left, int* right)
{
    /* 三分木の葉の要素の数を求める */
    int size = 1;
    while(size < right-left)size *= 3;
    return _cal_pivot(left,right,left,right,size); 
}

 ```

 # 最終的な実装

 ```c
#include <stdlib.h>
#define BOOL char

/* ソートは昇順にする return a<b にすると降順にソートができる */
BOOL cmp(int a,int b)
{
    return a>b;
}

/* 値を交換する */
void swap(int* i,int* j)
{
    int tmp = *i;
    *i = *j;
    *j = tmp;
    return;
}

/* [left,right)の中の[l,r)について再帰的に計算する */
int* _cal_pivot(int* left, int* right,int* l, int* r,int size)
{
    /* [l,r)が[left,right)の範囲外のときNULL(ポインタ)を返す */
    if(right <= l)return NULL;

    /* 範囲に要素が1個だけのときそのポインタを返す */
    if(size==1)return l;

    /* 三分木の子の要素の中央値を再帰的に求める */
    int *m1 = _cal_pivot(left,right, l,l+size/3,size/3);
    int *m2 = _cal_pivot(left,right,l+size/3,l+((size/3)*2),size/3);
    int *m3 = _cal_pivot(left,right,l+((size/3)*2),r,size/3);

    /* ３つあるこの要素の内、一つでもNULLポインタが返ってきた場合はm1を返す */
    if(m3 == NULL)return m1;

    /* m1 >= m2 >= m3 になるように並び替える */
    if(*m1 < *m2)swap(m1,m2);
    if(*m1 < *m3)swap(m1,m3);
    if(*m2 < *m3)swap(m2,m3);
    return m2;
}

/* 半開区間 [left,rigth)  */
int* cal_pivot(int* left, int* right)
{
    /* 三分木の葉の要素の数を求める */
    int size = 1;
    while(size < right-left)size *= 3;
    return _cal_pivot(left,right,left,right,size); 
}

/* 半開区間 [left,right) をソートする */
void quick_sort(int* left, int* right,BOOL (*cmp)(int a,int b))
{
    /* ソート範囲に値が１つしかない場合はreturn */
    if(right-left <= 1)return;
    
    /* cal_pivot関数でpivotを求める */
    int *pivot = cal_pivot(left,right);  
    
    /* 配列の両端から線形探査し、値を交換していく */
    int *l = left , *r = right-1;
    while(l<r)
    {
        while(cmp(*pivot,*l) && l < pivot)l++;
        while(cmp(*r,*pivot) && pivot < r)r--;
        swap(l,r);
        if(l ==  pivot) pivot = r;
        else if(r == pivot) pivot = l;
        if(l != pivot)l++;
        if(r != pivot)r--;
    }

    /* pivotを除き、残りの範囲をそれぞれソート */
    quick_sort(left,pivot,cmp);
    quick_sort(pivot+1,right,cmp);
}
 ```

 # 計算速度の実測

 以上の実装を実際に使用してランダムに生成した100万要素を持つ配列に対しての処理時間を測ってみました。比較のため以下の３つで図ります。
 - sort1 : pivotをソート範囲の中央の値にしたものでクイックソートするもの
 - sort2 : pivotをcal_pivot関数使って求めクイックソートするもの
 - sort3 : c++のsort関数を使ってソートするもの

# 結果
以下は結果です。
```
sort1
real    0m0.318s
user    0m0.306s
sys     0m0.012s
sort2
real    0m0.624s
user    0m0.617s
sys     0m0.007s
sort3
real    0m0.538s
user    0m0.527s
sys     0m0.011s
```
　結果は、そのままsort1が一番早くて、その次にsort3,sort2という順になりました。
　cal_pivot関数は計算量が$O(n$$log_3$$n)$程度になるので、sort2は大体その分計算速度がsort1より遅くなった感じです。しかし、最悪計算量は改善されるはずなので、配列の真ん中の値をpivotとした場合とcal_pivot関数を用いた場合で、999個の要素を持つランダム配列に対して1億回ほど試行し実際の中央値との誤差の量と回数をプロットしました。

### 配列の真ん中の値をpivotとした場合
 ![Figure_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3472759/9add2346-f7da-8fe1-5afc-2489c8c8faf3.png)

### cal_pivot関数を用いた場合
![Figure_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3472759/22a51b09-ec99-a112-ff9a-4c4fa7a6b4c1.png)

配列の真ん中の値をpivotとした場合は、ランダムに値を選ぶのと変わらないため、誤差が一様に現れました。その反対に、cal_pivot関数を用いた場合では、グラフのような分布となり、-181から+173の範囲で誤差が現れ、標準偏差を計算したところ31.4809となりました。
　ここからもわかるように最悪計算量について大きく削減されていると推測できます。（数学的証明については門外漢のためご容赦ください）

# 終わりに

　以上、簡単なクイックソートの実装と最悪計算量の改善をしてみました。何か間違いなどございましたら、ご指摘いただけると幸いです。