c/c++开发分享精确计算缩放互补误差函数,erfcx()

通常由erfcx指定的(指数)缩放互补误差函数在数学上定义为erfcx(x):= e x 2 erfc(x)。 它经常发生在物理学和化学中的扩散问题中。 虽然一些数学环境(如MATLAB和GNU Octave )提供此function,但它不存在于C标准数学库中,它只提供erf()erfc()

虽然可以直接在数学定义上实现自己的erfcx() ,但这仅适用于有限的输入域,因为在正半平面erfc()下溢中等幅度的参数,而exp()溢出,例如,正如这个问题所述 。

为了与C一起使用,可以调整一些erfcx()开源实现,例如Faadeeva包中的实现 ,正如对这个问题的回答所指出的那样。 但是,这些实现通常不能为给定的浮点格式提供完全准确性。 例如,使用2 32个测试向量的测试显示由Faadeeva包提供的erfcx()的最大误差在正半平面中为8.41ulps,在负半平面中为511.68ulps。

准确实现的合理界限是4 ulps,对应于英特尔矢量数学库的LA配置文件中数学函数的精度界限,我发现这对于需要两者的非平凡数学函数实现是一个合理的界限。准确性好,性能好。

如何只使用C标准数学库,并且不需要外部库, erfcxf() erfcx()和相应的单精度版本erfcxf()如何准确实现? 我们可以假设C的float nad double类型映射到IEEE 754-2008 binary32binary64浮点类型。 可以假设硬件支持融合乘法 – 加法运算(FMA),因为此时所有主要处理器架构都支持此function。

到目前为止我发现的erfcx()实现的最佳方法是基于以下文章:

MM Shepherd和JG Laframboise,“Chebyshev逼近(1 + 2 x)exp(x 2 )erfc x在0≤x<∞。” 计算数学 ,第36卷,第153期,1981年1月,第249-253页(在线)

本文提出了巧妙的转换,它将缩放的互补误差函数映射到一个紧密有界的辅助函数,该函数适用于简单的多项式逼近。 为了提高性能,我已经尝试了各种变换,但所有这些都对精度产生了负面影响。 变换中常数K的选择(x-K)/(x + K)与核近似的精度具有非显而易见的关系。 我凭经验确定了“最佳”值,这与论文不同。

核心近似的参数转换和中间结果转换为erfcx结果会导致额外的舍入误差。 为了减轻它们对准确性的影响,我们需要应用补偿步骤,我在之前关于erfcf问题和答案中详细概述了这些步骤。 FMA的可用性极大地简化了这项任务。

生成的单精度代码如下所示:

/* * Based on: MM Shepherd and JG Laframboise, \"Chebyshev Approximation of * (1+2x)exp(x^2)erfc x in 0 <= x [-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 0x1.fffffep127f) r = e; // 3.40282347e+38 // avoid creating NaN } return r; }

负半平面中此实现的最大误差将取决于标准数学库的expf()实现的准确性。 使用英特尔编译器,版本13.1.3.198,并使用/fp:strict编译我在正半平面中观察到最大误差为2.00450 ulps,在穷举测试中观察到负半平面中的最大误差为2.38412 ulps。 我现在能说的最好,一个忠实的expf()将导致最大误差小于2.5 ulps。

请注意,虽然代码包含两个可能很慢的操作,但它们以特殊的倒数forms出现,因此可以在许多平台上使用快速倒数近似。 只要倒数近似被忠实地舍入,基于实验,对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() ,但需要使用更多项的erfcxf()多项式近似。 这在优化核心近似时提出了挑战,因为当搜索空间非常大时,许多启发式算法将会崩溃。 下面的系数代表了我迄今为止的最佳解决方案,并且肯定有改进的余地。 使用英特尔编译器和/fp:strict构建,并使用2 32个随机测试向量,观察到的最大误差在正半平面中为2.83788 ulps,在负半平面中为2.77856 ulps。

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 0x1.fffffffffffffp1023) r = e; // avoid creating NaN } return r; }

以上就是c/c++开发分享精确计算缩放互补误差函数,erfcx()相关内容,想了解更多C/C++开发(异常处理)及C/C++游戏开发关注(猴子技术宅)。

本站无法对海量内容真伪性鉴别,请勿相信本站任何号码,邮件,站外网址等信息,如有需要,请自行甄别。版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容, 请发送邮件至net@163.com举报,一经查实,本站将立刻删除。
(0)
上一篇 2022年7月29日 上午8:43
下一篇 2022年7月29日 上午8:53

相关推荐

发表回复

登录后才能评论