缩放互补误差函数的精确计算，erfcx()
Accurate computation of scaled complementary error function, erfcx()
通常由 erfcx 指定的（指数）缩放互补误差函数在数学上定义为 erfcx(x) := ex2erfc(x)。它经常出现在物理和化学的扩散问题中。虽然某些数学环境（例如 MATLAB and GNU Octave）提供此函数，但 C 标准数学库中没有它，它仅提供 erf() 和 erfc().
虽然可以直接根据数学定义实现自己的erfcx()，但这只适用于有限的输入域，因为在正半平面erfc()下溢参数中等幅度，而 exp() 溢出，例如 this question 中所述。
为了与 C 一起使用，可以采用一些 erfcx() 开源实现，例如 Faadeeva package, as pointed to in responses to this question 中的实现。然而，这些实现通常不会为给定的浮点格式提供完全的准确性。例如，使用 232 测试向量的测试显示 Faadeeva 包提供的最大误差 erfcx() 在正半平面中为 8.41 ulp，在正半平面中为 511.68 ulp负半平面。
准确实现的合理界限是 4 ulps，对应于 Intel's Vector Math 库的 LA 配置文件中数学函数的准确界限，我发现这是非需要良好准确性和良好性能的简单数学函数实现。
erfcx() 和相应的单精度版本 erfcxf() 如何在仅使用 C 标准数学库而不需要外部库的情况下准确实现？我们可以假设 C 的 float nad double 类型映射到 IEEE 754-2008 binary32 和 binary64 浮点类型。可以假定硬件支持融合乘加运算 (FMA)，因为目前所有主要处理器体系结构都支持它。
目前我发现的 erfcx() 实施的最佳方法基于以下论文：
米。 M. Shepherd 和 J. G. Laframboise，"Chebyshev Approximation of (1 + 2 x) exp(x2) erfc x in 0 ≤ x < ∞." 计算数学，第 36 卷，第 153 期，1981 年 1 月，第 249-253 页 (online)
这篇论文提出了巧妙的变换，将缩放后的互补误差函数映射到适用于直接多项式逼近的紧界辅助函数。为了性能，我尝试了各种变换，但所有这些都对准确性产生了负面影响。变换 (x - K) / (x + K) 中常数 K 的选择与核心近似的精度具有不明显的关系。我根据经验确定了 "optimal" 值，这与论文不同。
将参数转换为核心近似值并将中间结果返回到 erfcx 结果会产生额外的舍入误差。为了减轻它们对准确性的影响，我们需要应用补偿步骤，我在之前的 中对此进行了一些详细概述。 FMA 的可用性大大简化了这项任务。
生成的单精度代码如下所示：
/*
* 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_erfcxf (float x)
{
float a, d, e, m, p, q, r, s, t;
a = fmaxf (x, 0.0f - x); // NaN-preserving absolute value computation
/* Compute q = (a-2)/(a+2) accurately. [0,INF) -> [-1,1] */
m = a - 2.0f;
p = a + 2.0f;
#if FAST_RCP_SSE
r = fast_recipf_sse (p);
#else
r = 1.0f / p;
#endif
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,1] */
p = 0x1.f10000p-15f; // 5.92470169e-5
p = fmaf (p, q, 0x1.521cc6p-13f); // 1.61224554e-4
p = fmaf (p, q, -0x1.6b4ffep-12f); // -3.46481771e-4
p = fmaf (p, q, -0x1.6e2a7cp-10f); // -1.39681227e-3
p = fmaf (p, q, 0x1.3c1d7ep-10f); // 1.20588380e-3
p = fmaf (p, q, 0x1.1cc236p-07f); // 8.69014394e-3
p = fmaf (p, q, -0x1.069940p-07f); // -8.01387429e-3
p = fmaf (p, q, -0x1.bc1b6cp-05f); // -5.42122945e-2
p = fmaf (p, q, 0x1.4ff8acp-03f); // 1.64048523e-1
p = fmaf (p, q, -0x1.54081ap-03f); // -1.66031078e-1
p = fmaf (p, q, -0x1.7bf5cep-04f); // -9.27637145e-2
p = fmaf (p, q, 0x1.1ba03ap-02f); // 2.76978403e-1
/* Divide (1+p) by (1+2*a) ==> exp(a*a)*erfc(a) */
d = a + 0.5f;
#if FAST_RCP_SSE
r = fast_recipf_sse (d);
#else
r = 1.0f / d;
#endif
r = r * 0.5f;
q = fmaf (p, r, r); // q = (p+1)/(1+2*a)
t = q + q;
e = (p - q) + fmaf (t, -a, 1.0f); // residual: (p+1)-q*(1+2*a)
r = fmaf (e, r, q);
if (a > 0x1.fffffep127f) r = 0.0f; // 3.40282347e+38 // handle INF argument
/* Handle negative arguments: erfcx(x) = 2*exp(x*x) - erfcx(|x|) */
if (x < 0.0f) {
s = x * x;
d = fmaf (x, x, -s);
e = expf (s);
r = e - r;
r = fmaf (e, d + d, r);
r = r + e;
if (e > 0x1.fffffep127f) r = e; // 3.40282347e+38 // avoid creating NaN
}
return r;
}
此实现在负半平面中的最大误差将取决于标准数学库实现 expf() 的准确性。使用 Intel 编译器 13.1.3.198 版并使用 /fp:strict 进行编译，我在详尽测试中观察到正半平面的最大误差为 2.00450 ulp，负半平面的最大误差为 2.38412 ulp。目前我能说的最好的是，expf() 的忠实圆整实现将导致最大误差小于 2.5 ulps。
请注意，虽然代码包含两个除法，这可能是缓慢的操作，但它们以倒数的特殊形式出现，因此适合在许多平台上使用快速倒数近似。根据实验，只要忠实地舍入倒数近似值，对 erfcxf() 精度的影响似乎可以忽略不计。即使是稍大的错误，例如在快速 SSE 版本中（最大错误 < 2.0 ulps）似乎也只有很小的影响。
/* Fast reciprocal approximation. HW approximation plus Newton iteration */
float fast_recipf_sse (float a)
{
__m128 t;
float e, r;
t = _mm_set_ss (a);
t = _mm_rcp_ss (t);
_mm_store_ss (&r, t);
e = fmaf (0.0f - a, r, 1.0f);
r = fmaf (e, r, r);
return r;
}
双精度版本 erfcx() 在结构上与单精度版本 erfcxf() 相同，但需要具有更多项的极小极大多项式逼近。这在优化核心近似时提出了挑战，因为当搜索 space 非常大时，许多启发式方法都会失效。下面的系数代表了我迄今为止最好的解决方案，而且肯定还有改进的余地。使用英特尔编译器和 /fp:strict 构建，并使用 232 个随机测试向量，观察到的最大误差在正半平面为 2.83788 ulp，在负半平面为 2.77856 ulp -飞机。
double my_erfcx (double x)
{
double a, d, e, m, p, q, r, s, t;
a = fmax (x, 0.0 - x); // NaN preserving absolute value computation
/* Compute q = (a-4)/(a+4) accurately. [0,INF) -> [-1,1] */
m = a - 4.0;
p = a + 4.0;
r = 1.0 / p;
q = m * r;
t = fma (q + 1.0, -4.0, a);
e = fma (q, -a, t);
q = fma (r, e, q);
/* Approximate (1+2*a)*exp(a*a)*erfc(a) as p(q)+1 for q in [-1,1] */
p = 0x1.edcad78fc8044p-31; // 8.9820305531190140e-10
p = fma (p, q, 0x1.b1548f14735d1p-30); // 1.5764464777959401e-09
p = fma (p, q, -0x1.a1ad2e6c4a7a8p-27); // -1.2155985739342269e-08
p = fma (p, q, -0x1.1985b48f08574p-26); // -1.6386753783877791e-08
p = fma (p, q, 0x1.c6a8093ac4f83p-24); // 1.0585794011876720e-07
p = fma (p, q, 0x1.31c2b2b44b731p-24); // 7.1190423171700940e-08
p = fma (p, q, -0x1.b87373facb29fp-21); // -8.2040389712752056e-07
p = fma (p, q, 0x1.3fef1358803b7p-22); // 2.9796165315625938e-07
p = fma (p, q, 0x1.7eec072bb0be3p-18); // 5.7059822144459833e-06
p = fma (p, q, -0x1.78a680a741c4ap-17); // -1.1225056665965572e-05
p = fma (p, q, -0x1.9951f39295cf4p-16); // -2.4397380523258482e-05
p = fma (p, q, 0x1.3be1255ce180bp-13); // 1.5062307184282616e-04
p = fma (p, q, -0x1.a1df71176b791p-13); // -1.9925728768782324e-04
p = fma (p, q, -0x1.8d4aaa0099bc8p-11); // -7.5777369791018515e-04
p = fma (p, q, 0x1.49c673066c831p-8); // 5.0319701025945277e-03
p = fma (p, q, -0x1.0962386ea02b7p-6); // -1.6197733983519948e-02
p = fma (p, q, 0x1.3079edf465cc3p-5); // 3.7167515521269866e-02
p = fma (p, q, -0x1.0fb06dfedc4ccp-4); // -6.6330365820039094e-02
p = fma (p, q, 0x1.7fee004e266dfp-4); // 9.3732834999538536e-02
p = fma (p, q, -0x1.9ddb23c3e14d2p-4); // -1.0103906603588378e-01
p = fma (p, q, 0x1.16ecefcfa4865p-4); // 6.8097054254651804e-02
p = fma (p, q, 0x1.f7f5df66fc349p-7); // 1.5379652102610957e-02
p = fma (p, q, -0x1.1df1ad154a27fp-3); // -1.3962111684056208e-01
p = fma (p, q, 0x1.dd2c8b74febf6p-3); // 2.3299511862555250e-01
/* Divide (1+p) by (1+2*a) ==> exp(a*a)*erfc(a) */
d = a + 0.5;
r = 1.0 / d;
r = r * 0.5;
q = fma (p, r, r); // q = (p+1)/(1+2*a)
t = q + q;
e = (p - q) + fma (t, -a, 1.0); // residual: (p+1)-q*(1+2*a)
r = fma (e, r, q);
/* Handle argument of infinity */
if (a > 0x1.fffffffffffffp1023) r = 0.0;
/* Handle negative arguments: erfcx(x) = 2*exp(x*x) - erfcx(|x|) */
if (x < 0.0) {
s = x * x;
d = fma (x, x, -s);
e = exp (s);
r = e - r;
r = fma (e, d + d, r);
r = r + e;
if (e > 0x1.fffffffffffffp1023) r = e; // avoid creating NaN
}
return r;
}
通常由 erfcx 指定的（指数）缩放互补误差函数在数学上定义为 erfcx(x) := ex2erfc(x)。它经常出现在物理和化学的扩散问题中。虽然某些数学环境（例如 MATLAB and GNU Octave）提供此函数，但 C 标准数学库中没有它，它仅提供 erf() 和 erfc().
虽然可以直接根据数学定义实现自己的erfcx()，但这只适用于有限的输入域，因为在正半平面erfc()下溢参数中等幅度，而 exp() 溢出，例如 this question 中所述。
为了与 C 一起使用，可以采用一些 erfcx() 开源实现，例如 Faadeeva package, as pointed to in responses to this question 中的实现。然而，这些实现通常不会为给定的浮点格式提供完全的准确性。例如，使用 232 测试向量的测试显示 Faadeeva 包提供的最大误差 erfcx() 在正半平面中为 8.41 ulp，在正半平面中为 511.68 ulp负半平面。
准确实现的合理界限是 4 ulps，对应于 Intel's Vector Math 库的 LA 配置文件中数学函数的准确界限，我发现这是非需要良好准确性和良好性能的简单数学函数实现。
erfcx() 和相应的单精度版本 erfcxf() 如何在仅使用 C 标准数学库而不需要外部库的情况下准确实现？我们可以假设 C 的 float nad double 类型映射到 IEEE 754-2008 binary32 和 binary64 浮点类型。可以假定硬件支持融合乘加运算 (FMA)，因为目前所有主要处理器体系结构都支持它。
目前我发现的 erfcx() 实施的最佳方法基于以下论文：
米。 M. Shepherd 和 J. G. Laframboise，"Chebyshev Approximation of (1 + 2 x) exp(x2) erfc x in 0 ≤ x < ∞." 计算数学，第 36 卷，第 153 期，1981 年 1 月，第 249-253 页 (online)
这篇论文提出了巧妙的变换，将缩放后的互补误差函数映射到适用于直接多项式逼近的紧界辅助函数。为了性能，我尝试了各种变换，但所有这些都对准确性产生了负面影响。变换 (x - K) / (x + K) 中常数 K 的选择与核心近似的精度具有不明显的关系。我根据经验确定了 "optimal" 值，这与论文不同。
将参数转换为核心近似值并将中间结果返回到 erfcx 结果会产生额外的舍入误差。为了减轻它们对准确性的影响，我们需要应用补偿步骤，我在之前的
生成的单精度代码如下所示：
/*
* 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_erfcxf (float x)
{
float a, d, e, m, p, q, r, s, t;
a = fmaxf (x, 0.0f - x); // NaN-preserving absolute value computation
/* Compute q = (a-2)/(a+2) accurately. [0,INF) -> [-1,1] */
m = a - 2.0f;
p = a + 2.0f;
#if FAST_RCP_SSE
r = fast_recipf_sse (p);
#else
r = 1.0f / p;
#endif
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,1] */
p = 0x1.f10000p-15f; // 5.92470169e-5
p = fmaf (p, q, 0x1.521cc6p-13f); // 1.61224554e-4
p = fmaf (p, q, -0x1.6b4ffep-12f); // -3.46481771e-4
p = fmaf (p, q, -0x1.6e2a7cp-10f); // -1.39681227e-3
p = fmaf (p, q, 0x1.3c1d7ep-10f); // 1.20588380e-3
p = fmaf (p, q, 0x1.1cc236p-07f); // 8.69014394e-3
p = fmaf (p, q, -0x1.069940p-07f); // -8.01387429e-3
p = fmaf (p, q, -0x1.bc1b6cp-05f); // -5.42122945e-2
p = fmaf (p, q, 0x1.4ff8acp-03f); // 1.64048523e-1
p = fmaf (p, q, -0x1.54081ap-03f); // -1.66031078e-1
p = fmaf (p, q, -0x1.7bf5cep-04f); // -9.27637145e-2
p = fmaf (p, q, 0x1.1ba03ap-02f); // 2.76978403e-1
/* Divide (1+p) by (1+2*a) ==> exp(a*a)*erfc(a) */
d = a + 0.5f;
#if FAST_RCP_SSE
r = fast_recipf_sse (d);
#else
r = 1.0f / d;
#endif
r = r * 0.5f;
q = fmaf (p, r, r); // q = (p+1)/(1+2*a)
t = q + q;
e = (p - q) + fmaf (t, -a, 1.0f); // residual: (p+1)-q*(1+2*a)
r = fmaf (e, r, q);
if (a > 0x1.fffffep127f) r = 0.0f; // 3.40282347e+38 // handle INF argument
/* Handle negative arguments: erfcx(x) = 2*exp(x*x) - erfcx(|x|) */
if (x < 0.0f) {
s = x * x;
d = fmaf (x, x, -s);
e = expf (s);
r = e - r;
r = fmaf (e, d + d, r);
r = r + e;
if (e > 0x1.fffffep127f) r = e; // 3.40282347e+38 // avoid creating NaN
}
return r;
}
此实现在负半平面中的最大误差将取决于标准数学库实现 expf() 的准确性。使用 Intel 编译器 13.1.3.198 版并使用 /fp:strict 进行编译，我在详尽测试中观察到正半平面的最大误差为 2.00450 ulp，负半平面的最大误差为 2.38412 ulp。目前我能说的最好的是，expf() 的忠实圆整实现将导致最大误差小于 2.5 ulps。
请注意，虽然代码包含两个除法，这可能是缓慢的操作，但它们以倒数的特殊形式出现，因此适合在许多平台上使用快速倒数近似。根据实验，只要忠实地舍入倒数近似值，对 erfcxf() 精度的影响似乎可以忽略不计。即使是稍大的错误，例如在快速 SSE 版本中（最大错误 < 2.0 ulps）似乎也只有很小的影响。
/* Fast reciprocal approximation. HW approximation plus Newton iteration */
float fast_recipf_sse (float a)
{
__m128 t;
float e, r;
t = _mm_set_ss (a);
t = _mm_rcp_ss (t);
_mm_store_ss (&r, t);
e = fmaf (0.0f - a, r, 1.0f);
r = fmaf (e, r, r);
return r;
}
双精度版本 erfcx() 在结构上与单精度版本 erfcxf() 相同，但需要具有更多项的极小极大多项式逼近。这在优化核心近似时提出了挑战，因为当搜索 space 非常大时，许多启发式方法都会失效。下面的系数代表了我迄今为止最好的解决方案，而且肯定还有改进的余地。使用英特尔编译器和 /fp:strict 构建，并使用 232 个随机测试向量，观察到的最大误差在正半平面为 2.83788 ulp，在负半平面为 2.77856 ulp -飞机。
double my_erfcx (double x)
{
double a, d, e, m, p, q, r, s, t;
a = fmax (x, 0.0 - x); // NaN preserving absolute value computation
/* Compute q = (a-4)/(a+4) accurately. [0,INF) -> [-1,1] */
m = a - 4.0;
p = a + 4.0;
r = 1.0 / p;
q = m * r;
t = fma (q + 1.0, -4.0, a);
e = fma (q, -a, t);
q = fma (r, e, q);
/* Approximate (1+2*a)*exp(a*a)*erfc(a) as p(q)+1 for q in [-1,1] */
p = 0x1.edcad78fc8044p-31; // 8.9820305531190140e-10
p = fma (p, q, 0x1.b1548f14735d1p-30); // 1.5764464777959401e-09
p = fma (p, q, -0x1.a1ad2e6c4a7a8p-27); // -1.2155985739342269e-08
p = fma (p, q, -0x1.1985b48f08574p-26); // -1.6386753783877791e-08
p = fma (p, q, 0x1.c6a8093ac4f83p-24); // 1.0585794011876720e-07
p = fma (p, q, 0x1.31c2b2b44b731p-24); // 7.1190423171700940e-08
p = fma (p, q, -0x1.b87373facb29fp-21); // -8.2040389712752056e-07
p = fma (p, q, 0x1.3fef1358803b7p-22); // 2.9796165315625938e-07
p = fma (p, q, 0x1.7eec072bb0be3p-18); // 5.7059822144459833e-06
p = fma (p, q, -0x1.78a680a741c4ap-17); // -1.1225056665965572e-05
p = fma (p, q, -0x1.9951f39295cf4p-16); // -2.4397380523258482e-05
p = fma (p, q, 0x1.3be1255ce180bp-13); // 1.5062307184282616e-04
p = fma (p, q, -0x1.a1df71176b791p-13); // -1.9925728768782324e-04
p = fma (p, q, -0x1.8d4aaa0099bc8p-11); // -7.5777369791018515e-04
p = fma (p, q, 0x1.49c673066c831p-8); // 5.0319701025945277e-03
p = fma (p, q, -0x1.0962386ea02b7p-6); // -1.6197733983519948e-02
p = fma (p, q, 0x1.3079edf465cc3p-5); // 3.7167515521269866e-02
p = fma (p, q, -0x1.0fb06dfedc4ccp-4); // -6.6330365820039094e-02
p = fma (p, q, 0x1.7fee004e266dfp-4); // 9.3732834999538536e-02
p = fma (p, q, -0x1.9ddb23c3e14d2p-4); // -1.0103906603588378e-01
p = fma (p, q, 0x1.16ecefcfa4865p-4); // 6.8097054254651804e-02
p = fma (p, q, 0x1.f7f5df66fc349p-7); // 1.5379652102610957e-02
p = fma (p, q, -0x1.1df1ad154a27fp-3); // -1.3962111684056208e-01
p = fma (p, q, 0x1.dd2c8b74febf6p-3); // 2.3299511862555250e-01
/* Divide (1+p) by (1+2*a) ==> exp(a*a)*erfc(a) */
d = a + 0.5;
r = 1.0 / d;
r = r * 0.5;
q = fma (p, r, r); // q = (p+1)/(1+2*a)
t = q + q;
e = (p - q) + fma (t, -a, 1.0); // residual: (p+1)-q*(1+2*a)
r = fma (e, r, q);
/* Handle argument of infinity */
if (a > 0x1.fffffffffffffp1023) r = 0.0;
/* Handle negative arguments: erfcx(x) = 2*exp(x*x) - erfcx(|x|) */
if (x < 0.0) {
s = x * x;
d = fma (x, x, -s);
e = exp (s);
r = e - r;
r = fma (e, d + d, r);
r = r + e;
if (e > 0x1.fffffffffffffp1023) r = e; // avoid creating NaN
}
return r;
}