从log1pf()计算asinhf()的最准确方法是什么？

反双曲函数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%,是正确的圆润.可以使用定向舍入的基本算术,而无需额外的性能成本.