C语言浮点数相等判定

等价关系

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

  1. 反身性: (x,x)R,xX,自己等于自己
  2. 对称性: 如果 (x,y)R,则 (y,x)R,我等于你,则你也等于我
  3. 传递性: 如果 (x,y)R(y,z)R,则 (x,z)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+253+1.0 不能被精确保存,得到的是 253,而 y=253+2 是被精确保存的。那么按照上述判定,它们不相等。

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

于是有如下判定的 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 标准值为 252
  2. 这里的 10 为经验值,用于误差累计后依然能判过去,另外 DBL_EPSILON * 10 一定会被编译器优化成一个常量
  3. 返回值用 int 不用 bool 是因为 C23 才定义 bool
  4. 上述判定存在一个问题,x=252+1, y=2529fIsEqual(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=4503599627370487y=4503599627370497fabs(x - y) = 10,此时 fIsEqual(x, y) 成立,但 fIsEqual(y, x) 不成立,这是难以让人接受的,详见下图

posted @   cuzperf  阅读(88)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 没有源码,如何修改代码逻辑?
· PowerShell开发游戏 · 打蜜蜂
· 在鹅厂做java开发是什么体验
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战
点击右上角即可分享
微信分享提示