Speed of native R function vs C equivalent-CodePudding

When I compare the speed of the native Gamma function gamma in R to the C equivalent std::tgamma I find that the latter is about 10-15x slower. Why? I expect some differences, but this is huge, isn't it?

My implementation:

library("Rcpp")
library("microbenchmark")

cppFunction("
     double gammacpp(double x) {
         return(std::tgamma(x));
     }
")

x = max(0, rnorm(1, 50, 25))
microbenchmark(gamma(x), gammacpp(x))

Results on my machine:

Unit: nanoseconds
expr         min   lq    mean median uq  max neval cld
gamma(x)    100   101  124.91  101  102 1001   100   a
gammacpp(x) 1000 1101 1302.98 1200 1300 9401   100   b

CodePudding user response：

Look at the code of gamma() (for example by typing gamma<Return>): it just calls a primitive.

Your Rcpp function is set up to be convenient. All it took was a two-liner. But it has some overhead is saving the state of the random-number generator, setting up the try/catch block for exception handling and more. We have documented how you can skip some of these steps, but in many cases that is ill-advised.

In a nutshell, my view is that you have a poorly chosen comparison: too little happens in the functions you benchmark. Also note the unit: nanoseconds. You have a lot of measurement error here.

But the morale is that with a 'naive' (and convenient) C function as the one you wrote as a one-liner, you will not beat somewhat optimised and tuned and already compiled code from inside R. That is actually a good thing because if you did, you would have to now rewrite large chunks of R .

Edit: For kicks, here is a third variant 'in the middle' where we use Rcpp to call the same C API of R.

Code

#include <Rcpp.h>

// [[Rcpp::export]]
double gammaStd(double x) { return (std::tgamma(x)); }

// [[Rcpp::export]]
Rcpp::NumericVector gamma_R(Rcpp::NumericVector x) { return (Rcpp::gamma(x)); }

/*** R
set.seed(123)
x <- max(0, rnorm(1, 50, 25))
res <- microbenchmark::microbenchmark(R = gamma(x), Cpp = gammaStd(x), Sugar = gamma_R(x) )
res
*/

Output

> Rcpp::sourceCpp("~/git/stackoverflow/72383007/answer.cpp")
> set.seed(123)
> x <- max(0, rnorm(1, 50, 25))
> res <- microbenchmark::microbenchmark(R = gamma(x), Cpp = gammaStd(x), Sugar = gamma_R(x) )
> 
res
Unit: nanoseconds
  expr  min     lq     mean median     uq     max neval cld
     R  102  112.0   136.95    124  134.5    1068   100   a
   Cpp 1111 1155.5 11930.02   1186 1260.0 1054813   100   a
 Sugar 1142 1201.0  6355.92   1246 1301.5  506628   100   a
>