从log1pf（）计算asinhf（）最准确的方法是什么？

2018-06-30 16:24:12

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

我的主要目标是最大限度地减少ulp错误，次要目标是最大限度地减少错误舍入结果的数量，因为改进的代码至多会比下面发布的版本慢得多。任何对准确度的改进，比如说0.2 ulp，都会受到欢迎。添加一对FMA（融合乘法）将会很好，另一方面，我希望有人能够确定一个采用快速rsqrtf() （倒数平方根）的解决方案。

由此产生的C99代码应该适用于矢量化，可能是通过一些简单直接的转换。所有中间计算都必须以函数参数和结果的精度出现，因为任何切换到更高精度都会对性能产生严重的负面影响。该代码必须在IEEE-754非正常支持和FTZ（清零）模式下正常工作。

到目前为止，我已经确定了以下两个候选实现。请注意，可以通过对log1pf()的单个调用将代码轻松转换为无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 log1pf()的特定log1pf()实现，当在所有232个可能的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的round-to-nearest-or-even模式正确舍入。此外，更快，几乎正确舍入的版本的倒数和rsqrtf()可能是可用的，其中“几乎正确舍入”意味着最大ulp误差将被限制为诸如0.53 ulps之类的事情，并且绝大多数结果，例如> 95％，正确圆整。具有定向舍入的基本算术可以在不增加性能成本的情况下提供。

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

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

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

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

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

由于我的问题是关于最小化除了log1pf()本身中的错误之外发生的参数转换的错误， 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;
}

使用我的问题中的第一个变体的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
}

使用全部232个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()所观察到的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()情况下的所有麻烦参数的幅度都很小，因此可以通过对原点周围的区域使用多项式极大极小近似来解决准确性问题。由于这会创建另一个代码分支，它可能会使向量化更困难。

链接地址: http://www.djcxy.com/p/85595.html

上一篇: Most accurate way to compute asinhf() from log1pf()?

下一篇: When is it more efficient to use CORDIC or a polynomial approximation?