glibcのexp(x)をいじめる

glibcのexp関数は、倍精度浮動小数点数の指数関数をソフトウェアで計算するものです。 glibcの実装（IBM Accurate Mathematical Libraryの実装を用いたもの）は、おそらく完全精度を達成しています。 完全精度とは、入力がどんな値であっても真の値に最も近い浮動小数点数に丸めた結果を返すことを言います。

しかし、これを実現するのは非常に困難であることが知られています（テーブルメーカーのジレンマ）。 glibcの実装では、以下の三ステップを踏むことで、完全精度を達成しています。

1. まず倍精度浮動小数点数演算とテーブル引きを駆使して近似値を求め、誤差を加味しても丸めた結果に影響がないと判断されれば、その値を返す
2. 144bitの多倍長演算を使って近似値を求め、誤差を加味しても丸めた結果に影響がないと判断されれば、その値を返す
3. 768bitの多倍長演算を使って近似値を求め、その値を返す

なお、倍精度指数関数については、triple-double（疑似六倍精度）で注意深く計算すれば、完全精度を達成できることが知られています（crlibmによる証明があります→https://hal-ens-lyon.archives-ouvertes.fr/ensl-01529804/file/crlibm.pdf）。

このような仕組みになっているため、多くの場合には高速であるものの、特定のワーストケース入力の場合には非常に遅くなります。

2.のステップに突入する数として、簡単な例に-0.1841があります。exp(-0.1841)を計算すると、普通の数に比べて90倍ほど低速です。

3.のステップに突入する数はそんなに簡単には見つかりません。 exp(x)の丸めがぎりぎりになる値を探索するという研究の結果（→https://hal.inria.fr/inria-00072594/file/RR2000-35.pdf）を参照すると、0x1.9e9cbbfd6080bp-31が最も丸めが難しいと記されていました。 これを入力してみると、普通の数に比べて1000倍ほど低速でした。

[Wandbox]三へ( へ՞ਊ ՞)へ ﾊｯﾊｯ

exp(0x1.9e9cbbfd6080bp-31)を正確に計算してみると、以下のようになり、丸める方向を正しく判断するために小数点以下112bit以上の計算が必要なことがわかります。

	0x1.0000000000000000000000000000000
	0x0.000000033D3977FAC10160000000000
	0x0.00000000000000053EFE9FFA55F3F7B...
	0x0.000000000000000000000005AA0EAD6...
	0x1.000000033D397800000000000002A51...