C標準数学ライブラリを使用した標準正規分布のCDFの正確な計算

Question

標準C数学ライブラリは、標準正規分布のCDFを計算する関数、normcdf()を提供しません。しかし、これは密接に関連する機能を提供します。誤差関数erf()と相補誤差関数erfc()です。 CDFを計算するための最速の方法は√½を表すために事前に定義された定数M_SQRT1_2を用いて、誤差関数を介してしばしばある：

double normcdf (double a) 
{ 
    return 0.5 + 0.5 * erf (M_SQRT1_2 * a); 
} 

明らかに、これは負の半平面における大規模な減法キャンセル苦しんでいるとは適していませんアプリケーションの大半。キャンセルの問題は簡単しかし、一般的erf()よりも幾分低いパフォーマンスのある、erfc()を使用することによって回避されているので、最も頻繁に推奨される計算は次のとおりです。

double normcdf (double a) 
{ 
    return 0.5 * erfc (-M_SQRT1_2 * a); 
} 

いくつかのテストでは、最大ULPエラーが負の半分に発生することを示していますしかし、飛行機はまだかなり大きいです。 erfc()の倍精度実装を0.51 ulpsに正確に使用すると、normcdf()で から1705.44 ulpまでの誤差を観測することができます。ここでの問題は、erfc()への入力における計算誤差が、erfc()に固有の指数スケーリングによって拡大されることです（累乗による誤差拡大の説明については、answer を参照してください）。

次の用紙が正しく、製品例えば√½として任意精度定数を有する浮動小数点オペランドを乗算丸め（ほぼ）一つは達成することができる方法を示しています正しく丸め」、

ニコラスBrisebarreとジャンミシェル・ミューラーを任意の精度定数による乗算 "、IEEE Transactions on Computers、Vol。この論文が主張する方法は、最近のすべての一般的なプロセッサアーキテクチャの実装で利用可能なヒューズの積和演算に依存しており、C言語で公開されています標準の数学関数fma()を使用してください。これは、次のバージョンにつながる：

double normcdf (double a) 
{ 
    double SQRT_HALF_HI = 0x1.6a09e667f3bcd0p-01; // 7.0710678118654757e-01 
    double SQRT_HALF_LO = -0x1.bdd3413b264560p-55; // -4.8336466567264567e-17 

    return 0.5 * erfc (fma (-SQRT_HALF_HI, a, -SQRT_HALF_LO * a)); 
} 

テストでは、これは、以前のバージョンに比べて半分におおよそ最大誤差をカットすることを示しています。前と同じ高精度のerfc()の実装を使用すると、観測される最大エラーは842.71 ulpsです。これは、せいぜい数ulpsの誤差を持つ基本的な数学関数を提供するという通常の目標とはまだまだ離れています。

normcdf()の正確な計算を可能にし、標準のC数学ライブラリで使用できる関数のみを使用する効率的な方法がありますか？

Answer 1

質問で概説されているアプローチの精度の制限に関する1つの方法は、ダブル・ダブル計算の制限された使用です。これはdouble変数hとlの頭/尾様式のペアとして-sqrt (0.5) * aの計算を含みます。製品の上位部分hはerfc()に渡され、一方、低位部分lは、hの相補誤差関数のローカルスロープに基づいて、結果を補間するために使用される。

erfc（x）の派生は-2 * exp（-x * x）/√πです。しかし、exp（-x * x）のかなり高価な計算を避けたい。したがって、漸近的には、 erfc（x * x + 4 /π）となるので、x> 0の場合、erfc（x）〜= 2 * exp（-x * x）/ | | | | h |に対して、erfc（h + 1）〜= erfc（h）-2 * h * l * erfc（x） 。H）、後者の用語の否定を容易lの計算に引き込むことができる1つの（）はIEEE-754 binary64を使用して倍精度のための次の実装に到達：使用

double my_normcdf (double a) 
{ 
    double h, l, r; 
    const double SQRT_HALF_HI = 0x1.6a09e667f3bcd0p-01; // 7.0710678118654757e-01 
    const double SQRT_HALF_LO = -0x1.bdd3413b264560p-55; // -4.8336466567264567e-17 

    /* clamp input as normcdf(x) is either 0 or 1 asymptotically */ 
    if (fabs (a) > 38.625) a = (a < 0.0) ? -38.625 : 38.625; 

    h = fma (-SQRT_HALF_HI, a, -SQRT_HALF_LO * a); 
    l = fma (SQRT_HALF_LO, a, fma (SQRT_HALF_HI, a, h)); 
    r = erfc (h); 
    if (h > 0.0) r = fma (2.0 * h * l, r, r); 
    return 0.5 * r; 
} 

最大観察誤差、前述の同じerfc()の実装は、1.96 ulpsです（IEEE-754 binary32を使用して）対応する単精度実装は、

0123です。

float my_normcdff (float a) 
{ 
    float h, l, r; 
    const float SQRT_HALF_HI = 0x1.6a09e6p-01f; // 7.07106769e-1 
    const float SQRT_HALF_LO = 0x1.9fcef4p-27f; // 1.21016175e-8 

    /* clamp input as normcdf(x) is either 0 or 1 asymptotically */ 
    if (fabsf (a) > 14.171875f) a = (a < 0.0f) ? -14.171875f : 14.171875f; 

    h = fmaf (-SQRT_HALF_HI, a, -SQRT_HALF_LO * a); 
    l = fmaf (SQRT_HALF_LO, a, fmaf (SQRT_HALF_HI, a, h)); 
    r = erfcf (h); 
    if (h > 0.0f) r = fmaf (2.0f * h * l, r, r); 
    return 0.5f * r; 
}