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

Question

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

使用精确到 < 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() 可能可用，其中 "nearly correctly rounded" 表示最大 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 的旧系统 V 错误代码（IEEE-754 异常仍然有效）。

Answer 2

经过各种额外的实验，我确信一个简单的参数转换 不使用比参数和结果更高的精度 无法实现比所实现的更严格的误差范围通过我发布的代码中的第一个变体。

因为我的问题是关于最小化除了 log1pf() 本身的错误之外的参数转换的错误，最直接的方法用于 experimentation 是利用该对数函数的正确舍入实现。请注意，correctly-rounded 实现极不可能存在于 high-performance 环境的上下文中。根据 J.-M. 的作品。穆勒等。等，为了产生准确的 single-precision 结果，x86 扩展精度计算应该足够了，例如

float accurate_log1pf (float a)
{
    float res;
    __asm fldln2;
    __asm fld     dword ptr [a];
    __asm fyl2xp1;
    __asm fst     dword ptr [res];
    __asm fcompp;
    return res;
}

asinhf() 使用我的问题的第一个变体的实现如下所示：

float my_asinhf (float a)
{
    float fa, s, t;
    fa = fabsf (a);
    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 = accurate_log1pf (t);
    }
    return copysignf (t, a); // restore sign
}

对所有 2³² IEEE-754 single-precision 操作数进行的测试表明，最大误差 1.49486070 ulp 出现在 ±0x1.ff5022p-9 处，并且有 353,521,140 个不正确四舍五入的结果。如果整个参数转换使用 double-precision 算术会怎样？代码改为

float my_asinhf (float a)
{
    float fa, s, t;
    fa = fabsf (a);
    if (fa >= 0x1.0p64f) { // prevent overflow in intermediate computation
        t = log1pf (fa) + 0x1.62e430p-1f; // log(2)
    } else {
        double tt = fa;
        tt = fma (tt / (1.0 + sqrt (fma (tt, tt, 1.0))), tt, tt);
        t = (float)tt;
        t = accurate_log1pf (t);
    }
    return copysignf (t, a); // restore sign
}

但是，此更改并没有改善错误范围！最大误差 1.49486070 ulp 仍然出现在 ±0x1.ff5022p-9 处，现在有 350,971,046 个不正确的四舍五入结果，比以前略少。问题似乎是 float 操作数无法将足够的信息传递给 log1pf() 以产生更准确的结果。计算 sinf() 和 cosf() 时也会出现类似的问题。如果将表示为正确舍入的 float 操作数的简化参数传递给核心多项式，则 sinf() 和 cosf() 中产生的误差略低于 1.5 ulp，正如我们my_asinhf().

在这里观察

一个解决方案是将转换后的参数计算为高于单精度，例如作为 double-float 操作数对（double-float 技术的有用简要概述可以在 this paper by Andrew Thall).在这种情况下，我们可以利用附加信息对结果进行线性插值，基于对数的导数是倒数的知识。这给了我们：

float my_asinhf (float a)
{
    float fa, s, t;
    fa = fabsf (a);
    if (fa >= 0x1.0p64f) { // prevent overflow in intermediate computation
        t = log1pf (fa) + 0x1.62e430p-1f; // log(2)
    } else {
        double tt = fa;
        tt = fma (tt / (1.0 + sqrt (fma (tt, tt, 1.0))), tt, tt);
        t = (float)tt;                // "head" of double-float
        s = (float)(tt - (double)t);  // "tail" of double-float
        t = fmaf (s, 1.0f / (1.0f + t), accurate_log1pf (t)); // interpolate
    }
    return copysignf (t, a); // restore sign
}

此版本的详尽测试表明最大误差已降至 0.99999948 ulp，发生在 ±0x1.deeea0p-22。有 349,653,534 个不正确的四舍五入结果。 faithfully-rounded 实现了 asinhf()。

不幸的是，这个结果的实际用途是有限的。根据硬件平台的不同，double算术运算的吞吐量可能只有float运算吞吐量的1/2到1/32。 double-precision 计算可以替换为 double-float 计算，但这也会产生非常大的成本。最后，我在这里的方法是使用 single-precision 实现作为后续 double-precision 工作的试验场，许多硬件平台（当然是我感兴趣的所有平台）不提供对数字的硬件支持精度高于 IEEE-754 binary64（双精度）的格式。因此，任何解决方案都不应在中间计算中需要 higher-precision 算法。

由于在 asinhf() 的情况下所有麻烦的论点都很小，因此可以 [部分？] 通过对原点周围区域使用多项式极小极大近似值来解决精度问题。由于这会创建另一个代码分支，因此可能会使向量化变得更加困难。

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

Most accurate way to compute asinhf() from log1pf()?

c

math

floating-point

floating-accuracy