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

Question

反双曲函数asinh()与自然对数密切相关。我试图确定从 C99 标准数学函数 log1p() 计算 asinh() 的最准确方法。为了便于实验，我现在限制自己使用 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
}

使用精确到 < 0.6 ulps 的特定 log1pf() 实现，我在对所有 2³² 可能的 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 ulp 之类的范围内，并且绝大多数结果(例如 > 95%)都被正确舍入。具有定向舍入的基本算术可能可用，不会对性能产生额外成本。

Answer 1

首先，您可能想查看 log1pf 函数的准确性和速度:这些在 libms 之间可能会有所不同(我发现 OS X 数学函数很快，glibc速度较慢但通常正确舍入的)。

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

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

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