双精度浮点数z=x-y能否保证IEEE 754浮点数中的z+y等于x?

4
我有一个问题,可以简化为以下问题陈述:
给定一系列双精度浮点数,每个数在范围[0, 1e7]之间,修改最后一个元素,使得这些数字的总和恰好等于目标数。这些双精度浮点数已经在epsilon(1e-7)内总和等于目标数,但它们不相等。
以下代码是有效的,但它是否保证对符合第一句中描述的所有输入都有效?
public static double[] FixIt(double[] input, double targetDouble)
{
    var result = new double[input.Length];
    if (input.Length == 0) return result;

    double sum = 0;
    for (int i = 0; i < input.Length - 1; i++)
    {
        sum += input[i];
        result[i] = input[i];
    }

    double remainder = targetDouble - sum;
    result[result.Length - 1] = remainder;
    return result;
}

var arr1 = Enumerable.Repeat(Math.PI / 13, 13).ToArray();
var arr2 = FixIt(arr1, Math.PI);

Debug.Print(Math.PI.ToString("R")); //3.1415926535897931
Debug.Print(arr1.Sum().ToString("R")); //3.1415926535897922
Debug.Print(arr2.Sum().ToString("R")); //3.1415926535897931

之前的问题是关于修改第一个元素的,但是修改最后一个元素可以将问题简化为已知总和和目标值,只需要回答是否 last = target-sum 意味着 sum+last == target。(当然不包括 NaN,并且范围限制也意味着对 last 有一些限制,这可能有所帮助。)
关于实际问题:我们在各种情况下多次遇到了这个问题,但是目前我们正在尝试减少线性规划求解器(Coin-OR CBC)中由于数值不稳定性而出现的浮点误差。例如,有6个变量都必须在 [0,X] 范围内,这些变量的和也必须为 X。由于数值不稳定性,求解器偶尔会返回略微负数和不完全等于 X 的值。我们已经解决了负数问题,现在只是尝试解决等于 X 的问题。(是的,我们可能违反了某些约束条件来更改结果,但确保这些数字加起来等于 X 是更重要的,而其他约束条件并不那么重要。)

1
@PanagiotisKanavos 是的,因为数组的目的是要加和到一个特定的已知数字。因此,如果==不成立,至少必须更改一个数字。 - MineR
我已经编辑了最后一个元素。我觉得编辑第一个元素不起作用,因为浮点数相加的顺序很重要。通过更改最后一个元素,我认为更有可能奏效。 - MineR
@PeterCordes 不,这不是一项家庭作业任务。这是一个真实世界的问题。这是否意味着通过更改最后一个元素而不是第一个元素来解决它? - MineR
@MineR 无论您选择更改哪个元素,算术运算(以及数据类型的后果)都是相同的。 - Marc Gravell
@PeterCordes 我承认通过更改最后一个元素,你可能可以更接近目标,因为任何中间的大舍入问题都已经“完成”。假设总和将始终按严格向前的顺序完成!但是:根本问题仍然存在;很难保证任何特定水平的“正确性”(至少在人们认为舍入应该如何工作方面)。 - Marc Gravell
显示剩余5条评论
3个回答

7

z = x-y;不能保证z+y == x,并且对于寻找一个使得z+y == xz并不总是有解的。以下是一个证明。

我们假设使用IEEE-754二进制浮点数算术,并采用最近舍入、平局向偶数的规则。基本的64位格式被使用,但结果适用于其他格式。请注意,64位格式使用53位有效数字,这意味着只能表示53个或更少的二进制有效数字。

考虑一个目标值x等于1+2−52。让y等于2−53。然后,在z = x-y;执行之后,z+y == x的结果为false。具体的算术细节如下:

  • z = x-y;z设置为1,然后z+y产生1,小于x
  • 如果我们将z增加到下一个可表示的数字1+2−52,那么z+y产生1+2−51,大于x
  • 因此,没有一个z的值使得z+y == x为真。

细节如下:

xy的数学结果为1+2−53。由于它具有54个有效位(从20到2−53),因此它无法表示,并且x-y的计算结果必须舍入。最接近的两个数字是1和1+2−52。平局向偶数规则将前者1作为其有效数字的最低位为0,而对于1+2−52,其低位为1。

因此,z = x-y;z设置为1。

然后,z+y的数学结果为1+2−53。与上面一样,这被四舍五入为1,因此z+y的计算结果为1。因此,z+y == x将1与1+2−52进行比较,并产生false。

此外,任何z的值都无法使比较为真。如果我们将z增加最小可用步长,从1到1+2−52,则z+y的数学和为1+2−52+2−53。这是介于两个可表示数字1+2−52和1+2−51之间的中间值。前者的低位为1,后者的低位为0,因此计算出的z+y结果为1+2−51,显然不等于1+2−52。浮点加法是弱单调的,因此没有任何z的值会产生z+y等于1+2−52的结果。

非常感谢 - 我测试了这个例子,但至少在我的硬件上它不起作用,特别是在C#中。 - MineR
3
@MineR: 使这个结果成立的情况非常罕见。若干其他条件可能会排除这种情况并解决你的问题。例如,如果y至少是x的一半,则x-y是精确的,且z+y==x必定为真。如果y比2^-52乘以x还要高至少一个二进制阶(指数值),那么有足够的余地使得z可以满足z+y==x即使x-y不精确。而且,如果x的低位是零,我认为在此答案中所述的情况就不会发生。因此,如果你的实际问题可以更具体地确定,可能会有解决方案。 - Eric Postpischil
是的,看起来我们可以通过像你建议的那样去掉低位来使修复工作,因为在我们的情况下,这种精度水平并不重要。 - MineR
1
@MineR:我认为这个问题只会在数量级极端不同的情况下发生,这可能不会在您的实际用例中发生。调整较小的数字应该总是能够达到所需的目标,因此也许可以使用if(last > sum_of_previous*(1ULL<<40)),调整倒数第二个元素,因为这将使您对最终总和的较小输入进行微调。(但要注意,这可能意味着微小元素的相对变化很大,甚至可能使其变为负数。也许您可以使用该条件来查找数值上有问题的情况...2^40只是我编造的一个数字) - Peter Cordes

3
不,它并不会。这里有一个具体的反例;虽然是用Python编写的,但你可以很容易地在C#中重复进行同样的实验:
>>> x = 0.24999916553497312
>>> y =  1.0000153779983518
>>> z = -0.7500162124633787
>>> z == x - y
True
>>> z + y == x
False

以下是一个小的反例,xyz都是正数:

>>> x = 0.4500000000000001
>>> y = 0.20000000000000004
>>> z = 0.2500000000000001
>>> z == x - y
True
>>> z + y == x
False

原帖中将事情限制为非负数(包括 z,因为目标总和必须在实际总和的小幅度范围内)。负数比Eric的情况更容易丢失精度的抵消,Eric的情况是最后一个数字 (z) 比前面数字(y) 的总和大几个数量级。 - Peter Cordes
@PeterCordes 没问题;添加了一个反例,其中包含“小”但是正值。这些值都足够接近OP所要求的。 - alias

1

浮点算术本质上是不精确的(除非你只涉及整数(更正:最大可达253,即9007199254740992)),你将总是会有舍入差异。如果你想要舍入结果与人类期望一致:请使用decimal而不是double。如果你使用decimal进行相同的操作,对于任何在十进制位上不具备病态特征的数字集合,它都可以正确工作。


1
将会添加的原因是decimal是基于10进制,而double是基于2进制。基于2进制的数字无法准确表示所有基于10进制的数字。 - Zer0
4
这并没有回答这个问题。原帖的问题是在问,对于IEEE-754 binary64 double,在限制条件下(最重要的是非负数,因此不会发生灾难性的抵消),这个算法是否总是有效的,而且它并没有说“期望”的目标不能小于其他元素的总和,需要一个负的第一个元素...例如,一个目标为1e-200的值将无法通过改变第一个元素来实现,因为它没有足够的有效数字位。 - Peter Cordes
6
这个回答与所问问题无关,它声称浮点运算总是有差异是错误的。 - Eric Postpischil
1
整数是二进制浮点数可以得到精确结果的一种情况。还有其他情况,包括微不足道的0.25 + 0.25 == 0.5。(二进制浮点格式可以精确表示分母为2的幂的分数)。总之,作为一般指导,这很好,但这个问题正在问一个具体的问题。 - Peter Cordes
5
这个答案提供了错误的陈述并不是卖弄学问。浮点数实践者在可用的情况下利用精确性和浮点数算术的其他属性来设计和证明代码。像这样误导性和错误的回答会对此造成损害。 - Eric Postpischil
显示剩余7条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接