浏览量:1,501

浮点数误差

转载请注明原文章链接:http://www.twinklingstar.cn/2015/1501/roundoff-error-of-floating-point-2/

在计算机中,小数有两种表示方式:定点数和浮点数。小数点隐含在某一个固定位置上的数据,就称为定点数。为了能正确的表示定点数,我们必需规定数据的位数和小数点的位置,以8位的数据为例,规定低3位是小数位,则一个定点数\({\left( {00010.110} \right)_2} = {2^1} + {2^{ - 1}} + {2^{ - 2}} = 2.75\)。定点表示法较为简单,但它表示的数据非常有限,无法表示数值很大的数据和数值很小的数据。为了表示更大范围的数据,数学上通常采用科学计数法,把数据表示成一个小数乘以一个以10为底的指数,例如,太阳的质量\(1.989 \times {10^{33}}\)克。科学计数法改写为\(0.1989 \times {10^{34}}\),把有效数字\(0.1989\)和以10底的指数\(34\)存储在机器中的固定单元内,就可以表示出数值如此大的数据,如果指数是一个负数,例如\( - 34\),则表示的数据的数值非常小。像这样,把一个数的有效数字和指数部分在计算机中分别予以表示,就是浮点表示法。这种表示法,相当于数的小数点位置随比例因子的不同而在一定范围内自由浮动,改变指数部分的数值相当于改变小数点的位置。

IEEE 754标准[4,5]是由1985年[1]引入的浮点数算术标准,是最广泛使用的浮点数运算标准,在2008年[2]对该标准进行了修改。

2015-1-30 15-18-57

 图1. 单精度和双精度浮点数的存储格式

IEEE 754标准规定的单精度和双精度浮点数的存储格式如图1所示,最高位是符号位\(S\),规定浮点数的正负;偏置指数\(E\),紧跟在符号位之后,占\(w\)位;有效数字在最后,占\(p\)位。IEEE 754规定偏置指数\(E = 0\)\(E = {2^w} - 1\)保留,用于表示\( \pm 0\)\( \pm \infty \)\(NaN\)。通过与偏置量的差来表示实际的指数值,规定偏置量\(bias = {2^{w - 1}} - 1\),有效数字\(F\)的二进制表示形式为\({d_0}{d_1} \cdots {d_{p - 1}},{d_k} \in \left\{ {0,1} \right\},k \in \left\{ {0,1, \cdots ,p - 1} \right\}\)\(emin = 2 - {2^{w - 1}}\),\(emax = {2^{w - 1}} - 1\)。对于单精度和双精度表示中,不同的参数如表1.所示:

1. 单精度和双精度浮点数格式的参数表

2015-1-30 15-19-14

浮点数的数值可以用三元组\(\left( {s,e,f} \right)\)表示,\(\left( {s,e,f} \right) = {\left( { - 1} \right)^s} \times {2^e} \times f\),不同数据的浮点数表示的数值分别如下所示:

1.              若\(1 \le E \le {2^w} - 2\),则\(value = \left( {S,E - bias,1.F} \right) = {\left( { - 1} \right)^S} \times {2^{E - bias}} \times 1.{d_0}{d_1} \cdots {d_{p - 1}}\),称该浮点数为正规数;

2.              若\(E = {2^w} - 1\)\(F \ne 0\),则\(value = NaN\)(无论\(S\)为0或1);

3.              若\(E = {2^w} - 1\)\(F = 0\),则\(value = {\left( { - 1} \right)^S} \times \left( { + \infty } \right)\)

4.              若\(E = 0\)\(F \ne 0\),则\(value = \left( {S,emin,0.F} \right) = {\left( { - 1} \right)^S} \times {2^{emin}} \times 0.{d_0}{d_1} \cdots {d_{p - 1}}\),称该浮点数为非正规数;

5.              若\(E = 0\)\(F = 0\),则\(value = {\left( { - 1} \right)^S} \times \left( { + 0} \right)\)

以单精度浮点数为例,单精度编码和对应的数值如表2.所示:

2. 单精度的位模式和对应的数值

2015-1-30 15-19-36

例如,将数“-7.5”用单精度浮点数的二进制表示,它的转化过程如下所示:

2015-1-30 20-33-20

数“-7.5”的二进制浮点数表示为\(\underline {{\rm{1100 0000}}} {\rm{ }}\underline {{\rm{1111 0000}}} {\rm{ }}\underline {{\rm{0000 0000}}} {\rm{ }}\underline {{\rm{0000 0000}}} {\rm{ = 0x40 f0 00 00}}\)

为了方便理解,设浮点数的底数为10,有3个有效数字,\(emin = - 98\)。存在两个浮点数\(x = 6.87 \times {10^{ - 97}}{\rm{,y}} = 6.81 \times {10^{ - 97}}\),则\(x - y = 0.06 \times {10^{ - 97}} = 6.0 \times {10^{ - 99}}\)太小,无法用一个正规数表示,如果不存在非正规数,那么\(x - y\)将被舍入为0。由于\(x \ne y\),但是\(x - y = 0\),这显然就是错误的,因此IEEE 754标准规定,当\(E = 0,F \ne 0\)时,浮点数的数值为\({\left( { - 1} \right)^S} \times {2^{emin}} \times 0.{d_0}{d_1} \cdots {d_{p - 1}}\),增加了非正规数。如图2.所示,非正规数弥补了数轴上区间\(\left( {0,{2^{emin}}} \right)\)内的数,显然,绝对值越小的数,能精确表示的概率越大。

2015-1-30 15-20-13

图2. 浮点数在数轴上的分布

浮点数表示的数据范围比定点数更广泛,但它并不能精确的表示在最大数和最小数之间的所有实数,例如无理数\(\sqrt 3 \),有理数\(0.1\)。IEEE 754定义了四种舍入模式(Rounding Mode):(1)向最接近的可表示的值舍入,当有两个最接近的可表示的值时首选“偶数”值;2)向负无穷大舍入;3)向正无穷大舍入;4)向0截断。以十进制数为例,不同的数采用四种舍入模式得到的结果如表3.所示。在最接近舍入模式中,当存在两个最接近的表示的值时,Reiser& Knuth[6,7]解释了为什么首选“偶数”而不是“奇数”。

3. 不同的浮点数采用不同的舍入模式

2015-1-30 15-20-23

浮点数不仅在实数的表示上存在误差,在数值运算上也存在误差。为了方便理解,设浮点数的底数为10,有3个有效数字。计算\(2.15 \times {10^{12}} - 1.25 \times {10^4}\),则先进行指数对齐,得到:

\(x = 2.15 \times {10^{12}}\)

\(y = 0.00000000125 \times {10^{12}}\)

\(x - y = 2.15 \times {10^{12}}\)

如果对两个绝对值差别很大的两个数作加减操作,就会出现上述所示的“大数吃小数”的现象。被消除的小数\(1.25 \times {10^4}\),就是本次运算的误差值。用尾数最后位(units in the last, ulp)来作为衡量误差的单位,例如,四位有效数的数值,\(ulp( - 1100) = 1,ulp(1) = 0.001\)。两个数值运算的误差表示为\(err(a\)\(\widetilde * b)\),符号\(\widetilde * \)表示对两个数值进行*(加、减、乘、除)运算得到的估计值,则有\(a\widetilde * b = a * b + err(a\widetilde * b)\),IEEE 754标准可以保证\(err(a\widetilde * b) \le \frac{1}{2}ulp(a\widetilde * b)\) [6,8,9]。为了降低浮点数的运算误差,引入保护位(guard bit)、舍入位(round bit)和粘贴位(sticky bit),额外添加的位能有效的降低运算误差[5,6]

此外,由于浮点数范围的局限,如果两个特别大(小)的数相乘或除以一个特别小(大)的数,就可能造成上(下)溢,这也是误差的一个来源。

在运算过程中,重复利用上一次浮点运算的结果,可能造成误差的累加。例如,计算多项式\({x^4}\),如果先计算\({x^2} = x \cdot x\),再计算\({x^3} = {x^2} \cdot x\),最后计算\({x^4} = {x^3} \cdot x\),共有3次运算的叠加,叠加的次数越多,可能造成的误差也将越大。

附件代码

IEEE_754_unions.h

//IEEE_754_unions.h
//Modified by twinklingstar.

// Copyright (c) 1999  
// Utrecht University (The Netherlands),
// ETH Zurich (Switzerland),
// INRIA Sophia-Antipolis (France),
// Max-Planck-Institute Saarbruecken (Germany),
// and Tel-Aviv University (Israel).  All rights reserved. 
//
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 3 of the License,
// or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/releases/CGAL-4.1-branch/Number_types/include/CGAL/IEEE_754_unions.h $
// $Id: IEEE_754_unions.h 67093 2012-01-13 11:22:39Z lrineau $
//
//
// Author(s)     : Stefan Schirra

#ifndef SR_IEEE_754_UNIONS_H_
#define SR_IEEE_754_UNIONS_H_

#include <iomanip>

#if defined (__GLIBC__)
#  include <endian.h>
#  if (__BYTE_ORDER == __LITTLE_ENDIAN)
#    define SR_LITTLE_ENDIAN
#  elif (__BYTE_ORDER == __BIG_ENDIAN)
#    define SR_BIG_ENDIAN
#  else
#    error Unknown endianness
#  endif
#elif defined(__sparc) || defined(__sparc__) \
	|| defined(_POWER) || defined(__powerpc__) \
	|| defined(__ppc__) || defined(__hppa) \
	|| defined(_MIPSEB) || defined(_POWER) \
	|| defined(__s390__)
#  define SR_BIG_ENDIAN
#elif defined(__i386__) || defined(__alpha__) \
	|| defined(__x86_64) || defined(__x86_64__) \
	|| defined(__ia64) || defined(__ia64__) \
	|| defined(_M_IX86) || defined(_M_IA64) \
	|| defined(_M_ALPHA) || defined(_WIN64)
#  define SR_LITTLE_ENDIAN
#else
#  error Unknown endianness
#endif

union IEEE_754_double
{
  double   a;
#ifdef SR_BIG_ENDIAN
  struct { unsigned sign : 1;
           unsigned exp  :11;
           unsigned high :20;
           unsigned low  :32;
         } b;
  struct { unsigned H    :32;
           unsigned L    :32;
         } c;
#else
  struct { unsigned low  :32;
           unsigned sign : 1;
           unsigned exp  :11;
           unsigned high :20;
         } b;
  struct { unsigned L    :32;
           unsigned H    :32;
         } c;
#endif
};

union IEEE_754_float
{
  float    a;
  struct { unsigned sign : 1;
           unsigned exp  : 8;
           unsigned high :23;
         } b;
  unsigned c;
};

inline
void
show( IEEE_754_double* p)
{
	std::cout << std::endl;
	std::cout << std::hex << std::setw(8) << std::setfill('0') << p->c.H;
	std::cout << ' ';
	std::cout << std::hex << std::setw(8) << std::setfill('0') << p->c.L;
	std::cout << std::endl;
}

inline
void
show( IEEE_754_float* p)
{
	std::cout << std::endl;
	std::cout << std::hex << std::setw(8) << std::setfill('0') << p->c;
	std::cout << std::endl;
}

#endif // SR_IEEE_754_UNIONS_H

Example.cpp

/************************************************************************/
/*					    copyright by twinklingstar                       */
/*					       www.twinlingstar.cn                           */
/*					           2015.01.30					             */
/************************************************************************/
#include <iostream>
#include "IEEE_754_unions.h"

int main()
{
	float f1 = 100.1f;
	IEEE_754_float* p1 = reinterpret_cast<IEEE_754_float*>(&f1);
	show(p1);

	float f2 = 0.0f;
	IEEE_754_float* p2 = reinterpret_cast<IEEE_754_float*>(&f2);
	show(p2);

	float f3 = -0.0f;
	IEEE_754_float* p3 = reinterpret_cast<IEEE_754_float*>(&f3);
	show(p3);

	return 0;
}

参考

[1]      IEEE. "IEEE Standard for Binary Floating-Point Arithmetic." website <http://homepages.math.uic.edu/~jan/mcs471/Lec3/ieee754.pdf>.

[2]      Dan Zuras, et al. "IEEE standard for floating-point arithmetic." IEEE Std 754-2008, pp.1-70, 2008.

[3]      Christer Ericson. Real-time collision detection. Amsterdam/Boston: Elsevier, 2005.

[4]      James M. Van Verth, and Lars M. Bishop. Essential Mathematics for Games and Interactive Applications: A Programmer's Guide. CRC Press, 2008.

[5]      David A. Patterson, and John L. Hennessy. Computer organization and design: the hardware/software interface. Newnes, 2013.

[6]      David Goldberg. "What every computer scientist should know about floating-point arithmetic." ACM Computing Surveys (CSUR), vol.23, no.1. pp.5-48, 1991.

[7]      John F. Reiser, and Donald E. Knuth. "Evading the drift in floating-point addition." Information Processing Letters, vol.3, no.3, pp.84-87, 1975.

[8]      Jonathan Richard Shewchuk. "Adaptive precision floating-point arithmetic and fast robust geometric predicates." Discrete & Computational Geometry, vol.18, no.3, pp.305-363, 1997.

[9]      Johnathan Richard Shewchuk. "Robust adaptive floating-point geometric predicates." Proceedings of the twelfth annual symposium on Computational geometry, ACM, 1996.

spacer

One comment on “浮点数误差

Leave a reply