从 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
}

通过精确到 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 ulps 和绝大多数结果，比如 > 95%，都是正确四舍五入的。具有定向舍入的基本算术可能可以在不增加性能成本的情况下使用。

Answer 1

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

Openlibm，基于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 异常仍然有效）。

Answer 2

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

由于我的问题是关于最小化参数转换的错误，除了log1pf() 本身的错误之外，用于实验的最直接的方法是使用正确的四舍五入该对数函数的实现。请注意，在高性能环境的上下文中不太可能存在正确舍入的实现。根据 J.-M. 的作品。穆勒等。例如，要产生准确的单精度结果，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 单精度操作数进行测试表明，最大误差 1.49486070 ulp 发生在 ±0x1.ff5022p-9 处，并且有 353,521,140 个不正确的舍入结果。如果整个参数转换使用双精度算术会发生什么？代码更改为

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，正如我们在此处使用@987654336 观察到的那样@。

一种解决方案是将转换后的参数计算为高于单精度，例如作为双浮点操作数对（可以在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 个不正确舍入的结果。 asinhf() 的完整实现已经实现。

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

由于asinhf() 案例中的所有麻烦参数的量级都很小，因此可以[部分地？]通过对原点周围的区域使用多项式极小极大近似来解决准确性问题。由于这会创建另一个代码分支，因此可能会使矢量化更加困难。