C++中非常快速的对数(自然对数)函数?
我们发现代替std::sqrt
(Timing Square Root)和std::exp
(Using Faster Exponential Approximation)的一些技巧,但我找不到替代std::log
的东西。它是我程序中循环的一部分,它被多次调用,而exp和sqrt被优化,因此英特尔VTune现在建议我优化std::log
,之后似乎只有我的设计选择会受到限制。C++中非常快速的对数(自然对数)函数?
很多谢谢。
这取决于你需要的准确度。通常会通过调用log来了解数字的大小,通过检查浮点数的指数字段,您可以基本上免费进行此操作。这也是你的第一个近似值。我将为我的“基本算法”一书提供一个插件,它解释了如何从第一原理实现标准库数学函数。
我正在寻找真正的数学应用的自然对数,不需要双精度,浮点精度甚至10-3,10-4会很好 – user3091460
链接或书籍参考没有引用的相关部分不是答案 – BeyelerStudios
在开始设计和部署性能超越函数的自定义实现之前,强烈建议在算法级别以及通过工具链进行优化。不幸的是,我们没有任何有关此处优化的代码的信息,我们也没有关于工具链的信息。
在算法级别,检查是否所有对超越函数的调用都是真正必需的。也许有一个数学转换需要较少的函数调用,或将超越函数转换为代数运算。任何超越函数调用都可能是多余的,例如因为计算是不必要地切入和切出对数空间?如果精度要求适中,整个计算是否可以单精度执行,全部使用float
而不是double
?在大多数硬件平台上,避免使用double
计算可以显着提高性能。
编译器倾向于提供影响数字密集代码性能的各种开关。除了将通用优化级别提高到-O3
之外,通常还有一种方法可以禁用非正常支持,即打开清零或FTZ模式。这在各种硬件平台上具有性能优势。此外,通常会有一个“快速数学”标志,其使用会导致准确度略有下降,并消除了处理特殊情况(如NaN和无穷大)的开销,以及处理errno
。一些编译器还支持自动向量化代码,并附带一个SIMD数学库,例如英特尔编译器。
一个对数函数通常涉及分离二进制浮点参数x
成指数e
和尾数m
的自定义实现,使得x = m * 2
e
,因此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;
}
Thks,但是我不能找到在win10上使用Msvc 15的int_as_float和float_as_int。我发现它是cuda的一部分,但没有下载完整的软件包。 – user3091460
@ user3091460这些功能是机器特定功能的*抽象*。作为第一步,你可以简单地使用'memcpy()',例如'float __int_as_float(int32_t a){float r; memcpy(&r,&a,sizeof(r)); return r;}'一个好的编译器可能会适当地优化它,但取决于你所针对的硬件(你没有透露),可能有更好的方法,可能涉及到内在函数或内联汇编。 – njuffa
@ user3091460和njuffa:由于XMM寄存器用于标量/矢量浮点数和矢量整数,x86的最佳asm可能会使用SSE2整数指令对浮点数进行整数操作。所以你应该'_mm_set_ss(your_float)'和'_mm_castps_si128'来获得你可以操作的'__m128i'。 (这可能会浪费一条使xmm寄存器的高位为零的指令[由于内部函数的设计限制。](http://*.com/q/39318496/224132))。从一个整数寄存器获取浮点数的MOVD也可能很好。 –
两个upvotes –
稀释是 - 在精确度和性能的问题 - 但没有说明什么精度是可以接受的,或者是被审判我不要以为你会得到'答案' – UKMonkey
浮点精度就足够了。我试图从log2开始,然后转回,但非常快的log2只是输出一个int,导致很差的近似值。也尝试使用ln(x)在t = 0时是t-> x^t的导数的事实,但它的好处不在于计算。 – user3091460