C语言浮点数相等判定

等价关系

按照离散数学的等价关系叙述，集合 \(X\) 上的关系 \(R \subset (X, X)\) 如果满足

1. 反身性: \((x, x) \in R, \forall x \in X\)，自己等于自己
2. 对称性: 如果 \((x, y) \in R\)，则 \((y, x) \in R\)，我等于你，则你也等于我
3. 传递性: 如果 \((x, y) \in R\) 且 \((y, z) \in R\)，则 \((x, z) \in R\)，我等于你且你等于他，则我等于他

这里考虑 IEEE754 中双精度浮点数（即 C 语言中的 double）相等判定，按照 C 语言中 == 判定，则它确实是等价关系，注意 0.0 == -0.0 成立(虽然它们的二进制中符号位不同）

由于机器舍入误差的存在，浮点数使用 == 作为相等判定条件无法满足现实中的要求。例如 0.1 + 0.2 == 0.3 是不成立的（它们在二进制中均为无限循环小数）。

因此一般用绝对误差判定：fabs(x - y) < eps 时认为“相等”，其中 eps 是一个指定的小数。显然它满足反身性和对称性，但不满足传递性，因此此判定不是等价关系：例如取 x = eps/2, y = 0, z = -eps/2。

但绝对误差判定存在另一个问题：\(x = 1.0 + 2^{53} + 1.0\) 不能被精确保存，得到的是 \(2^{53}\)，而 \(y = 2^{53} + 2\) 是被精确保存的。那么按照上述判定，它们不相等。

因此再引入相对误差的判定：fabs(x - y) < eps1 || fabs(x - y) < eps2 * fabs(x) || fabs(x - y) < eps2 * fabs(y)，其中 eps1 和 eps2 均为指定的小数且一般相等。显然它满足反身性和对称性，但不满足传递性：例如 \(x = (1 - eps2) \cdot 2^n\), \(y = 2^n\), \(z = (1 + eps2) \cdot 2^n\)。

于是有如下判定的 C 代码（仅满足反身性和对称性），全篇最终代码

#include <math.h>
#include <float.h>

int fIsEqual(double x, double y)
{
	double z = fabs(x - y);
	return z < DBL_EPSILON * 10 || z < DBL_EPSILON * 10 * fabs(x) || z < DBL_EPSILON * 10 * fabs(y);
}

1. 其中 DBL_EPSILON 是 1 与比 1 大的最小浮点数之间的差值，按照 IEEE754 标准值为 \(2^{-52}\)
2. 这里的 10 为经验值，用于误差累计后依然能判过去，另外 DBL_EPSILON * 10 一定会被编译器优化成一个常量
3. 返回值用 int 不用 bool 是因为 C23 才定义 bool
4. 上述判定存在一个问题，\(x = 2^{52} + 1\), \(y = 2^{52} - 9\) 时 fIsEqual(x, y) = 1 然而 fIsEqual(x - 1, y - 1) = 0, 即对线性运算不封闭，这是可以接受的，因为浮点数本来就不满足结合率，加法本身就有大数吃小数的现象。再者仅用绝对误差作为判据也无法满足对线性运算封闭的性质。

牺牲对称性换来所谓性能提升是没有意义的

上述代码中需要做三次判定才能确定不相等，考虑到如果 fabs(x - y) < DBL_EPSILON * 10 * fabs(y) 时, fabs(x) 与 fabs(y) 差距较小，再判断 fabs(x - y) < DBL_EPSILON * 10 * fabs(x) 有损性能且意义不大。于是删掉中间的判定条件，并用 x 的值覆盖 z，就推出下面所谓的"优化版"（仅满足反身性，不满足对称性和传递性，假优化请勿使用！！！）

#include <math.h>
#include <float.h>

int fIsEqual(double x, double y)
{
	x = fabs(x - y);
	return x < DBL_EPSILON * 10 || x < DBL_EPSILON * 10 * fabs(y);
}

但此时当 \(x = 4503599627370487\) 且 \(y = 4503599627370497\) 时 fabs(x - y) = 10，此时 fIsEqual(x, y) 成立，但 fIsEqual(y, x) 不成立，这是难以让人接受的，详见下图