浮点计算中的精度和准确性

原始 KB 数: 125056

概要

在很多情况下,浮点计算的精度、舍入和准确性都可用于生成对程序员来说令人惊讶的结果。 它们应遵循四个常规规则:

  1. 在涉及单精度和双精度的计算中,结果通常不会比单精度更准确。 如果需要双精度,请确保计算中的所有项(包括常量)都指定为双精度。

  2. 从不假定计算机中准确表示简单的数值。 大多数浮点值不能精确表示为有限的二进制值。 例如, .1.0001100110011... 进制中(无限重复),因此不能使用二进制算术在计算机上完全准确地表示(包括所有电脑)。

  3. 从不假定结果准确到最后一个小数位。 所谓“真实”答案与任何浮点运算单元有限精度所能计算的结果之间,总是存在很小的差异。

  4. 永远不要比较两个浮点值来判断它们是否相等或不相等。 这是规则 3 的必然结果。 “应该”等于“的数字之间几乎总是会有一些小的区别。 相反,请始终检查数字是否几乎相等。 换句话说,检查它们之间的差异是否很小或微不足道。

详细信息

通常,上述规则适用于所有语言,包括 C、C++ 和汇编程序。 下面的示例演示了使用 FORTRAN PowerStation 的一些规则。 所有示例都是使用 FORTRAN PowerStation 32 编译的,没有任何选项,但最后一个示例是用 C 编写的。

示例 1

第一个示例演示了两个事项:

  • 默认情况下,FORTRAN 常量是单精度(默认情况下 C 常量为双精度)。
  • 包含任何单精度术语的计算并不比所有字词都是单精度的计算更准确。

使用 1.1(单精度常量)初始化后,y 与单个精度变量一样不准确。

x = 1.100000000000000  y = 1.100000023841858

将单个精度值乘以准确的双精度值的结果几乎与乘以两个单精度值一样差。 这两个计算的误差是两个双精度值相乘的数千倍。

true = 1.320000000000000 (multiplying 2 double precision values)
y    = 1.320000052452087 (multiplying a double and a single)
z    = 1.320000081062318 (multiplying 2 single precision values)

示例代码

C Compile options: none

       real*8 x,y,z
       x = 1.1D0
       y = 1.1
       print *, 'x =',x, 'y =', y
       y = 1.2 * x
       z = 1.2 * 1.1
       print *, x, y, z
       end

示例 2

示例 2 使用二次公式。 它表明,即使是双精度计算也不完美,当小错误可能导致重大结果时,应该在依赖计算结果之前对其进行测试。 示例 2 中平方根函数的输入仅略为负,但它仍然无效。 如果双精度计算没有轻微的错误,则结果为:

Root =   -1.1500000000

而是生成以下错误:

运行时错误 M6201:MATH

  • sqrt:域错误

示例代码

C Compile options: none

       real*8 a,b,c,x,y
       a=1.0D0
       b=2.3D0
       c=1.322D0
       x = b**2
       y = 4*a*c
       print *,x,y,x-y
       print "(' Root =',F16.10)",(-b+dsqrt(x-y))/(2*a)
       end

示例 3

示例 3 演示即使未启动优化,仍由于某些优化过程,值可能会暂时保持比预期更高的精度,因此测试两个浮点数是否相等是不明智的。

在此示例中,两个值在某些情况下既相等又不相等。 在第一个 IF 中,Z 的值仍位于协处理器的堆栈上,其精度与 Y 相同。因此 X 不等于 Y,第一条消息打印出来。在第二个 IF 时,必须从内存中加载 Z,因此精度和值与 X 相同,并且还会打印第二条消息。

示例代码

C Compile options: none

       real*8 y
       y=27.1024D0
       x=27.1024
       z=y
       if (x.ne.z) then
         print *,'X does not equal Z'
       end if
       if (x.eq.z) then
         print *,'X equals Z'
       end if
       end

示例 4

示例代码 4 的第一部分计算两个数字接近 1.0 之间的最小可能差值。 它通过在 1.0 的二进制表示中添加一个比特来实现这一点。

x   = 1.00000000000000000  (one bit more than 1.0)
y   = 1.00000000000000000  (exactly 1.0)
x-y =  .00000000000000022  (smallest possible difference)

某些版本的 FORTRAN 在显示数字时将其四舍五入,以使得固有的数字不精确性不那么明显。 这就是为什么 x 和 y 在显示时看起来相同的原因。

示例代码 4 的第二部分计算两个接近 10.0 的数字之间的最小可能差值。 同样,它通过在十进制数 10.0 的二进制表示中添加一个位来实现。 请注意,接近 10 的数字之间的差异大于接近 1 的差值。 这演示了一般原则,即一个数字的绝对值越大,它就越不精确地存储在给定的位数中。

x   = 10.00000000000000000  (one bit more than 10.0)
y   = 10.00000000000000000  (exactly 10.0)
x-y =   .00000000000000178

还会显示这些数字的二进制表示形式,以显示它们只存在 1 位差异。

x = 4024000000000001 Hex
y = 4024000000000000 Hex

示例代码 4 的最后一部分显示,简单的非重复十进制值通常只能通过重复分数以二进制形式表示。 在本例中,x=1.05 需要在尾数中使用重复因子 CCCCCCCC....(十六进制)。 在 FORTRAN 中,最后一个数字“C”向上舍入为“D”,以保持最高的准确性:

x = 3FF0CCCCCCCCCCCD (Hex representation of 1.05D0)

即使在舍入后,结果也不完全准确。 在最不重要数字之后会出现一些错误,通过删除第一个数字可以看到。

x-1 = .05000000000000004

示例代码

C Compile options: none

       IMPLICIT real*8 (A-Z)
       integer*4 i(2)
       real*8 x,y
       equivalence (i(1),x)

       x=1.
       y=x
       i(1)=i(1)+1
       print "(1x,'x  =',F20.17,'  y=',f20.17)", x,y
       print "(1x,'x-y=',F20.17)", x-y
       print *

       x=10.
       y=x
       i(1)=i(1)+1
       print "(1x,'x  =',F20.17,'  y=',f20.17)", x,y
       print "(1x,'x-y=',F20.17)", x-y
       print *
       print "(1x,'x  =',Z16,' Hex  y=',Z16,' Hex')", x,y
       print *

       x=1.05D0
       print "(1x,'x  =',F20.17)", x
       print "(1x,'x  =',Z16,' Hex')", x
       x=x-1
       print "(1x,'x-1=',F20.17)", x
       print *
       end

示例 5

在 C 中,浮点常量默认为双精度。 使用“f”指示浮点值,如“89.95f”。

/* Compile options needed: none
*/ 

#include <stdio.h>

void main()
   {
      float floatvar;
      double doublevar;

/* Print double constant. */ 
      printf("89.95 = %f\n", 89.95);      // 89.95 = 89.950000

/* Printf float constant */ 
      printf("89.95 = %f\n", 89.95F);     // 89.95 = 89.949997

/*** Use double constant. ***/ 
      floatvar = 89.95;
      doublevar = 89.95;

printf("89.95 = %f\n", floatvar);   // 89.95 = 89.949997
      printf("89.95 = %lf\n", doublevar); // 89.95 = 89.950000

/*** Use float constant. ***/ 
      floatvar = 89.95f;
      doublevar = 89.95f;

printf("89.95 = %f\n", floatvar);   // 89.95 = 89.949997
      printf("89.95 = %lf\n", doublevar); // 89.95 = 89.949997
   }