通过使用 R 和 Rcpp 进行递归的阶乘函数答案

【问题标题】：Factorial function through recursion using R with Rcpp通过使用 R 和 Rcpp 进行递归的阶乘函数
【发布时间】：2014-04-17 18:52:08
【问题描述】：

我的基本问题是为什么阶乘函数的这四种实现的结果不同，更具体地说，为什么函数在 n=13 时开始不同？

library(Rcpp)
cppFunction('   int facCpp(int n)
                {
                    if (n==0) return 1;
                    if (n==1) return 1;
                    return n*facCpp(n-1);
                }
            ')



cppFunction('   double fac2Cpp(int n)
                {
                    if (n==0) return 1;
                    if (n==1) return 1;
                    return n*fac2Cpp(n-1);
                }
            ')

cppFunction('   long int fac3Cpp(long int n)
                {
                    if (n==0) return 1;
                    if (n==1) return 1;
                    return n*fac3Cpp(n-1);
                }
            ')

c(factorial(12),prod(1:12),facCpp(12),fac2Cpp(12),fac3Cpp(12))
c(factorial(13),prod(1:13),facCpp(13),fac2Cpp(13),fac3Cpp(13))
c(factorial(20),prod(1:20),facCpp(20),fac2Cpp(20),fac3Cpp(20))
c(factorial(40),prod(1:40),facCpp(40),fac2Cpp(40),fac3Cpp(40))

我意识到这个问题可能是重复的，因为这里可能会建议一个答案 Rcpp, creating a dataframe with a vector of long long 这也说明了为什么 f(13) 的函数开始不同

2^31-1>facCpp(12)
#> [1] TRUE
2^31-1>13*facCpp(12)
#> [1] FALSE


c(factorial(12),prod(1:12),facCpp(12),fac2Cpp(12),fac3Cpp(12))
#>[1] 479001600 479001600 479001600 479001600 479001600
c(factorial(13),prod(1:13),facCpp(13),fac2Cpp(13),fac3Cpp(13))
#> [1] 6227020800 6227020800 1932053504 6227020800 1932053504
c(factorial(20),prod(1:20),facCpp(20),fac2Cpp(20),fac3Cpp(20))
#> [1]  2.432902e+18  2.432902e+18 -2.102133e+09  2.432902e+18 -2.102133e+09

【问题讨论】：

12！既不会溢出int，也不会溢出long int。 13！溢出 int 但不是 long int。 40！两者都溢出。您看到的是 C++ 代码中的溢出。
@josilber thx 用于清除 C++ 相关部分的答案。但是如果你运行代码，你会注意到 fac3Cpp(13) 也返回了错误的结果。所以我的猜测是，这与不支持long int 的 Rcpp 包有关，因此 fac3Cpp 成为 facCpp。
@user2055639 fac3Cpp(13) 在我的机器上返回正确的结果。不过，C++ 标准允许这种差异。
与@Dason 类似，fac3Cpp(13) 在我的计算机上给出了正确的结果。请使用您看到的输出更新您的帖子，以获得有关您看到的差异的帮助。
@Dason ad "fac3Cpp(13) 在我的机器上返回正确的结果" 这让我很困惑？ “不过，C++ 标准允许这种差异。”不知道我明白你的目标是什么？我想最后我正在寻找一个关于何时不使用int 的通用工作规则。所以我想一般来说不能确定long int 会解决这个问题。 ...我将用输出更新我的帖子。

标签： r recursion rcpp factorial

【解决方案1】：

你基本上做错了。请参阅 R 帮助页面了解阶乘：

‘factorial(x)’（x！对于非负整数‘x’）被定义为 ‘gamma(x+1)’ 和 ‘lfactorial’ 为 ‘lgamma(x+1)’。

您不应该以这种方式计算它。为什么？看看这个：

R> evalCpp("INT_MAX")
[1] 2147483647
R>

您将遇到数字溢出。因此，例如在 R 的 factorial() 函数中实现的不同算法只执行 gamma(x+1)。你也可以在 C++ 中做到这一点：

R> cppFunction('double myFac(int x) { return(R::gammafn(x+1.0)); }')
R> myFac(4)
[1] 24
R> myFac(12)
[1] 479001600
R> myFac(13)
[1] 6227020800
R> myFac(20)
[1] 2.4329e+18
R> myFac(40)
[1] 8.15915e+47
R>

【讨论】：

@Eddelbuettel 感谢您的回答，重要的是要了解您应该采取的方式。我并没有真正瞄准阶乘的正确实现，而只是为了使用 Rcpp 得到 goint。如果您对为什么在 fac3Cpp 中使用 long int 会返回错误答案有任何想法，请不要犹豫？
没有时间详细查看，并且会怀疑那里存在一些操作系统和编译器版本依赖性，因为我们都知道正确答案应该是什么。我会设置算法，使您不会被不足或溢出咬伤。 [顺便说一句：这是有两个 D 的 Eddelbuettel ]
@user2055639 比较 evalCpp("LONG_MAX") 和 factorial(21) -- 超过 20 点仍然会溢出。double 能够存储更大范围的值。