特定舍入模式的意外结果

Question

首先,抱歉标题不好:/

我试图在计算三对角对称矩阵的特征值时重现论文的结果.我通过分别舍入到正负无限来确定一些值'上限和下限'.

我不是每次都改变舍入模式,而是使用'技巧':fl - (y)=-fl⁺(-y),其中fl - (y)是使用负无限舍入模式时的y值fl⁺(y)是使用舍入模式加上无穷大时的y值.所以,我在C中有以下代码:

fesetround(FE_UPWARD);
first =  - (-d[i] + x);
second =  ( - (( e[i-1]*e[i-1] ) / a_inf ));
a_inf = first + second;

first = d[i] - x;
second = - ( ( e[i-1]*e[i-1] ) / a_sup );
a_sup = first + second;

除了一个例子,其中a_inf给了我正确的结果,但是a_sup给出了错误的结果,尽管第一个和第二个变量似乎都具有相同的值,但它工作正常.

但是,如果我喜欢这样:

  fesetround(FE_UPWARD);
  first =  - (-d[i] + x);
  second =  ( - (( e[i-1]*e[i-1] ) / a_inf ));

  fesetround(FE_DOWNWARD);
  first =  - (-d[i] + x);
  second =  ( - (( e[i-1]*e[i-1] ) / a_sup ));

我得到了正确的结果.所以,如果我使用技巧fl - (y)=-fl⁺(-y),我得到正确的结果,如果我改变舍入模式并使用原始表达式我得到错误的结果.知道为什么吗？

在这两种情况下,变量的第一个和第二个值如下:

first  1.031250000000000e+07,  second -1.031250000000000e+07
first  1.031250000000000e+07,  second -1.031250000000000e+07

并且a_inf和a_sup的正确值分别为-1.862645149230957e-09和+ 1.862645149230957e-09,但在第一种情况下a_sup = 0,这是错误的

我猜它正在发生的是某种灾难性的取消,但我不知道在这种情况下如何解决它...

提前致谢!

Answer 1

您遇到的第一个问题是您仅使用“技巧”来四舍五入first到 -inf，而不是四舍五入second，因此它仍然四舍五入到 +inf。

第二个问题是，C 并没有就它如何进行浮点计算提供任何形式的保证。特别是，编译器可以自由地对操作进行重新排序和重新关联，因为它认为适合性能，即使这种重新排列可能会改变程序的舍入行为。所以当你说：

first =  - (-d[i] + x);

编译器可能会重新排列它并执行一次减法，而不是求反、加法、求反，这会反转舍入的方向。现在，有时您可以通过禁用所有优化来使事情按您期望的方式工作，但即使如此，也不能保证。