浮点数的比较
首先，这个不算原创，原文是洋文的，我翻译了一下
写这个文章的人绝对是个大师，虽然知识并不是很深奥，不过想法真的很不错，值得学习
两个月前忽然看见的这篇文章，昨天仔细读了一遍，翻译了一下，原文在此
http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
这是我翻译的，里面有些地方翻译的不太对，有些地方我也不懂，你们谁要是看懂了就告诉我
浮点数的比较
Bruce Dawson
两数完全相等
整型数是精确的，浮点数并不是那样。举个简单的例子，0.2这个数无法用二进制的浮点数精确表示，并且有限的精度意味着微小的误差在经过多次操作后会改变结果。不同的编译器和CPU架构使用不同的精度存储临时数据，所以计算结果取决于你的运行环境。所以当你计算出结果并与你预期的结果进行比较时，会发现它们也许大不一样。
换句话说，如果你在计算中做如下比较：
if (result == expectedResult)
这个比较的结果几乎不可能为真。如果结果确实为真，那说明这个结果是不稳定的——输入数据、编译器、CPU的微小改变也许会使程序产生截然不同的结果。
比较e-绝对误差
既然浮点数的计算引入了一点点不确定性，那么我们忽略它，方法就是判断两个数是否足够近。不管是通过什么方法（误差分析、测试、甚至是瞎蒙），如果你认为两个数只要相差不到0.00001就算相等的话，你可以将程序改成这样：
if (fabs(result - expectedResult) < 0.00001)
这个最大的误差值称为（epsilon）
绝对误差可能有点用处，不过大家不常用它。大部分时候，误差是指的一个百分比。大家不用绝对误差的原因是，一个1.0的误差，你不知道结果是一百万还是0.1。
对于有限精度的浮点数还有一个更杯具的问题，就是如果绝对误差太小的话，有限的精度可能根本无法表示这个误差。
举个例子，假如你的预期计算结果是10000，由于有限精度的浮点数天生的缺陷，你也许算出来的数不是严格的10000，说不定二进制的浮点数只有一位偏差，如果用4字节表示一个浮点数的话，假设你得到了这个数+10000.000977，写成程序如下：
float expectedResult = 10000;
float result = +10000.000977; // 最接近10000的 4字节浮点数
float diff = fabs(result - expectedResult);
误差是0.000977，这是我们预定的（epsilon）的97.7倍。所以我们的系统说：这俩数不够接近，尽管它们的二进制码只有最后一位不同！使用0.00001作为（epsilon）没有任何用处，它跟判断完全相等是一样的，唯一的不同是这么干更麻烦。
绝对误差有一个确定的值，如果预期的结果的范围是确定的，那么检查绝对误差简单而有效。只需要确定绝对误差值比浮点数能表达的最小数字大就行，当然，这也取决于选用什么类型的float。
比较e-相对误差
0.00001的误差对于一个与1差不多大的数而言，比较合适，对于0.00001这么大的数，则显得太大，对于10000这种数量级的数，又显得太小。一个更好的方法，就是不依赖于结果的大小，采用相对误差。
相对误差就是用绝对误差与预期结果进行比较。一个简单的方法就是如下：
relativeError = fabs((result - expectedResult) / expectedResult);
如果结果是99.5，而预期结果是100，那么相对误差就是0.005。
有时候我们并没有一个“预期”的结果，我们只有两个数让我比较来看它们是相等。比如下面的这个函数：
// 非最优 AlmostEqual function – 不推荐。
bool AlmostEqualRelative(float A, float B, float maxRelativeError)
{
if (A == B)
return true;
float relativeError = fabs((A - B) / B);
if (relativeError <= maxRelativeError)
return true;
return false;
}
这个maxRelativeError变量限定了我们所能容忍的最大误差。如果我们想要99.999%的精确度，那么maxRelativeError就应该是0.00001。
最初的A == B比较看起来有些奇怪：如果A == B那么relativeError就应该等于0？显然这么干是不行的。如果A、B都等于0那么relativeError是0.0 / 0.0。这个是个无效数NAN。在<=比较中，NAN的返回值永远不为真，所以当A、B都等于0时函数会返回假。（在一些平台上NAN的比较也许不是这么处理的，函数也许会返回真，但是这时所有NAN的输入都会导致函数返回真，这会带来更多的麻烦）
这个函数还有个问题就是AlmostEqualRelative(x1, x2, epsilon)的返回值未必等于AlmostEqualRelative(x2, x1, epsilon)的返回值，因为它只用第二个变量作为除数。下面给出一个升级版，它用绝对值大的数作为除数：
// 这个函数比AlmostEqual稍微强点 – 不过依然不推荐
bool AlmostEqualRelative2(float A, float B, float maxRelativeError)
{
if (A == B)
return true;
float relativeError;
if (fabs(B) > fabs(A))
relativeError = fabs((A - B) / B);
else
relativeError = fabs((A - B) / A);
if (relativeError <= maxRelativeError)
return true;
return false;
}
知道现在，我们的函数依然不完美。这个函数当输入参数为0附近的数时，效果并不好。一个接近0的正数和一个接近0的负数也许可以认为相等，但是说不定会得到2.0的相对误差（误差比绝对值都大）。如果你想判断两个符号相反但是都接近0的数字是否相等，可以再引入一个变量。当absoluteError和AbsoluteError任何一个小于各自的规定值时，就会返回真。一个标准的maxAbsoluteError应该非常小- FLT_MAX或者更小（这句话谁看懂了？），取决于你的系统是否支持非正常情况。
// 比AlmostEqual稍微好点的函数 – 依然不推荐（外国人写文章真拐弯，都不推荐你写这么多干什么）
bool AlmostEqualRelativeOrAbsolute(float A, float B,
float maxRelativeError, float maxAbsoluteError)
{
if (fabs(A - B) < maxAbsoluteError)
return true;
float relativeError;
if (fabs(B) > fabs(A))
relativeError = fabs((A - B) / B);
else
relativeError = fabs((A - B) / A);
if (relativeError <= maxRelativeError)
return true;
return false;
}
使用整型数字比较
下面才是这篇文章的NB之处，就是使用整型数来对两个比较相近的浮点数进行比较。回想一下前面那个处理绝对误差的问题，那个函数并没有考虑到检查范围内的所有的值。对于一个合法的绝对误差为0.00001并且预期结果为10000的问题，应该能够接受9,999.99999 到10,000.00001之间的所有数据，使用4字节的浮点数，在这个范围内的数字只有一个合法的，就是10000。有没有一种方法可以让我们方便的照出来，在规定范围内，有多少个合法的数字呢？或者说，是否有一种方便的方法能满足：“我觉得结果是10000，但是鉴于浮点数天生有缺陷，我可以接受10000前面5个或者后面5个数。”
好吧，确实有个好方法。
IEEE的浮点数和双精度浮点数在设计时，是字典顺序向量排列的，IEEE设计师William Kahan说“如果两个具有相同格式的浮点数是排列好的话，那么他们将它们的二进制码表示为符号-幅值整型时仍然满足这种关系”（这里有个Sign-Magnitude integer，是指的二进制源码，要区别于二进制补码。计算机中使用的是二进制补码，不满足这个关系）。
这就意味着内存中的两个浮点数，如果把它们的二进制码映射为整型并比较，就可以知道哪个大，完全可以绕过浮点比较。使用C/C++实现就是：
if (*(int*)&f1 < *(int*)&f2)
这个美妙的算法将f1的地址看成一个整型指针，然后废除它。这些指针操作看起来很费工夫，不过它抵消了这些并且意味着“f1是一个整型数”。f2也做了同样的操作，整行语句的意思就是“将它们在内存中的二进制码表示为整型，并且比较大小”。
Kahan曰：把它们表示成符号-幅值整型，就可以比较它们的大小了。不过杯具的是，当代大部分处理器使用二进制补码表示整型数。所以这种方法只在一个或者多个浮点数为正的条件下才能作比较。如果两个数都是负数的话，结果就错了-结果会正好相反。下面我们来搞定这个问题。
由于浮点数是字典排序的，那么如果我们把一个浮点数表示为一个整型数。换句话说，执行下面的代码：
(*(int*)&f1) += 1;
会在这个浮点数上面增加一个值，从而得到下一个浮点数。对于正数，得到下一个更大的数，对于负数，得到下一个更小的数。总之，得到了下一个更远离0的数。
从相反的角度讲，如果将两个浮点数的整型表示相减，就可以知道它们有多接近。如果相减的结果为0，它们就是完全相同。如果差距为1，它们就是相邻的浮点数。总之，如果相差为n，则它们之间相邻n-1个数。
下面的表是一些浮点数和整型数在内存中的表示。这里是离2.0最近的5个浮点数以及它们的整型表示。这个表表示了减法整型表示（subtracting integer representations是什么意思？谁能告诉我？），并且告诉我们在1.99999988和 2.0之间没有别的浮点数了。
|
|
Representation |
|
|
Float value |
Hexadecimal |
Decimal |
|
+1.99999976 |
0x3FFFFFFE |
1073741822 |
|
+1.99999988 |
0x3FFFFFFF |
1073741823 |
|
+2.00000000 |
0x40000000 |
1073741824 |
|
+2.00000024 |
0x40000001 |
1073741825 |
|
+2.00000048 |
0x40000002 |
1073741826 |
考虑到上面的知识，我们继续改写AlmostEqual函数
// 最初的AlmostEqualULPs 版 – 快速简单，但是有些限制
bool AlmostEqualUlps(float A, float B, int maxUlps)
{
assert(sizeof(float) == sizeof(int));
if (A == B)
return true;
int intDiff = abs(*(int*)&A - *(int*)&B);
if (intDiff <= maxUlps)
return true;
return false;
}
这样当然简单，尤其是你发现这里没有浮点除法和浮点数绝对值计算fabs()!（作者在这里用了感叹号！看样子这是个革命性的改进。估计浮点数除法和浮点数绝对值计算是个费工夫的活儿）
上面这个方程的最后一个参数跟前面的AlmostEqual不大一样。我们在最后一个参数使用了一个最大的误差来代替前面的maxRelativeError。它表示了在一个浮点数的末尾，我们能接受的最大的误差是多少。maxUlps也可以表示为在A、B之间可以接受多少个浮点数。这个方程可以接受maxUlps-1个浮点数。
如果两个数几乎一样，只有最后的尾数差一位，那么这个函数就会计算出intDiff=1
如果一个数字是某个指数所能表示的最大数，比如1.99999988，并且另一个数是下一个指数所能表示的最小的数，那么根据这个方程，计算出来的intDiff是1
在这两个例子下两数是最接近的浮点数。
给一个maxRelativeError 和 maxUlps不完全解释。对于一个正常的浮点数maxUlps等于1意味着maxRelativeError处于1/8,000,000 和1/16,000,000之间。这个差距是由于一个浮点数的精准度取决于它是接近它的指数所能表达的浮点数范围的上限还是下限。看看2.0附近的几个数，大于2.0的数的两数差距是小于2.0的两数差距的两倍。
AlmostEqualUlps函数一开头先检查A、B是否相等- 就像AlmostEqualRelative一样，但是有个不一样的原因，欲知后事，请继续看。
有关编译器
在AlmostEqualUlps的最新版本里，我们用了指针，并且告诉编译器将内存中的浮点数看作是整型数。这么干会带来很多隐患，一个问题就是int和float可能不一样长。Float也许是32bits，但是int可以是几乎任何长度。这个问题必须考虑清楚，不过现在的编译器都有32-bits的int，所以这个问题就好解决了，如果你的编译器有其他位数的int，找个32-bits的int代替就行了。
另一个问题来自于混淆现象优化（？aliasing optimizations，谁来翻译一下）。严格来说C/C++标准编译器可以假设不同的类型在内存中并不重叠（有些许异常情况比如char指针）。比如，可以允许一个指向int的指针和一个指向float的指针不指向重叠的内存。这个可以打开很多有价值的优化，但是对于代码，这违反了规定（虽然很平常），它产生了不确定的结果。尤其是，一些版本的g++默认了非常严格的混淆规定，不想AlmostEqualUlps中使用的技术那样。
幸运的是g++知道这么干会出问题，所以给了个如下的warning：
warning: dereferencing type-punned pointer will break strict-aliasing rules
对于这个问题，有两个解决方案。使用-fno-strict-aliasing选项关闭aliasing设置，或者使用一个int和float的结合来将float重新解释为int。欲知详情，请看-fstrict-aliasing的文档。这里给一个参考链接http://blog.csdn.net/quickbasic411/archive/2010/10/04/5921417.aspx
一些更恶心的注意事项
浮点数一点都不简单。AlmostEqualUlps没有很好的解决浮点数的全部稀奇古怪的问题。事发后将这些问题都解决取决于你准备用它干什么，当然，我们需要一个改进版。
IEEE的浮点数分为以下几类：（自己查去吧）
- Zeroes
- Subnormals
- Normal numbers
- Infinities
- NANs
零的处理
AlmostEqual是为了处理正常数字的（废话），并且它将这个功能实现的最好。第一个问题是如何处理0。IEEE浮点数可以有正0和负0。如果你按照浮点数比较，它们是相等的，但是如果你按照整型来比较，它们差远了-正的0.0是整型的0，而负0是0x80000000!（十进制是-2147483648）。下面的表是正0和负0附近的浮点数，以及它们的整型表示
|
|
Representation |
|
|
Float value |
Hexadecimal |
Decimal |
|
+4.2038954e-045 |
0x00000003 |
3 |
|
+2.8025969e-045 |
0x00000002 |
2 |
|
+1.4012985e-045 |
0x00000001 |
1 |
|
+0.00000000 |
0x00000000 |
0 |
|
-0.00000000 |
0x80000000 |
-2147483648 |
|
-1.4012985e-045 |
0x80000001 |
-2147483647 |
|
-2.8025969e-045 |
0x80000002 |
-2147483646 |
|
-4.2038954e-045 |
0x80000003 |
-2147483645 |
在AlmostEqualUlps我解决了这个问题，方法是检查A、B是否完全一致，当输入是一个正0和一个负0的情况，这个问题算是解决了。当然，这还不够，在这种情况下正0和最小的正浮点数会被认为是一个ulp数，并且会被认为是相等。但是负0和最小的负浮点数却在那遥远的地方，这就违反了正0和负0相等的基本原则。
解决负数的一个更好的方法就是把他们换成二进制补码，具体做法就是当检测到负数时，用0x80000000减去这个数。这种方法把负0转化为整型的0（这样就跟正0相同了），并且这么干使得最小的负数对应了整型的-1，并且向下也是一样。这样一转换，处理0附近的数就方便多了。
|
Remapping for twos complement |
||
|
|
Representation |
|
|
Float value |
Hexadecimal |
Decimal |
|
+4.2038954e-045 |
0x00000003 |
3 |
|
+2.8025969e-045 |
0x00000002 |
2 |
|
+1.4012985e-045 |
0x00000001 |
1 |
|
+0.00000000 |
0x00000000 |
0 |
|
-0.00000000 |
0x00000000 |
0 |
|
-1.4012985e-045 |
0xFFFFFFFF |
-1 |
|
-2.8025969e-045 |
0xFFFFFFFE |
-2 |
|
-4.2038954e-045 |
0xFFFFFFFD |
-3 |
但是，这么转换了之后，就不能用IEEE浮点数的标准去处理他们了。（负数的值都变了），但是我们可以将它们看做是整型来比较大小，看下面。
这个方法有了一个新的优点，就是现在两个数字之间的距离可以通过与0之间的距离来测量了。换句话说，最小的负数和最小的正数现在非常近，只有几个ulps的距离，看起来很酷，这等同于在相对误差检查中又引入了一个绝对误差检查，一举两得一箭双雕一石双鸟，代码如下：
// 可以使用的AlmostEqual函数（最终版）
bool AlmostEqual2sComplement(float A, float B, int maxUlps)
{
// Make sure maxUlps is non-negative and small enough that the
// default NAN won't compare as equal to anything.
assert(maxUlps > 0 && maxUlps < 4 * 1024 * 1024);
int aInt = *(int*)&A;
// Make aInt lexicographically ordered as a twos-complement int
if (aInt < 0)
aInt = 0x80000000 - aInt;
// Make bInt lexicographically ordered as a twos-complement int
int bInt = *(int*)&B;
if (bInt < 0)
bInt = 0x80000000 - bInt;
int intDiff = abs(aInt - bInt);
if (intDiff <= maxUlps)
return true;
return false;
}
下一个问题是极小数问题（subnormal，我翻译为极小数，不知是否合适），极小数就是指那些非常接近0的数，它们跟正常数不大一样。缺点就是它们的精度很小（越接近0，精度越小）。当比较两个极小数的时候，一个ulp的误差可能会产生很大的相对误差，甚至可能100%。但是ulps作为测量两个数字的误差的工具仍然有意义。因此，在相对误差描述中的这个变化说不定是件好事-比较浮点数的比较的另一个优点。
无穷大
IEEE浮点数对无穷大有一个特殊的描述。无穷大被用来表示溢出或者除零。无穷大是与能表示的最大的浮点数相连的。所以AlmostEqualUlps函数中FLT_MAX与无穷大是一样的。从某种意义上讲，这是对的，因为它们之间没有可表示的浮点数了，当然从另外的角度来说也不对，毕竟，任何有限的数都不能接近于无穷大
如果你认为有限数不能接近FLT_MAX的话你就加个额外的检查吧。
非法
IEEE浮点数有对于NAN的描述。NAN的指数部分与一些合法数是共享的，但是在尾数上与合法数有差别，所以，当作为int进行比较时，它们就与合法数相邻（也许会带来不必要的麻烦）。所以，有可能一个有限数或者FLT_MAX会与NAN相邻。所以，一旦你的代码输出一个NAN的结果那就杯具了。有两个东西来保护这种情况，一是大部分浮点代码不会轻易产生NAN，事实上大多数浮点数代码会将NAN看做是除零错误并报出异常。另一个东西是通常只有一个NAN值，在x87协处理器中这个值是0xFFC00000，跟最近的有效数有四百万的距离。这个值之所以设成这样是有特殊目的的，因为在NAN的比较中它们会卷曲（？卷曲？缠绕？什么东东？可能是指浮点数是转圈轮回的意思吧）。一个值为0xFFFFFFFF的NAN会被认为是接近于0。在AlmostEqual2sComplement中使用的转换负数的方法避免了这个问题，是因为各种NAN之间只能互相转换，是个封闭集合。但是值为0xFFC00000的NAN也可以避免这个问题，因为它将NAN放在离转换的距离四百万个ulp之外。（我送你离开千里之外你是否还在）
还有一个问题是从常理上来讲出现了NAN之后，程序应该返回假，但是AlmostEqual2sComplement在出现两个整型表示的NAN相等时，会返回相等。如果你想解决这个问题，自己动手吧。
极限
maxUlps也不能说随便多大都行，如果maxUlps大于等于四百万的话就会有跑入NAN的风险。如果maxUlps大于等于一千六百万，就会有最大的负数等于最大的正数的风险。（看来这个函数还是暗藏杀机啊）
一般的用于实际使用的情况下这么大的maxUlps是没用的。一个一千六百万的maxUlps意味着+100%和-50%的差距都是可以认为是相等的。一个四百万的maxUlps意味着+25%和-12.5%的差距可以认为是相等的。如果真要用这么大的maxUlps，就要加上检查是否存在大于NAN的数，还有相反符号的数。为了防止大maxUlps带来的各种问题，最好把maxUlps限制在一个安全的范围。
AlmostEqual2sComplement对IEEE浮点数格式的依赖性很高，并且有个前提就是二进制补码的整型是同样大小的。在大多数机器上，尤其是家庭和商业用机，这些限制都是正常的，但是还是有些个特殊情况的。由于这个原因，再加上这个技术很乱并且不那么显而易见，所以有必要写一个有文档、有判断、有条件检查的函数。
总结
AlmostEqual2sComplement是一个浮点数比较的有效方法，它的表现没有AlmostEqualRelative这么完美，但是在许多方面，它比AlmostEqualRelative更好。总结一下，AlmostEqual2sComplement有如下特性：
- 衡量两个浮点数是否足够接近，接近的程度由ulps定义，并且算出两个数之间有多少个浮点数。
- 将无穷大看做是接近FLT_MAX。
- 将NAN扔到了千里之外，与正常数之间有四百万的距离（NAN是x87默认值）。
- 接受更大的相对误差逐渐地下溢到无穷小（这句话有语法错误，我反正看不懂）。
- 将极小的负数看做与极小的正数相邻。
有一些AlmostEqual2sComplement的特殊问题在需要的情况下可以另行考虑。如果你愿意，你可以#ifdefed一些简单的控制，来生成一个带有必要检查的新版本。
AlmostEqual2sComplement在能将浮点数快速转化为整型的机子上能够很好的工作。因为它需要过内存（指针操作）所以会比较慢。有些机器的向量处理可以在同一个寄存器里做浮点和整型的操作，这个函数就能发挥作用了。因为这么做不需要走内存。（不需要走内存也得走缓存吧？）
这个技术也能用于双精度浮点数，对应到64-bit整型。双精度浮点有53位尾数，所以每个ulp对应的相对误差在1/4,000,000,000,000,000 和 1/8,000,000,000,000,000之间。
参考文献：
IEEE Standard 754 Floating Point Numbers by Steve Hollasch
Lecture Notes on the Status of IEEE Standard 754 for Binary Floating-Point Arithmetic by William Kahan
Source code for compare functions and tests
Riding the STP on a unicycle, again
Ride the Lobster -- riding Nova Scotia on a unicycle
Comparing Floats - as Integers
Game Scripting in Python from GDC 2002
What Happened to my Colours!?! Displaying Console Computer Graphics on a TV - available here.
Game Programming Gems 2 - Micro-Threads for Game Object AI
Game Programming Gems 2 - Game Input Recording and Playback
Debugging release executables with VC++
Structured Exception Handling, Jan 1999 Game Developer Magazine - code available
Playing Perfect Pikmin - a guide to the nine day, zero deaths Pikmin challenge
Printing Digital Pictures - a guide to photo printing sites
Calculating PI - fourteen decimal places with pen and paper
