nju*_*ffa 7 c math floating-point floating-accuracy
反双曲函数asinh()与自然对数密切相关.我试图asinh()从C99标准数学函数中确定最准确的计算方法log1p().为了便于实验,我现在限制自己使用IEEE-754单精度计算,我正在考虑asinhf()和log1pf().我打算重新使用双精度计算,即完全相同的算法asinh()和log1p(),后来.
我的主要目标是最小化ulp错误,次要目标是最小化错误舍入结果的数量,在改进代码最多比下面发布的版本最低的约束下.任何提高准确度的改进,比如说0.2 ulp,都是值得欢迎的.添加几个FMA(融合乘法 - 加法)会很好,另一方面我希望有人可以找到一个采用快速rsqrtf()(倒数平方根)的解决方案.
生成的C99代码应该适用于矢量化,可能是通过一些简单的直接转换.所有中间计算必须以函数参数和结果的精度发生,因为任何切换到更高精度可能会产生严重的负面性能影响.代码必须在IEEE-754非正常支持和FTZ(刷新到零)模式下正常工作.
到目前为止,我已经确定了以下两个候选实现.请注意,可以通过单次调用将代码轻松转换为无分支矢量化版本log1pf(),但在此阶段我还没有这样做以避免不必要的混淆.
/* for a >= 0, asinh(a) = log (a + sqrt (a*a+1))
= log1p (a + (sqrt (a*a+1) - 1))
= log1p (a + sqrt1pm1 (a*a))
= log1p (a + (a*a / (1 + sqrt(a*a + 1))))
= log1p (a + a * (a / (1 + sqrt(a*a + 1))))
= log1p (fma (a / (1 + sqrt(a*a + 1)), a, a)
= log1p (fma (1 / (1/a + sqrt(1/a*a + 1)), a, a)
*/
float my_asinhf (float a)
{
float fa, t;
fa = fabsf (a);
#if !USE_RECIPROCAL
if (fa >= 0x1.0p64f) { // prevent overflow in intermediate computation
t = log1pf (fa) + 0x1.62e430p-1f; // log(2)
} else {
t = fmaf (fa / (1.0f + sqrtf (fmaf (fa, fa, 1.0f))), fa, fa);
t = log1pf (t);
}
#else // USE_RECIPROCAL
if (fa > 0x1.0p126f) { // prevent underflow in intermediate computation
t = log1pf (fa) + 0x1.62e430p-1f; // log(2)
} else {
t = 1.0f / fa;
t = fmaf (1.0f / (t + sqrtf (fmaf (t, t, 1.0f))), fa, fa);
t = log1pf (t);
}
#endif // USE_RECIPROCAL
return copysignf (t, a); // restore sign
}
对于log1pf()精确到<0.6 ulps 的特定实现,我在所有2 32个可能的IEEE-754单精度输入中进行详尽测试时观察到以下错误统计.何时USE_RECIPROCAL = 0,最大误差为1.49486 ulp,并且有353,587,822个错误的舍入结果.使用时USE_RECIPROCAL = 1,最大误差为1.50805 ulp,并且只有77,569,390个错误的舍入结果.
在性能方面,USE_RECIPROCAL = 0如果倒数和完全划分花费的时间大致相同,则变量USE_RECIPROCAL = 1会更快,但如果可以获得非常快速的互惠支持,变量可能会更快.
答案可以假设所有基本算术,包括FMA(融合乘法 - 加法)都是根据IEEE-754舍入到最接近或偶数模式正确舍入的.此外,速度更快,近正确的舍入,互惠版本,并rsqrtf() 可能是可用的,其中"近正确舍入"是指最大ULP错误将被限制在类似0.53 ULPS最广大的结果,说> 95%,是正确的圆润.可以使用定向舍入的基本算术,而无需额外的性能成本.
经过各种额外的实验后，我确信，不使用比参数和结果更高的精度的简单参数转换无法实现比我发布的代码中的第一个变体所实现的更严格的误差范围。
\n\n由于我的问题是关于最小化除了log1pf()本身的错误之外所产生的参数转换的错误，因此用于实验的最直接的方法是利用该对数函数的正确舍入实现。请注意，在高性能环境中不太可能存在正确舍入的实现。根据 J.-M. 的著作。穆勒等。al.，为了产生准确的单精度结果，x86 扩展精度计算应该足够了，例如
float accurate_log1pf (float a)\n{\n float res;\n __asm fldln2;\n __asm fld dword ptr [a];\n __asm fyl2xp1;\n __asm fst dword ptr [res];\n __asm fcompp;\n return res;\n}\n使用我的问题中的第一个变体的实现asinhf()如下所示：
float my_asinhf (float a)\n{\n float fa, s, t;\n fa = fabsf (a);\n if (fa >= 0x1.0p64f) { // prevent overflow in intermediate computation\n t = log1pf (fa) + 0x1.62e430p-1f; // log(2)\n } else {\n t = fmaf (fa / (1.0f + sqrtf (fmaf (fa, fa, 1.0f))), fa, fa);\n t = accurate_log1pf (t);\n }\n return copysignf (t, a); // restore sign\n}\n对所有 2 32 个IEEE-754 单精度操作数进行测试表明，最大错误 1.49486070 ulp 出现在 \xc2\xb1 处0x1.ff5022p-9，并且有 353,521,140 个不正确舍入的结果。如果整个参数转换使用双精度算术会发生什么？代码更改为
float my_asinhf (float a)\n{\n float fa, s, t;\n fa = fabsf (a);\n if (fa >= 0x1.0p64f) { // prevent overflow in intermediate computation\n t = log1pf (fa) + 0x1.62e430p-1f; // log(2)\n } else {\n double tt = fa;\n tt = fma (tt / (1.0 + sqrt (fma (tt, tt, 1.0))), tt, tt);\n t = (float)tt;\n t = accurate_log1pf (t);\n }\n return copysignf (t, a); // restore sign\n}\n然而，错误界限并没有随着此更改而改善！最大错误 1.49486070 ulp 仍然出现在 \xc2\xb1 处0x1.ff5022p-9，并且现在有 350,971,046 个错误舍入结果，比以前略有减少。问题似乎是float操作数无法传达足够的信息来log1pf()产生更准确的结果。sinf()计算和时也会出现类似的问题cosf()。如果将简化的参数（表示为正确舍入的float操作数）传递给核心多项式，则 和 中产生的误差sinf()仅cosf()略低于 1.5 ulp，正如我们在此处观察到的那样my_asinhf()。
一种解决方案是将转换后的参数计算为高于单精度，例如作为双浮点操作数对（Andrew Thall 的论文中可以找到双浮点技术的有用简要概述）。在这种情况下，我们可以基于对数导数是倒数的知识，使用附加信息对结果执行线性插值。这给了我们：
\n\nfloat my_asinhf (float a)\n{\n float fa, s, t;\n fa = fabsf (a);\n if (fa >= 0x1.0p64f) { // prevent overflow in intermediate computation\n t = log1pf (fa) + 0x1.62e430p-1f; // log(2)\n } else {\n double tt = fa;\n tt = fma (tt / (1.0 + sqrt (fma (tt, tt, 1.0))), tt, tt);\n t = (float)tt; // "head" of double-float\n s = (float)(tt - (double)t); // "tail" of double-float\n t = fmaf (s, 1.0f / (1.0f + t), accurate_log1pf (t)); // interpolate\n }\n return copysignf (t, a); // restore sign\n}\n该版本的详尽测试表明最大错误已减少到 0.99999948 ulp，它发生在 \xc2\xb1 处0x1.deeea0p-22。有 349,653,534 个错误舍入结果。asinhf()已经实现了忠实全面的实施。
不幸的是，这个结果的实际用途是有限的。根据硬件平台的不同，算术运算的吞吐量double可能仅为运算吞吐量的 1/2 至 1/32 float。双精度计算可以用双浮点计算代替，但这也会产生非常大的成本。最后，我在这里的方法是使用单精度实现作为后续双精度工作的试验场，并且许多硬件平台（当然是我感兴趣的所有硬件平台）不提供对更高精度数字格式的硬件支持比 IEEE-754binary64（双精度）。因此任何解决方案都不应该在中间计算中要求更高精度的算术。
由于在这种情况下所有麻烦的争论asinhf()量级都很小，因此可以通过对原点周围的区域使用多项式极小极大近似来[部分？]解决准确性问题。由于这会创建另一个代码分支，因此可能会使矢量化变得更加困难。