C++中非常快速的近似对数(自然对数)函数?

inn6fuwd  于 2022-12-01  发布在  其他
关注(0)|答案(7)|浏览(384)

我们找到了各种各样的技巧来替换std::sqrtTiming Square Root),也找到了一些技巧来替换std::expUsing Faster Exponential Approximation),但我找不到任何技巧来替换std::log
它是我的程序中循环的一部分,它被多次调用,虽然exp和sqrt已经过优化,但英特尔VTune现在建议我优化std::log,在此之后,似乎只有我的设计选择会受到限制。
现在我使用ln(1+x)的三阶泰勒近似,其中x-0.5+0.5之间(90%的情况下最大误差为4%),否则回到std::log。这给了我15%的加速。

ghg1uchk

ghg1uchk1#

在开始设计和部署超越函数的自定义实现以提高性能之前,最好在算法级别以及通过工具链进行优化。遗憾的是,我们没有关于要优化的代码的任何信息,也没有关于工具链的信息。
在算法层面,检查是否所有超越函数的调用都是真正必要的。也许有一个数学变换需要更少的函数调用,或者将超越函数转换为代数运算。是否有任何超越函数调用可能是多余的,例如,因为计算不必要地在对数空间内外切换?如果精度要求不高,整个计算是否可以使用float而不是double以单精度执行?在大多数硬件平台上,避免double计算可以显著提高性能。
编译器倾向于提供各种开关,这些开关会影响数字密集型代码的性能。除了将常规优化级别提高到-O3之外,通常还有一种方法可以关闭非规格化支持,即打开刷新为零或FTZ模式。这在各种硬件平台上都有性能优势。此外,通常会有一个“快速数学”标志,使用该标志会略微降低精度,并消除处理特殊情况(如NaN和无穷大)以及处理errno的开销。一些编译器还支持代码的自动矢量化,并附带SIMD数学库,例如Intel编译器。
对数函数的定制实现通常涉及将二进制浮点自变量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()Map到软件仿真的平台上,此代码将非常慢。

#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;
}
wsewodh2

wsewodh22#

看一下this的讨论,公认的答案是基于Zeckendorf分解计算对数的函数的implementation
在实现文件的注解中有一个关于复杂度的讨论和一些达到O(1)的技巧。
希望这对你有帮助!

vom3gejh

vom3gejh3#

#include <math.h>
#include <iostream>

constexpr int LogPrecisionLevel = 14;
constexpr int LogTableSize = 1 << LogPrecisionLevel;

double log_table[LogTableSize];

void init_log_table() {
    for (int i = 0; i < LogTableSize; i++) {
        log_table[i] = log2(1 + (double)i / LogTableSize);
    }
}

double fast_log2(double x) { // x>0
    long long t = *(long long*)&x;
    int exp = (t >> 52) - 0x3ff;
    int mantissa = (t >> (52 - LogPrecisionLevel)) & (LogTableSize - 1);
    return exp + log_table[mantissa];
}

int main() {
    init_log_table();

    double d1 = log2(100); //6.6438561897747244
    double d2 = fast_log2(100); //6.6438561897747244
    double d3 = log2(0.01); //-6.6438561897747244
    double d4 = fast_log2(0.01); //-6.6438919626096089
}
yacmzcpb

yacmzcpb4#

我矢量化了@njuffa的答案。自然对数,适用于AVX 2:

inline __m256 mm256_fmaf(__m256 a, __m256 b, __m256 c){
    return _mm256_add_ps(_mm256_mul_ps(a, b), c);
}

//https://stackoverflow.com/a/39822314/9007125
//https://stackoverflow.com/a/65537754/9007125
// vectorized version of the answer by njuffa
/* natural log on [0x1.f7a5ecp-127, 0x1.fffffep127]. Maximum relative error 9.4529e-5 */
inline __m256 fast_log_sse(__m256 a){

    __m256i aInt = *(__m256i*)(&a);
    __m256i e =    _mm256_sub_epi32( aInt,  _mm256_set1_epi32(0x3f2aaaab));
            e =    _mm256_and_si256( e,  _mm256_set1_epi32(0xff800000) );
        
    __m256i subtr =  _mm256_sub_epi32(aInt, e);
    __m256 m =  *(__m256*)&subtr;

    __m256 i =  _mm256_mul_ps( _mm256_cvtepi32_ps(e), _mm256_set1_ps(1.19209290e-7f));// 0x1.0p-23
    /* m in [2/3, 4/3] */
    __m256 f =  _mm256_sub_ps( m,  _mm256_set1_ps(1.0f) );
    __m256 s =  _mm256_mul_ps(f, f); 
    /* Compute log1p(f) for f in [-1/3, 1/3] */
    __m256 r =  mm256_fmaf( _mm256_set1_ps(0.230836749f),  f,  _mm256_set1_ps(-0.279208571f) );// 0x1.d8c0f0p-3, -0x1.1de8dap-2
    __m256 t =  mm256_fmaf( _mm256_set1_ps(0.331826031f),  f,  _mm256_set1_ps(-0.498910338f) );// 0x1.53ca34p-2, -0x1.fee25ap-2

           r =  mm256_fmaf(r, s, t);
           r =  mm256_fmaf(r, s, f);
           r =  mm256_fmaf(i, _mm256_set1_ps(0.693147182f),  r);  // 0x1.62e430p-1 // log(2)
    return r;
}
xwbd5t1u

xwbd5t1u5#

我也需要一个快速的对数近似,到目前为止,最好的一个似乎是一个基于Ankerls算法。
说明:http://martin.ankerl.com/2007/10/04/optimized-pow-approximation-for-java-and-c-c/
代码(复制自):

double log_fast_ankerl(double a)
{
  static_assert(__BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__, "Little endian is required!");
  union { double d; int x[2]; } u = { a };
  return (u.x[1] - 1072632447) * 6.610368362777016e-7;
}

只是一个减法和乘法。它是令人惊讶的好和无与伦比的快。

zqry0prt

zqry0prt6#

精度和性能略有提高:

#include <bit>   // C++20

//fast_log abs(rel) : avgError = 2.85911e-06(3.32628e-08), MSE = 4.67298e-06(5.31012e-08), maxError = 1.52588e-05(1.7611e-07)
const float s_log_C0 = -19.645704f;
const float s_log_C1 = 0.767002f;
const float s_log_C2 = 0.3717479f;
const float s_log_C3 = 5.2653985f;
const float s_log_C4 = -(1.0f + s_log_C0) * (1.0f + s_log_C1) / ((1.0f + s_log_C2) * (1.0f + s_log_C3)); //ensures that log(1) == 0
const float s_log_2 = 0.6931472f;

// assumes x > 0 and that it's not a subnormal.
// Results for 0 or negative x won't be -Infinity or NaN
inline float fast_log(float x)
{
    unsigned int ux = std::bit_cast<unsigned int>(x);
    int e = static_cast<int>(ux - 0x3f800000) >> 23; //e = exponent part can be negative
    ux |= 0x3f800000;
    ux &= 0x3fffffff; // 1 <= x < 2  after replacing the exponent field
    x = std::bit_cast<float>(ux);
    float a = (x + s_log_C0) * (x + s_log_C1);
    float b = (x + s_log_C2) * (x + s_log_C3);
    float c = (float(e) + s_log_C4);
    float d = a / b;
    return (c + d) * s_log_2;
}

虽然它有被认为是慢操作的除法,但有许多操作可以由处理器并行执行。

mfpqipee

mfpqipee7#

这取决于你需要的精确度。通常,log被调用来得到数字的大小,你可以通过检查浮点数的指数字段来免费地得到。这也是你的第一近似。我会在我的书《基本算法》中插入一个插件,这本书解释了如何从第一原理实现标准库数学函数。

相关问题