This article is a Private article. Only a writer and users who know the URL can access it.
Please change open range to public in publish setting if you want to share this article with other users.

クイックソート徹底解説

Last updated at 2023-12-04Posted at 2023-12-04

はじめに

どうも。aqualengthです。¹

この記事はアルゴリズムを解説する企画でハブられた腹いせに書いた記事(Part2)です。
ソートというのは基本的に標準のソートが最適化されており、特殊な状況(取り得る値が有限,要素数が極めて少ない,整列がほとんど完了している,処理がGPU上で完結する等²)でもない限りはそれを使うべきです。
しかし、C言語など、言語によっては標準のソートが遅い場合があります。
その場合はCPU上で動作する汎用性の高いソートの中でも高速なクイックソートが選択肢に入ると思います。
クイックソートはあらゆるアルゴリズムの中でも上位に入る実装難易度だと思う³ので、実装のレパートリーや実装上での注意点などを交えて解説していきます。
なお、言語はソートが遅いC言語ですが、末尾再帰最適化のない言語で実装する必要が無いとも言い切れないので、末尾再帰を使わない実装についても解説します。⁴

はじめに
目次
クイックソートの特徴について
末尾再帰を使った実装
(末尾再帰でない)再帰を使った実装
ピボットについて
Lomutoのパーティション
Hoareのパーティション
非再帰クイックソート
クイックソート+挿入ソート
乱択クイックソートvs決定的クイックソートvs std::sort

クイックソートの特徴について

クイックソートは分割統治法を用いたソートアルゴリズムで、平均計算量はO($N\log{N}$)です。
一般的に、ランダムなデータに対しては(GPU上で実行するソートや比較を行わないソートを除き)最も高速とされています。

メモリの消費も少なく、最悪O($\log{N}$)のメモリで済みます。

反面、ソート前後で元の順番が保持される保証がない,データによっては非常に遅くO($N^2$)かかる場合があるなどの欠点があります。
後者に関してはイントロソートと呼ばれるクイックソートの改良版を使うことで解決できますが、ここでは解説しません。

なお、ここでは昇順ソートについて解説します。

末尾再帰を使った実装

簡単のため、ここでいくつかの関数を定義します。
具体的な実装については後の章で解説します。

pivot(a[],l,r) … 配列aの半開区間[l,r)から適切なピボットのインデックスを返す
partition(a,l,r,pivot) … 配列aの半開区間[l,r)をpivot以下の半開区間[l,m),pivotより大きい半開区間[m,r)に分割し,その境界mを返す

これらの関数さえ定義してしまえば、末尾再帰を使った実装は非常にシンプルです。

ピボットを求める
パーティションを実行する
分割後の半開区間について、小さい方の半開区間でクイックソートを行う
分割後の半開区間について、大きい方の半開区間でクイックソートを行う

ポイントとしては、大きい方の半開区間を末尾で実行することです。
これにより、大きい方の半開区間は末尾再帰最適化によってループに展開されます。
よって、再帰によって半開区間のサイズは常に半分以下となり、最悪でも$log_2 N$しかメモリを消費しません。

void quickSort(type a[],unsigned l,unsigned r){
	unsigned p=pivot(a,l,r); //ピボットを選択
	unsigned m=partition(a,l,r,p); //パーティション
	if(m-l<r-m){ //前半の方が小さい
		quickSort(a,l,m); //前半を再帰
		quickSort(a,m,r); //後半を再帰(末尾再帰)
	}else{ //後半の方が小さい
		quickSort(a,m,r); //後半を再帰
		quickSort(a,l,m); //前半を再帰(末尾再帰)
	}
}

(末尾再帰でない)再帰を使った実装

全ての言語に末尾再帰最適化があるわけではありません。
例えば、Pythonという言語は(末尾再帰を多用するにも関わらず⁵)つい最近まで末尾再帰最適化がありませんでした。
末尾再帰はループに変換できるという特性があるため、手動でループに変換すれば良いです。

void quicksort(type a[],unsigned l,unsigned r){
	while(r-l>1){
		unsigned p=pivot(a,l,r);
		unsigned m=partition(a,l,r,p);
		if(m-l<r-m){
			quicksort(a,l,m);
			l=m;
		}else{
			quicksort(a,m,r);
			r=m;
		}
	}
}

ピボットについて

ソートの性能を決める重要な要素です。
パーティションが正しく実装されている限り、配列内の任意の値をピボットにすればソートが完了します。
しかし、例えばピボットが常に最大,最小値を返すような実装の場合、ソートの計算量は常にO($N^2$)になります。
逆にピボットが常に中央値を返すような実装の場合、ソートの計算量は常にO($N\log{N}$)になります。
実際にそのような実装は可能(中央値の中央値)ですが、ピボットの選択に時間がかかるため、実際には中央値に近い値が期待される実装を行います。
典型的には、配列内の3つの値を抜き出し、その中央値をピボットするような実装です。

抜き出す3つの値に関しては、先頭3要素や末尾3要素,先頭,末尾,中央の3要素といった決定的な選び方とランダムな3要素を選ぶ確率的な実装があります。
決定的な実装だと人為的に最悪となるケースを入力される可能性があるため、そのような可能性のある状況(競技だとCodeforcesなど)では確率的な実装が望ましいです。
なお、確率的な実装を乱択クイックソートと呼びます。

個人的には配列外参照を気にせずに済む乱数か、先頭中央末尾の実装がおすすめです。

先頭中央末尾の実装は以下の通りです。

unsigned pivot(type a[],unsigned l,unsigned r){
	unsigned m=(l+r)>>1;
	r--;
	if(a[l]<a[m]){
		if(a[m]<a[r]){ //l<m<r
			return m;
		}else{ //l,r <m
			if(a[l]<a[r]){ //l<r<m
				return r;
			}else{ //r<l<m
				return l;
			}
		}
	}else{ //l>m
		if(a[m]<a[r]){ //l,r >m 
			if(a[l]<a[r]){ //r>l>m
				return l;
			}else{ //l>r>m
				return r;
			}
		}else{ //l>m>r
			return m;
		}
	}
}

乱数を使った実装は以下の通りです。
乱数は実装が簡単で、速度と品質のバランスが良いxorshift32を使います。

unsigned _seed=(unsigned)(time(NULL)%4294967295+1); //シード値:0でない任意の値
unsigned rnd(unsigned n){
	_seed^=_seed<<13;
	_seed^=_seed>>17;
	_seed^=_seed<<5;
	return _seed%n;
}

unsigned pivot(type a[],unsigned l,unsigned r){
	unsigned x,y,z;
	x=l+rnd(r-l);y=l+rnd(r-l);z=l+rnd(r-l);
	if(a[x]<a[y]){
		if(a[y]<a[z]){ //x<y<z
			return y;
		}else{ //x,z <y
			if(a[x]<a[z]){ //x<z<y
				return z;
			}else{ //z<x<y
				return x;
			}
		}
	}else{ //x>y
		if(a[y]<a[z]){ //x,z >y 
			if(a[x]<a[z]){ //z>x>y
				return x;
			}else{ //x>z>y
				return z;
			}
		}else{ //x>y>z
			return y;
		}
	}
}

パーティションについて

半開区間[l,r)についてクイックソートがバグらない条件として
戻り値mが[l+1,r-1)の範囲である必要があります。

このためには、ピボットと等しい値はどちらの半開区間にも含まれ得るように分割する必要があります。

パーティションの実装についてはいくつかバリエーションがありますが、代表的な実装としてLomutoのパーティションとHoareのパーティションが有名です。

Lomutoのパーティション

Lomutoのパーティションは、配列を先頭から見ていき、ピボットより大きな値が見つかれば前半に移動させるようなパーティションです。

具体例を見てみましょう。

※ 具体例は後ほど追加予定です

実はstd::uniqueとほとんど同じ処理です。

私は普段std::uniqueに相当する操作をこのように書いています。

for(l=r=0;r<n;l++){
	register type reg=a[r];
	do{r++;}while(r<n&&reg==a[r]);
	a[l]=reg;
}
n=l;

範囲を半開区間に変更し、while内の条件式をrpivotに変更して、値の上書きを値の交換に変更し、do～whileをwhileに書き換えればLomutoのパーティションです。(ほとんど同じと言いつつ、実際に書いてみたら結構変わりました)

一つ注意点があって、Lomutoのパーティションはピボットが末尾の要素である必要があります。
処理の最初に末尾とピボットがある要素を交換してしまえば大丈夫です。

unsigned partition(type a,unsigned left,unsigned right,unsigned pivot){
	unsigned l,r;
	register type p=a[pivot];
	a[pivot]=a[right];
	a[right]=p;
	
	for(l=r=left;;r<right;l++){
		register type reg=a[r];
		while(r<right&&a[r]>p) r++;
		if(r<right){ 
			a[r]=a[l];
			a[l]=reg;
		}
	}
	if(right==l) return l-1;
	return l;
}

Hoareのパーティション

※ この章は書きかけです

クイックソートの最難関であり、本記事のメインコンテンツです。

まずはループの終了条件をlとrの大小が入れ換わる…つまり、left>=rightに設定します。
半開区間ではインクリメントは最後,デクリメントは最初に行うのがセオリーです。⁶
ループの最初と最後のデクリメントとインクリメントを追加します。

具体例を考えます。
出来るだけバグらないよう以下の3つのケースについて考えます。

ピボットが最小値かつ先頭のケース
ピボットが最大値かつ末尾のケース
全ての値が同じケース

ピボットが最小値かつ先頭の場合、正しく動作しないことが分かったので、回避するif文を入れます(入れなくても動作するスマートな方法があれば教えてください)

unsigned partition(unsigned a[],unsigned left,unsigned right,unsigned pivot){
	register unsigned p=a[pivot];
	while(left<right){
		right--;
		while(a[left]<p) left++;
		while(a[right]>p) right--;
		if(left<right){
			unsigned tmp=a[left];
			a[left]=a[right];
			a[right]=tmp;
			left++;
		}
	}
	return left;
}

非再帰クイックソート

ここからは、更なる速度を求めて改良を行います。
まず、再帰呼び出しにはコストがかかるので再帰を無くします。

再帰を非再帰にするには大きく分けて2つの方法があります。
1つはロジックを工夫して再帰を不要にする方法、もう一つはスタックを実装して再帰呼び出しをシミュレーションする方法です。
前者はマージソートや平衡二分探索木が該当し、非常に大きな効果を発揮します。
クイックソートは後者なわけですが、正直高速化としてはあまり効果はありません。
しかし、慣れれば機械的に変換できますし言語によっては効果も大きいため、非再帰化してみます。

#define med3(x,y,z) ((x)<(y)?(y)<(z)?(y):(x)<(z)?(z):(x):(y)<(z)?(x)<(z)?(x):(z):(y))
unsigned partition(unsigned a[],unsigned left,unsigned right,unsigned pivot){
	register unsigned p=pivot;
	while(left<right){
		right--;
		while(a[left]<p) left++;
		while(a[right]>p) right--;
		if(left<right){
			unsigned tmp=a[left];
			a[left]=a[right];
			a[right]=tmp;
			left++;
		}
	}
	if(0==left) return left+1;
	return left;
}
void quicksort(unsigned a[],unsigned l,unsigned r){
	unsigned stack[33<<1]; //最大で32要素(l,r合わせて64要素)
	unsigned slast=0; //スタック内の要素数
	if(r-l>1){
		do{
			if(r-l<=1){
				l=stack[--slast];
				r=stack[--slast]; //取り出すときはl→rの順
			}
			unsigned p=med3(a[l],a[r-1],a[(l+r)>>1]);
			unsigned m=partition(a,l,r,p);
			if(m-l<r-m){
				if(r-m>1){
					stack[slast++]=r; //大きい方をスタックに追加(小さい方を再帰に回す)
					stack[slast++]=m; //追加するときはr→lの順
				}
				r=m;
			}else{
				if(m-l>1){
					stack[slast++]=m; //大きい方をスタックに追加(小さい方を再帰に回す)
					stack[slast++]=l; //追加するときはr→lの順
				}
				l=m;
			}
		}while(slast);
	}
}

クイックソート+挿入ソート

クイックソートは一般的に高速とされていますが、要素数の少ない配列や整列がほとんど完了した配列に対しては遅いです。

逆に要素数の少ない配列や整列がほとんど完了した配列に対しては挿入ソートやノームソートが高速なため、クイックソートである程度ソートした後に挿入ソートを適用することで高速化が見込めます。

打ち切る要素数ですが、およそ50まで要素では挿入ソートが高速と噂されているため、ある程度キリの良い48にします。

挿入ソートはやねうらおさんのものを使います。

#define med3(x,y,z) ((x)<(y)?(y)<(z)?(y):(x)<(z)?(z):(x):(y)<(z)?(x)<(z)?(x):(z):(y))
unsigned partition(unsigned a[],unsigned left,unsigned right,unsigned pivot){
	register unsigned p=pivot;
	while(left<right){
		right--;
		while(a[left]<p) left++;
		while(a[right]>p) right--;
		if(left<right){
			unsigned tmp=a[left];
			a[left]=a[right];
			a[right]=tmp;
			left++;
		}
	}
	if(0==left) return left+1;
	return left;
}
void quicksort(unsigned a[],unsigned left,unsigned right){
	unsigned stack[32<<1]; //最大で32要素(l,r合わせて64要素)
	unsigned slast=0; //スタック内の要素数
	unsigned l=left,r=right;
	if(r-l>47){
		do{
			if(r-l<=47){
				l=stack[--slast];
				r=stack[--slast]; //取り出すときはl→rの順
			}
			unsigned p=med3(a[l],a[r-1],a[(l+r)>>1]);
			unsigned m=partition(a,l,r,p);
			if(m-l<r-m){
				if(r-m>47){
					stack[slast++]=r; //大きい方をスタックに追加(小さい方を再帰に回す)
					stack[slast++]=m; //追加するときはr→lの順
				}
				r=m;
			}else{
				if(m-l>47){
					stack[slast++]=m; //大きい方をスタックに追加(小さい方を再帰に回す)
					stack[slast++]=l; //追加するときはr→lの順
				}
				l=m;
			}
		}while(slast);
	}
	// yaneurao's insertion sort
	for (int i = left+1; i < right; i++) 
	{
		register unsigned tmp = a[i];
		if (a[i-1] > tmp)
		{
			int j = i;
			do {
				a[j] = a[j-1];
				--j;
			} while ( j > 0 && a[j-1] > tmp);
			a[j] = tmp;
		}
	}
}

ポインタを使った実装

std::sortより遅かったので、更なる高速化を試みます。
先人の知恵を借りて、C言語らしくポインタを用いた実装に書き換えます。

ついでに、今までサブルーチンとして実装していたpartitionをクイックソート内に直接埋め込みます。

#define med3(x,y,z) ((x)<(y)?(y)<(z)?(y):(x)<(z)?(z):(x):(y)<(z)?(x)<(z)?(x):(z):(y))
void quicksort(unsigned *left,unsigned *right){
	unsigned *stack[32<<2]; //最大で32要素(l,r合わせて64要素)
	unsigned **slast=stack; //スタック内の要素数
	unsigned *l=left,*r=right;
	if(r-l>47){
		do{
			if(r-l<=47){
				l=*(--slast);
				r=*(--slast); //取り出すときはl→rの順
			}
			unsigned p=med3(*l,*(r-1),*(l+((r-l)>>1)));
			unsigned* m,*rp=r,*lp=l;
			do{
				rp--;
				while(*lp<p) lp++;
				while(*rp>p) rp--;
				if(lp<rp){
					unsigned tmp=*lp;
					*lp=*rp;
					*rp=tmp;
					lp++;
				}
			}while(lp<rp);
			m=lp+(r==lp);
			if(m-l<r-m){
				if(r-m>47){
					*(slast++)=r; //大きい方をスタックに追加(小さい方を再帰に回す)
					*(slast++)=m; //追加するときはr→lの順
				}
				r=m;
			}else{
				if(m-l>47){
					*(slast++)=m; //大きい方をスタックに追加(小さい方を再帰に回す)
					*(slast++)=l; //追加するときはr→lの順
				}
				l=m;
			}
		}while(slast>stack);
	}
	// insertion sort
	for(unsigned *i=left+1;i<right;i++){
		register unsigned tmp=*i;
		if(*(i-1)>tmp){
			unsigned *j=i;
			do{
				*j=*(j-1);
				--j;
			}while(j>left&&*(j-1)>tmp);
			*j=tmp;
		}
	}
}

乱択クイックソートvs決定的クイックソートvs std::sort