C语言 从 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()版本，其中"几乎正确舍入"意味着最大 ulp 误差将限制在 0.53 ulps 左右，绝大多数结果(例如 95%(>正确舍入。带有有向四舍五入的基本算术可能无需额外的性能成本。

Answer 1

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

Openlibm，基于BSD libm，而BSD libm又基于Sun的fdlibm，按范围使用多种方法，但主要的一点是关系：

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

您可能还想尝试使用 -fno-math-errno 选项进行编译，该选项禁用sqrt的旧系统 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^{32 个 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，就像我们在这里用 my_asinhf() 观察到的那样。

一种解决方案是将转换后的参数计算到高于单精度，例如作为双浮点操作数对(双浮点数技术的有用简要介绍可以在 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()情况下所有麻烦的参数在量级上都很小，因此可以通过对原点周围的区域使用多项式最小最大值近似来解决精度问题。由于这会创建另一个代码分支，因此可能会使矢量化更加困难。