CODEWITHSUNDEEP
rrecursionrcppfactorial
Jesper Hybel Pedersen
Jesper Hybel Pedersen

Reputation: 268

Factorial function through recursion using R with Rcpp

My basic question is why do the results differ for these four implementations of the factorial function and more specifically why do the functions start to differ for 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))

I realize that the question is perhaps a duplicate since an answers is probably suggested here Rcpp, creating a dataframe with a vector of long long which also shows suggests why the functions start to differ for 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

Upvotes: 1

Views: 718

Answers (1)

Dirk is no longer here
Dirk is no longer here

Reputation: 368261

You are essentially doing this wrong. See the R help page for factorial:

‘factorial(x)’ (x! for non-negative integer ‘x’) is defined to be ‘gamma(x+1)’ and ‘lfactorial’ to be ‘lgamma(x+1)’.

You are not supposed to compute it this way. Why? Well look at this:

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

You will hit numerical overflow. Hence the different algorithm as implemented eg in R's factorial() function which just does gamma(x+1). And you can do that in C++ too:

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>

Upvotes: 3

Related Questions