さまざまなアプローチを検討した結果、 （1 + 2 X）EXP（X ）ERFCの
MMシェパードとJG Laframboise、「チェビシェフ近似X 0≤Xで：アルゴリズムは、次の論文で提案されています<∞。計算の数学、36巻、第153号、1981年1月、頁249-253（online copy）
紙の基本的な考え方は、（1 + 2 X）EXPの近似を生成することです（ X ）をEXPと乗算することにより、我々は、単に（1 + 2 X）で割ることによってerfcx（X）を計算することができ、そこからERFC（X）、およびERFC（X）（ - x ）。おおよそ[1、1.3]の関数値とその一般的な「平坦度」を持つ関数の密接に境界がある範囲は、多項式近似によく役立ちます。このアプローチの数値特性をさらに近似間隔狭くすることにより改善される：元の引数Xを =（X- K）qで形質転換されたKが適切に一定に選択される/（X+ K） p（q）を計算する。ここで、pは多項式である。 ERFCので
からX = 2 - X ERFC 、我々は、この変換により、区間[-1,1]にマッピングされている区間[0、∞]を考慮する必要があります。 IEEE-754単精度の場合、x> 10.0546875の場合、erfcf()はゼロになりますので、x∈[0、10.0546875]のみを考慮する必要があります。この範囲のKの「最適な」値は何ですか？私は解を提供する数学的解析がないことを知っています。この論文では、実験に基づいてK = 3.75を示唆しています。
単精度計算、 Kが1.5と4の間で1/16のステップで変化するようなRemezアルゴリズムを用いてこのような近似を体系的に生成すると、最も近い近似誤差がK = （x -K）/（x + K）の非常に正確な計算に役立ち、this questionに示すように、K = 2が最も有利な選択である。
値K = 2とxの入力ドメインは、my answerのバリアント4を使用する必要があることを示唆していますが、安価なバリアント5がここで同じ正確さを達成したことを経験的に実証することができます。 q> -0.5の近似関数の非常に浅い勾配であり、これは引数qの誤差をおよそ10倍に減少させます。
erfc()の計算には、初期近似に加えて後処理ステップが必要であるため、十分に正確な最終結果を得るためには、これらの計算の両方の精度が高い必要があることは明らかです。エラー訂正技術を使用する必要があります。
一つはEXP（X 1 + 2 ）（X ）ERFC（X）の多項式近似で最も重要な係数は、フォーム（1 + S）であることが観察しますs < 0.5。これは、1を分割し、多項式で、のみを使用して、先頭の係数をより正確に表すことができることを意味します。従って、多項式p（q）を計算して逆数を掛ける代わりに、 = 1 /（1 + 2 x）と数学的には同等ですが、コア近似をp（q） +1、pを使用してfma (p, r, r)を計算します。Q *（1 + 2 X -
分割の精度は、残留E = P +1を計算し、相互Rから初期商Qを計算することによって向上させることができます）FMAの助けを借りて、その後、補正Q = Q +（E * R）を適用するEを使用し、再びFMAを使用して。 X慎重に行わなければならない -
累乗したがってEの計算は、誤差倍率特性を有します。 S 低い：ダブルfloat S 高いとしてX - FMAの利用可能性は自明の計算を可能にします。一つはE S 高いを計算することができるように、E Xは、それ自体の誘導体である：S 低い E Sとして高い + E S 高い * S 低い。この計算は、Rを生成する前の中間結果の乗算と組み合わせることができるR = R * E S 高い + R * E S 高い * S 低い 。 FMAを使用することにより、最も重要な用語r * e s highが可能な限り正確に計算されます。例外的な場合と負の引数を処理するためのいくつかの簡単な選択で上記の手順を組み合わせる
、1は、次のCコードに到着：
float my_expf (float);
/*
* Based on: M. M. Shepherd and J. G. Laframboise, "Chebyshev Approximation of
* (1+2x)exp(x^2)erfc x in 0 <= x < INF", Mathematics of Computation, Vol. 36,
* No. 153, January 1981, pp. 249-253.
*/
float my_erfcf (float x)
{
float a, d, e, m, p, q, r, s, t;
a = fabsf (x);
/* Compute q = (a-2)/(a+2) accurately. [0, 10.0546875] -> [-1, 0.66818] */
m = a - 2.0f;
p = a + 2.0f;
r = 1.0f/p;
q = m * r;
t = fmaf (q + 1.0f, -2.0f, a);
e = fmaf (q, -a, t);
q = fmaf (r, e, q);
/* Approximate (1+2*a)*exp(a*a)*erfc(a) as p(q)+1 for q in [-1, 0.66818] */
p = -0x1.a48024p-12f; // -4.01020574e-4
p = fmaf (p, q, -0x1.42a172p-10f); // -1.23073824e-3
p = fmaf (p, q, 0x1.585784p-10f); // 1.31355994e-3
p = fmaf (p, q, 0x1.1ade24p-07f); // 8.63243826e-3
p = fmaf (p, q, -0x1.081b72p-07f); // -8.05991236e-3
p = fmaf (p, q, -0x1.bc0b94p-05f); // -5.42047396e-2
p = fmaf (p, q, 0x1.4ffc40p-03f); // 1.64055347e-1
p = fmaf (p, q, -0x1.540840p-03f); // -1.66031361e-1
p = fmaf (p, q, -0x1.7bf612p-04f); // -9.27639678e-2
p = fmaf (p, q, 0x1.1ba03ap-02f); // 2.76978403e-1
/* Divide (1+p) by (1+2*a) ==> exp(a*a)*erfc(a) */
t = a + a;
d = t + 1.0f;
r = 1.0f/d;
q = fmaf (p, r, r); // q = (p+1)/(1+2*a)
e = (p - q) + fmaf (q, -t, 1.0f); // (p+1) - q*(1+2*a)
r = fmaf (e, r, q);
/* Multiply by exp(-a*a) ==> erfc(a) */
s = a * a;
e = my_expf (-s);
t = fmaf (a, -a, s);
r = fmaf (r, e, r * e * t);
/* Handle NaN arguments to erfc() */
if (!(a <= 0x1.fffffep127f)) r = x + x;
/* Clamp result for large arguments */
if (a > 10.0546875f) r = 0.0f;
/* Handle negative arguments to erfc() */
if (x < 0.0f) r = 2.0f - r;
return r;
}
/* Compute exponential base e. Maximum ulp error = 0.87161 */
float my_expf (float a)
{
float c, f, r;
int i;
// exp(a) = exp(i + f); i = rint (a/log(2))
c = 0x1.800000p+23f; // 1.25829120e+7
r = fmaf (0x1.715476p+0f, a, c) - c; // 1.44269502e+0
f = fmaf (r, -0x1.62e400p-01f, a); // -6.93145752e-1 // log_2_hi
f = fmaf (r, -0x1.7f7d1cp-20f, f); // -1.42860677e-6 // log_2_lo
i = (int)r;
// approximate r = exp(f) on interval [-log(2)/2,+log(2)/2]
r = 0x1.6a98dap-10f; // 1.38319808e-3
r = fmaf (r, f, 0x1.1272cap-07f); // 8.37550033e-3
r = fmaf (r, f, 0x1.555a20p-05f); // 4.16689515e-2
r = fmaf (r, f, 0x1.55542ep-03f); // 1.66664466e-1
r = fmaf (r, f, 0x1.fffff6p-02f); // 4.99999851e-1
r = fmaf (r, f, 0x1.000000p+00f); // 1.00000000e+0
r = fmaf (r, f, 0x1.000000p+00f); // 1.00000000e+0
// exp(a) = 2**i * exp(f);
r = ldexpf (r, i);
// handle special cases
if (!(fabsf (a) < 104.0f)) {
r = a + a; // handle NaNs
if (a < 0.0f) r = 0.0f;
if (a > 0.0f) r = 1e38f * 1e38f; // + INF
}
return r;
}
私は自分を隔離するために上記のコードでexpf()の私の独自の実装を使用しましたさまざまなコンピューティングプラットフォームで実装されたexpf()の違いから作業します。しかし、最大誤差が0.5 ulpに近いexpf()の実装はうまくいくはずです。上記のように、すなわち、my_expf()を使用する場合、my_erfcf()は、2.65712 ulpsの最大誤差を持ちます。ベクター化可能なexpf()が提供されていれば、上記のコードは問題なくベクトル化する必要があります。私はIntelコンパイラ13.1.3.198を使って素早くチェックしました。インテルコンパイラはループがベクトル化されていたことを報告し
/Qstd=c99 /O3 /QxCORE-AVX2 /fp:precise /Qfma /Qimf-precision:high:expf /Qvec_report=2
、：スイッチ私は、ループ内でmy_erfcf()に電話を入れ#include <mathimf.h>を追加し、その後、これらのコマンドラインを使用してコンパイルexpf()への呼び出し、とmy_expf()への呼び出しを置き換えますこれは分解されたバイナリコードの検査によって二重にチェックされます。
my_erfcf()は、完全な除算ではなく逆数のみを使用するため、ほぼ正確に丸められた結果が得られれば、高速な相互実装を使用することができます。ハードウェアで高速の単精度の逆数近似を提供するプロセッサの場合、これを3次収束のハレー反復と組み合わせることで簡単に達成できます。 x86プロセッサのこのアプローチの（スカラー）例は、次のとおりです。
/* Compute 1.0f/a almost correctly rounded. Halley iteration with cubic convergence */
float fast_recipf (float a)
{
__m128 t;
float e, r;
t = _mm_set_ss (a);
t = _mm_rcp_ss (t);
_mm_store_ss (&r, t);
e = fmaf (r, -a, 1.0f);
e = fmaf (e, e, e);
r = fmaf (e, r, r);
return r;
}