C中非常快速的近似对数（自然对数）函数？

Question

我们找到了各种技巧来替换 std::sqrt （ 时序平方根）和一些 std::exp （ 使用更快的指数近似），但我找不到可以替换的东西 std::log

它是我程序中循环的一部分，它被多次调用，虽然 exp 和 sqrt 被优化，但英特尔 VTune 现在建议我优化 std::log ，之后似乎只有我的设计选择会受到限制。

For now I use a 3rd order taylor approximation of ln(1+x) with x between -0.5 and +0.5 (90% of the case for max error 4%）并回退到 std::log 否则。这给了我 15% 的加速。

原文由 qwark 发布，翻译遵循 CC BY-SA 4.0 许可协议

Answer 1

在着手设计和部署性能超越函数的定制实现之前，强烈建议在算法级别以及通过工具链进行优化。不幸的是，我们没有任何关于这里要优化的代码的信息，也没有关于工具链的信息。

在算法层面，检查所有对超越函数的调用是否真的有必要。也许有一个数学变换需要更少的函数调用，或者将超越函数转换为代数运算。是否有任何超越函数调用可能是多余的，例如因为计算不必要地进出对数空间？如果精度要求适中，整个计算能否以单精度执行，始终使用 float 而不是 double ？在大多数硬件平台上，避免 double 计算可以显着提升性能。

编译器倾向于提供各种影响数字密集型代码性能的开关。除了将一般优化级别提高到 -O3 ，通常还有一种方法可以关闭非正规支持，即打开清零或 FTZ 模式。这在各种硬件平台上具有性能优势。此外，通常有一个“快速数学”标志，其使用会导致准确性略有降低，并消除了处理特殊情况（如 NaN 和无穷大）以及处理 errno 的开销。一些编译器还支持代码的自动矢量化并附带 SIMD 数学库，例如英特尔编译器。

A custom implementation of a logarithm function typically involves separating a binary floating-point argument x into an exponent e and a m , such that x = m * 2 e ，因此 log(x) = log(2) * e + log(m) 。选择 m 使其接近于一，因为这提供了有效的近似，例如 log(m) = log(1+f) = log1p(f) 通过 极小多项式近似。

C++ 提供了 frexp() 函数将浮点操作数分离为尾数和指数，但在实践中，通常使用更快的机器特定方法，通过将浮点数据重新解释为相同大小的整数。下面的单精度对数代码 logf() 演示了这两种变体。 Functions __int_as_float() and __float_as_int() provide for the reinterpretation of an int32_t into an IEEE-754 binary32 floating-point number and vice-versa.此代码严重依赖于大多数当前处理器、CPU 或 GPU 的硬件中直接支持的融合乘加操作 FMA。在 fmaf() 映射到软件仿真的平台上，此代码的速度会慢得令人无法接受。

 #include <cmath>
#include <cstdint>
#include <cstring>

float __int_as_float (int32_t a) { float r; memcpy (&r, &a, sizeof r); return r;}
int32_t __float_as_int (float a) { int32_t r; memcpy (&r, &a, sizeof r); return r;}

/* compute natural logarithm, maximum error 0.85089 ulps */
float my_logf (float a)
{
    float i, m, r, s, t;
    int e;

#if PORTABLE
    m = frexpf (a, &e);
    if (m < 0.666666667f) {
        m = m + m;
        e = e - 1;
    }
    i = (float)e;
#else // PORTABLE
    i = 0.0f;
    if (a < 1.175494351e-38f){ // 0x1.0p-126
        a = a * 8388608.0f; // 0x1.0p+23
        i = -23.0f;
    }
    e = (__float_as_int (a) - __float_as_int (0.666666667f)) & 0xff800000;
    m = __int_as_float (__float_as_int (a) - e);
    i = fmaf ((float)e, 1.19209290e-7f, i); // 0x1.0p-23
#endif // PORTABLE
    /* m in [2/3, 4/3] */
    m = m - 1.0f;
    s = m * m;
    /* Compute log1p(m) for m in [-1/3, 1/3] */
    r =             -0.130310059f;  // -0x1.0ae000p-3
    t =              0.140869141f;  //  0x1.208000p-3
    r = fmaf (r, s, -0.121483512f); // -0x1.f198b2p-4
    t = fmaf (t, s,  0.139814854f); //  0x1.1e5740p-3
    r = fmaf (r, s, -0.166846126f); // -0x1.55b36cp-3
    t = fmaf (t, s,  0.200120345f); //  0x1.99d8b2p-3
    r = fmaf (r, s, -0.249996200f); // -0x1.fffe02p-3
    r = fmaf (t, m, r);
    r = fmaf (r, m,  0.333331972f); //  0x1.5554fap-2
    r = fmaf (r, m, -0.500000000f); // -0x1.000000p-1
    r = fmaf (r, s, m);
    r = fmaf (i,  0.693147182f, r); //  0x1.62e430p-1 // log(2)
    if (!((a > 0.0f) && (a < INFINITY))) {
        r = a + a;  // silence NaNs if necessary
        if (a  < 0.0f) r = INFINITY - INFINITY; //  NaN
        if (a == 0.0f) r = -INFINITY;
    }
    return r;
}

如代码注释中所述，上面的实现提供了忠实四舍五入的单精度结果，并且它处理与 IEEE-754 浮点标准一致的例外情况。通过消除特殊情况支持、消除对非规范参数的支持以及降低准确性，可以进一步提高性能。这导致以下示例性变体：

 /* natural log on [0x1.f7a5ecp-127, 0x1.fffffep127]. Maximum relative error 9.4529e-5 */
float my_faster_logf (float a)
{
    float m, r, s, t, i, f;
    int32_t e;

    e = (__float_as_int (a) - 0x3f2aaaab) & 0xff800000;
    m = __int_as_float (__float_as_int (a) - e);
    i = (float)e * 1.19209290e-7f; // 0x1.0p-23
    /* m in [2/3, 4/3] */
    f = m - 1.0f;
    s = f * f;
    /* Compute log1p(f) for f in [-1/3, 1/3] */
    r = fmaf (0.230836749f, f, -0.279208571f); // 0x1.d8c0f0p-3, -0x1.1de8dap-2
    t = fmaf (0.331826031f, f, -0.498910338f); // 0x1.53ca34p-2, -0x1.fee25ap-2
    r = fmaf (r, s, t);
    r = fmaf (r, s, f);
    r = fmaf (i, 0.693147182f, r); // 0x1.62e430p-1 // log(2)
    return r;
}

原文由 njuffa 发布，翻译遵循 CC BY-SA 4.0 许可协议

Answer 2

看看 这个 讨论，接受的答案是指基于 Zeckendorf 分解计算对数的函数的 实现。

在实现文件的注释中，讨论了复杂性和达到 O(1) 的一些技巧。

希望这可以帮助！

原文由 Mo Abdul-Hameed 发布，翻译遵循 CC BY-SA 3.0 许可协议