浮点数据的精度分析

浮点数据的精度分析

    如果浮点算术计算的结果不如预期的精确，可能是由于计算机硬件的限制所致。由于硬件缺乏足够的位而无法呈现具有完美精度的结果，计算机可能会将结果值截断，使得结果值不够准确（本文分析基于MATLAB进行）。

1.双精度数的精度

    由于双精度数的数量有限，因此您无法在双精度存储中表示所有数值。在任何计算机上，每个双精度数和下一个最小的更大双精度数之间都存在一个较小的间隔。您可以使用 eps 函数确定此间隔的大小，该大小限制了您的结果的精度。例如，要计算 5 和下一个更大的双精度数之间的间距，请输入

format long

eps(5)

ans =
     8.881784197001252e-16

    这表明，5 和 5 + eps(5) 之间不存在任何双精度数。如果某个双精度计算返回答案 5，则结果仅精确到 eps(5) 之内。

eps(x) 的值取决于 x。以下示例显示，当 x 变得更大时，eps(x) 也会变得更大：

eps(50)
ans =
     7.105427357601002e-15

    如果您输入不带输入参数的 eps，MATLAB 将返回 eps(1) 的值（从 1 到 下一个最小的更大双精度数之间的间距）。

eps()

    eps(5)和eps(50)的精度间距为什么不一致呢？我们只要分析eps(5)和eps(50)答案是怎么得到即可？下面按照IEEE 754 标准来分析(IEEE 754 标准介绍见整型与浮点型在计算机中的表示与储存)：
浮点型在计算机中的表示：

V = (−1)^S * M ∗ 2^E

双精度5.0的表示：
写成二进制是101.0，相当于1.01×2^2。那么，按照上面V的格式，可以得出S=0，M=1.01，E=2。

5.0 = (−1)⁰ * 1.0100…000 ∗ 2² (1.01后边有50个0)
则，基于此表示规则，要表示比5.0仅大一丢丢（这个“一丢丢”尽可能的小）的数，改变的部分只能是 1.0100…001部分，因为如果改变此部分 2² ，数值至少发生2倍的改变，不满足“一丢丢”。5.0+“一丢丢”表示为： (−1)⁰ * 1.0100…001 ∗ 2² (1.01后边有49个“0” + 1个“1”)
由此可得： “一丢丢” = 0.00…001 ∗ 2² = 2^-52+2 = 2^-50 = 8.8817841970012523233890533447266e^-16
通过windows自带计算器计算得到，与eps(5)一致。

同理，

50.0 = (−1)⁰ * 1.100100…000 ∗ 2⁵ (1.1001后边有49个0)
eps(50)的“一丢丢” = 0.00…001 ∗ 2⁵ = 2^-52+5 = 2^-47 = 7.1054273576010018587112426757813e^-15
通过windows自带计算器计算得到，与eps(50)一致。

2.单精度数的精度

    同样，两个单精度数之间也存在间隔。如果 x 的类型为 single，则 eps(x) 返回 x 和下一个更大的单精度数之间的间距。例如，

x = single(5);
eps(x)

返回

ans =

  single

  4.7684e-07

单精度5.0的表示：写成二进制是101.0，相当于1.01×2^2。那么，按照上面V的格式，可以得出S=0，M=1.01，E=2。

5.0 = (−1)⁰ * 1.0100…000 ∗ 2² (1.01后边有21个0)

同理：

“一丢丢” = 0.00…001 ∗ 2² = 2^-23+2 = 2^-21 = 476837158203125e^-7
通过windows自带计算器计算得到，与eps(5)一致。

    请注意，此结果大于 eps(5)。由于单精度数的数量少于双精度数的数量，因此单精度数之间的间隔也大于双精度数之间的间隔。这意味着，单精度算术运算的结果精度要低于双精度算术运算的结果精度。

    对于类型为 double 的双精度数 x，eps(single(x)) 的计算值即为将 x 从 double 转换为 single 时的舍入上限。
例如，当您将双精度数 3.14 转换为 single 时，它会通过以下方式进行舍入

double(single(3.14) - 3.14)

ans =
   1.0490e-07

3.14 舍入的数量小于

eps(single(3.14))
ans =

  single

  2.3842e-07

3.避免浮点算术运算出现常见问题

    MATLAB 中几乎所有运算都是通过符合 IEEE 754 标准的双精度算术运算执行的。由于计算机仅将数值表示为有限精度（双精度需要 52 个尾数位），因此计算有时会生成数学上的非预期结果。务必注意，这些结果并非 MATLAB 中的错误。

使用以下示例来帮助您确定这些情况：

示例3.1 - 舍入或您所获得的不是所期望的

    十进制数 4/3 不能精确表示为二进制分数。为此，以下计算的结果不是零，而是显示数量 eps。

e = 1 - 3*(4/3 - 1)

e =
   2.2204e-16

    同样，0.1 也不能精确表示为二进制数（整型与浮点型在计算机中的表示与储存：小数的二进制表示）。因此，您会发现以下非预期行为：

a = 0.0;
for i = 1:10
  a = a + 0.1;
end
a == 1
ans =

  logical

   0

    请注意，计算中运算的顺序会很重要：

b = 1e-16 + 1 - 1e-16;
c = 1e-16 - 1e-16 + 1;
b == c
ans =

  logical

   0

    浮点数之间存在间隔。当数值变得越大时，间隔也会变得越大，如以下所示：

(2^53 + 1) - 2^53

ans =
     0

    由于 pi 实际上不是 π，因此 sin(pi) 不精确为零并不足为奇：

sin(pi)

ans =
     1.224646799147353e-16

示例 3.2 - 具有灾难性后果的取消操作

    对几乎相等的操作数执行减法时，有时可能会发生意外取消。以下是因淹没（使加法没有意义的精度损失）导致的取消的示例。

sqrt(1e-16 + 1) - 1

ans =
     0

    MATLAB 中的某些函数（例如 expm1 和 log1p）可用于弥补这种灾难性取消所造成的影响。

示例 3.3 - 浮点运算和线性代数

    解决线性代数的问题时，舍入、取消和浮点算术运算的其他一些特性结合起来时，可能产生令人意想不到的运算。MATLAB 会警告下面的矩阵 A 是病态的，因此即便是细微的扰动，都可能对方程组 Ax = b 产生很大的影响：

A = diag([2 eps]); 
b = [2; eps]; 
y = A\b; 
Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 1.110223e-16.

    这些只是几个例子，用于说明 IEEE 浮点算术运算如何影响 MATLAB 中的计算。请注意，IEEE 754 算术运算中执行的所有计算都会受到影响，这包括用 C 或 FORTRAN 编写的应用程序以及 MATLAB。
参考资料：mathwork浮点数据的精度