r - 在 R 中，为什么 factorial(100) 与 prod(1 :100)?) 的显示不同

Question

在 R 中，我发现了一些我无法解释的奇怪行为，我希望这里有人可以。我相信100的值(value)!这是 big number 吗？

控制台中的几行显示了预期的行为......

>factorial( 10 )
[1] 3628800
>prod( 1:10 )
[1] 3628800
> prod( as.double(1:10) )
[1] 3628800
> cumprod( 1:10 )
[1]       1       2       6      24     120     720    5040   40320  362880 3628800

但是，当我尝试 100 时!我得到(注意结果数字如何开始相差大约 14 位数字):

> options(scipen=200) #set so the whole number shows in the output
> factorial(100)
[1] 93326215443942248650123855988187884417589065162466533279019703073787172439798159584162769794613566466294295348586598751018383869128892469242002299597101203456
> prod(1:100)
[1] 93326215443944102188325606108575267240944254854960571509166910400407995064242937148632694030450512898042989296944474898258737204311236641477561877016501813248
> prod( as.double(1:100) )
[1] 93326215443944150965646704795953882578400970373184098831012889540582227238570431295066113089288327277825849664006524270554535976289719382852181865895959724032
> all.equal( prod(1:100) , factorial(100) , prod( as.double(1:100) ) )
[1] TRUE

如果我对设置为“已知”数字 100 的变量进行一些测试!然后我看到以下内容:

# This is (as far as I know) the 'true' value of 100!
> n<- as.double(93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000)
> factorial(100) - n
[1] -1902315522848807765998160811905210717565551993186466795054798772271710903343294674760811531554315419925519536152107160826913610179566298858520576
> prod(1:100) - n
[1] -48777321098687378615337456715518223527321845979140174232174327494146433419058837814379782860367062049372295798771978482741374619988879457910784
> prod(as.double(1:100)) - n
[1] 0

最终结果计算为零，但为 prod( as.double( 1:100 ) ) 返回的数字没有按我预期的那样显示，即使它正确计算了 prod( as.double( 1:100 ) ) - n，其中 n 是一个设置为 100! 的变量。

任何人都可以向我解释这种行为吗？据我所知，它不应该与溢出等有关，因为我使用的是 x64 系统。版本和机器信息如下:

> .Machine$double.xmax
[1] 1.798e+308
> str( R.Version() )
List of 14
 $ platform      : chr "x86_64-apple-darwin9.8.0"
 $ arch          : chr "x86_64"
 $ os            : chr "darwin9.8.0"
 $ system        : chr "x86_64, darwin9.8.0"
 $ status        : chr ""
 $ major         : chr "2"
 $ minor         : chr "15.2"
 $ year          : chr "2012"
 $ month         : chr "10"
 $ day           : chr "26"
 $ svn rev       : chr "61015"
 $ language      : chr "R"
 $ version.string: chr "R version 2.15.2 (2012-10-26)"
 $ nickname      : chr "Trick or Treat"

任何人都可以向我解释这一点吗？我不怀疑 R 做的一切都是正确的，这很可能与 useR 相关。您可能会指出，由于 prod( as.double( 1:100 ) ) - n 正确评估了我所担心的问题，但我正在执行 Project Euler Problem 20 所以我需要显示正确的数字。

谢谢

Answer 1

您使用 all.equal 进行的测试不会产生您期望的结果。 all.equal 只能比较 两个 值。第三个参数在位置上与 tolerance 匹配，它给出了比较操作的容差。在您对 all.equal 的调用中，您给了它一个 100! 的容差，这肯定会导致比较对于非常不同的值是正确的:

> all.equal( 0, 1000000000, prod(as.double(1:100)) )
[1] TRUE

但即使你只给它两个参数，例如

all.equal( prod(1:100), factorial(100) )

它仍然会产生 TRUE 因为默认容差是 .Machine$double.eps ^ 0.5 ，例如这两个操作数必须匹配到大约 8 位数字，这绝对是这种情况。另一方面，如果您将容差设置为 0 ，那么从比较中，这三种可能的组合都不相等:

> all.equal( prod(1:100), factorial(100), tolerance=0.0 )
[1] "Mean relative difference: 1.986085e-14"
> all.equal( prod(1:100), prod( as.double(1:100) ), tolerance=0.0 )
[1] "Mean relative difference: 5.22654e-16"
> all.equal( prod(as.double(1:100)), factorial(100), tolerance=0.0 )
[1] "Mean relative difference: 2.038351e-14"

另请注意，仅仅因为您告诉 R 打印 200 个有效数字并不意味着它们都是正确的。事实上，1/2^53 大约有 53 个十进制数字，但只有前 16 个被认为是有意义的。

这也使您与“真实”值的比较存在缺陷。观察这一点。 R 为您提供的 factorial(100) 的结尾数字是:

...01203456

您从中减去 n，其中 n 是 100 的“真”值!所以它最后应该有 24 个零，因此差异也应该以与 factorial(100) 相同的数字结尾。但它以:

...58520576

这仅表明所有这些数字都是不重要的，人们不应该真正研究它们的值(value)。

需要 525 位二进制精度才能准确表示 100! - 这是 double 精度的 10 倍。