std::pow 不同指数的行为非常不同

Question

我目前正在尝试优化一些代码，其中 50% 的时间花在 std::pow() 上。我知道指数将始终为正整数，而底数将始终为区间 (0, 1) 中的双精度数。为了好玩，我写了一个函数：

inline double int_pow(double base, int exponent)
{
    double out = 1.0;
    for(int i = 0; i < exponent; i++)
    {
        out *= base;
    }

    return out;
}

我正在编译：

> g++ fast-pow.cpp -O3 --std=c++11

我在 (0, 1) 之间生成了 1 亿个双精度，并比较了 (1) std::pow (2) 我自制的 int_pow 函数和 (3) 直接乘法的时序。这是我的计时程序的草图（这是一个非常快速的组合测试）：

void time_me(int exp, size_t reps)
{
    volatile double foo = 0.0;
    double base = 0.0;

    size_t i;
    for (i = 0; i < reps; ++i)
    {
        base = ((double) rand() / (RAND_MAX)) + 1;
        foo = pow(base, exp);
        // foo = int_pow(base, exp);
        // foo = base * base * base;
    }

    // check that the loop made it to the end
    std::cout << foo << "  " << i <<  std::endl;
}

int main()
{
    std::clock_t start;

    start = std::clock();
    time_me(3, 1e8);
    std::cout << "Time: " << (std::clock() - start) / (double)(CLOCKS_PER_SEC / 1000) << std::endl;

    return 0;
}

以下是我观察到的各种指数的时间安排：

• 0： std::pow 0.71s，int_pow 0.77s
• 2： std::pow 1.31s，int_pow 0.80s，直接乘法 0.86s
• 3： std::pow 6.9s (!!)，int_pow 0.84s，直接多 0.76s
• 5：类似于3：

我的问题

因此，我的问题是：

1. 为什么std::pow 的性能在大于 2 的幂时会出现如此严重的下降？
2. 当提前知道基数或指数类型时，是否存在现有的更快的幂函数？
3. 我忽略了一些完全明显的东西吗？对于已知整数指数的情况，我将通过直觉std::pow，并且不愿意错过一些完全微不足道的事情。

谢谢！！

Answer 1

std::pow() 是一个通用函数，旨在接受任何一对浮点值。它执行昂贵的计算，应该被认为是一个慢速函数。然而，显然，很多人滥用它来进行平方，因此 IBM Accurate Mathematical Library（由 glibc 使用）中的 pow() 实现针对该特定情况进行了优化：

sysdeps/ieee754/dbl-64/e_pow.c:

double
__ieee754_pow (double x, double y)
{
  ...
  ...
  if (y == 1.0)
    return x;
  if (y == 2.0)
    return x * x;
  if (y == -1.0)
    return 1.0 / x;
  if (y == 0)
    return 1.0;

如您所见，指数值 0、1 和 -1 也经过特殊处理，但这些至少在数学上是重要的特殊情况，而平方只是统计上重要的情况，否则不应该进行特殊处理）。 编辑：指数值0、1、2 和-1 是唯一允许用（快得多）算术运算表达std::pow(x,n) 而不会损失任何准确性的值。有关详细信息，请参阅this answer。因此，2 的指数值不仅仅是一个具有统计意义的案例。 结束编辑

如果您想快速替代std::pow() 以获得指数的非负整数值并且不关心轻微的精度损失，那么

1. 对于足够小的指数值，请使用您的 int_pow() 实现；
2. 否则，请使用exponentiation by squaring approach。

必须通过仔细的基准测试找到选择第一种和第二种方法的指数边界值。

Answer 2

switch (n)
{
case 0:
  return 1;
case 1:
  return x;
case 8:
  x*= x;
case 4:
  x*= x;
case 2:
  return x * x;
case 6:
  x*= x;
case 3:
  return x * x * x;
case 5:
  y= x * x; return x * y * y;
case 7:
  y= x * x * x; return x * y * y;
...
};