doubleの基数ソートができないといつから錯覚していた？

Fastestメモ

Fastest道を志す方なら、整数のソートには基数ソート(radix sort)が有力だという事を知っていると思います。実際、32bit整数のソートなら普通にstd:sortに任せるよりかなり速いです。今回の記事は整数のソートに基数ソートを用いるのは当然として受け入れてもらった上で、(正の)doubleについても基数ソートができるよという話をします。それに加えて、競技プログラミングでありがちな正整数a,bが1<=a,b<=1,000,000,000という制約下での(long double)b/aのデータを64bitに無理やり詰め込んで基数ソートができそうだという報告をします。

ではまずdoubleで基数ソートができるという話です。doubleというデータ型のbitごとの内部表現を見てみます。64bitの内訳は以下になっています。
・1 bit (63番目のbit) 正負の表記。正だと0、負だと1。
・11 bit (52-62番目のbit) 指数部
・52 bit (0-51番目のbit) 仮数部
このうち指数部の11bitの中身はまんま整数と同じルールで書かれています。仮数部の52bitも小数を2進法で表したものになっており、この52bitの大小関係も上位bitから確認すれば大丈夫です。そして上位が指数部、下位が仮数部になっているので、最上位の正負を表す1bitを除いた部分は内部表現としては整数と同様に上位のbitから確認すれば大小比較できることが分かります。
で、どうやってdoubleを基数ソートするの？と思うかもしれませんがこれは単にdoubleのbitの内部情報を保ったままlong longに入れてしまうだけでとても簡単です。

typedef long long ll;
int main() {
	int N;
	cin >> N;
	ll A[200001];
	for (int i = 0; i < N; i++) {
		int a, b;
		cin >> a >> b;
		double d = (double)b / a;
		A[i] = *(ll*)&d;
	}
}

「*(ll*)&d;」の部分でdoubleのbitの内部情報を保ったままlong longに入れています。&dでdoubleのポインタを持ってきて、(ll*)でそれをlong longのポインタに変換して、*でポインタの中身をlong longとして取り出しています。これによりdoubleが正であれば大小関係を保ったままlong long A[]にデータが入っているので基数ソートができます。ソートした後、逆に*(double*)(A + i)といった具合に書けばA[i]をdoubleとして値を取り出せると思います。
doubleが負の場合は残念ながら場合分けをしてやるしかありません。long longの場合は負でも最上位の正負bitだけを反転させれば全体の大小関係が正しくなりますが、負のdoubleは指数部と仮数部で絶対値を表しているようなのでうまくいきません。

続いて、正整数a,bが1<=a,b<=1,000,000,000という制約下での(long double)b/aのデータを64bitにうまく詰め込んで基数ソートができそうという話です。
1,000,000,000は2^30弱なので、b/aという割り算の結果の大小関係を表現するには60bitの精度が必要です。再度doubleの内部表現を見ると、仮数部が52bitなので精度が足りません。単純な対策として知られるのがlong doubleを使う方法で、これは15bitの指数部と64bit(実質63bit)の仮数部なので、60bitの精度という要求をクリアします。今回は(long double)b/aを基数ソートしたいのですが80bitのlong doubleは64bitのlong longに収まらず、80 bit分を2つの整数にまたがって管理・比較すると基数ソートに時間がかかるようになり、他の方法に対し速度的な優位性がなくなっていくうえ煩雑です。(他の方法はa,bを整数ペアのまま保持して比較関数を工夫してstd:sortにかける方法を想定しています。)
そこでまずは64bitに収めるため仕分けをします。まず正負bitは仕分けします、これで1bit稼ぎました。次にdoubleの指数部は11bitが使われていますが、今回は1e-9から1e9までの小数なので2^-30から2^30までの範囲の指数ということで整数60=30+30が表現できるよう64=2^6で、指数部が6bitあれば十分です。この工夫により仮数部に58bitが使えることになりましたが、今回の要求精度は60bitと見積もっており、実際に実験してみるとやっぱり最悪60bit必要で不足でした。
そこであきらめず訳が分からないコードを書くのがFastest道ということで、場合分けをします。精度が60bit必要な最悪ケースというのはAとBがともに1e9の近辺のときに限られるので、(long double)B/Aはほぼ1です。つまり、最悪ケースを考えて精度を高める必要があるのは指数部の絶対値が小さい時であり、指数部の絶対値が大きくなったら精度は下げてもよい(仮数部のbit数を減らしてもよい)という事になります。60bitを仮数部に捧げると指数部は4bitしかないため16種しか値をとれませんが、16種の指数部のbitパターンのうち最小の0000と最大の1111のときは精度を下げ(仮数部のbit数を減らし)、6bitの長さの第2指数部を持つようにします。
つまり、

〇最上位4bit(第1指数部)が0000でも1111でもないとき、
・4 bit (60-63番目のbit) 第1指数部 (0001から1110までの14種の数字で絶対値の小さい指数を表現(-6~+7))
・60 bit (0-59番目のbit) 仮数部
〇最上位4bit(第1指数部)が0000もしくは1111のとき、
・4 bit (60-63番目のbit) 第1指数部 (0000もしくは1111)
・6 bit (54-59番目のbit) 第2指数部 (本来の-30から+30までの指数を表現できる)
・54 bit (0-53番目のbit) 仮数部

これにより「上位bitから比較すれば大小が分かる」という条件を保ったまま、指数部の絶対値が小さいとき60bitの精度を持ち、指数部の絶対値が大きくなると54bit精度に切り替えるのが実現できるという寸法です。第一指数部の絶対値が6とか7とかまでしか表現できませんが指数部で2^6=64倍動けば要求精度はそのまま6bitほど減りますのでちょうどうまい具合に収まります。

正のlong doubleを無理やりunsigned long longに入れ込むコードを示します。

typedef unsigned long long ull;
ull ld2ull(long double LD) {
	ull A = *(ull*)&LD ^ 1ull << 63;
	ull B = *((ll*)&LD + 1) - 16352 & 63;
	if (B >= 39) return 15ull << 60 | B << 54 | A >> 9;
	else if (B <= 24) return B << 54 | A >> 9;
	else return B - 24 << 60 | A >> 3;
}
int main() {
	int N;
	cin >> N;
	ull A[200001];
	for (int i = 0; i < N; i++) {
		int a, b;
		cin >> a >> b;
		long double LD = (long double)b / a;
		A[i] = ld2ull(LD);
	}
}

こちらではビットシフトの演算子を使うのでAにlong longでなくunsigned long longを使っています。Aは内部表現を変えてしまったのでこれを小数に戻すのは面倒です。Fastest上、いつか必要に迫られたらやるかもしれませんが今回はA自体を浮動小数点数に戻すコードは考えません。どちらかというと、このAと必要な情報(元配列のindexや元のa,bなど)をペア等にしておいて、Aのbit情報に従い基数ソートすることを応用として想定します。速度の面で考えると、long doubleの割り算や64bitの基数ソートのせいで爆速というわけではありませんが、それでもstd:sortに比べれて速くなる可能性を秘めております。
少しコードの説明をします。*(ull*)&LDは先ほどと同様の処理ですが、80bitの情報を64bitに入れるのでlong doubleの下位bitである仮数部64bitだけがそのままでてきます。long doubleの仮数部はdoubleと違って最上位bitは常に1になっていてその一つ下のbitからが目的のbitなので注意です。*((ll*)&LD + 1)で取り出しているのが指数部です。0-63までの値に収めるために引き算などをしています。