我有 R 版本 4.1.2 (2021-11-01)。
当输入数字具有大量十进制值时trunc()
函数似乎不一致。
trunc(3.99999999999999977799999999999999999999900)
[1] 4
trunc(3.999999999999999777999999999999999999999000)
[1] 3
或
trunc(3.9999999999999997778888888888880)
[1] 4
trunc(3.999999999999999777888888888888)
[1] 3
我不确定是什么导致了这种不一致?
这里有两个问题:正确的答案是什么,为什么R在不同的情况下得到不同的答案?
这个数字 3.999999999999999779999999999999999999999 显然非常接近 4。 事实上,它比任何其他 IEEE-754 双精度浮点数更接近 4。下一个较低的可表示数字约为 3.9999999999999995,距离稍远一点。 所以,严格来说,trunc(3.999999999999999777999999999999999999999)
应该是trunc(4.0)
,这显然是4。 也就是说,当 R 采用输入 3.9999999999999997799999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999 这看起来"错误",因为你和我可以清楚地看到截断 3.999...应该给出 3,但并非每个实数都可以在有限精度浮点表示中表示这一事实偶尔会导致这样的异常。 (另请参阅这三个问题,它们收集了SO对这些二进制浮点异常的规范答案。
对于本答案的其余部分,我们将离开使用 R 的领域,深入研究实现R 的世界。 (我是一个C程序员,而不是R用户,这个答案很可能背叛了这种偏见,因为它对R的任何细微差别一无所知。 对此表示歉意。 但是,无论如何,R是用C编写的,并且使用C的double
类型进行大部分算术。 在当今绝大多数流行的通用计算机上,Cdouble
都是使用 IEEE-754 双精度实现的,这就是为什么我通过提到该标准来引导这个答案。
但是为什么R会根据有多少尾随0得到不同的答案呢? 答案在于 R 解释器深处的一个函数,该函数将用户键入的字符实际转换为内部 R 数据结构。
我们如何将像"123.456"这样的字符串转换为其内部浮点表示形式? 一种方法是暂时忽略小数点并将其转换为整数,得到数字123456
,然后计算小数点后的位数,然后除以十的幂。 事实上,123456 ÷ 10³
是123.456。
但是使用该策略,转换
3.99999999999999997779999999999999999999900 将涉及取一个 42 位数字并将其除以 10 41,而转换 3.9999999999999999977999999999999999999999000 将涉及取一个 43 位数字并将其除以 10 42。这些数字都不能在二进制浮点中精确表示。 他们会有点偏差,这有时会导致差异。 特别是,当数字如此之大时,不能保证÷b会给你与10a÷10 b完全相同的答案。
对于当前示例,差异在于一个除法导致一个接近 4 的数字,一个除法导致一个接近 3.9999999999999995 的数字。 (而且,请记住,我说的是 R 解释器深处在 C 代码中发生的划分,而不是您认为在 R 中所做的任何划分。
这里还涉及几个其他因素。(特别是,R使用"二进制幂"来计算 10N,这最终也会有所作为。 我现在没有时间写这些细节;也许以后。 有兴趣的读者可以查阅 R 源代码分发中的文件src/main/util.c
,特别是函数R_strtod5
。
但带回家的教训是,在二进制浮点数和人类可读的 10 进制表示之间准确来回转换是很困难的。 除此之外,获得正确的舍入结果通常需要以一些更高精度的表示形式进行计算,这样你才能得到一些勉强足够精确的东西,可以在最后四舍五入(也就是说,产生所需的"正确舍入结果")。具有讽刺意味的是,R的实现试图在这方面做正确的事情,使用C的long double
类型计算两个数字(即要除的两个数字)。 我本以为这足以避免这样的异常情况,但显然不是。
这也值得报告为 R 中的错误。 一个真正高质量的strtod
实现不会有这样的异常,并且在走上了实现自己的路线之后,R(我会说)需要重新发明任何必要的轮子,以便在所有情况下获得适当的全面结果。
为了补充@SteveSummit的出色答案,让我们将这两个数字存储在它自己的变量中,并在应用之前查看它们的外观trunc()
,将它们打印到最大可用精度:
x1 <- 3.99999999999999977799999999999999999999900
x2 <- 3.999999999999999777999999999999999999999000
print(x1, digits = 22)
## [1] 4
print(x2, digits = 22)
## [1] 3.999999999999999555911
如果您想查看@SteveSummit所指的确切代码(通过将前一个值的 10* 连续添加到下一个数字来找到 n 位数字,然后将适当的次数除以 10),它在这里......
我不确定是什么导致了这种不一致?
执行薄弱。
3.999999999999997799999999999999999999900(0)有什么特别之处?
使用浮点数的常见编码,4.0 及其后续和前面的值正好是:
// As decimal As hexadecimal
4.00000000000000088817841970012523233890533447265625 0x4.0000000000002
4.0 0x4.0
3.999999999999999555910790149937383830547332763671875 0x3.ffffffffffffe
值之间的值以及 OP 的常量(使用更广泛的数学运算)是:
// As decimal As hexadecimal
4.000000000000000444089209850062616169452667236328125 0x4.0000000000001
3.9999999999999997779553950749686919152736663818359375 0x3.fffffffffffff
3.99999999999999977799999999999999999999900 Let us call this C1
3.999999999999999777999999999999999999999000 Let us call this C2
1 23456789 123456789 123456789 123456789 123 Significant digit count
很明显,选择OP的两个常数(C1
42位和C2
43位)来测试OPR
的文本到浮点值的转换。
在完美的文本到浮点值转换中,文本C1
,C2
将转换为更接近的 4.0。 由于C1
转换为较小的 3.9999999999999995559...只是反映了 R 实现质量的弱点。
任何文本"3.999999999999999555910790149937383830547332763671875"或更多(以及≤"4.0000000000000000444089209850062616169452667236328125")应变为 4.0。