C++，如何优化浮点算术运算？

Question

在 x86 架构上测试极限情况下的简单算术运算时，我观察到一个令人惊讶的行为：

const double max = 9.9e307; // Near std::numeric_limits<double>::max()
const double init[] = { max, max, max };

const valarray<double> myvalarray(init, 3);
const double mysum = myvalarray.sum();
cout << "Sum is " << mysum << endl;             // Sum is 1.#INF
const double myavg1 = mysum/myvalarray.size();
cout << "Average (1) is " << myavg1 << endl;    // Average (1) is 1.#INF
const double myavg2 = myvalarray.sum()/myvalarray.size();
cout << "Average (2) is " << myavg2 << endl;    // Average (2) is 9.9e+307

（在 release 模式下使用 MSVC 以及通过 Codepad.org 的 gcc 进行测试。MSVC 的调试模式将平均值 (2) 设置为 #INF。）

我预计平均值 (2) 等于平均值​​ (1)，但在我看来，C++ 内置除法运算符已被编译器优化，并以某种方式阻止了累积达到 #INF。
简而言之：大数字的平均值不会产生#INF。

我在 MSVC 上使用 std 算法观察到相同的行为：

const double mysum = accumulate(init, init+3, 0.);
cout << "Sum is " << mysum << endl;             // Sum is 1.#INF
const double myavg1 = mysum/static_cast<size_t>(3);
cout << "Average (1) is " << myavg1 << endl;    // Average (1) is 1.#INF
const double myavg2 = accumulate(init, init+3, 0.)/static_cast<size_t>(3);
cout << "Average (2) is " << myavg2 << endl;    // Average (2) is 9.9e+307

（不过，这次 gcc 将平均值 (2) 设置为 #INF: http://codepad.org/C5CTEYHj。）

1. 有人愿意解释一下这种“效果”是如何实现的吗？
2. 这是一个“功能”吗？或者我可以认为这是一种“意外行为”而不是简单的“令人惊讶”吗？

谢谢

Answer 1

只是一个猜测，但是：Average (2) 可能是直接在浮点寄存器中计算的，浮点寄存器的宽度为 80 位，并且比内存中双精度数的 64 位存储更晚溢出。你应该检查你的代码的反汇编，看看是否确实如此。

Answer 2

这是一种功能，或者至少是有意的。 基本上，x86 上的浮点寄存器有更多 精度和范围比双精度（15 位指数，而不是 11、64 位 matissa，而不是 52）。 C++ 标准允许 对中间值使用更高的精度和范围，以及 几乎任何英特尔的编译器都会在某些情况下这样做 情况;性能差异显着。 是否获得扩展精度取决于何时 以及编译器是否溢出到内存。 （将结果保存在 命名变量需要编译器将其转换为实际的 双精度，至少根据标准。） 我见过的更糟糕的情况是一些基本上可以做到的代码：

return p1->average() < p2->average()

，average() 在内部表上做你期望的事情 在数据中。在某些情况下，p1 和 p2 实际上会指向 到同一个元素，但返回值仍然为真； 函数调用之一的结果将溢出到 内存（并截断为double），其他的结果 保留在浮点寄存器中。

(该函数被用作sort的排序函数， 结果代码崩溃了，因为由于这种影响，它 没有定义足够严格的排序标准，并且 sort 超出传递给它的范围时的代码。）

Answer 3

在某些情况下，编译器可以使用比声明类型所隐含的类型更宽的类型，但 AFAIK，这不是其中之一。

因此我认为我们的效果类似于Gcc bug 323 的效果，在不应该使用额外精度的情况下。

x86 有 80 位 FP 内部寄存器。虽然 gcc 倾向于以最大精度使用它们（因此是错误 323），但我的理解是 MSVC 将精度设置为 53 位，即 64 位中的两倍。但延长的显着性并不是 80 位 FP 的唯一区别，指数的范围也增加了。而 IIRC，x86 中没有强制使用 64 位范围的设置。

键盘现在似乎无法访问，否则我会在没有 80 位长 double 的架构上测试您的代码。

Answer 4

g++ -O0 -g -S test.cpp  -o test.s0
g++ -O3 -g -S test.cpp  -o test.s3

比较 test.s[03] 表明确实 valarray::sum 甚至没有被再次调用。我没看很久，但下面的 sn-ps 似乎是定义片段：

    .loc 3 16 0 ; test.s0

    leal    -72(%ebp), %eax
    movl    %eax, (%esp)
    call    __ZNKSt8valarrayIdE3sumEv
    fstpl   -96(%ebp)
    leal    -72(%ebp), %eax
    movl    %eax, (%esp)
    call    __ZNKSt8valarrayIdE4sizeEv
    movl    $0, %edx
    pushl   %edx
    pushl   %eax
    fildq   (%esp)
    leal    8(%esp), %esp
    fdivrl  -96(%ebp)
    fstpl   -24(%ebp)

    .loc 3 17 0

对

    .loc 1 16 0 ; test.s3
    faddl   16(%eax)
    fdivs   LC3
    fstpl   -336(%ebp)
LVL6:
LBB449:
LBB450:
    .loc 4 514 0