5

首先，抱歉标题不好：/

我试图重现一篇论文关于计算三对角对称矩阵的特征值的结果。我分别通过四舍五入到正无穷和负无穷来确定一些值“上限和下限”。

我不是每次都更改舍入模式，而是使用“技巧”：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 和 second 的值如下：

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，这是错误的

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

提前致谢！

4

1 回答 1

1

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

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

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

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

于 2012-06-17T01:43:07.240 回答