【问题标题】:Very fast approximate Logarithm (natural log) function in C++?C ++中非常快速的近似对数(自然对数)函数?
【发布时间】:2017-02-10 19:03:00
【问题描述】:

我们找到了各种技巧来替换 std::sqrt (Timing Square Root) 和一些 std::exp (Using Faster Exponential Approximation) ,但我找不到可以替换 std::log

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

现在我使用ln(1+x) 的三阶泰勒近似值,x-0.5+0.5 之间(90% 的情况下最大误差为 4%),否则回退到 std::log。这给了我 15% 的加速。

【问题讨论】:

  • 在现代 CPU 上 std::sqrt 编译为一条指令。很难相信您可以以类似的准确性更快地完成任何事情。
  • @user3091460 如果float 精度足够,为什么不从cmath 调用logf()?还是您需要double 的完整输入域,但计算结果的精度仅相当于float(约6 个十进制数字)?
  • @user3091460 那么该站点上的错误计算不正确。 sqrtss 精确到完全精确,而 rsqrtss * x 后跟单个 Newton-Raphson 步骤仍然不能完全精确。
  • 是什么让您认为您正在实施的std::log 尚未使用可用于您的系统的最有效算法?如果你愿意为了速度而牺牲准确性(我可能会说快速得到错误答案),你需要在你的问题中这样说。
  • 现在我使用 ln(1+x) 的三阶泰勒近似,其中 x 在 -0.5 和 +0.5 之间(最大误差为 4% 的情况的 90%)并回退到 std: :log 否则。给了我 15% 的加速。

标签: c++ math logarithm micro-optimization sqrt


【解决方案1】:

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

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

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

对数函数的自定义实现通常包括将二进制浮点参数x 分离为指数e 和尾数m,这样x = m * 2e ,因此log(x) = log(2) * e + log(m)。选择m 使其接近于一,因为这提供了有效的近似值,例如log(m) = log(1+f) = log1p(f) by minimax polynomial approximation

C++ 提供frexp() 函数将浮点操作数分离为尾数和指数,但在实践中,通常使用更快的机器特定方法,通过将浮点数据重新解释为相同的方式在位级别操作浮点数据-size 整数。下面的单精度对数代码logf() 演示了这两种变体。函数__int_as_float()__float_as_int() 提供将int32_t 重新解释为IEEE-754 binary32 浮点数,反之亦然。此代码严重依赖于大多数当前处理器、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;
}

【讨论】:

  • 感谢,但是我无法在 win10 上使用 Msvc 15 找到 int_as_float 和 float_as_int。我发现它是 cuda 的一部分,但没有下载完整的包。
  • @user3091460 这些函数是机器特定功能的抽象。作为第一步,您可以简单地使用memcpy(),例如float __int_as_float(int32_t a) { float r; memcpy (&amp;r, &amp;a, sizeof(r)); return r;} 一个好的编译器可能会对此进行适当的优化,但根据您所针对的硬件(您尚未披露),可能有更好的方法,可能涉及内在函数或内联汇编。
  • @user3091460 和 njuffa:x86 的最佳 asm 可能会使用 SSE2 整数指令将浮点数作为整数进行任何操作,因为 XMM 寄存器用于标量/向量浮点数和向量整数。所以你可能应该_mm_set_ss(your_float)_mm_castps_si128 得到一个你可以操纵的__m128i。 (这可能会浪费一条指令将 xmm 寄存器的高位归零,because of design limitations in intrinsics.)。一个 MOVD 来获取/从整数寄存器中获取浮点位也可能很好。
  • @PeterCordes 明白了。我不打算投入交钥匙 SIMD 内在解决方案所需的大量工作,特别是考虑到仍然不清楚询问者的硬件上有哪些 ISA 扩展。考虑使用 SIMD 内在函数发布您自己的版本,我很乐意投票。
  • 我将尽可能使用一个高效的 float_to_int,它使用联合来进行类型双关,并使用 clang 和 gcc for x86 编译为单个 movd eax, xmm0godbolt.org/g/UCePpA。就像您希望的那样简单,@user3091460 :) 将整数操作为uint32_t 甚至可能更有效,因为整数指令更短,并且在 Haswell 上可以在 port6 上运行(它没有任何向量 ALU)。但可能留在 XMM 寄存器中会更好,因为您不会使用整数。
【解决方案2】:

查看this 讨论,接受的答案是指基于 Zeckendorf 分解计算对数的函数的implementation

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

希望这会有所帮助!

【讨论】:

  • 我去看看,谢谢
【解决方案3】:
#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
}

【讨论】:

  • 您的类型双关语违反了严格别名。 (使用memcpy 而不是指针转换。另外,您可能应该使用unsigned long long,因为您不需要算术移位。对于2 的补码机器上的正确性无关紧要,但仍然如此。)这也需要整数endian 匹配 float endian,就像在 x86 上一样,所以你至少应该记录一下。
  • 一些文本来解释表格查找策略和在某些输入范围内的实际相对/绝对精度,以及诸如 0 或负输入会发生的限制,将是一个好主意。
  • 您的表格只需为float。这会将您的缓存占用空间减少一半。 (但是 2^14 * 4 字节的表仍然是 64kiB。在大多数用例中,您会遇到很多缓存未命中,这就是为什么大多数快速日志实现在现代 CPU 上使用多项式近似,而不是表查找。尤其是什么时候可以使用 SIMD:Efficient implementation of log2(__m256d) in AVX2)
  • 彼得,很抱歉评论了一个非常古老的答案,但是严格的别名规则真的在这里被打破了吗?我猜你指的是 fast_log2 函数的第一行。我假设这里真的没有别名,并且 x 的值被复制,重新解释为 long long(与 memcpy 非常相似的行为)。除非我遗漏了什么,否则 t 不是别名,对吧?
【解决方案4】:

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

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

【讨论】:

  • 请注意,您的 mm256_fmaf 可以编译为单独的 mul 和 add 操作,中间产品的舍入。它保证是 FMA。 (当目标像大多数 AVX2 机器一样支持 FMA 时(并非全部:一种 VIA 设计),只有一些编译器,如 GCC,会为您“收缩”它为 FMA 指令。可能最好只针对 AVX2+FMA3 和使用_mm256_fmadd_ps,如果你愿意,可以使用可选的后备,但不是一个误导性的名称,默认情况下可能更慢的fma函数。
【解决方案5】:

这取决于您需要多准确。通常调用 log 来了解数字的大小,您可以通过检查浮点数的指数字段基本上免费做到这一点。这也是你的第一个近似值。我将为我的书“基本算法”插入一个插件,它解释了如何从第一原理实现标准库数学函数。

【讨论】:

  • 我正在寻找真正数学应用的自然对数,不需要双精度、浮点精度甚至 10-3、10-4 都很好
  • 没有引用相关部分的链接或参考书目不是答案
猜你喜欢
  • 1970-01-01
  • 2012-10-15
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2019-09-13
  • 1970-01-01
  • 2011-01-21
  • 1970-01-01
相关资源
最近更新 更多