格物致知-Floating Point

image

作者 | 曲健
类型 | 译文
全文 | 14791字, 案例术语代码较多,半小时读不完

之前陆陆续续写了很多架构、设计、思想、组织方向的文字,突然感觉到有些厌烦。因为笔者不断看到有些程序员“高谈阔论、指点江山”之余,各种定律、原则、思想似乎都能信手拈来侃侃而谈,辩论的场合就更喜欢扯这些大旗来佐证自己的”金身”。殊不知,这些人的底座脆弱到不堪一击,那些“拿来”的东西都是空中楼阁罢了。优秀程序员区别于其他的一项重要指标,就是基础知识的底蕴足够强大。靠看靠学靠实战靠日积月累,绝无捷径。

image

而今天要讲的”浮点运算”绝对是不少程序员知其表不知其内的一个知识点,甚至在某些人的编程世界里已经慢慢沦为鸡肋。比如在涉及精度必须准确的货币类财务运算中,要不更换单位*100转为Integer,要不使用BigDecimal,稍微有点”经验”的人这里都不会敢用Float/Double。

这里的“经验”指的是那种有一定的编程经验,但是理论知识结构有死角,在FP这个场景就是知道FP有精度丢失问题,但不知道为什么的那一类。

有人会问既然精度不准那浮点数还有什么价值, 直接用精度不会丢失的Integer或者BigDecimal就好了呗?这其实就是在数据精度**、数值范围、**运算速度之间的一个权衡,浮点数可表示的数值范围远大于整型,浮点数的运算速度也快于BigDecimal(现代的CPU都会集成FPU浮点运算单元)。在很多软件系统的实现里不用浮点数确实没问题,但是在一些游戏引擎、工程处理甚至天体物理方面的超大型数据运算,为了快速得到结果也防止数据范围溢出,只要误差在可接受范围内,浮点数当然可以是首选。

Lua 在 5.3 之前没有整型只有浮点数,你还敢做数值计算么?

与浮点数相对的就是定点数,Java里的整型, BigDecimal,Mysql的Numeric/Decimal类型本质就是定点数实现,用他们不会有精度问题(当然前提是别越界)。

https://dev.mysql.com/doc/refman/8.0/en/fixed-point-types.html

(知识点:十进制的定点数。

可研究下BigDecimal源码一目了然)

到现在我没解释过任何东西,有些看官可能越看越糊涂,没错,我还要再码一些代码,让大家更糊涂一些。文末还有大量的问答和练习题,可以自己去做,一次性把浮点基础全部夯实了。

下方截图的代码都在这里了:

https://github.com/NicholasQu/snippets/blob/master/src/main/java/cc/xiaoquer/data/types/FloatingPointDemo.java

FP知识点1

无穷的二进制表示及转换

浮点数除以0抛出什么异常?

image

FP知识点2

浮点运算的精度丢失

十进制到二进制的转换发生了什么?

image

FP知识点3

非规格化浮点数的处理性能损耗

什么是非规格化浮点数? 为什么它的运算会慢?

image

FP知识点4

不同精度的运算速度对比

小数的运算,为什么浮点数相比定点数快?

image

FP知识点5

浮点的取值范围 VS 可表示的个数

取值范围好理解,大家看一下Float.MAX_VALUE(3.4*10^38) 和 MIN_VALUE 即可。那么32 bits的Float可以表示多少个浮点数呢? (2^32次方而已),再次表达了离散不连续的数字世界,也可看出有限的浮点数在其近乎无限的取值范围内是有多么的稀疏。

翻译正文

现在开始

传统计算机科学与科学计算的一个显著区别就是离散数学(0和1)取代了连续数学(微积分为代表)的使用。将整数转为实数绝对不只是一种外表的变化。数字计算机无法精确表达每一个实数,所以当我们设计实数的计算机算法时不得不面临新的挑战。现在,除了分析运行时间和内存占用,我们必须更加关注最终解决方案的“正确性”。事实上这个极具挑战性的问题正在进一步恶化,因为为了适应离散结构的计算机,很多重要的科学算法不得不妥协做出很多额外的近似。正如我们发现一些离散算法自身实在是太慢了(算法时间复杂度:多项式 VS 指数型),我们将会看到一些浮点算法太不准确(稳定 VS 不稳定)。有时可以设计一个更聪明的算法来纠正这个问题。离散问题有时候时常需要面临的困难正是来自于它本身的痼疾(比如NP完全问题/NPC问题)。浮点数问题的困难同样来自于本身(病态性),例如,准确的长期天气预报。要成为一位有效的计算机科学家,我们必须有能力归类这些算法和相应的问题。

译者注
『有效的科学家』听起来很拗口,但实在没找到更合适的词汇来对应。
实际上真的有本书就是叫 The Effective Scientist:
https://www.cambridge.org/core/books/effective-scientist/EC65B924792BB2C4DAAADAF21FF0AC0F#fndtn-information

浮点数

如果数字计算机没有浮点运算能力,20世纪的很多伟大成就是不可能实现的。然而,大多数程序员都没有很好的理解这个源头上经常容易引起混乱的课题。在1998年2月一个名为”Java数值计算扩展”的主题演讲中,James Gosling宣称“95%的人对浮点数一无所知”。但其实浮点背后主要的思想并不是那么难, 接下来我们将一一揭秘那些困扰着大多数新手的疑惑。

IEEE 754二进制浮点表示法。首先,我们将描述浮点数是如何表示的。Java使用IEEE 754二进制浮点标准的一个子集来表示浮点数并定义算术运算的结果。几乎所有的现代计算机都符合这个标准。一个浮点数用32位表示,每个可能的位组合代表一个实数。这意味着即使有无穷多个实数(甚至在0和1之间),最多也能精确地表示232个可能的实数。这涵盖了从±1.40129846432481707±3.40282346638528860 e+38 e-45。具有6或7位有效的小数,包括正无穷、负无穷和NaN(非数字)。这个数字包含一个符号位s(表示正负),8位表示指数e, 23位表示尾数m。

IEEE 754

IEEE 754二进制浮点数表示/标准(后文统一用IEEE代指)。首先, 我们将描述一下浮点数应该如何表示。Java使用IEEE标准的一个子集来表示浮点数和定义算术运算的结果。几乎所有的现代电脑都符合这个标准。浮点数使用32位表示, 比特位的每一种可能的组合都代表一个实数。这意味着可以被精确表达最多2^32种实数, 但其实实数是无限多的(甚至在0和1之间)。

IEEE标准使用一种类似于科学记数法的内部表示, 只是用的是二进制不是十进制罢了。它可以覆盖从±1.40129846432481707e-45 到 ±3.40282346638528860e+38的实数区间,可确保6到7位的有效小数位。

译者注
这里的底数e并不是那个自然常数e,它其实就是10,这也叫E-表示法。
这种方法更多的是为了区分前后两部分数字,方便书写并理解直观,不然你换成10试试,是不是都拧巴到一块了…
https://en.wikipedia.org/wiki/Scientific_notation#E_notation

数字的精度是一般指的就是所有数字的个数,小数位数是小数点右边的数字个数。例如123.45 的精度是 5,小数位数是2。*

浮点的精度由尾数的位数决定,单精度浮点数尾数是23位,取值范围为0~2^23,而2^23=8388608≈10^6.92369,所以float的精度为6~7位小数,6位保证是绝对精确,7位一般也是正确的,8位再往后就完全不靠谱了。

这其中也包括三个特殊值:正无穷, 负无穷, NaN(Not-a-Number 不是一个数字)。它由1位符号位(可理解为正负), 8位的指数e,和23位的尾数m组成。 十进制表达公式如下:

image

感谢提供的在线Latex功能https://www.codecogs.com/latex/eqneditor.php*

image

  • 符号位 s(31位). 这一位最显著的符号位表示数值的正负,1为负数,0为正数。

  • 指数位 e (30 - 23位). 接下来的8位就是指数位啦(译者注:无符号整数0~255)。按照约定指数需做个127的偏差。也就是,要表示二进制指数5,偏差127加上5等于132(译者注:这个值才是e),132的二进制编码是10000100。要表示二进制指数-5,偏差127加上-5等于122,对应二进制编码是01111010。我们知道通常负数采用的是补码表示法,而这个约定(127偏差)可视作补码的一种替代方案。

译者注
建议大家再复习一下三个概念:原码、反码(1’s / ones’ complement)、补码(2’s / two’s complement)的知识。要展开讲透的话不重开一篇基本不现实,这里只提几个要点先:

1. 数值是有正负的,也就意味着有符号位的存在。而计算机为了简化加减乘除运算的复杂度,需要想办法让符号位也参与运算,让减法转变为最基本的加法运算。

2. 反码可以视作补码的一个中间步骤,因为补码就等于反码+1。现而今大家都知道大部分主流计算机都是使用补码方式来存储数字的,而80年代之前的计算机其实用的是反码。另外,原本的英文名称比中文翻译似乎更能体现两者关系。(Tips: 反码的英文维基上是ones’,很多其他英文文献用的one’s,从Redirected from这个表述来看,两者应该是一致的,但为何选了ones’不是one’s呢?尝试各种搜索,原因仍不可得)

https://en.wikipedia.org/wiki/Ones%27_complement
https://en.wikipedia.org/wiki/Two%27s_complement
image

3. 补码运算有个关键的“模”的概念,尝试以钟表的指针转动就是不断做加法来类比,钟表的模就是12,超过12小时就从0重新开始计数。补码的进位(Carry Bit)一旦溢出就会被丢弃。

4. 对于8 bits表示的数字来说,补码求补运算使用的偏移量是 2^8=256,而浮点数指数位e的偏移量使用了中间数127, 也就是IEEE在这里没有直接采用补码方式。这里也潜藏着更加深刻的数学意义,用于区分浮点和整型还是区分浮点和定点,大家可以看下方wiki理解.

https://en.wikipedia.org/wiki/Exponent_bias

  • 尾数 m(22 - 0位). 剩下的23位表示尾数,规格化为0.5到1之间。这种规格化依靠相应地调整二进制指数总是可能实现的。二进制的小数与十进制小数原理并没有什么不同:

译者注
后文为了方便直观理解,都会用后缀B和D区别表示二进制和十进制数。

0.1101B=1/2 + 1/4 + 1/16 = 13/16 = 0.8125D
当然也不是所有的十进制数都能被转换为二进制小数的,比如
1/10 = 1/16 + 1/32 + 1/256 + 1/512 + 1/4096 + 1/8192 + …
在这种情况下,十进制的0.1用最接近的二进制近似值来表示:即23位的二进制小数0.000110011001100110011…
进一步优化,既然尾数永远是以1开头的,就没有必要显性的存储这个隐藏位了。
举例说明:十进制小数0.085会被存储为二进制浮点数
00111101101011100001010001111011
image

套用前文十进制的表达公式:
image

译者注
乍看这个公式并不是完全套用上文的那个,这里简单解释下标红的两部分:
1. m需要+1,这在尾数部分说明过,因为尾数默认都是1.XXXX的小数,为了节省一位存储就把1省略了,转换的时候要记得加回来;
2. m实际上就是0.XXXX的小数,而直接转十进制是按照23位的整型来算的,所以这个整型需要再除以2^23用来移小数点位;
译者再提供一种对m计算方式可能理解起来更直观一些:
二进制的m实际应该是 “1”+”.”+”上图22-0位”
1.01011100001010001111011
image
红线标注部分是完全的二进制运算,而后2的负数次方相加和正文的整型除以2^23本质相同就是做个变换罢了。

双精度浮点double与单精度float相似,只不过double内部使用64位展示,其中11位的指数,52位的尾数,指数偏移量不再是127而是1023。double覆盖范围从±4.94065645841246544e-324±1.79769313486231570e+308,可保证14或15位(注)的有效精度。

译者注
按照前文单精度浮点的精度推算方法2^52=10^15.65356,这里应该是15-16位。大家理解如何推论就OK, 不需要死记硬背。

精确度Precision vs 准确度Accuracy

精确度=规格的紧密性,准确度=正确性(译者注:完全是只可意会的概念解释…请往下直接看例子会清晰的多),千万不要把精确度与准确度混为一谈。指定精度为10的数字3.133333333是对数学常量 π 的估算值,但它只有两位数字3.1是准确的。正如约翰•冯•诺依曼曾经说过的:“当你甚至不知道自己在说什么时,精确到什么程度都没有意义。” Java通常会以16或17位精度来打印浮点数,但不要盲目认为这么多位数字都是准确的!计算器通常显示10位数字,但计算精度为13位。哈勃太空望远镜的镜片以超高的精确度打磨的,但使用了错误的规格。因此,既然它不能像预期那样拍摄高分辨率的图像,那么它就是一个巨大的失败。尽管如此,它的精度使宇航员能够安装一个校正透镜来抵消误差。货币计算通常是根据给定的精度来定义的,例如,欧元汇率的报价必须是6位数。

译者注
哈勃空间望远镜镜片的精度问题造成巨大事故的著名案例
https://zh.wikipedia.org/wiki/哈勃空间望远镜
Hubble’s Mirror Flaw
https://www.nasa.gov/content/hubbles-mirror-flaw

舍入误差

对于不熟悉浮点数的人来说,浮点数编程是一个令人困惑和危机四伏的过程。整数运算是精确的,除非运算结果不在整数的范围内(溢出)。相比之下,浮点运算并不精确,因为一些实数的表示需要无限的数字,例如数学常量e、π和1/3。然而,大多数初级Java程序员惊讶地发现,在标准二进制浮点数中,1/10也不是完全可以表示的。这些舍入误差会以非直观的方式在计算过程中不断传播。

小数的四舍五入

例如,下面来自于FloatingPoint.java

(https://introcs.cs.princeton.edu/java/91float/FloatingPoint.java.html)

的代码片段,第一个打印为false,而第二个打印true。

image

牛顿迭代法

来个更实际的例子,下面的代码采用牛顿迭代法来实现计算c的平方根。数学上, 在 t*t - c > 0 的前提下基于上一轮数据不断迭代可最终收敛于√c。然而, 浮点数只有有限的几个比特的准确度,所以最终我们可能期望t的平方等于c,直至达到机器的精度。

image

对于某些值的c,该方法”确实”是有效的。

image

这可能会让我们相信我们的代码实现是正确的。但当我们试着计算20的平方根时,一个令人惊讶的事情发生了。我们的程序将陷入无限死循环! 此类错误即称为**舍入错误机器准确度**可以这样理解:若存在一个最小的数ε可以使(1.0 +ε! = 1.0),那么数ε就是机器准度。将容错误差ε改为很小的正数会有所帮助,但没从根源上解决这个问题。

我们必须满足于近似平方根。一个可靠的计算方法是选择某个容错误差ε,比如1E-15,然后再试着找到一个值t可以使 |t - c/t| <εt。我们用相对误差代替绝对误差,否则程序可能进入无限死循环。

image

调和级数求和

译者注
这部分翻的我想死,太难了…所以我决定放弃这段!
H(n)表示调和级数,是发散序列,极限值是无穷,也就是不存在。但是调和级数的拉马努金求和(Ramanujan summation)是存在的。 这块内容超纲了…

(总结一下) 每次执行算术运算, 你至少会引入一个额外的误差ε。Kernighan和Plauger说:“浮点数就像一堆沙子; 每次你移动, 就会失去一些沙子, 捡起一些泥土”。如果错误是随机发生的, 那我们可以预见对N开平方根的误差会不断累加。然而, 如果我们在设计数值算法的时候不保持警惕, 这些错误将以非常不好和不直观的方式传播, 导致错误累积甚至更糟。

财务计算

我们举例来说明可能破坏财务计算的那些舍入误差。涉及美元和美分的财务计算以10为基数。下面的示例演示了使用二进制浮点系统(如IEEE 754)的一些风险。

销售税

源码Financial.java说明了这种危险

(https://introcs.cs.princeton.edu/java/91float/Financial.java.html)

计算一通50美分的电话需要征收9%的销售税(或0.70美元的通话征收5%的销售税)。使用IEEE 754标准, 套入公式 1.09 * 50 = xxx。这个结果四舍五入到0.xx而实际上的确切答案却是0.yy, 因为按规定电话公司必须使用银行家舍入法(Banker’s rounding)。

译者注
Banker’s rounding银行家舍入法
不同于我们传统意义的那种四舍五入,银行家舍入法又有个直白叫法四舍六入五奇偶。具体规则有个口诀:四舍六入五考虑,五后非零就进一,五后为零看奇偶,五前为偶应舍去,五前为奇要进一
下面的例子感受一下:
round(0.5) = 0
round(1.5) = 2
round(1.4) = 1
round(1.6) = 2
round(2.4) = 2
round(2.5) = 2
round(2.6) = 3

Round up 向上舍入 vs Round down 向下舍入
这两者相对比较好理解,大家可以参见一下Excel里的这两个公式,算法相同。这个翻译我就复用Excel里的说法的,免得引起无谓的理解门槛,其实向上应该叫直接进位,向下就是直接舍位
roundup(0.14, 1) = 0.2
roundup(0.15, 1) = 0.2
roundup(0.16, 1) = 0.2
rounddown(0.14, 1) = 0.1
rounddown(0.15, 1) = 0.1
rounddown(0.16, 1) = 0.1

相比之下,对一个75美分的电话征收14%的税向上舍入到得到x.xx,尽管按规定电话公司必须使用银行家舍入法将其向下**舍入。这种便士(美分)级别的差异可能看起来不太明显,您也可能指望这些影响最终会相互抵消。然而,放大到数亿交易量级,这种香肠切片可能导致数百万美元。参考电影:《超人3》(1983年)、《黑客》(1995年)和《上班一条虫》(1999年)。在《上班一条虫》里,三个人在银行会计系统里植入一种电脑病毒,会将每笔交易做向下舍入,*也就是每笔交易里的小数零头部分转移到他们自己的账户中。(日积月累也发财了…)*

译者注
Salami Slicing 萨拉米香肠切片
https://en.wikipedia.org/wiki/Salami_slicing
指的是一系列许多小动作,通常由秘密手段进行,随着不断的累积整体会产生更大的动作或结果,这些动作或结果通常会长期实施或非法实施。该术语通常用于贬义。尽管萨拉米香肠切片经常被用来进行非法活动,但它只是一种通过以较小的增量累积它来获得优势的策略,因此它也可以以合法的方式使用。

这个定义有些晕,上文在《上班一条虫》中提到的极少量资金的反复盗窃,不断累积达到巨大危害的就是一个典型案例。在信息安全、政治、学术界都有相应的场景案例对应。

出于这个原因,一些程序员认为应该始终使用整型来存储金融值,而不是浮点类型。下一个示例将告诉你使用int类型存储财务值的风险。

复利

此示例介绍舍入误差的危险。假设你以5%的利息投资1000美元,每日复利。一年后你会得到多少钱?如果银行正确计算了价值,那么你应该得到 $1051.27,因为确切的公式是

image

此公式运算后的实际结果是1051.2674964674473…,假设银行将您的余额存储为整数(以美分为度量)。在每天结束时,银行会计算你的余额乘以1.05,并将结果四舍五入到最接近的美分。如此你最终只会得到1051.10美元,并被骗了17美分。再假设银行在每天结束时向下舍入到最近的美分,那么你将只有1049.40美元,被骗走1.87美元。不存储美分小数位的微小错误会累积并最终明显放大,甚至是成为虚假的。

CompoundInterest.java源码

(https://introcs.cs.princeton.edu/java/91float/Financial.java.html)

你应该使用的是Java里的BigDecimal库,而不是整数或浮点。该库有两个主要优点。首先,它可以准确地表示十进制小数,这可以防止使用二进制浮点的销售税问题。其次,它可以以任意精度存储数字,这使程序员能够控制舍入误差对计算的影响程度。

其他错误来源

除了使用浮点算法时固有的舍入误差之外,在科学应用中还经常出现很多不同类型的近似误差问题。

测量误差

在日常的计算过程中使用的原始数据本身就是不准确的。这种情况经常源于实际(不准确或不精确的测量仪器)和理论(海森堡测不准原理)的考虑。此种场景下,我们更关注的是模型和解决方案,而它们的结论对原始数据并不是非常敏感。

离散误差

另一个不精确的来源来自于离散连续系统,例如,通过截断其泰勒展开以逼近超越函数,使用矩形的有限和近似积分,找到微分方程的近似解的有限差分方法,或基于格估计连续函数。

即使使用再精确的算法,离散误差依然存在,这通常是不可避免的,但我们可以通过使用更精细的离散化来减少截断误差。当然,这需要以使用更多资源为代价的,无论是内存还是时间。在实际应用中,离散误差往往比舍入误差更重要。

统计误差

没有足够的随机样本。

灾难性消除

当通过加法或减法从大的数计算小的数时,精确度损失很大。 举个例子,如果x和y两个大数字除了最后几位数字不同其他都完全一样,那么如果我们计算z = x-y,z可能只有几位精度(译者注:只有最后几位不同,所以可能无限接近于0)。 如果我们在随后的计算中使用z,那么结果同样只可能有几位精度。

绘制函数曲线

试着绘制函数

f (x) = (1 - cos (x) / (x^2)

x 区间[410^8, 410^8]

在这个区域,数学函数f(x)近似为常数0.5。然而在IEEE浮点中,它绝对不是! 代码**Catastrophic.java*(https://introcs.cs.princeton.edu/java/91float/Catastrophic.java.html*)向大家展示了一些令人意想不到的结果。

image

指数函数

现在我们考虑灾难性消除的一个更微妙更有害的后果。代码Exponential.java(https://introcs.cs.princeton.edu/java/91float/Exponential.java.htm)使用泰勒级数计算e^x。

e^x = 1 + x + x^2/2! + x^3/3! + x^4/4! + …

这个级数对于所有的x值都是均匀且绝对地收敛,但是对于x的许多负值(如-25),无论级数中有多少项相加,程序都得不到正确的数字。例如,当x = -25时,该级数以浮点形式收敛到-7.129780403672078E-7(一个负数!),但真正的答案是1.3887943864964021E-11。要知道为什么,观察总和中的第24项是 25^24/24! 第25项是 -25^25/25! 原则上它们应该完全抵消掉,而在实践它们会灾难性地相互抵消。这些项(5.7*10^9)的大小比正确答案大20个数量级,消除的任何误差在计算的结果中都会被放大。幸运的是,在这个场景里,问题很容易得到纠正。

数值分析

Lloyd Trefethen (1992) 定义了“数值分析是对连续数学问题算法的研究”。

此章节从略…有兴趣的自己读读吧…

数值灾难真实案例

我们上面讨论的例子相当简单。然而,这些问题在实际应用中会经常出现。如果处理不当,灾难会迅速降临。

阿丽亚娜5型火箭

阿丽亚娜5型火箭在欧洲航天局发射40秒后爆炸。经过10年和70亿美元的研发方才首航。根据传感器报告由于加速度过大导致程序中负责重新校准惯性制导的部分溢出。64位浮点数被转换为16位有符号整数,而该数字大于32,767,故而转换失败。通用诊断系统捕获了这个意外的溢出,并将调试数据转储到用于引导火箭发动机的内存区域。控制权被切换到备份计算机上,但这台备机存储了一模一样的数据。这导致了系统剧烈运行试图纠正一个不存在的问题,将发动机与他们的安装件分离,最终导致阿里亚娜5型的悲剧(http://www.around.com/ariane.html)。

爱国者导弹意外

1991年2月25日,一枚美国爱国者导弹未能追踪并摧毁一枚伊拉克飞毛腿导弹。相反它击中了一个军营,造成26人死亡。后来确定原因是由于用十分之一秒来度量时间导致计算的不准确,因为24位浮点不能准确地表示1/10。修复问题的软件于2月26日抵达达卡兰。(http://www.ima.umn.edu/~arnold/disasters/patriot.html)

奔腾电路中的英特尔FDIV BUG

英特尔于1994年7月发现,1994年9月由数学教授重新发现并公布。1994年12月英特尔的召回价值3亿美元。 1997年发现了另一个浮点错误。

斯莱普纳钻井平台沉没

1991年8月,价值7亿美元的石油和天然气开采平台发生泄漏,沉入北海。在不准确的有限元近似中,误差低估了剪应力达47%。

温哥华证券交易所

经过22个月的累计四舍五入误差,温哥华证券交易所指数被低估了50%以上。一个显而易见的算法是将所有股票的价格加起来,而“聪明”的分析师认为,在每笔交易后加上一只股票的净变化来重新计算指数会更加有效率。这个计算是使用四位小数和截断(而不是四舍五入)结果到三位来完成的。(http://www.ualberta.ca/~carolina/CHE374/notes/U01_errors/flerrors/fl_and_errors.html)

教训

  • 使用double代替float来提高准确性。

  • 仅在您确实需要节省内存时才使用float,并且准确地知晓相关风险。通常它不会使事情变得更快,反而偶尔会使事情变得更慢。

  • 小心计算两个非常相似的值的差异,已经在随后的计算中使用两者的中间结果。

  • 合计两种不同量级的数据时要小心。

  • 多次重复略微不准确的计算时要小心。例如,计算行星随时间的位置变化。

  • 设计稳定的浮点算法非常重要。可能的话尽量使用现成函数库。

问答时间

问:浮点数有什么好的参考物么?

答:这里有两篇关于浮点精度的文章:

David Goldberg的《每个计算机科学家应该知道的浮点运算》http://docs.sun.com/source/806-3568/ncg_goldberg.html

以及由图灵奖得主William Kahn共同撰写的《Java的浮点如何做到四处伤害每个人》http://www.cs.berkeley.edu/~wkahan/JAVAhurt.pdf

这是维基百科关于数值分析的条目http://en.wikipedia.org/wiki/List_of_numerical_analysis_topics

问:如何将IEEE位表示转换为双精度?

答:要在Java中获取双变量x的IEEE 754位表示,请使用Double.doubleToLongBits(x)。根据代码项目,要获得大于x的最小双精度数(假设x是正的和有限的)是Double.longBitsToDouble(Double.doubleToLongBits(x)+ 1)。

这里提供一个不错的二进制与Double转换的在线网站:

http://www.binaryconvert.com/result_double.html

image

问:有没有直接的方法来检查整数类型的溢出?

答:不可以。整数类型不以任何方式表示溢出。只有当分母为零时,整数做除法和取余会抛出异常。

问:如果我输入一个太大的数字,例如1E400,会发生什么?

答:Java返回错误消息“浮点数太大”。

image.png

问:(上面两个问题换成)浮点类型会怎么样呢?

答:溢出操作会得到正无穷或负无穷,下溢操作会导致正负零,数学上没有明确定义值的操作设置为NaN(不是数字),例如0/0,sqrt(-3),acos(3.0),log(-3.0)。

问:如何测试我的变量是否具有NaN值?

答:使用方法Double.isNaN()。请注意,NaN是无序的,因此涉及一个或两个NaN的比较操作<,>和==始终计算为false。任何涉及NaN的不等于的比较都为真,甚至NaN != NaN。

问:-0.0和0.0之间有什么区别?

答:两者都是数字零的表示。0.0==-0.0返回true。然而,1/0.0返回正无穷,而1/-0.0输出负无穷。代码NegativeZero.java(https://introcs.cs.princeton.edu/java/91float/NegativeZero.java.html)

给出了一个令人惊讶的反例来反驳这个误解:

如果x == y,那么f(x) == f(y)

log 0 = -infinity vs log(-1)= NaN。 log(ε)vs log(-ε)。

问:听说程序员不应该用完全的等式来比较两个实数,而应该总是使用测试的方法,比如

|a-b| < ε 或**者**

|a-b| < ε * min(|a|,|b|)

这是正确的吗?

答:这取决于具体情况。通常首选相对误差方法,但如果数字非常接近零那就失效了。绝对值方法可能会产生意想不到的后果,而ε使用什么值也不是非常明确。如果a=+ε/4,b =-ε/ 4,我们应该认为它们是相等的。从传递性角度却不同:如果a和b是“相等的”,b和c也是“相等的”,却无法证明a和c一定是“相等的”。

问:Java如何打印双精度数?

答:通过将所有指数位设置为1。参见java.lang.Double API

(https://docs.oracle.com/javase/6/docs/api/java/lang/Double.html#toString(double))。它始终在小数点后打印至少一位数。之后,它根据需要使用尽可能多的数字(但不会很多)来区别最接近的可表示双精度数。

问:使用IEEE 754如何表示零,无穷和NaN?

答:通过将所有指数位设置为1。正无穷=0x7ff0000000000000(所有指数位1,符号位0和所有尾数位0),负无穷=0xfff0000000000000(所有指数位1,符号位1和所有尾数位0),NaN=0x7ff8000000000000(所有指数位1,至少设置一个尾数位)。正零=所有位0,负零=所有位0,除了符号位1。

问:使用双精度IEEE标准存储0.1时所表示的确切值是多少。

答:它是0.1000000000000000055511151231257827021181583404541015625。您可以使用System.out.println(new BigDecimal(0.1));亲眼看看。

问:什么是StrictMath?

答:Java的Math库保证其计算结果距离准确答案在1到2个ulps内。而StrictMath库可以保证所有计算结果都精确到到1/2个ulp以内。这无非是经典的速度与准确性两者的权衡。

译者注
ULP(Unit in the Last Place)
[Java Glossary]**http://mindprod.com/jgloss/ulp.html
[WIKI]**https://zh.wikipedia.org/wiki/最后一位上的单位值

在计算机科学与数值分析中,最后一位上的单位值或称最小精度单位,缩写为ULP,是毗邻的浮点数值之间的距离,也即浮点数在保持指数部分的时候最低有效数字为1所对应的值。ULP用于度量数值计算的精度。例如:圆周率位于毗邻的双精度浮点数3.1415926535897927与3.1415926535897936之间。

1.0和2.0之间有多少个数?在数学中是无限的,但是在计算机中只能是有限的(资源有限),所以就有了ulp。java.lang.Math.ulp(double d) 返回的double ulp是该浮点值与下一个数值较大浮点值之间的正距离。

对于非NaN x, ulp(-x) == ulp(x)

问:什么是strictfp修饰符?

答:在声明类或方法时可以使用strictfp修饰符(译者注:strict float point精确浮点数的意思)。这可确保浮点结果在不同JVM之间逐位精确。如果结果溢出,IEEE标准允许处理器使用更高的精度执行中间计算。 strictfp要求将每个中间结果截断为64位双精度格式。由于英特尔奔腾处理器寄存器使用IEEE 754的80位双扩展格式运行,因此这可能是一个重大的性能影响。当然只有极少数人才用得到这种模式。

问:编译器标记javac -ffast-math有什么作用?

答:它放宽了IEEE的一些舍入要求, 使一些浮点计算更快。

问:整数是否总是可以使用IEEE浮点精确表示?

答:是的,除非52位尾数溢出。使用double不能精确表示的最小正整数是2^53 + 1 = 9,007,199,254,740,993。

问:当a和b是浮点数时,(a + b)总是等于(b + a)吗?

答:是的。

问:x/y 总是等于相同的值,与平台无关吗?

答:是的。IEEE要求精确执行操作加减乘除,然后四舍五入到最近的浮点数(如果存在平局,则使用银行家舍入法: 舍入到最接近的偶数)。这以牺牲效率为代价提高了可移植性。

问:(x-y) 总是等于 (x+(-y)) ?

答:是的。由IEEE 754保证。

问:-x总是与0-x相同吗?

答:大部分情况是的,除非x=0则-x=-0,但0-x=0。

问:这三个等式成立么?

x+x==2*x

1*x==x

0.5*x==x/2

答:是的,由IEEE 754保证。

问:x/1000.0总是等于0.001*x吗?

答:不可以。例如,如果x = 1138,它们是不同的。

问:如果x!=y,则z=1/(x-y) 保证始终不会发生除零的情况。

答:是的,由IEEE 754保证(使用非规范化数字)。

问:(x>=y) 和 !(x<y) 相同意思吗?

答:不相同,在x和y一个或两者为NaN的时候。

问:为什么不使用十进制浮点代替二进制浮点?

答:十进制浮点比二进制浮点提供了几个优点,特别是对于财务计算。尽管如此,它通常需要超出大约20%的存储空间(假设它使用二进制硬件存储)并且代码运行出结果会有些慢。

问:为什么不使用定点表示代替浮点?

答:定点数在小数点后有固定的位数,可以使用整数运算表示。浮点有一个基于精度的滑动窗口,提供了大的动态范围和高精度。固定点数用于一些没有FPU(省钱)的嵌入式硬件设备,例如音频解码或3D图形。适用于数据受限于某个范围的情况。

问:Java中数字的任何好资源?

答:Java numerics(https://math.nist.gov/javanumerics/)为Java中的数值计算提供了信息的焦点。 JSci(https://sourceforge.net/projects/jsci/):Java中用于数值计算的免费科学API。

问:为什么Java的Math.sin和Math.cos函数比它们的C函数慢?

A.在 -pi/4 到 pi/4 的区间内,Java使用硬件级sin和cos函数,因此它们与C的时间大致相同。超出此区间的参数必须通过模数pi转换到此区间内。正如James Gosling在x87平台上的博客一样,硬件sin和cos使用近似值,这使得计算速度很快,但它并达不到IEEE要求的准确度。C实现会对所有参数使用硬件级sin和cos,牺牲速度以换取IEEE的一致性。

尾语

原文文末还有大量的练习题(可点击原文查看),有兴趣的同学一定要亲自去做做,对于理解上面大段理论有极大益处,顺带不要忘记篇首的几个知识点哟。当然如果对进制转换和精度舍入等等原理都已吃透,这些问答和练习题都不在话下。

古人云“格物致知”, 就是要我们穷究事物的真理达到知其然知其所以然的通透感,想来必是极好的。

- 致“天上飘着”的程序员们

(读者若是同时了解大学、朱熹、王阳明各种思想的神人,可能对这个词有不同见解,咱不在这辩论。反正现代汉语词典解释的是:研究事物原理而获得知识。语出《礼记·大学》:欲诚其意者,先致其知,致知在格物)

You forgot to set the qrcode for Alipay. Please set it in _config.yml.
You forgot to set the qrcode for Wechat. Please set it in _config.yml.
You forgot to set the business and currency_code for Paypal. Please set it in _config.yml.
You forgot to set the url Patreon. Please set it in _config.yml.
Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×