完全精度 expf 関数の作り方

C言語では、expf関数などの超越関数の精度は一切保証されないことになっています*1。 そうは言っても、数学的な結果に近い結果が返ってくることが、常識的に信じられています。

数学的な結果に最も近いfloatやdoubleが返ってくることを、完全精度（丸め誤差0.5ulp保証*2）といいます。また、数学的な結果に最も近い、または二番目に近い値が返ってくることを、丸め誤差1ulp保証といいます。

一般に、超越関数の結果は、丸め境界に非常に近い値になり得るため、完全精度を実現するのは困難です（テーブルメーカーのジレンマ）。一方、1ulp保証ならば現実的な計算量で実現可能です。

ちょっと調べてみたところ、手元でビルドしたg++8.2では、double std::exp(double)は完全精度っぽいですが、float std::exp(float)は単にdouble std::exp(double)関数の結果を丸めた結果*3を返しているような挙動を示しており、完全精度になっていません。 そこで、自分で作ってみます。

完全精度expf関数の作り方

ここでは、高速化については考えず、とにかく完全精度を追求するものとします。また、テーブルに頼ることはできるだけ避けます。

正しい方針

実は、expfであれば完全精度を達成することは意外と簡単です。

最も丸め境界に近いのはexpf(-0x1.d2259ap+3) = 0x1.fa6636p-22 であり、真の値は0x1.fa66350000001p-22付近にあります。この場合であっても、double精度で正確に計算し最後に一回だけ丸めればぎりぎり有効数字が足ります。

しかし、そのためには簡易的なdouble-double演算*4が必要になってくるため、やや計算量が高くつきます。

どのくらい雑にやってよいのか

実際には丸め境界に近い値はそれほど多くありません。double精度の下2bitを正確に求める必要があるのは、上記の例に加えて

  • expf(-0x1.e1dbe2p-8) = 0x1.fc3fd2p-1 真の値は0x1.fc3fd10000002p-1付近
  • expf(-0x1.c1c4b8p-10) = 0x1.ff1f4ep-1 真の値は0x1.ff1f4effffffcp-1付近
  • expf(0x1.fdff02p-17) = 0x1.000100p+0 真の値は0x1.0000ff0000003p+0付近
  • expf(0x1.cd3982p-14) = 0x1.000734p+0 真の値は0x1.000734ffffffcp+0付近
  • expf(0x1.8d7cb6p-12) = 0x1.0018dap+0 真の値は0x1.0018d90000004p+0付近

くらいしか存在しないため、運が悪くなければ慎重にdoubleで計算するだけで丸め誤差の問題を発生させずに済むかもしれません。

形式的な証拠とはなりえませんが、定義域のすべての値について正しく丸められていることを確認すれば、完全精度を達成していること自体は確かめられます。 expfの定義域は-104.0fから88.0fほどであり、22.4億通り（231よりちょっと多いくらい）で正しく丸められているかを確認すれば目的を達成できます。 これは、最近の*5パソコンであれば、特に工夫をしなくても1時間程度で終わります。

実際にとった方針

基本的にはテーラー展開を使って求めることになりますが、引数の値が大きい場合、収束に必要な項数が多くなります。これは誤差の増大も意味するため、避けるべきです。

方針としては、 \exp(s+t) = \exp(s) \times \exp(t) をうまく活用していきます。 まず、 \exp(s)が簡単に計算できるように、かつtの絶対値が小さくなるように、うまくsを選びます。 こうすると、 \exp(t)をテーラー展開で求めるときに必要な項数が少なくなります。

\exp(s) の計算法

\exp(s) を簡単に計算する手法として、メジャーなのは以下の二つです。

  • s を \ln(2) の整数倍になるようにとる。すると単に二のべき乗となり、std::ldexp関数で効率よく計算できる。tの絶対値は0.5以下となる。
  • \exp(s) のテーブルを持つ。tの絶対値を小さくすることができるため、上の方法より収束が速い。収束の速さはテーブルサイズとのトレードオフとなる。

しかし、今回は難しいことを考えず、sとして整数をとりました。 \exp(s)は繰り返し二乗法により求めました。これによりいくらか誤差が乗りますが、最終的な結果への影響はないようです。

\exp(t) の計算法

tの絶対値が0.5以下の時、| t15 / 15! | < 2^{-55} であるため、double精度を出すためには14次までのテーラー近似を利用すれば十分です。

なお、テーラー展開を13次までとした場合は最終的な結果に誤差が発生しました。

テーラー近似は展開中心点から離れるほど近似度が悪くなることが知られています。この問題を解決する方法として、minimax近似というものあります。 minimax近似は、特定の区間内での最大誤差が最も小さくなる近似多項式です。しかしその係数は複雑になります。

今回はややこしいことを考えず、単純なテーラー近似を利用しました。

ここで、テーラー近似を求める段階で注意が必要です。

一般に多項式の値を求める場合、乗算回数が少なくなる*6ホーナー法が用いられます。

ホーナー法とは、 a_0 + x \times ( a_1 + x \times ( a_2 + x \times a_3 ) ) のように計算する多項式計算法のことです。

この方法は、係数の誤差が蓄積する性質があるため、限られた精度で正確な値を求めるのに向いていないようです（むしろ乗算回数が少なくなることから、多倍長演算向けだと思われます）。

今回は愚直に、 \sum_{i=0}^{14} c_i t^{i} を求めることにしました。ここで、絶対値の小さいほうから足していくことで計算誤差を最小化します。

実際のコード

constexpr std::array<double, 15> make_exp_coeff() {
        std::array<double, 15> arr {};
        double fact = 1.0;
        for( int i = 1; i < 15; ++i ) {
                fact *= i;
                arr[i] = 1. / fact;
        }
        return arr;
}

double power( double x, unsigned int n ) {
        double ret = 1.0;
        for( ; n; n >>= 1 ) {
                if( n & 1 ) { ret *= x; }
                x *= x;
        }
        return ret;
}

// |t| <= 0.5  --> | t^15/15! | < 2^-55
double exp_taylor14( double t1 ) {
        constexpr auto coeff = make_exp_coeff();
        const double t2 = t1 * t1;
        const double t4 = t2 * t2;
        const double t8 = t4 * t4;
        double ret = 0.0;
        ret += t8 * t4 * t2      * coeff[14];
        ret += t8 * t4      * t1 * coeff[13];
        ret += t8 * t4           * coeff[12];
        ret += t8      * t2 * t1 * coeff[11];
        ret += t8      * t2      * coeff[10];
        ret += t8           * t1 * coeff[9];
        ret += t8                * coeff[8];
        ret +=      t4 * t2 * t1 * coeff[7];
        ret +=      t4 * t2      * coeff[6];
        ret +=      t4      * t1 * coeff[5];
        ret +=      t4           * coeff[4];
        ret +=           t2 * t1 * coeff[3];
        ret +=           t2      * coeff[2];
        ret +=                t1 * coeff[1];
        return ret + 1.0;
}

float myExp( float x ) {
        constexpr double exp_plus_one  = 0x1.5bf0a8b145769p+1;
        constexpr double exp_minus_one = 0x1.78b56362cef38p-2;
        const int s = std::round( x ); // round to nearest(ties to max magnitude)
        const double t = x - s; // |t| <= 0.5

        const double exp_s = s > 0 ? power( exp_plus_one, s ) : power( exp_minus_one, -s );
        const double exp_t = exp_taylor14( t );
        return static_cast<float>( exp_s * exp_t );
}

検証コード

#include <iostream>
#include <cstdint>
#include <cstring>
#include <cassert>

#include <kv/interval.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

using idd = kv::interval<kv::dd>;

template<class To, class From>
To bit_cast(const From& from) noexcept {
        To to;
        static_assert( sizeof to == sizeof from );
        std::memcpy(&to, &from, sizeof to);
        return to;
}


float next( float x ) {
        assert( !std::signbit(x) );
        return bit_cast<float>( bit_cast<uint32_t>(x) + 1 );
}

float before( float x ) {
        assert( !std::signbit(x) );
        if( x == 0 ) { return -0.0f; }
        return bit_cast<float>( bit_cast<uint32_t>(x) - 1 );
}

idd intervalFloat( float x ) {
        double y = next(x);
        double z = before(x);
        return idd( (x+z)/2, (x+y)/2 );
}

float nearestFloat( idd ia ) {
        // 二重丸め
        float x = (float)(double)mid(ia);
        float y = next(x);
        float z = before(x);
        idd ix = intervalFloat(x);
        idd iy = intervalFloat(y);
        idd iz = intervalFloat(z);
        if( subset(ia, ix) ) { return x; }
        if( subset(ia, iy) ) { return y; }
        if( subset(ia, iz) ) { return z; }
        std::cerr.precision(50);
        std::cerr << x << std::endl;
        std::cerr << y << std::endl;
        std::cerr << z << std::endl;
        std::cerr << ia.lower() << std::endl;
        std::cerr << ia.upper() << std::endl;
        throw nullptr;
}

#include "my_exact_expf.hpp"

void check( float x ) {
        idd ix(x);
        auto exact_exp_range = exp(ix);
        float nearest = nearestFloat(exact_exp_range);
        float my = myExp(x);

        if( nearest != my ) {
                std::cerr << "x     " << x << std::endl;
                std::cerr << "x     " << std::hexfloat << x << std::defaultfloat << std::endl;
                std::cerr << "my    " << my << std::endl;
                std::cerr << "exact " << exact_exp_range << std::endl;
                std::cerr << abs(mid(exact_exp_range) - my) / (nearest - before(nearest)) << " ulp" << std::endl;
        }
}

int main() {
        std::cerr.precision(50);
        const uint32_t minimum = bit_cast<uint32_t>(-104.0f);
        const uint32_t maximum = bit_cast<uint32_t>(88.0f);
        for( uint32_t i = bit_cast<uint32_t>(-0.0f); i < minimum; ++i ) {
                check( bit_cast<float>(i) );
        }
        for( uint32_t i = bit_cast<uint32_t>(0.0f); i < maximum; ++i ) {
                check( bit_cast<float>(i) );
        }
        return 0;
}

この検証コードのコンパイルには、区間演算ライブラリkvが必要です。

github.com

なお、そのままコンパイルするにはboostが必要ですが、このライブラリはC++03対応を維持するためにC++11で導入された標準関数を使っていないようです。 そのため、以下の強引な置換により、boostへの依存は解消されます。

sed -i s/boost/std/g
sed -i s/enable_if_c/enable_if/g

コンパイルオプションは、高速化のため、g++-8 -std=c++17 -march=native -O3 -DKV_FASTROUND -DKV_USE_TPFMA -I ../../kv/を指定しました。使ったCPUはIntel(R) Core(TM) i5-7200U CPU @ 2.50GHzです。かかった時間はほぼ一時間ぴったりでした。

区間演算によりexp(x)の真値を含む区間が求まります。正しくfloatに丸めた場合の可能性が二つ以上ある場合はthrow nullptrされるはずですが、そのようなことは起こりませんでした。 また、(float)(double)mid(ia)の部分はdoubleに丸めてからfloatに丸めているので正しく丸められていない可能性がありますが、依然1ulp保証ではあるので前後のfloat値を確認しています。 本当は正しく求める方法もあるのですが、手を抜きました。

このプログラムが何も出力せずに終了したことにより、実装したmyExp関数は完全精度であることを確かめました。0.0付近の非正規化数であっても0.5ulp保証されています。

ただし、inf付近でどちらに丸めるべきかについては、よくわからなかったので、その部分について0.5ulp（？？）保証されているかはわかりません。

g++のよくわからない挙動

g++では、なぜかstd::expf関数が使えません。::expf関数はあるらしいです。::expf関数はおしいですが完全精度になっていません（この記事で書いたコードが意味のないものとならなくてよかったです）。また、type detectionイディオムを使った確認により、std::exp(float)オーバーロード自体は存在し、返り値の型はfloatであることを確認しました。

#include <cmath>

template<class T>
struct TD;

int main() {
  TD<decltype(std::exp(1.0f))> td;
}

なぜstd::exp(double)と::expf(float)は完全精度であるのにstd::exp(float)だけ精度が低いのでしょうか？

::expf 関数の精度

以下の98個の入力について0.5ulpを超える誤差があることを確認しました。最大の誤差は0.5000031997 ulp 未満でした。

  • -0x1.c1c4b8p-10 （上で紹介した入力のうちの一つ）
  • -0x1.abb574p-9
  • -0x1.f5723ep-9
  • -0x1.f9c67ep-9
  • -0x1.fd2ddcp-9
  • -0x1.021d74p-8
  • -0x1.2e4ab4p-8
  • -0x1.400432p-8
  • -0x1.4b0f34p-8
  • -0x1.4be042p-8
  • -0x1.4d8260p-8
  • -0x1.4dcdc2p-8
  • -0x1.58ef7ap-8
  • -0x1.db8546p-7
  • -0x1.f691dap-7
  • -0x1.f99ab6p-7
  • -0x1.02f812p-6
  • -0x1.074b54p-6 （最大誤差となる入力。真値は0x1.f7d67affff94ap-1付近ですが、0x1.f7d67cp-1が返ってきます）
  • -0x1.a59e4ep-6
  • -0x1.b17ee0p-6
  • -0x1.b1db4ap-6
  • -0x1.8e02c8p-5
  • -0x1.cd9502p-4
  • -0x1.d0a5b6p-4
  • -0x1.fdcbf6p-4
  • -0x1.2b6a50p-3
  • -0x1.0faa34p-2
  • -0x1.19662ep-2
  • -0x1.d66d74p-2
  • -0x1.074b34p-1
  • -0x1.c2aaeap-1
  • -0x1.df1c00p-1
  • -0x1.efdceap-1
  • -0x1.08abe2p+0
  • -0x1.74a396p+0
  • -0x1.7f4296p+0
  • -0x1.b75242p+0
  • -0x1.ff2602p+0
  • -0x1.7dd9aap+1
  • -0x1.d2445ap+1
  • -0x1.f66118p+1
  • -0x1.2b0b04p+2
  • -0x1.90ced0p+2
  • -0x1.bc933ep+2
  • -0x1.e4c6c0p+2
  • -0x1.e8ebf0p+2
  • -0x1.1089b0p+3
  • -0x1.287d08p+3
  • -0x1.93e788p+3
  • -0x1.d2259ap+3 （上で紹介した、最も丸め境界に近くなる入力）
  • -0x1.de6f7cp+3
  • -0x1.1f8dc4p+4
  • -0x1.433aaep+4
  • -0x1.b43edcp+4
  • -0x1.c95efap+4
  • -0x1.23134ap+5
  • -0x1.444328p+5
  • -0x1.2219dep+6
  • -0x1.2c4726p+6
  • -0x1.5ce26ap+6
  • 0x1.6438b6p-8
  • 0x1.9785f6p-8
  • 0x1.1c2c32p-6
  • 0x1.2dc3f0p-6
  • 0x1.2e3554p-6
  • 0x1.c649eap-6
  • 0x1.3a7784p-5
  • 0x1.980654p-5
  • 0x1.eaddb4p-5
  • 0x1.eff8cep-5
  • 0x1.4feb94p-4
  • 0x1.599d64p-3
  • 0x1.72793ap-3
  • 0x1.b28da4p-3
  • 0x1.f3518cp-3
  • 0x1.54f076p-2
  • 0x1.0d5f92p-1
  • 0x1.39a300p-1
  • 0x1.70f4e4p-1
  • 0x1.b94234p-1
  • 0x1.bf4b50p-1
  • 0x1.13f37ap+0
  • 0x1.c8264ap+0
  • 0x1.4c1e18p+1
  • 0x1.69512ep+1
  • 0x1.172096p+2
  • 0x1.24b172p+3
  • 0x1.3983a8p+3
  • 0x1.56f768p+3
  • 0x1.ddc228p+3
  • 0x1.5d4bd4p+4
  • 0x1.6255f4p+5
  • 0x1.91fceap+5
  • 0x1.d8afbep+5
  • 0x1.097f6ap+6
  • 0x1.15298ap+6
  • 0x1.1ddcb8p+6
  • 0x1.3829eep+6

*1:超越関数だけでなく、sqrt関数以外の無理関数は精度が保証されません。

*2:ulpはunit in the last placeの略で、最後の桁の重みのことです。数学的な結果に最も近い値に丸めた場合、その丸め誤差は0.5ulpになります。

*3:これでは二回丸めが発生し、正しくありません。なぜ二回丸めると正しくないかというと、0.495→0.5→1.0のような間違いが発生するためです。

*4:doubleを二個使う多倍長表現のことをdouble-doubleと呼びます。疑似四倍精度とも呼ばれます。

*5:FMA命令に対応しているHaswell以降のCPUでは、double-double演算が高速化します。FMA命令は計算結果を一回だけ丸めるため、演算誤差の補正が簡単になるためです。

*6:必ずしも乗算回数が最小となるわけではありません。さらに乗算回数が少ない方法として、Knuthのadaptation of coefficientsという方法があるようです。