C（++）中双打的排序/比较不稳定？答案

【问题标题】：Sorting/comparison of doubles in C(++) not stable?C（++）中双打的排序/比较不稳定？
【发布时间】：2012-07-05 01:02:42
【问题描述】：

我遇到了一个非常奇怪的双打问题。我有一个按降序排序的浮点数（双精度）列表。后来在我的程序中，我发现它们不再完全排序。例如：

0.65801139819
0.6545651031    <-- a
0.65456513001   <-- b
0.64422968678

中间的两个数字翻转。有人可能会认为这个问题在于数字的表示，它们只是打印错误。但是我使用我用来对它们进行排序的相同运算符将每个数字与前一个数字进行比较 - 没有转换为基数 10 或类似的情况：

double last_pt = 0;
for (int i = 0; i < npoints; i++) {
  if (last_pt && last_pt < track[i]->Pt()) {
    cout << "ERROR: new value " << track[i]->Pt()
         << " is higher than previous value " << last_pt << endl;
  }
  last_pt = track[i]->Pt();
}

在排序过程中比较值

bool moreThan(const Track& a, const Track& b) {
  return a.Pt() > b.Pt();
}

我确保它们始终是双精度的，而不是转换为浮点数。 Pt() 返回一个双精度值。列表中没有NaN，排序后我不碰列表。

为什么会这样，这些数字有什么问题，以及（如何）我可以对数字进行排序以使其保持排序状态？

【问题讨论】：

AFAIK，float 几乎总是 32 字节，而 double 是 64 字节，在两种架构上。
谢谢，这就是我的假设。所以至少它可能不是那样......
会不会是排序不对？你用什么来排序？
这是一个错字:-)。我删除了问题的那一部分，因为我认为它太令人困惑并且与原始问题没有太大关系。
请添加一个可以重现您的错误的工作程序

标签： c++ floating-point

【解决方案1】：

您正在比较函数的返回值。浮点数返回值在浮点寄存器中返回，该寄存器具有更高的精度比双倍。当比较两个这样的值（例如a.Pt() > b.Pt()）时，编译器将调用其中一个函数，存储返回 double 类型的未命名临时值中的值（因此四舍五入结果到double)，然后调用另一个函数，并比较它的结果（仍然在浮点寄存器中，而不是四舍五入到 double) 与存储的值。这意味着您最终可以得到 a.Pt() > b.Pt() 和 b.Pt() > a.Pt() 或 a.Pt() > a.Pt() 的情况。这将导致sort 变得有点困惑。（正式地，如果我们在这里谈论std::sort，这会导致未定义的行为，我听说过它确实导致核心的情况转储。）

另一方面，您说Pt()“只返回一个双字段”。如果Pt() 不进行任何计算；如果只是：

double Pt() const { return someDouble; }

，那么这应该不是问题（前提是 someDouble 具有类型 double)。扩展精度可以表示所有可能的双精度准确的值。

【讨论】：

有趣的是，您对 Pt() 函数的看法是正确的 (1.usa.gov/MJbljh)。我刚刚再次检查，它确实执行了计算（但是，我确实希望调用 Pt() 总是返回相同的值）。
这仅适用于使用 x87 FPU 的代码。 x64 ABI 要求在 XMM 寄存器中传递浮点值。
@jdm 调用本身将始终返回相同的值。但是该值具有扩展精度，并且编译器将在进行第二次函数调用之前将其中一个返回值“溢出”到内存中。溢出到内存的值将四舍五入为双精度。解决方案是通过将它们分配给局部变量并比较局部变量，例如，将两个返回值强制为双精度。 double t1 = a.Pt(); double t2 = b.Pt(); return t1 > t2;.
返回值的扩展精度不会导致询问者报告的行为。这两个无序浮点值的双精度不同，因此它们的扩展精度也会不同（如果一个是双精度的，一个是扩展的，一个被转换为另一个的类型然后进行比较，它们也会有所不同）。这种差异会导致它们被正确排序。
@EricPostpischil 问题是将扩展精度值与相同的值进行比较，但该值已存储在双精度值中。我实际上已经看到了 a.Pt() > a.Pt() 的等价物会返回 true 的情况。这足以混淆std::sort。

【解决方案2】：

您确定没有将double 转换为float 吗？让我们看一下这两个数字的二进制表示：

0 01111111110 0100111100100011001010000011110111010101101100010101
0 01111111110 0100111100100011001010010010010011111101011010001001

double 有 1 位符号、11 位指数和 53 位尾数，而在float 中有 1 位符号、8 位指数和 23 位尾数。请注意，两个数字的尾数在前 23 位是相同的。

根据舍入方法，会有不同的行为。如果只是修剪 >23 位，这两个数字 float 是相同的：

0 011111110 01001111001000110010100 (trim: 00011110111010101101100010101)
0 011111110 01001111001000110010100 (trim: 10010010011111101011010001001)

【讨论】：

好观察！这应该可以帮助我缩小范围。这是一个巨大的代码库，不同的模块通过配置脚本连接：-/。这就是为什么我不能轻易地想出一个最小的例子。我有 90% 的把握，但不能排除这种可能性。
请注意，除非舍入模式已明确更改，否则它在 x86 和 x64 上是舍入到最近的，并且当转换为 float 时，第二个数字的有效位的 lsb 中应该有 1 （使用 gcc 和 icc 对 32 位和 64 位 Linux 目标进行测试，以及使用 VS2010 对 32 位和 64 位目标进行测试，均使用 x87 和 SSE 指令）。