很高兴和你相遇
这里正在记录我的所思所学
首页 标签 归档 播客 关于

眼见不为实:为什么29不等于29

浮点数计算和比较的坑

29怎么会不等于29

今天有一个小伙伴给我发来了一个截图如下。我下意识的说了一句:0.58*50不等于29。

感概,有时候,眼见也不一定为实。这就是一个浮点运算导致的「算不准」问题。

在历史上,其实已经有过类似因为浮点运算导致的事故。

1994 年,英特尔的奔腾微处理器芯片的浮点计算单元出现了一个 Bug。在某些情况下,浮点除法并不能返回正确的结果。例如 4195835.0/3145727.0 产生的是 1.33374 而不是 1.33382。虽然这个 Bug 可能对绝大多数人没有影响,但是市面上 500 万左右流通的芯片都存在该缺陷,随着整个事件的发酵,最终这个 Bug 给 intel 造成了 4.75 亿美元的损失。这次事件也被称为「Pentium FDIV bug」而载入了 Bug 史册。

R 语言中的实际示例

我们分别在 R 中计算 0.58*50,0.29*100,然后让他们和 29 互相比一比。

(a=29)
(b=0.58*50)
(c=0.29*100)
a==c
a==b
b==c

结果如下所示:

> (a=29)
[1] 29
> (b=0.58*50)
[1] 29
> (c=0.29*100)
[1] 29
> a==c
[1] FALSE
> b==c
[1] TRUE
> a==b
[1] FALSE

我们看到屏幕上显示的 a , b和 c 都是29,但是如果用==进行一下比较,发现0.58*500.29*100竟然都不等于29。嗯,很迷。

浮点数的精度和计算问题,在任何编程语言里都类似。

浮点数在计算机中不能以任意精读存储,都会有一个准确性的限制,通过有限的连续字符来保存浮点数。

目前常见的浮点格式包括:单精度,双精度和扩展双精度。他们的准确性从前到后依次提高,其中单精度可以用于一般计算,双精度用于科学计算。在 R 中就是采用了双精度来对浮点数进行保存。如下所示:

> (d <- 10)
[1] 10
> class(d)
[1] "numeric"
> typeof(d)
[1] "double"

我们让 d 等于 10,通过 class 可以看出它是一个数值型变量,如果用 typeof 查看,就会发现他的存储模式(storage mode)是一个双精度浮点型数字。

这里又涉及到了存储模式和类型的问题。在R中class返回的是一个对象的高级类,而typeof返回的是对象的内部低级类。如果用calss查看一个数据框,其返回结果就是他的类(calss)为数据框,如果使用typeof查看,返回的结果就是内部类型为列表,如下所示。(嗯,想起了果子老师课上讲的,data.frame 就是一种特殊的 list)

> class(iris)
[1] "data.frame"
> typeof(iris)
[1] "list"

继续回来说浮点数的比较。因为精度问题,在R语言中使用 == 进行浮点数的比较是非常危险的。

就像文章开头的例子,此时假设你有一个 if 语句,想进行一下结果的判断。如果判断条件是 0.29*100 == 29,就会返回 F。因为 0.29*100 的实际值在计算机中比 29 要小一些。当你使用 as.interger去处理的时候,结果为28。

> as.integer(0.29*100)
[1] 28

这也就是为什么rep(1,29)有29个数字,但是rep(1,0.29*100)有28个数字。

如何避免浮点数计算的坑

如何解决这个问题呢?在 R 语言中有一个函数叫做 all.equal,虽然名字是 equal,但是他其实是 almost equal。

这让我不禁想到了一款果子和洲更(当然还有我)都在用全面屏 Almost 手机。

all.equal() 的功能是 Test If Two Objects Are (Nearly) Equal。也就是它可以略微忍受两个对象之间有一丢丢差别。如下所示:

> (a=29)
[1] 29
> (b=0.58*50)
[1] 29
> (c=0.29*100)
[1] 29
> a==b
[1] FALSE
> all.equal(a,b)
[1] TRUE
> a==c
[1] FALSE
> all.equal(a,c)
[1] TRUE

你可能会问,这个一丢丢究竟是多少呢?

通过帮助文档,我们可以发现这个一丢丢实际是tolerance = sqrt(.Machine$double.eps)。如果实际运行一下会发现.Machine$double.eps 等于2.220446e-16,那么sqrt(.Machine$double.eps)就是 1.490116e-08。

这里的.Machine是一个存储了 R 正在运行机器的数值特性的变量,其中包括了最小正浮点数(double.eps)。因此,以后想进行类似的比较,记得使用all.equal()。另外,多说一句,在 if 判断中不要直接使用它,而是写成 isTRUE(all.equal(....))

回到最开始的问题,如果确实需要 rep(1,rep(1,0.29*100))返回 29 个数字,就可能需要用到取整。

在 R 中有3种常见的取整,分别是:

  • 向下取整:floor()
  • 向上取整:ceiling()
  • 四舍五入:round()

此外还有一种不常用的「向0看齐」取整方法trunc()

他们的差别如下:

> ( x1 <- seq(-2, 4, by = .5) )
 [1] -2.0 -1.5 -1.0 -0.5  0.0  0.5  1.0  1.5  2.0  2.5  3.0  3.5  4.0
> floor(x1)
 [1] -2 -2 -1 -1  0  0  1  1  2  2  3  3  4
> ceiling(x1)
 [1] -2 -1 -1  0  0  1  1  2  2  3  3  4  4
> round(x1)
 [1] -2 -2 -1  0  0  0  1  2  2  2  3  4  4
> trunc(x1)
 [1] -2 -1 -1  0  0  0  1  1  2  2  3  3  4

因此通过rep(1,round(0.58*50))就可以得到和rep(1,29)一致的结果了。

写到这里并没有结束,因为如果用 python 算一下,就会发现 0.58*100 真的不等于 58。感觉R语言的对于浮点运算结果的展示还是有点问题?同时,0.58 和 0.59 还是不一样。

>>> 0.58*100
57.99999999999999
>>> 0.59*100
59.0
>>> 0.57*100
56.99999999999999
>>> 0.56*100
56.00000000000001

看来,要知其所以然,还得再仔细学习一下浮点精度和计算的实质。


参考资料:
R FAQ 7.31


本文作者:思考问题的熊

版权声明:本博客所有文章除特别声明外,均采用 知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议 (CC BY-NC-ND 4.0) 进行许可。

扫码即刻沟通

· 分享链接 https://kaopubear.top/blog/2019-12-19-rfloat/