glibcのexp関数を読んでみる

最近のglibcは数学関数の高速化に取り組んでいるため確実なことは言えませんが、多くの環境ではC/C++の数学関数はこのglibcの実装を使うことになっているはずです。

glibcのexp関数の実装は、IBM Accurate Mathematical Libraryを使っているようです。

参考：glibc/e_exp.c at 3737a0d70cac3ba09aa8f32709c4be8d70f41207 · lattera/glibc · GitHub

このライブラリは、（おそらく）（IEEE754 binary64の最近接丸めに対して）完全精度を達成するライブラリ実装です。完全精度とは、数学的な結果に最も近い、doubleで表せる値に丸めた値が返ってくることを言います。

一般に、超越関数を完全精度で計算するライブラリ実装は非常に計算コストが高いものとなってしまいます。このライブラリでは、雑な計算で大丈夫な場合は雑な計算ですまし、どちらに丸めてよいかきわどい時のみ高精度な計算を行うという仕組みで極力コストを削減しています。

どちらに丸めてよいかきわどい時は計算が遅くなります。たとえば、std::exp(1.0000000001412739)の計算は、他の引数に比べて約90倍ほど低速です。

ライブラリ実装の中身を見ていきます。

このライブラリ実装では、union { int32_t i[2]; double x; }が濫用されているため読みづらいです。しかし、大枠でやっていることは非常に単純で、

入力xの絶対値が小さすぎず大きくない値なら、まじめに計算する
入力xの絶対値が小さすぎる場合は、1.0を返す
入力xの絶対値が大きすぎる場合なら、0.0かinfを返す
入力xがNaNなら、NaNを返す
入力xが、exp(x)がオーバーフロー/アンダーフローギリギリの値の場合、気を付けながら真面目に計算する

という仕組みになっています。

「真面目に計算する」の部分は、きれいに書くと以下のような感じです。

#include <cstring>
#include <cstdint>
#include <cmath>

template<class T, class U>
T bit_cast( U u ) {
        T t;
        std::memcpy( &t, &u, sizeof t );
        return t;
}

#include "tbl.h"

double expm1_poly3( double x ) {
        static constexpr double p2 = 0x1.00000000003dcp-1;
        static constexpr double p3 = 0x1.5555555555a0fp-3;
        return x + x*x*( p2 + x * p3 );
}

constexpr double log2e   = 0x1.71547652b82fep+0;
constexpr double ln_two1 = 0x1.62e42fefa3800p-1;
constexpr double ln_two2 = 0x1.ef35793c76730p-45;

constexpr double three33 = 0x3.p+33;
constexpr double three51 = 0x3.p+51;

double e_exp( double x ) {
        if( std::abs(x) > 0x1.fffffp-55 && std::abs(x) < 0x1.62002p+9 ) {
                const double tmp1 = x * log2e + three51;
                const double bexp = tmp1 - three51;

                const double eps1 = bexp * ln_two2;
                const double t    = x - bexp * ln_two1;

                const double tmp2 = t + three33;
                const double base = tmp2 - three33;

                const double del = (t-base) - eps1;
                const double eps2 = expm1_poly3(del);

                const double binexp = bit_cast<double>((bit_cast<uint64_t>(tmp1) + 1023) << 52);
                const uint64_t base_i = bit_cast<uint64_t>(tmp2) - bit_cast<uint64_t>(three33);
                const uint64_t i = (base_i +91136) >> 9;
                const uint64_t j = base_i & 511;

                const double al = coar.x[i*2] * fine.x[j*2];
                const double bet = ((coar.x[i*2]*fine.x[j*2+1] + coar.x[i*2+1] * fine.x[j*2]) + coar.x[i*2+1] * fine.x[j*2+1]);

                const double rem = (bet + bet*eps2) + al * eps2;
                const double res = al + rem;
                const double cor = (al - res) + rem;

                const double retval = res * binexp;
                if( res == (res + cor * 1.000014) ) {
                        return retval;
                } else {
                        return __slow_exp(x);
                }
        }
        // 以下略
}

一つずつ説明していきます。

log2e は $\log_2e$ の近似値です。
bexpは、がに近くなるような整数です。「大きな数を足すことで丸める手法」を用いています。参考：高速な倍精度指数関数expの実装の17ページ目
- binexpは、 $2^{bexp}$ です。ビット演算テクニックによって作成しています。
ln_two1とln_two2は、二つ足し合わせるとほぼ正確にになるような二数で、かつln_two1は仮数部の下11bitが0であるような数です。
- bexp * ln_two1が丸め誤差なく求まるようにするため、そのように選んでいます。bexpは11bitに収まることを利用しています。
baseは、tに最も近い $2^{-18}$ の倍数です。ここでも「大きな数を足すことで丸める手法」が用いられています。
delは、に非常に近い値です。
- eps2は $e^{del}-1$ を多項式近似で求めたものです。
base_iは、に相当する整数で、ビット演算テクニックによって作成しています。
- $e^{base}$ は、表引きで求めます。と言っても10万エントリ以上もあるテーブルを作るわけにはいかないため、上位9bit（i-178）と下位9bit（j）に分けて表を引き、その値を掛け合わせることで求めます。上位9bitを使う側の表がcoarであり、下位9bitを使う側の表がfineです。
- 表引きの結果得られるのは二つのdouble値であり、二つ足すとほぼ正確にやになるような二数で、かつ片方の下26bit（仮数部の下半分）が0であるような数です。
  - alが丸め誤差なく求まるようにするため、そのように選んでいます。
betは $e^{base} - al$ （の近似値）です。
remは $e^{base+del} - al$ （の近似値）です。
resは $e^x / binexp = e^{base+del}$ の近似値です。
corは、res = al + rem;と計算した時の丸め誤差（の近似値）です。倍倍精度演算で用いられるtwo_sumと同じことをしています。

最後のif文では、res = (res + cor * 1.000014)かを判定しています。1.000014の根拠は不明ですが、主にeps2を多項式近似で求めた時に生じた誤差の評価分が含まれているものと思われます（確認していないので自信なしですが、その他の部分は十進八桁精度はあるはずなので多項式近似部分が主たる誤差要因だと思います）。つまり、丸め誤差corが十分に小さく、途中計算で生じた丸め誤差を加味（1.000014倍）しても丸める方向が変わらないと結論付けられる場合はres * binexpを返し、丸め誤差corが大きいため途中計算で生じた丸め誤差を加味するとどちらに丸めるべきか怪しい場合はさらに精度の高い計算（__slow_exp(x)）を行うという仕組みになっています。