c - 从log1pf()计算asinhf()的最准确方法？

反双曲函数asinh()与自然对数密切相关。我正在尝试确定从C99标准数学函数asinh()计算log1p()的最准确方法。为了便于实验，我现在将自己限制为IEEE-754单精度计算，也就是说，我正在查看asinhf()和log1pf()。我打算稍后再使用完全相同的算法进行 double 计算，即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()实现，在对所有232种可能的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％)是正确舍入。可以直接进行带舍入的基本算术运算，而不会增加性能。

最佳答案

首先，您可能需要研究log1pf函数的准确性和速度:这些函数在libms之间可能会有所不同(我发现OS X数学函数较快，而glibc函数较慢，但通常会四舍五入)。

基于BSD libm(基于Sun的fdlibm)的Openlibm按范围使用了多种方法，但是主要位是the relation:

t = x*x;
w = log1pf(fabsf(x)+t/(one+sqrtf(one+t)));

您可能还想尝试使用-fno-math-errno选项进行编译，该选项将禁用sqrt的旧系统V错误代码(IEEE-754异常仍将起作用)。

ULP

c - 从log1pf()计算asinhf()的最准确方法？