在 R 中，为啥 factorial(100) 与 prod(1:100) 的显示不同？

Posted 2023-02-14

【中文标题】在 R 中，为啥 factorial(100) 与 prod(1:100) 的显示不同？

【英文标题】：In R why is factorial(100) displayed differently to prod(1:100)?在 R 中，为什么 factorial(100) 与 prod(1:100) 的显示不同？

【发布时间】：2012-12-28 07:56:23

【问题描述】：

在 R 中，我发现了一些我无法解释的奇怪行为，我希望这里有人可以。我相信100的价值！这是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 做的一切都是正确的，这很可能与 userR 相关。您可能会指出，由于 prod( as.double( 1:100 ) ) - n 正确评估了我在乎的东西，但我正在做 Project Euler Problem 20 所以我需要显示正确的数字。

谢谢

【问题讨论】：

要使用 R 计算 100! 的精确值，请执行：library(gmp); factorialZ(100)

@JoshO'Brien 非常感谢！

感谢所有关于这个问题的海报。我想我现在对 R 中的大整数有了更好的理解。使用library(gmp)，正如一些人所建议的那样，我注意到我可以做identical(factorialZ(100) , prod(as.bigz(1:100)))，这将返回[1]TRUE。

【参考方案1】：

嗯，你可以从factorial 的正文中看出它调用了gamma，它调用了.Primitive("gamma")。 .Primitive("gamma") 是什么样的？ Like this.

对于大型输入，.Primitive("gamma") 的行为在该代码的 line 198 上。它在召唤

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

这是just an approximation。

---

顺便说一句，the article on Rmpfr 以factorial 为例。因此，如果您想解决问题，“只需使用 Rmpfr 库”。

【参考方案2】：

我将添加第三个答案，以图形方式描述您遇到的行为。本质上，阶乘计算的双精度足以达到 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))

【参考方案3】：

您使用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

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

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

【讨论】：

非常感谢以 R 为中心的解释。我会将这个标记为正确答案，而不是蒂姆的（对不起！），因为我觉得它更好地回答了我原来的问题，给出了 [r] 标签。谢谢。

【参考方案4】：

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

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

顺便说一句，100! 实际上正如您所怀疑的那样，

93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

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

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

【讨论】：

请仔细检查我上面n的值。它与您发布的内容完全相同。在我将此问题标记为已回答之前，我将等待一天左右，看看是否有任何其他意见。谢谢。

@SimonO101：你是对的；在进行比较时，我忽略了这一点。我已经相应地编辑了我的答案。

谢谢蒂姆。如果 R 首先没有正确计算数字，您能否解释为什么 R 正确评估 prod( as.double( 1:100 ) - 100!（即 =0）？此外，如果 158 位整数在第 14 位之后开始出现差异，这似乎是大量的舍入错误？

一个可能的解决方案是 gmp 包，它为 R 中的 GNU 多精度库提供接口。

@TimPietzcker，这是正确的 - R 确实 not 像 Python 那样具有用于大整数的内置类型。它的目标是统计数据处理，而不是解决数论中的挑战。