今更だけど爆速な基数ソート(Radix Sort)をC++(ほぼC言語)で実装してみた

なぜ最速といわれるクイックソートがあるのに基数ソートを使うのか

・数値かつ要素数が多い場合はクイックソートより圧倒的に速い
・比較を行わないので非常に高速
・元データの状態に関わらず安定した時間で完了する
・個人的に好き

実装しつつテストデータでシミュレーションしてみる

unsigned int arr[10] = {8, 4, 3, 9, 9, 6, 7, 1, 0, 2};

わかりやすい例としてこのようなテストデータを使用

i	0	1	2	3	4	5	6	7	8	9
arr[i]	8	4	3	9	9	6	7	1	0	2

/* ０初期化 */
inline void _clear(unsigned int* arr, int len) {
    long long* p = (long long*)arr;
    len >>= 1;
	while (len-- > 0) {
		*p++ = 0;
	}
}

/* 累積和書き込み */
inline void _satulate(int* arr, int len) {
	for (int i = 1; i < len; i++) {
		arr[i] += arr[i - 1];
	}
}

内部で使う関数を先に定義しておき本体を書いていく

/* 基数ソート(配列, 要素数) */
void radix_sort(unsigned int arr[], int len){
    const int BITS = 2;
    const int RADIX = 1 << BITS;    /* 基数 = 4 (0100) */
    const int BIT_MASK = RADIX - 1; /* 桁マスク = 3 (0011) */

    int group;
    int max_value = 0;
    unsigned int* buf;
    unsigned int* pDest;
    unsigned int* pSrc;
    int counts[RADIX];
    
    _clear(counts, RADIX);
    pSrc = arr;
    /* 数値のBIT数分だけを見て振り分け */
	for (i = 0; i < len; i++) {
		group = pSrc[i] & BIT_MASK;
		counts[group]++;

        /* 初回振り分けついでに最大値を求めておく */
		if (pSrc[i] > max_value) {
			max_value = pSrc[i];
		}
	}

counts に基数グループ毎の数を入れた結果

i	counts[i]	カウント対象
0 ( 00 )	3	8 ( 10 00 )	4 ( 01 00 )	0 ( 00 00
1 ( 01 )	3	9 ( 10 01)	9 ( 10 01)	1 ( 00 01)
2 ( 10 )	2	6 ( 01 10)	2 ( 00 10)
3 ( 11 )	2	3 ( 00 11)	7 ( 01 11)

max_value	9

	/* 最大値が基数未満なら結果を入れて戻る */
	if (max_value < RADIX) {
		/* この場合は単純なバケツソート（不安定） */
		for (i = 0; i < RADIX; i++) {
			while(counts[i] > 0){
				*arr++ = i;
				counts[i]--;
			}
		}
		return 0;
	}

最大値が基数未満ならカウントされた分を順番に入れるだけでソート終了
今回の例では最大値が９なので次の桁のソートをしないといけない

	/* 基数以上の値がある場合は作業領域が必要 */
	buf = new unsigned int[len];
	if (!buf) {
		return -1;
	}

	pDest = buf;
    _satulate(counts, RADIX);

出力先を設定して書き込みの前準備をする

累積和を取ることでその「基数グループの終端位置がわかる」ようになる

counts

i	0	1	2	3
counts[i]	3	6	8	10

	for (i = len - 1; i >= 0; --i) { /* 後ろから見ることによって安定 */
		group = pSrc[i] & BIT_MASK;
		pDest[--counts[group]] = pSrc[i];	/* 属する基数グループの終端から配置していく */
	}
	max_value >>= BITS;	/* 見終わった桁数を破棄 */

pSrc

i	0	1	2	3	4	5	6	7	8	9
pSrc[i]	8	4	3	9	9	6	7	1	0	2
group(2進)	10 00	01 00	00 11	10 01	10 01	01 10	01 11	00 01	00 00	00 10
group(10進)	0	0	3	1	1	2	3	1	0	2

pDest と counts の変化

i	0	1	2	3	4	5	6	7	8	9
pDest[i]	8	4	0	9	9	1	6	2	3	7
group	0			1			2		3
counts[group] : 後 ← 前	0 ← 3			3 ← 6			6 ← 8		8 ← 10

max_value	9 (1001) → 2 (0010)

	/* 全桁見るまでループ */
	for (shift = BITS; max_value > 0; shift += BITS, max_value >>= BITS) {
		/* 元配列と作業領域を交互に使う */
		unsigned int *w = pSrc;
		pSrc = pDest;
		pDest = w;

		/* 初回と同じ流れ（グループ毎カウント -> 累積和 -> インデックスでマッピング） */
		_clear(counts, RADIX);
		for (i = 0; i < len; i++) {
			group = (pSrc[i] >> shift) & BIT_MASK;
			counts[group]++;
		}
		_satulate(counts, RADIX);

		for (i = len - 1; i >= 0; --i) {
			group = (pSrc[i] >> shift) & BIT_MASK;
			pDest[ --counts[group] ] = pSrc[i];
		}
	}

次の桁も追ってみる

shift	2

counts

i	counts[i]	カウント対象
0 ( 00 )	4	0 ( 00 00 )	1 ( 00 01 )	2 ( 00 10 )	3 ( 00 11 )
1 ( 01 )	3	4 ( 01 00 )	6 ( 01 10 )	7 ( 01 11 )
2 ( 10 )	3	8 ( 10 00 )	9 ( 10 01 )	9 ( 10 01 )
3 ( 11 )	0

counts 累積和

i	0	1	2	3
counts[i]	4	7	10	10

pSrc

i	0	1	2	3	4	5	6	7	8	9
pSrc[i]	8	0	9	9	1	4	6	2	3	7
group(2進)	10 00	00 00	10 01	10 01	00 01	01 00	01 10	00 01	00 11	01 11
group(10進)	2	0	2	2	0	1	1	0	0	1

pDest と counts の変化

i	0	1	2	3	4	5	6	7	8	9
pDest[i]	0	1	2	3	4	6	7	8	9	9
group	0				1			2			3
counts[group] : 後 ← 前	0 ← 4				4 ← 8			8 ← 10			10 ← 10

max_value	2 (0010) → 0 (0000)

shift	4

最大の桁まで見終わったのでここでループ終了

	/* 最終書き込み先が作業領域の方なら元配列に反映 */
	if (pDest == buf) {
		for (i = 0; i < len; i++) {
			arr[i] = buf[i];
		}
	}

	delete buf; /* 作業領域の後始末 */
	return 0;
}

最後に元の配列に結果が入るように調整したら完了

最後に

データに適したソートを使えるのが一番いい
ソートアルゴリズムとその特性は知識として持っておきたい

因みに簡易測定結果は以下の通り
・要素数100000
・値重複ありランダム配置
・上記入力に対してclock()関数を使用して大雑把に計測

概要	関数	clock() 差分
C標準	qsort	24
非再帰型	heep_sort	21
再帰型	quick_sort	8
今回実装したもの	radix_sort_int RADIX=256	3

おまけ：負数への対応

bit数に注意が必要だけど負の数対応もとても簡単
イメージ的には -10～10 を 0～20 にしてソートした後もとに戻す感じ

/* int 型が 32bit 前提 */
void radix_sort_int(int arr[], int len) {
	int i;
	unsigned int* uarr = (unsigned int*)arr;
	for (i = len - 1; i >= 0; i--) {
		uarr[i] ^= 0x80000000;
	}

	radix_sort(uarr, len);

	for (i = len - 1; i >= 0; i--) {
		uarr[i] ^= 0x80000000;
	}
}