C++ 中非常快速的近似对数（自然对数）函数？
Very fast approximate Logarithm (natural log) function in C++?
我们找到了各种技巧来替换 std::sqrt (Timing Square Root) and some for std::exp (Using Faster Exponential Approximation) ，但我找不到任何东西可以替换 std::log.
它是我程序中循环的一部分，它被多次调用，虽然 exp 和 sqrt 已优化，Intel VTune 现在建议我优化 std::log，之后似乎只有我的设计选择会受到限制.
现在我使用 ln(1+x) 的三阶泰勒近似，x 在 -0.5 和 +0.5 之间（90% 的最大误差为 4% ) 否则退回到 std::log 。这给了我 15% 的加速。
这取决于您需要的准确度。通常会调用 log 来了解数字的大小，这基本上可以通过检查浮点数的指数字段来免费完成。这也是您的第一个近似值。我将为我的书 "Basic Algorithms" 添加一个插件，它解释了如何从第一原理实现标准库数学函数。
看一下 this discussion, the accepted answer refers to an implementation 基于 Zeckendorf 分解计算对数的函数。
在实现文件的注释中讨论了复杂性和达到 O(1) 的一些技巧。
希望对您有所帮助！
在着手设计和部署超越性能的自定义实现之前，强烈建议在算法级别以及通过工具链进行优化。不幸的是，我们没有关于这里要优化的代码的任何信息，也没有关于工具链的信息。
在算法层面，检查是否所有对超越函数的调用都是真正必要的。也许有一种数学变换需要更少的函数调用，或者将超越函数转换为代数运算。是否有任何超越函数调用可能是多余的，例如因为计算不必要地切换进出对数space？如果精度要求适中，整个计算是否可以单精度执行，始终使用 float 而不是 double？在大多数硬件平台上，避免 double 计算可以显着提高性能。
编译器往往会提供各种影响数字密集型代码性能的开关。除了将一般优化级别提高到 -O3 之外，通常还有一种方法可以关闭非正规支持，即打开清零或 FTZ 模式。这在各种硬件平台上都有性能优势。此外，通常有一个“快速数学”标志，其使用会导致准确性略有降低，并消除处理特殊情况（例如 NaN 和无穷大）以及 errno 处理的开销。一些编译器还支持代码的自动矢量化并附带 SIMD 数学库，例如 Intel 编译器。
对数函数的自定义实现通常涉及将二进制浮点参数 x 分成指数 e 和尾数 m，这样 x = m * 2 e，因此 log(x) = log(2) * e + log(m)。选择 m 使其接近统一，因为这提供了有效的近似值，例如 log(m) = log(1+f) = log1p(f) by 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>
#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;
}
#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
}
我矢量化了@njuffa 的答案。自然对数，适用于 AVX2：
inline __m256 mm256_fmaf(__m256 a, __m256 b, __m256 c){
return _mm256_add_ps(_mm256_mul_ps(a, b), c);
}
//
//
// 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;
}
我们找到了各种技巧来替换 std::sqrt (Timing Square Root) and some for std::exp (Using Faster Exponential Approximation) ，但我找不到任何东西可以替换 std::log.
它是我程序中循环的一部分，它被多次调用，虽然 exp 和 sqrt 已优化，Intel VTune 现在建议我优化 std::log，之后似乎只有我的设计选择会受到限制.
现在我使用 ln(1+x) 的三阶泰勒近似，x 在 -0.5 和 +0.5 之间（90% 的最大误差为 4% ) 否则退回到 std::log 。这给了我 15% 的加速。
这取决于您需要的准确度。通常会调用 log 来了解数字的大小，这基本上可以通过检查浮点数的指数字段来免费完成。这也是您的第一个近似值。我将为我的书 "Basic Algorithms" 添加一个插件，它解释了如何从第一原理实现标准库数学函数。
看一下 this discussion, the accepted answer refers to an implementation 基于 Zeckendorf 分解计算对数的函数。
在实现文件的注释中讨论了复杂性和达到 O(1) 的一些技巧。
希望对您有所帮助！
在着手设计和部署超越性能的自定义实现之前，强烈建议在算法级别以及通过工具链进行优化。不幸的是，我们没有关于这里要优化的代码的任何信息，也没有关于工具链的信息。
在算法层面，检查是否所有对超越函数的调用都是真正必要的。也许有一种数学变换需要更少的函数调用，或者将超越函数转换为代数运算。是否有任何超越函数调用可能是多余的，例如因为计算不必要地切换进出对数space？如果精度要求适中，整个计算是否可以单精度执行，始终使用 float 而不是 double？在大多数硬件平台上，避免 double 计算可以显着提高性能。
编译器往往会提供各种影响数字密集型代码性能的开关。除了将一般优化级别提高到 -O3 之外，通常还有一种方法可以关闭非正规支持，即打开清零或 FTZ 模式。这在各种硬件平台上都有性能优势。此外，通常有一个“快速数学”标志，其使用会导致准确性略有降低，并消除处理特殊情况（例如 NaN 和无穷大）以及 errno 处理的开销。一些编译器还支持代码的自动矢量化并附带 SIMD 数学库，例如 Intel 编译器。
对数函数的自定义实现通常涉及将二进制浮点参数 x 分成指数 e 和尾数 m，这样 x = m * 2 e，因此 log(x) = log(2) * e + log(m)。选择 m 使其接近统一，因为这提供了有效的近似值，例如 log(m) = log(1+f) = log1p(f) by 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>
#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;
}
#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
}
我矢量化了@njuffa 的答案。自然对数，适用于 AVX2：
inline __m256 mm256_fmaf(__m256 a, __m256 b, __m256 c){
return _mm256_add_ps(_mm256_mul_ps(a, b), c);
}
//
//
// 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;
}