从 log1pf（） 计算阿辛夫（）的最准确方法？

Question

反双曲函数< code>asinh()与自然对数密切相关。我试图确定从C99标准数学函数< code>log1p()计算< code>asinh()的最准确方法。为了便于实验，我现在将自己限制在IEEE-754单精度计算，也就是说，我正在研究< code>asinhf()和< code>log1pf()。我打算稍后重新使用完全相同的算法进行双精度计算，即< code>asinh()和< code>log1p()。

我的主要目标是最小化ulp错误，次要目标是在改进的代码最多比下面发布的版本慢一点的约束下，最小化不正确舍入的结果的数量。精度的任何增量改进，比如0.2 ulp，都是受欢迎的。添加几个fma(融合乘加)会很好，另一方面，我希望有人能找到一个采用快速< code>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
}

使用特定的< code>log1pf()实现，可以精确地

就性能而言，如果倒数除法和全除法花费的时间大致相同，则变量USE_RECIPROCAL=0将更快，但如果有非常快速的倒数支持，则变量USE_RECIPROCAL=1可能更快。

答案可以假设所有基本算法，包括FMA（融合乘法加法），都根据IEEE-754正确四舍五入到最接近或偶数模式。此外，可能会提供更快、几乎正确舍入的倒数和rsqrtf（）版本，其中“近似正确舍入”意味着最大ulp误差将限制在0.53 ulps左右，并且绝大多数结果都是如此

Answer 1

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

由于我的问题是关于最小化参数转换的错误，除了log1pf（）本身中的错误之外，最直接的实验方法是使用对数函数的正确舍入实现。请注意，在高性能环境中不太可能存在正确舍入的实现。根据J.-M. Muller等人的工作，要产生准确的单精度结果，例如，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;
}

使用我的问题中的第一个变体的< code>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的最大误差出现在< code>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 个错误舍入的结果，比以前略少。问题似乎是浮点操作数无法向 log1pf（） 传达足够的信息来产生更准确的结果。在计算辛夫（）和坐标（）时，也会出现类似的问题。如果将简化的参数（表示为正确舍入的浮点操作数）传递给核心多项式，则 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，它出现在< code > 0x 1 . dee EAP-22 处。有349，653，534个不正确的舍入结果。< code>asinhf()的完整实现已经实现。

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

由于在asinhf（）的情况下，所有麻烦的参数都很小，因此可以[部分地？]通过对原点周围区域使用多项式minimax近似来解决精度问题。由于这将创建另一个代码分支，因此可能会使矢量化更加困难。

Answer 2

首先，您可能需要查看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的旧System V错误代码（IEEE-754异常仍然有效）。