在研究了各种方法之后，似乎最合适的方法是以下论文中提出的算法：

MM Shepherd和JG Laframboise，“ 在 0≤x <∞中，（1 + 2 x ）exp（ x 2）erfc x的
切比雪夫逼近”。 计算数学
，第36卷，第153号，1981年1月，第249-253页（在线副本）

本文的基本思想是创建（1 + 2 x ）exp（ x 2）erfc（ x ）的近似值，从中我们可以通过简单地除以（1 + 2 x
）和erfc 来计算erfcx（ x ）（x）然后乘以exp（ -x
2）。函数的紧密边界范围（函数值大致在[1，1.3]中）及其一般的“平坦度”非常适合多项式逼近。通过缩小近似间隔，可以进一步改善此方法的数值属性：原始参数
x 通过 q =（ x -K）/（ x + K），其中K是适当选择的常数，然后计算 p （ q
），其中 p 是多项式。

由于erfc- x = 2-erfc x
，我们只需要考虑通过此变换映射到间隔[-1，1]的间隔[0，∞]。对于IEEE-754单精度，erfcf()消失（变为零）为 X >10.0546875，所以一个仅需要考虑_X_∈[0，10.0546875）。该范围内K的“最佳”值是多少？我不知道会提供答案的数学分析，该论文根据实验表明K = 3.75。

可以很容易地确定，对于单精度计算，对于那个大致附近的K的各个值，次数9的极小极大多项式近似就足够了。使用Remez算法系统地生成此类近似值，其中K在1.5和4之间以1/16的步长变化，对于K
= {2，2.625，3.3125}，观察到最低的近似误差。其中，K=2是最有利的选择，因为它可以非常精确地计算（ x -K）/（ x +K）

值K = 2和 x
的输入域将表明有必要使用我的答案中的变体4，但是一旦可以通过实验证明较便宜的变体5在此处达到相同的精度，这可能是由于非常浅 q>-0.5时近似函数的斜率，这会导致参数_q的_任何误差都减少大约十倍。

由于erfc()除了初始近似值外，还需要进行后处理步骤，因此很明显，为了获得足够准确的最终结果，这两个计算的精度都必须很高。必须使用纠错技术。

人们注意到，在（1 + 2 x ）exp（ x 2）erfc（ x ）的多项式近似中，最高有效系数的形式为（1 + s），其中 s
<0.5。这意味着我们可以通过除以1并仅在多项式中使用 s 来更准确地表示前导系数。因此，与其计算多项式p（q）而不是乘以倒数 r = 1 /（1
+ 2 x ），在数学上是等价的，但在数值上有利的是，将核近似值计算为p（ q ）+ 1并使用 p 计算fma (p, r, r)。

可以通过从倒数 r 计算初始商 q ，通过FMA 计算残差 e = p +1- q *（1 + 2 x ），然后使用 e
进行校正来提高除法的精度。 q = q +（ e * r ），再次使用FMA。

幂运算具有误差放大性质，因此必须小心执行 e - x 2的计算。FMA的可用性可让 -x 2计算为double float-s
high：s low。e x是其自身的导数，因此可以将e s high：s low计算为e s high + e s high * s
low。该计算可以与先前的中间结果 r 的相乘相结合，得出 r = r * es 高 + r * e s 高 * s
低。通过使用FMA，可以确保尽可能精确地计算出最高有效位 r * e s high。

将上面的步骤与几个简单的选择结合起来以处理异常情况和否定参数，可以得出以下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()不同计算平台上实现的差异区分开。但是，任何expf()最大误差接近0.5ulp的实现都应该可以正常工作。如上所示，即使用时my_expf()，my_erfcf()最大误差为2.65712ulps。如果提供了vectorizable的可用性expf()，则上面的代码应该可以正确地进行矢量化。我对Intel编译器13.1.3.198进行了快速检查。我将调用my_erfcf()循环放置，添加#include<mathimf.h>，将调用替换my_expf()为expf()，然后使用以下命令行开关进行编译：

/Qstd=c99 /O3 /QxCORE-AVX2 /fp:precise /Qfma /Qimf-precision:high:expf /Qvec_report=2

英特尔编译器报告说，循环已被矢量化，我通过检查反汇编的二进制代码对其进行了仔细检查。

由于my_erfcf()仅使用倒数而不是完整的除法，因此，只要它们提供几乎正确的舍入结果，就可以使用快速倒数实现。对于在硬件中提供快速单精度倒数逼近的处理器，可以通过将其与具有三次收敛性的Halley迭代耦合来轻松实现。针对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;
}