四舍五入以避免在随后的求和中四舍五入答案

【问题标题】：Rounding a double to avoid round off in subsequent summation四舍五入以避免在随后的求和中四舍五入
【发布时间】：2014-01-31 16:07:11
【问题描述】：

我该如何实现？

// round to the nearest double so that x + ref doesn't cause round off error
double round(double x, double ref) { }

这样

double x = ....;
double y = ....;

double x_new = round(x, y);
return x_new + y; // NO ROUND OFF!

换句话说 (y + x_new) - x_new 严格等于 y

【问题讨论】：

这通常不可能以您想要的方式进行。文本“到最近的双精度”建议您要在 x 附近返回一个 x_new。然而，如果x 大于y，那么任何在x 附近的某个数字对y 的加法都会强制将y 的一些低位排除在总和之外。例如，假设 x 是 264，y 是 1。满足要求的最接近的 double（使用 IEEE-754 64 位二进制）是 253-1，并且根本不在 2**64 附近。
@eric ，这是正确的。假设 x 归一化指数（由 frexp 返回）小于或等于 y/ref 归一化指数
我改了标题，因为“双舍入”已经有了特定的含义，“双”用作形容词。
请注意，检查(y + x_new) - x_new == y 与检查x_new + y 是否准确不同。后者暗示了前者，但不是必须的。示例：y=1e300, x_new=1。 (y + x_new) - x_new == 1e300 == y 但添加 y + x_new 不准确。

标签： c floating-point double ieee-754

【解决方案1】：

让我们假设x 和y 都是正数。

设S 为双精度和x + y。

有两种情况：

如果x ≤ y，则S - y 完全符合Sterbenz 引理。由此可见，加法(S - y) + y 是精确的（它精确地产生S，它是一个双精度数）。因此，您可以为x_new 选择S - y。不仅y + x_new 是准确的，而且它产生的结果S 与y + x 相同。
如果x > y，那么根据y 的有效位中设置的位数，您可能会遇到问题。例如，如果设置了y 的有效数字的最后一位，那么在y 的binade 之后的binade 中的任何数字z 都不能具有z + y 精确的属性。
李>

这个答案与that answer有模糊的关系。

【讨论】：

0b10000011000100100110111010010111100011010100111111100
@paolo_losi 我太慢地得出结论，我的解决方案略有偏差（很高兴你找到了一个例子）。准确地说，我想知道将我的解决方案更改为 2*y - (2*y - x) 是否会使其正确（仍然假设 0 ≤ x ≤ y）。至少现在我可以在您的示例上进行尝试，如果还没有证明，这将是一个线索。
@paolo_losi 我已经对我的答案进行了相当多的修改。它现在包含一个明显正确的解决方案，适用于我之前的答案想要的相同假设 (0 ≤ x ≤ y)。在x 大于y 的情况下，可能无法提供令人满意的x_new，例如，如果y 的有效位以一个设置位结束。

【解决方案2】：

可以直接翻译成 C 的 Python 中可能的解决方案

import math
from decimal import Decimal


def round2(x, ref):
    x_n, x_exp     = math.frexp(x)
    ref_n, ref_exp = math.frexp(ref)
    assert x_exp <= ref_exp
    diff_exp = ref_exp - x_exp

    factor = 2. ** (53 - diff_exp)

    x_new_as_int = int(round(x_n * factor))
    x_new_as_norm_float = float(x_new_as_int) / factor
    return math.ldexp(x_new_as_norm_float, x_exp)


x = 0.001
y = 1.0

assert (y + x) - x != y

x_new = round2(x, y)

assert (y + x_new) - x_new == y

print "x:", Decimal(x)
print "x_new:", Decimal(x_new)

print "relative difference:", (x_new/x - 1.)

【讨论】：

看起来这个解决方案依赖于 IEEE-754 64 位二进制文件。 C 没有将double 指定为IEEE-754 64-bit。您是否正在寻找在 C 中通用的解决方案？
@chux 它绝对应该是 C99 可移植的。我不知道 ieee-754 64 bit != c99 double。