浮点编程的五个技巧
4.91/5 (108投票s)
本文介绍了在处理浮点数时，了解的五件最重要的事情。
引言
即使是非常有经验的程序员，在编写依赖浮点数运算的代码时，也可能陷入几个陷阱。本文解释了在处理浮点数（即float和double数据类型）时需要牢记的五件事。
不要进行相等性测试
你几乎从不想写出如下代码：
double x;
double y;
...
if (x == y) {...}
大多数浮点数运算都至少会涉及微小的精度损失，因此即使两个数字在实践中相等，它们在最后一位二进制位上可能并不完全相等，所以相等性测试很可能会失败。例如，以下代码片段打印出-1.778636e-015。虽然理论上，平方可以抵消平方根，但往返操作的精度略有不准确。
double x = 10;
double y = sqrt(x);
y *= y;
if (x == y)
cout << "Square root is exact\n";
else
cout << x-y << "\n";
在大多数情况下，上面的相等性测试应该写成如下形式：
double tolerance = ...
if (fabs(x - y) < tolerance) {...}
这里的tolerance是某个阈值，它定义了什么叫做“足够接近”的相等。这引出了一个问题：什么样的距离才算足够近？这个问题不能笼统地回答；你必须了解你的具体问题才能知道在你的上下文中，“足够近”是多少。
比乘除法更要担心加减法
乘法和除法的相对误差总是很小的。另一方面，加法和减法可能导致精度完全丢失。实际上问题在于减法；加法只有在相加的两个数字符号相反时才会成为问题，所以你可以将其视为减法。尽管如此，代码可能用一个“+”号写成，但实际上是减法。
当相减的两个数字非常接近时，减法就会成为一个问题。数字越接近，精度损失的可能性越大。具体来说，如果两个数字在 n 位上相同，那么在减法中可能会损失 n 位的精度。这在极端情况下最容易看清楚：如果两个数字理论上不相等，但在其机器表示中相等，那么它们的差将被计算为零，精度损失 100%。
这里有一个例子，说明了这种精度损失经常出现。函数f在点x处的导数定义为当h趋近于零时，(f(x+h) - f(x))/h的极限。因此，计算函数导数的一种自然方法是为某个小的h值计算(f(x+h) - f(x))/h。理论上，h越小，这个分数越接近导数。实际上，精度会先提高一段时间，但超过某个点后，更小的h值会导致对导数的近似变差。随着h变小，近似误差变小，但数值误差增加。这是因为f(x+h) - f(x)的减法变得有问题。如果你取的h足够小（毕竟，理论上，越小越好），那么f(x+h)将等于f(x)到机器精度。这意味着所有导数都将被计算为零，不管函数是什么，只要你取的h足够小。这里有一个计算sin(x)在x = 1处的导数的例子。
cout << std::setprecision(15);
for (int i = 1; i < 20; ++i)
{
double h = pow(10.0, -i);
cout << (sin(1.0+h) - sin(1.0))/h << "\n";
}
cout << "True result: " << cos(1.0) << "\n";
这是上面代码的输出。为了使输出更容易理解，第一个错误数字之后的数字已替换为句点。
0.4...........
0.53..........
0.53..........
0.5402........
0.5402........
0.540301......
0.5403022.....
0.540302302...
0.54030235....
0.5403022.....
0.540301......
0.54034.......
0.53..........
0.544.........
0.55..........
0
0
0
0
True result: 0.54030230586814
当h = 10-8时，精度有所提高。之后，由于减法中的精度损失，精度开始下降。当h = 10-16或更小时，输出恰好为零，因为sin(1.0+h)到机器精度等于sin(1.0)。（事实上，到机器精度，1+h等于1。稍后会详细介绍。）
（以上结果是用 Visual C++ 2008 计算的。当在 Linux 上用 gcc 4.2.3 编译时，结果相同，除了最后四个数字。VC++ 产生零，而 gcc 产生负数：-0.017...、-0.17...、-1.7... 和 17....）
当你的问题需要减法并且会导致精度损失时，你会怎么做？有时精度损失不是问题；double类型一开始就有大量的备用精度。当精度很重要时，通常可以使用一些技巧来改变问题，使其不需要减法，或者不需要你最初需要的减法。
请参阅 CodeProject 文章 避免溢出、下溢和精度损失，其中有一个关于使用代数技巧将二次公式转换为更适合保留精度的形式的例子。另请参阅 比较计算标准差的三种方法，其中有一个关于代数等价方法如何表现出巨大差异的例子。
浮点数有有限的范围
每个人都知道浮点数有有限的范围，但这个限制可能会以意想不到的方式出现。例如，你可能会对以下代码行的输出感到惊讶。
float f = 16777216;
cout << f << " " << f+1 << "\n";
这段代码打印两次16777216。发生了什么？根据浮点数算术的 IEEE 规范，float类型宽度为 32 位。其中 24 位用于尾数（以前称为**基数**），其余用于指数。数字16777216是224，因此float变量f没有剩余的精度来表示f+1。如果f是double类型，则会发生类似的现象，因为 64 位double有 53 位用于尾数。以下代码打印0而不是1。
x = 9007199254740992; // 2^53
cout << ((x+1) - x) << "\n";
在将小数字加到中等大小的数字时，我们也可能遇到精度问题。例如，以下代码打印“Sorry!”，因为DBL_EPSILON（在*float.h*中定义）是最小的正数 e，使得在使用double类型时 1 + e != 1。
x = 1.0;
y = x + 0.5*DBL_EPSILON;
if (x == y)
cout << "Sorry!\n";
类似地，常量FLT_EPSILON是最小的正数 e，使得在使用float类型时 1 + e 不等于 1。
使用对数避免溢出和下溢
上一节描述的浮点数限制源于尾数中的位数有限。溢出和下溢是指数位数也有限的结果。有些数字太大或太小而无法存储在浮点数中。
许多问题似乎需要计算一个中等大小的数字，作为两个巨大数字的比值。最终结果可以表示为浮点数，即使中间结果不能。在这种情况下，对数提供了一种解决方案。如果你想计算大数M和N的M/N，计算log(M) - log(N)并将结果应用exp()。例如，概率通常涉及阶乘的比值，而阶乘会迅速变得天文数字。对于N > 170，N!大于DBL_MAX，这是double（不带扩展精度）可以表示的最大数字。但可以通过以下方式评估200!/(190! 10!)这样的表达式，而不会发生溢出：
x = exp( logFactorial(200)
- logFactorial(190)
- logFactorial(10) );
一个简单但效率低下的logFactorial函数可以这样写：
double logFactorial(int n)
{
double sum = 0.0;
for (int i = 2; i <= n; ++i)
sum += log((double)i);
return sum;
}
更好的方法是使用对数伽马函数（log gamma function），如果有一个可用的话。有关更多信息，请参阅 如何计算二项式概率。
数值运算并不总是返回数字
由于浮点数存在局限性，有时浮点运算会返回“无穷大”，以表示“结果超出我的处理能力”。例如，以下代码在 Windows 上打印1.#INF，在 Linux 上打印inf。
x = DBL_MAX;
cout << 2*x << "\n";
有时，返回有意义结果的障碍与逻辑有关，而不是有限精度。浮点数据类型表示实数（而不是复数），并且没有一个实数的平方是-1。这意味着当代码请求sqrt(-2)时，即使在无限精度下，也没有有意义的数字可以返回。在这种情况下，浮点运算返回NaN。这些是浮点值，它们代表错误代码而不是数字。NaN值在 Windows 上显示为1.#IND，在 Linux 上显示为nan。
一旦一系列运算遇到NaN，之后的所有结果都是NaN。例如，假设你有一些代码，其作用类似于以下内容：
if (x - x == 0)
// do something
什么可能阻止if语句之后的代码执行？如果 x 是NaN，那么x - x也是NaN，而NaN不等于任何东西。事实上，NaN甚至不等于它自己。这意味着表达式x == x可以用来测试x是否是（可能是无限的）数字。有关无穷大和NaN的更多信息，请参阅 C++ 中的 IEEE 浮点异常。
更多信息
文章 每个计算机科学家都应该了解的浮点数算术 详细介绍了浮点数算术。它可能是每个计算机科学家理想中应该知道的，但很少有人能完全理解其中的所有内容。
历史
- 2008 年 9 月 24 日：原始帖子
- 2008 年 10 月 29 日：添加了引用，修改了代码使其也能与 gcc 编译，报告了一个示例中的 VC++ 与 gcc 的差异