23

在 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 ) ) - nwhere nis a variable set to value 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 做的一切都是正确的，这很可能与 userR 相关。您可能会指出，因为prod( as.double( 1:100 ) ) - n正确评估了我在乎的东西，但我正在做Project Euler Problem 20，所以我需要显示正确的数字。

谢谢

4

4 回答 4

15

这与 a 的最大值无关，double而与它的精度有关。

100!有 158 个有效（十进制）数字。IEEE doubles（64 位）有 52 位的尾数存储空间，因此在超过大约 16 位十进制精度后会出现舍入错误。

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

93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

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

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

于 2013-01-14T11:00:40.910 回答
13

您的测试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

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

为了准确地表示 100，它需要 525 位的二进制精度！- 这是 10 倍的精度double。

于 2013-01-14T13:24:29.603 回答
7

我将添加第三个答案，以图形方式描述您遇到的行为。本质上，阶乘计算的双精度足以达到 22!，然后它开始越来越偏离实际值。

在 50! 左右，factorial(x) 和 prod(1:x) 这两种方法之间存在进一步的区别，正如您所指出的，后者产生的值更类似于“真实”因子。

R中的阶乘计算精度

附上代码：

# 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))
于 2013-07-11T23:29:46.783 回答
1

好吧，您可以从它调用的主体中factorial看出它gamma调用了.Primitive("gamma")。.Primitive("gamma")看起来像什么？像这样。

对于大输入，.Primitive("gamma")的行为在该代码的第 198 行。它在呼唤

exp((y - 0.5) * log(y) - y + M_LN_SQRT_2PI +
            ((2*y == (int)2*y)? stirlerr(y) : lgammacor(y)));

这只是一个近似值。

顺便说一下，文章Rmpfr以factorial它为例。因此，如果您想解决问题，“只需使用Rmpfr库”。

于 2015-09-03T17:06:04.317 回答