Christian Robert, whose blog I

commented-on here once before, had
followed up on a recent set of posts by

Radford Neal which had appeared both on Radford's blog and on the
r-devel mailing list.
Now, let me prefix this by saying that I really enjoyed Radford's posts. He obviously put a lot of time into finding a number of (all somewhat
small in isolation) inefficiencies in

R which, when taken together, can make a difference in
performance. I already spotted one commit by Duncan in the SVN logs for R so this is being looked at.
Yet Christian, on the other hand, goes a little overboard in bemoaning performance differences somewhere between ten and fifteen percent -- the
difference between curly and straight braces (as noticed in Radford's first post). Maybe he spent too much time waiting for his MCMC runs to
finish to realize the obvious: compiled code is evidently much faster.
And before everybody goes and moans and groans that that is

*hard*, allow me to just interject and note that it is not. It really
doesn't have to be. Here is a quick
cleaned up version of Christian's example code, with proper assigment operators and a second variable

`x`

. We then get to the
meat and potatoes and load our

Rcpp package as well as

inline to define the same little test function in C++. Throw in

rbenchmark which I am becoming increasingly fond of for these little timing tests,

*et voila*, we have ourselves a horserace:

# Xian's code, using <- for assignments and passing x down
f <- function(n, x=1) for (i in 1:n) x=1/(1+x)
g <- function(n, x=1) for (i in 1:n) x=(1/(1+x))
h <- function(n, x=1) for (i in 1:n) x=(1+x)^(-1)
j <- function(n, x=1) for (i in 1:n) x= 1/ 1+x
k <- function(n, x=1) for (i in 1:n) x=1/ 1+x
# now load some tools
library(Rcpp)
library(inline)
# and define our version in C++
l <- cxxfunction(signature(ns="integer", xs="numeric"),
'int n = as<int>(ns); double x=as<double>(xs);
for (int i=0; i<n; i++) x=1/(1+x);
return wrap(x); ',
plugin="Rcpp")
# more tools
library(rbenchmark)
# now run the benchmark
N <- 1e6
benchmark(f(N, 1), g(N, 1), h(N, 1), j(N, 1), k(N, 1), l(N, 1),
columns=c("test", "replications", "elapsed", "relative"),
order="relative", replications=10)

And how does it do? Well, glad you asked. On my i7, which the other three cores standing around and watching, we get an
eighty-fold increase relative to the best interpreted version:

/tmp$ Rscript xian.R
Loading required package: methods
test replications elapsed relative
6 l(N, 1) 10 0.122 1.000
5 k(N, 1) 10 9.880 80.984
1 f(N, 1) 10 9.978 81.787
4 j(N, 1) 10 11.293 92.566
2 g(N, 1) 10 12.027 98.582
3 h(N, 1) 10 15.372 126.000
/tmp$

So do we really want to spend time arguing about the ten and fifteen percent differences? Moore's law gets you
those gains in a couple of weeks anyway. I'd much rather have a conversation about how we can get people speed increases that are orders of
magnitude, not fractions.

Rcpp is one such tool. Let's get more of them.