23

我的程序在函数中花费了 90% 的 CPU 时间std::pow(double,int)。准确性不是这里的主要关注点，所以我想知道是否有更快的替代方案。我想尝试的一件事是强制浮动，执行操作，然后返回双倍（还没有尝试过）；我担心这不是一种提高性能的可移植方式（大多数 CPU 本质上不是在双打上运行吗？）

干杯

4

4 回答 4

17

看起来Martin Ankerl有几篇关于这方面的文章，Optimized Approximative pow() in C / C++是一篇，它有两个快速版本，一个如下：

inline double fastPow(double a, double b) {
  union {
    double d;
    int x[2];
  } u = { a };
  u.x[1] = (int)(b * (u.x[1] - 1072632447) + 1072632447);
  u.x[0] = 0;
  return u.d;
}

它依赖于通过联合的类型双关语，这是 C++ 中未定义的行为，来自标准草案9.5 [class.union]：

在一个union中，任何时候最多只能有一个非静态数据成员处于活动状态，即任何时候最多可以有一个非静态数据成员的值存储在一个union中。[...]

但是包括 gcc 在内的大多数编译器都以明确的行为支持这一点：

从不同的工会成员那里阅读而不是最近写入的成员（称为“类型双关语”）的做法很常见。即使使用 -fstrict-aliasing，也允许类型双关，前提是通过联合类型访问内存

但这不是通用的，正如本文所指出的那样，正如我在我的回答中指出的那样，使用memcpy应该生成相同的代码并且不会调用未定义的行为。

他还链接到Java、C/C++ 和 C# 的第二个优化 pow() 近似值。

第一篇文章还在这里链接到他的微基准测试

于 2013-05-28T02:03:47.440 回答
12

根据您需要做什么，在对数域中操作可能会起作用——也就是说，您将所有值替换为它们的对数；乘法变成加法，除法变成减法，取幂变成乘法。但是现在加法和减法变得昂贵并且有些容易出错的操作。

于 2013-05-28T01:59:06.283 回答
6

你的整数有多大？它们在编译时已知吗？与 . 相比x^2，计算要好得多。注意：几乎所有整数幂的应用都涉及将某个数字提高到二次或三次幂（或负指数情况下的乘法逆）。在这种情况下使用是多余的。为这些小整数幂使用模板，或者只使用.x*xpow(x,2)pow()pow()x*x

如果整数很小，但在编译时不知道，比如在 -12 和 +12 之间，乘法仍然会击败pow()并且不会失去准确性。您不需要 11 次乘法来计算 x^12。四个就可以了。使用 x^(2n) = (x^n)^2 和 x^(2n+1) = x*((x^n)^2) 的事实。例如，x^12 是 ((x*x*x)^2)^2。两次乘法计算 x^3 (x*x*x)，再乘一次计算 x^6，最后一次乘法计算 x^12。

于 2013-05-28T02:10:43.507 回答
3

是的！如果您只需要 'y'/'n' 作为 long/int，则速度非常快，这样您就可以避免缓慢的 FPU FSCALE 函数。如果您只需要 'y'/'n' 作为 INT 的结果，这是 Agner Fog 的 x86 手动优化版本。我将它升级到 __fastcall/__declspec(naked) 以获得速度/大小，利用 ECX 传递“n”（对于 32 位 MSVC++，浮点数总是在堆栈中传递），所以我的调整非常小，主要是 Agner 的工作. 它是在 MS Visual VC++ 2005 Express/Pro 上测试/调试/编译的，所以应该可以滑入较新的版本。针对通用 CRT pow() 函数的准确性非常好。

extern double __fastcall fs_power(double x, long n);

// Raise 'x' to the power 'n' (INT-only) in ASM by the great Agner Fog!

__declspec(naked) double __fastcall fs_power(double x, long n) { __asm {
    MOV  EAX, ECX     ;// Move 'n' to eax
;// abs(n) is calculated by inverting all bits and adding 1 if n < 0:
    CDQ               ;// Get sign bit into all bits of edx
    XOR  EAX, EDX     ;// Invert bits if negative
    SUB  EAX, EDX     ;// Add 1 if negative. Now eax = abs(n)
    JZ   RETZERO      ;// End if n = 0
    FLD1              ;// ST(0) = 1.0 (FPU push1)
    FLD  QWORD PTR [ESP+4] ;// Load 'x' : ST(0) = 'x', ST(1) = 1.0 (FPU push2)
    JMP  L2           ;// Jump into loop
L1: ;// Top of loop
    FMUL ST(0), ST(0) ;// Square x
L2: ;// Loop entered here
    SHR  EAX, 1       ;// Get each bit of n into carry flag
    JNC  L1           ;// No carry. Skip multiplication, goto next
    FMUL ST(1), ST(0) ;// Multiply by x squared i times for bit # i
    JNZ  L1           ;// End of loop. Stop when nn = 0
    FSTP ST(0)        ;// Discard ST(0) (FPU Pop1)
    TEST EDX, EDX     ;// Test if 'n' was negative
    JNS  RETPOS       ;// Finish if 'n' was positive
    FLD1              ;// ST(0) = 1.0, ST(1) = x^abs(n)
    FDIVR             ;// Reciprocal
RETPOS:    ;// Finish, success!
    RET  4 ;//(FPU Pop2 occurs by compiler on assignment

RETZERO:
    FLDZ   ;// Ret 0.0, fail, if n was 0
    RET  4
}}
于 2019-07-28T22:30:18.290 回答