在R中为什么阶乘(100)与prod(1:100)显示不同？

Question

在RI我发现了一些我无法解释的奇怪行为,我希望有人在这里.我相信100的价值!这是一个很大的数字.

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

>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

这必须不是a的最大值,double而是精度.

100!有158位有效数字(十进制).IEEE double(64位)具有52位的尾数存储空间,因此在超过大约16位十进制数字后会出现舍入错误.

顺便说一下,100!事实上,正如你所怀疑的,

93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

所以R计算的所有值都是不正确的.

现在我不知道R,但似乎在比较之前all.equal()将所有这三个值转换为floats,因此它们的差异就会丢失.

Answer 2

您的测试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它,n100的"真实"值在哪里!所以它最后应该有24个零,因此差异也应该以相同的数字结束factorial(100).但它的结尾是:

...58520576

这只表明所有这些数字都是非重要数字,而且不应该真正查看它们的值.

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

Answer 3

我将添加第三个答案,以图形方式描述您遇到的行为.从本质上讲,因子计算的双精度足够高达22 !,然后它开始越来越偏离实际值.

在50!左右,两个方法factorial(x)和prod(1:x)之间还有一个区别,后者产生,如你所示,更接近于"真实"因子.

附加代码:

# Precision of factorial calculation (very important for the Fisher's Exact Test)
library(gmp)
perfectprecision<-list()
singleprecision<-c()
doubleprecision<-c()
for (x in 1:100){
    perfectprecision[x][[1]]<-factorialZ(x)
    singleprecision<-c(singleprecision,factorial(x))
    doubleprecision<-c(doubleprecision,prod(1:x))
}


plot(0,col="white",xlim=c(1,100),ylim=c(0,log10(abs(doubleprecision[100]-singleprecision[100])+1)),
        ,ylab="Log10 Absolute Difference from Big Integer",xlab="x!")
for(x in 1:100) {
    points(x,log10(abs(perfectprecision[x][[1]]-singleprecision[x])+1),pch=16,col="blue")
    points(x,log10(abs(perfectprecision[x][[1]]-doubleprecision[x])+1),pch=20,col="red")
}
legend("topleft",col=c("blue","red"),legend=c("factorial(x)","prod(1:x)"),pch=c(16,20))