C语言 5位精度互补误差函数的快速实现

t1qtbnec  于 2024-01-06  发布在  其他
关注(0)|答案(2)|浏览(183)

[虽然这是一个自我回答的问题，但我很乐意投票并接受任何替代答案，无论是在相同的计算工作量下提供上级准确度，还是在保持相同准确度的情况下减少计算工作量。]
我以前有demonstrated如何计算互补误差函数erfcf()，其最大误差小于3个ulps。这可以作为其他函数的构建块，例如标准正态分布的CDF Φ（x）= 1/2 erfc（-λ 1/2 x）或高斯Q函数，Q（x）= 1-Φ（x）= 1/2 erfc（λ 1/2 x）。然而，对于某些使用情况，不需要完全精确到单精度的计算，而erfc()评估对总运行时间的贡献不可忽略。
文献提供了对互补误差函数的各种低精度近似，但它们要么限于完整输入域的子集，要么针对绝对误差进行了优化，要么计算过于复杂，例如，需要多次调用超越函数。如何才能实现erfcf()具有高性能和大约5位十进制数的相对精度 * 在整个输入域中 *？

c

2条答案

按热度按时间
6yt4nkrj

6yt4nkrj1#

以下假设平台符合IEEE-754（2008）浮点标准，在该平台上float被Map到IEEE-754 binary32，并且在32位整数和float之间使用相同的字节序。（如果需要，通过设置适当的命令行开关）保留IEEE-754语义。我使用的是带有开关-march=skylake-avx152 -O3 -fp-model=precise的Intel C/C++编译器。
由于互补误差函数关于（0，1）对称，因此可以关注正半平面中的函数输入。这里函数大致像exp（-x2）一样衰减，并且对于x > 10.5的参数，float计算下溢到零。如果将erfc（x）/ exp（-x2）绘制在[0，10.5]形状表明，用多项式近似有点困难，但应该很容易用有理函数近似，即，两个多项式的比率。一些初始实验表明，两个多项式的次数各为3应该足以达到五位数的精度。
虽然有许多工具可以生成多项式近似，但不幸的是，有理近似并非如此。我使用Remez算法的修改来生成初始极大极小近似R（x）= P（x）/Q（x）到erfc（x）/ exp（-x2），但必须进行相当广泛的启发式搜索，以获得一个近似值，该近似值提供 * 接近 * 误差峰值的等振荡，实现了10-5的相对误差，其余的差异对于我的需要来说可以忽略不计。
通过计算erfc（x）= exp（-x2）R（x），所获得的精度显然取决于平台expf()实现的精度。此功能的忠实全面的实现（最大误差<= 1ulp）是常见的。虽然英特尔编译器附带了一个高度精确的数学库，可以提供近乎正确的四舍五入实现（最大误差非常接近0.5 ulps），我也尝试了我自己的忠实四舍五入的替代my_expf()，误差更大，只观察到对fast_erfcf()精度的影响非常小。

  1. #include <stdlib.h>
  2. #include <stdio.h>
  3. #include <stdint.h>
  4. #include <string.h>
  5. #include <math.h>
  7. #define USE_FMA          (1)
  8. #define USE_BUILTIN_EXP  (0)
  10. #if !USE_BUILTIN_EXP
  11. float my_expf (float a);
  12. #endif // USE_BUILTIN_EXP
  14. /* Fast computation of the complementary error function. For argument x > 0
  15.    erfc(x) = exp(-x*x) * P(x) / Q(x), where P(x) and Q(x) are polynomials. 
  16.    If expf() is faithfully rounded, the following error bounds should hold:
  17.    Maximum relative error < 1.065e-5, maximum absolute error < 9.50e-6, and 
  18.    maximum ulp error < 176.5
  19. */
  20. float fast_erfcf (float x)
  21. {
  22.     float a, c, e, p, q, r, s;
  23.     a = fabsf (x);
  24.     c = fminf (a, 10.5f);
  25.     s = -c * c;
  26. #if USE_BUILTIN_EXP
  27.     e = expf (s);
  28. #else // USE_BUILTIN_EXP
  29.     e = my_expf (s);
  30. #endif // USE_BUILTIN_EXP
  31. #if USE_FMA
  32.     q =             3.82346243e-1f;  //  0x1.8785c6p-2
  33.     p =            -4.38094139e-5f;  // -0x1.6f8000p-15
  34.     q = fmaf (q, c, 1.30382288e+0f); //  0x1.4dc756p+0
  35.     p = fmaf (p, c, 2.16852024e-1f); //  0x1.bc1ceap-3
  36.     q = fmaf (q, c, 1.85278833e+0f); //  0x1.da5056p+0
  37.     p = fmaf (p, c, 7.23953605e-1f); //  0x1.72aa0cp-1
  38.     q = fmaf (q, c, 9.99991655e-1f); //  0x1.fffee8p-1
  39.     p = fmaf (p, c, 1.00000000e+0f); //  0x1.000000p+0
  40. #else // USE_FMA
  41.     q =         3.82346272e-1f; //  0x1.8785c8p-2f
  42.     p =        -4.38243151e-5f; // -0x1.6fa000p-15
  43.     q = q * c + 1.30382371e+0f; //  0x1.4dc764p+0
  44.     p = p * c + 2.16852218e-1f; //  0x1.bc1d04p-3
  45.     q = q * c + 1.85278797e+0f; //  0x1.da5050p+0
  46.     p = p * c + 7.23953605e-1f; //  0x1.72aa0cp-1
  47.     q = q * c + 9.99991596e-1f; //  0x1.fffee6p-1
  48.     p = p * c + 1.00000000e+0f; //  0x1.000000p+0
  49. #endif // USE_FMA
  50.     r = e / q;
  51.     r = r * p;
  52.     if (x < 0.0f) r = 2.0f - r;
  53.     if (isnan(x)) r = x + x;
  54.     return r;
  55. }
  57. float uint32_as_float (uint32_t a)
  58. {
  59.     float r;
  60.     memcpy (&r, &a, sizeof r);
  61.     return r;
  62. }
  64. /* Exponential function. Maximum error 0.86565 ulps */
  65. float my_expf (float a)
  66. {
  67.     float f, r, j, s, t;
  68.     int i;
  69.     unsigned int ia;
  71.     // exp(a) = 2**i * exp(f); i = rintf (a / log(2))
  72.     j = fmaf (1.442695f, a, 12582912.f); // 0x1.715476p0 // log2(e)
  73.     j = j - 12582912.f; // 0x1.8p23 // 2**23+2**22
  74.     f = fmaf (j, -6.93145752e-1f, a); // -0x1.62e400p-1  // log_2_hi 
  75.     f = fmaf (j, -1.42860677e-6f, f); // -0x1.7f7d1cp-20 // log_2_lo 
  76.     i = (int)j;
  77.     // approximate r = exp(f) on interval [-log(2)/2, +log(2)/2]
  78.     r =             1.37805939e-3f;  // 0x1.694000p-10
  79.     r = fmaf (r, f, 8.37312452e-3f); // 0x1.125edcp-7
  80.     r = fmaf (r, f, 4.16695364e-2f); // 0x1.555b5ap-5
  81.     r = fmaf (r, f, 1.66664720e-1f); // 0x1.555450p-3
  82.     r = fmaf (r, f, 4.99999851e-1f); // 0x1.fffff6p-2
  83.     r = fmaf (r, f, 1.00000000e+0f); // 0x1.000000p+0
  84.     r = fmaf (r, f, 1.00000000e+0f); // 0x1.000000p+0
  85.     // exp(a) = 2**i * r
  86.     ia = (i > 0) ? 0 : 0x83000000;
  87.     s = uint32_as_float (0x7f000000 + ia);
  88.     t = uint32_as_float ((i << 23) - ia);
  89.     r = r * s;
  90.     r = r * t;
  91.     // handle special cases: severe overflow / underflow
  92.     if (fabsf (a) >= 104.0f) r = (a < 0) ? 0.0f : INFINITY;
  93.     return r;
  94. }
  96. uint32_t float_as_uint32 (float a)
  97. {
  98.     uint32_t r;
  99.     memcpy (&r, &a, sizeof r);
  100.     return r;
  101. }
  103. uint64_t double_as_uint64 (double a)
  104. {
  105.     uint64_t r;
  106.     memcpy (&r, &a, sizeof r);
  107.     return r;
  108. }
  110. double floatUlpErr (float res, double ref)
  111. {
  112.     uint64_t refi, i, j, err;
  113.     int expoRef;
  114.     
  115.     /* ulp error cannot be computed if either operand is NaN, infinity, zero */
  116.     if (isnan (res) || isnan (ref) || isinf(res) || isinf (ref) ||
  117.         (res == 0.0f) || (ref == 0.0)) {
  118.         return 0.0;
  119.     }
  120.     i = ((int64_t)float_as_uint32 (res)) << 32;
  121.     expoRef = (int)(((double_as_uint64 (ref) >> 52) & 0x7ff) - 1023);
  122.     refi = double_as_uint64 (ref);
  123.     if (expoRef >= 129) {
  124.         j = (refi & 0x8000000000000000ULL) | 0x7fffffffffffffffULL;
  125.     } else if (expoRef < -126) {
  126.         j = ((refi << 11) | 0x8000000000000000ULL) >> 8;
  127.         j = j >> (-(expoRef + 126));
  128.         j = j | (refi & 0x8000000000000000ULL);
  129.     } else {
  130.         j = ((refi << 11) & 0x7fffffffffffffffULL) >> 8;
  131.         j = j | ((uint64_t)(expoRef + 127) << 55);
  132.         j = j | (refi & 0x8000000000000000ULL);
  133.     }
  134.     err = (i < j) ? (j - i) : (i - j);
  135.     return err / 4294967296.0;
  136. }
  138. int main (void)
  139. {
  140.     uint32_t argi, resi, refi, diff;
  141.     float arg, res, reff, abserrloc = NAN, relerrloc = NAN, ulperrloc = NAN;
  142.     double ref, relerr, abserr, ulperr;
  143.     double maxabserr = 0, maxrelerr = 0, maxulperr = 0;
  145.     argi = 0;
  146.     do {
  147.         arg = uint32_as_float (argi);
  148.         ref = erfc ((double)arg);
  149.         res = fast_erfcf (arg);
  150.         reff = (float)ref;
  151.         resi = float_as_uint32 (res);
  152.         refi = float_as_uint32 (reff);
  153.         ulperr = floatUlpErr (res, ref);
  154.         if (ulperr > maxulperr) {
  155.             ulperrloc = arg;
  156.             maxulperr = ulperr;
  157.         }
  158.         abserr = fabs (res - ref);
  159.         if (abserr > maxabserr) {
  160.             abserrloc = arg;
  161.             maxabserr = abserr;
  162.         }
  163.         if (fabs (ref) >= 0x1.0p-126) {
  164.             relerr = fabs ((res - ref) / ref);
  165.             if (relerr > maxrelerr) {
  166.                 relerrloc = arg;
  167.                 maxrelerr = relerr;
  168.             }
  169.         }
  170.         diff = (resi > refi) ? (resi - refi) : (refi - resi);
  171.         if (diff > 200) {
  172.             printf ("diff=%u @ %15.8e : res=% 15.8e  ref=% 15.8e\n", 
  173.                     diff, arg, res, ref);
  174.             return EXIT_FAILURE;
  175.         }
  176.         argi++;
  177.     } while (argi);
  179.     printf ("max rel err = %.6e @ % 15.8e\n"
  180.             "max abs err = %.6e @ % 15.8e\n"
  181.             "max ulp err = %.6e @ % 15.8e\n",
  182.             maxrelerr, relerrloc, 
  183.             maxabserr, abserrloc,
  184.             maxulperr, ulperrloc);
  185.     return EXIT_SUCCESS;
  186. }

字符串

展开查看全部
赞(0）分享  回复(0）举报  2024-01-06
exdqitrt

exdqitrt2#

OP提供了用于评估性能测试的范围的注解：1 ULP步骤中的0.0至10.5。
由于所有float [-1.6e-08...+1.6e-08]中约有40%返回值约为1.0，代码可以使用以下内容。对我来说，整体性能至少 * 翻了一番 *。

  1. float fast_erfcf(float x) {
  2.   if (fabsf(x) <= 1.60e-8) {  // or maybe a bit more.
  3.     return 1.0;
  4.   }
  5.   ...
  6. }

字符串
由于erfc(x)中float的大范围可以返回1.0，因此我建议进行这种预先测试。
优点包括：

  • 对于所有float的40%来说，速度要快得多。
  • erfc(x_near_0)预期返回1.0。OP的fast_erfcf(0.0)返回1.00000834。考虑到可以容忍的大ULP，可以接受。1.0仍然很好。
  • OP的fast_erfcf(small negative value)返回0.999991536。可以接受，因为可以容忍较大的ULP，但在可能的情况下，希望有erfc(x) >= erfc(next x towards + infinity)。

此外：建议在[0.1.... 10.0]范围内对最感兴趣的值进行性能测试。
这个建议对OP没有多大帮助，所以把它移到这个答案的底部。
与其用my_expf(float a)执行 ex，不如用my_exp2f(float a)执行 2x（这会影响有理函数的常数）。
j = fmaf(1.442695f, a, 12582912.f);被删除。
使用base 2比base e 执行取幂稍微容易一些。
我估计代码会快几个百分点（~3%），而ULP结果不会有太大变化。
确实，这是一个微优化，但它实际上是免费的，因为它不会给链接答案的方法增加任何真实的复杂性。

展开查看全部
赞(0）分享  回复(0）举报  2024-01-06
我来回答

相关问题

    查看更多

    热门标签

    更多
    JavaquerypythonNode开发语言requestUtil数据库Table后端算法LoggerMessageElementParser

    最新问答

    更多
    • 蜀ICP备13028337号-1 大数据知识库 https://www.saoniuhuo.com © All rights reserved
    • 本站内容来源互联网,如果侵犯您的权益请联系我们删除, 联系方式：448109455@qq.com