C++中非常快速的對數（自然對數）函數？

Question

我們發現代替std::sqrt（Timing Square Root）和std::exp（Using Faster Exponential Approximation）的一些技巧，但我找不到替代std::log的東西。它是我程序中循環的一部分，它被多次調用，而exp和sqrt被優化，因此英特爾VTune現在建議我優化std::log，之後似乎只有我的設計選擇會受到限制。

很多謝謝。

Answer 1

這取決於你需要的準確度。通常會通過調用log來了解數字的大小，通過檢查浮點數的指數字段，您可以基本上免費進行此操作。這也是你的第一個近似值。我將爲我的「基本算法」一書提供一個插件，它解釋瞭如何從第一原理實現標準庫數學函數。

Answer 2

看看this的討論，接受的答案是指基於Zeckendorf分解的計算對數函數implementation。

在實現文件的註釋中，討論了複雜性和一些技巧以達到O（1）。

希望這會有所幫助！

Answer 3

在開始設計和部署性能超越函數的自定義實現之前，強烈建議在算法級別以及通過工具鏈進行優化。不幸的是，我們沒有任何有關此處優化的代碼的信息，我們也沒有關於工具鏈的信息。

在算法級別，檢查是否所有對超越函數的調用都是真正必需的。也許有一個數學轉換需要較少的函數調用，或將超越函數轉換爲代數運算。任何超越函數調用都可能是多餘的，例如因爲計算是不必要地切入和切出對數空間？如果精度要求適中，整個計算是否可以單精度執行，全部使用float而不是double？在大多數硬件平臺上，避免使用double計算可以顯着提高性能。

編譯器傾向於提供影響數字密集代碼性能的各種開關。除了將通用優化級別提高到-O3之外，通常還有一種方法可以禁用非正常支持，即打開清零或FTZ模式。這在各種硬件平臺上具有性能優勢。此外，通常會有一個「快速數學」標誌，其使用會導致準確度略有下降，並消除了處理特殊情況（如NaN和無窮大）的開銷，以及處理errno。一些編譯器還支持自動向量化代碼，並附帶一個SIMD數學庫，例如英特爾編譯器。

一個對數函數通常涉及分離二進制浮點參數x成指數e和尾數m的自定義實現，使得x = m * 2e，因此log(x) = log(2) * e + log(m)。選擇m以使其接近於1，因爲這提供了有效的近似值，例如log(m) = log(1+f) = log1p(f),minimax polynomial approximation。

C++提供了frexp()函數將浮點操作數分離爲尾數和指數，但實際上通常使用更快的機器特定方法，通過將它們重新解釋爲相同位數來在位級操作浮點數據，大小整數。下面的單精度對數代碼logf()演示了兩種變體。功能__int_as_float()和__float_as_int()規定將int32_t重新解釋爲IEEE-754 binary32浮點數，反之亦然。該代碼嚴重依賴於大多數當前處理器，CPU或GPU上的硬件中直接支持的融合乘加操作FMA。在fmaf()映射到軟件仿真的平臺上，此代碼的速度會慢得令人無法接受。

#include <cmath> 
#include <cstdint> 

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

    if ((a > 0.0f) && (a <= 3.40282347e+38f)) { // 0x1.fffffep+127 
#if PORTABLE 
     m = frexpf (a, &e); 
     if (m < 0.666666667f) { 
      m = m + m; 
      e = e - 1; 
     } 
     i = (float)e; 
#else // PORTABLE 
     i = 0.0f; 
     /* fix up denormal inputs */ 
     if (a < 1.175494351e-38f){ // 0x1.0p-126 
      a = a * 8388608.0f; // 0x1.0p+23 
      i = -23.0f; 
     } 
     e = (__float_as_int (a) - 0x3f2aaaab) & 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] */ 
     f = m - 1.0f; 
     s = f * f; 
     /* Compute log1p(f) for f in [-1/3, 1/3] */ 
     r = fmaf (-0.130187988f, f, 0.140889585f); // -0x1.0aa000p-3, 0x1.208ab8p-3 
     t = fmaf (-0.121489584f, f, 0.139809534f); // -0x1.f19f10p-4, 0x1.1e5476p-3 
     r = fmaf (r, s, t); 
     r = fmaf (r, f, -0.166845024f); // -0x1.55b2d8p-3 
     r = fmaf (r, f, 0.200121149f); // 0x1.99d91ep-3 
     r = fmaf (r, f, -0.249996364f); // -0x1.fffe18p-3 
     r = fmaf (r, f, 0.333331943f); // 0x1.5554f8p-2 
     r = fmaf (r, f, -0.500000000f); // -0x1.000000p-1 
     r = fmaf (r, s, f); 
     r = fmaf (i, 0.693147182f, r); // 0x1.62e430p-1 // log(2) 
    } else { 
     r = a + a; // silence NaNs if necessary 
     if (a < 0.0f) r = 0.0f/0.0f; // NaN 
     if (a == 0.0f) r = -1.0f/0.0f; // -Inf 
    } 
    return r; 
} 

如代碼註釋指出，實施上述規定如實全面的單精度結果，並將其與與IEEE-754浮點標準一致的特殊情況下的交易。通過消除特殊情況支持，可以進一步提高性能，消除對非正常參數的支持，並降低準確性。這導致以下示例變體：8分鐘招搖題外話問題後

/* 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; 
}