为什么有些算术运算比平常花费更多时间？

Question

我发现在使用小精度的浮点数进行算术运算时存在不寻常的计算时间。以下简单的代码展示了这种行为：

#include <time.h>
#include <stdlib.h>
#include <stdio.h>

const int MAX_ITER = 100000000;

int main(int argc, char *argv[]){
    double x = 1.0, y;
    int i;
    clock_t t1, t2;
    scanf("%lf", &y);
    t1 = clock();
    for (i = 0; i < MAX_ITER; i++)
        x *= y;
    t2 = clock();
    printf("x = %lf\n", x);
    printf("Time: %.5lfsegs\n", ((double) (t2 - t1)) / CLOCKS_PER_SEC);
    return 0;
}

下面是程序的两个不同运行结果：

• y = 0.5时：

  x = 0.000000
  时间：1.32000秒

• y = 0.9时：

  x = 0.000000
  时间：19.99000秒

我使用以下规格的笔记本电脑测试代码：

• CPU: Intel® Core™2 Duo CPU T5800 @ 2.00GHz × 2
• 内存: 4 GB
• 操作系统: Ubuntu 12.04 (64 bits)
• 型号: Dell Studio 1535

有人可以详细解释一下这种行为发生的原因吗？我知道当y = 0.9时，x值下降得比y = 0.5慢，所以我认为问题直接与此有关。

- jace

---

在第二种情况下，您可能会在达到0之前获得更多的非规范化数。 - Paul R

请参考此答案，了解denormals的性能影响。 - Daniel Fischer

Answer 1

Denormal（或者说是subnormal）数字通常会影响性能。根据第二个示例，缓慢地收敛到0将生成更多的subnormals。可以在这里和这里阅读更多信息。对于更严肃的阅读，请查看经常引用的（非常深奥的）计算机科学家应该了解的浮点运算知识。

从第二个来源：

在IEEE-754标准下，浮点数以二进制形式表示为：Number = signbit * mantissa * 2exponent。同一个数可能有多种不同的表示方法，以十进制为例，数字0.1可以表示为1*10-1或0.1*100或0.01*10。该标准规定数字总是以第一位为1存储。在十进制中，这对应于1*10-1的示例。现在假设可以表示的最低指数为-100。因此，正常形式下可表示的最小数字为1*10-100。但是，如果我们放宽第一位必须是1的限制，则可以在相同的空间内实际上表示更小的数字。以十进制为例，我们可以表示0.1*10-100。这被称为次规格化数。拥有次规格化数的目的是平滑最小正常数和零之间的差距。非常重要的是要意识到，次规格化数的精度比正常数低。事实上，它们为了更小的尺寸而交换了精度。因此，使用次规格化数进行计算将不会具有与正常数进行计算时相同的精度。因此，对次规格化数进行大量计算的应用程序可能值得调查，以查看重新缩放（即通过某些缩放因子乘以数字）是否会产生更少的次规格化数和更准确的结果。

我本来想自己解释一下，但是上面的解释已经写得非常清晰简洁了。

- David Titarenco

Answer 2

您获得可测量的差异不是因为在数学上0.9^n收敛比0.5^n慢，而是因为在IEEE-754浮点实现中，它根本不会收敛到0。

IEEE-754表示中最小的正double是2^-1074，最小的正常值是2^-1021，因此使用y = 0.5时，循环会遇到53个次正常数。一旦到达最小的正次正常数，下一个乘积将是2^-1075，但由于默认舍入模式“向最后一位舍入为零”，它会被舍入为0。（IEEE-754浮点数的表示和默认的舍入模式“向最后一位舍入为零”几乎普及了标准消费硬件，即使标准没有完全实现。）从那时起，您就有了一个乘法0*y，这是一种普通的浮点乘法（即使y是一个次正常数，这也很快）。

对于0.5 < y < 1，一旦达到（正）次正常范围的较低端，x*y的结果再次舍入为x的值（对于y = 0.9，迭代的固定点是5*2^-1074）。由于在几千次迭代后就达到了这一点（0.9^7 < 0.5），因此您基本上是将一个次正常数与非零数相乘整个循环。在许多处理器上，这样的乘法无法直接处理，必须在微码中处理，这会慢得多。

如果速度比IEEE-754语义更重要（或者由于其他原因不希望使用该语义），许多编译器提供选项以禁用该行为，并在硬件支持时将次正常数刷新为0。我在gcc的man页面中找不到明确的选项，但-ffast-math在这里起了作用。

- Daniel Fischer

---

感谢您证明了通过精确计算分析问题是可能的，并提供了编译器选项进行优化。您能否再解释一下为什么在子规范范围内，当0.5 <y <1时，乘积x * y会四舍五入到x？ - jace

如果 x 是最小的正次规范数，且 y > 0.5，那么 x*y 的确切值大于 x/2，因此它不是一个平局，并且会被舍入到更接近 x 和 0 中的一个，即 x。如果 y 足够大于 0.5，那么在 x 达到最小的次规范数 s 之前就会发生这种情况，例如对于 y = 0.75，你有 (2*s)*y = 1.5*s 在数学上，这产生了一个平局 - 它恰好处于 s 和 2*s 之间，因此它再次被舍入为 2*s，其最后一位为 0。对于 y = 0.9 和 x = 5*s，你将得到 4.5*s，但最接近的 double 稍微大于 0.9，所以它会向上舍入。 - Daniel Fischer

谢谢，我之前对这些底层细节不太了解。 - jace