Dirk Eddelbuettel — Dec 20, 2012 | source

This example takes the standard simulation for the value of pi and in implements it in both R and C++. We will see that the C++ version is very similar to the R version.

We start with an R version.

The basic idea is that for a point (x,y), we compute the distance to the origin using Pythargoras’ well-known expression, or in this context a standard distance norm. We do this repeatedly for a number of points, the ratio of those below one (“inside the unit circle”) to the number N of simulation will approach pi/4 – as we were filling the area of one quarter of the unit circle by limiting ourselves the first quadrant (as we forced x and y to be between 0 and 1).

The key here is that we do this in a *vectorised* way: by drawing
all N x and y coordinates, then computing all distances and comparing to 1.0. After
scaling up to the full circle, we have an estimate of pi.

```
piR <- function(N) {
x <- runif(N)
y <- runif(N)
d <- sqrt(x^2 + y^2)
return(4 * sum(d < 1.0) / N)
}
set.seed(5)
c(piR(1000), piR(10000), piR(100000), piR(1000000))
```

[1] 3.156 3.155 3.139 3.141

The neat thing about Rcpp sugar enables us to write C++ code that looks almost as compact.

```
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
double piSugar(const int N) {
RNGScope scope; // ensure RNG gets set/reset
NumericVector x = runif(N);
NumericVector y = runif(N);
NumericVector d = sqrt(x*x + y*y);
return 4.0 * sum(d < 1.0) / N;
}
```

Apart from using types (hey, this is C++) and assuring the RNG gets set and reset, the code is essentially identical.

And by using the same RNG, so are the results.

```
set.seed(5)
c(piSugar(1000), piSugar(10000), piSugar(100000), piSugar(1000000))
```

[1] 3.156 3.155 3.139 3.141

**tags:**
sugar

- Call Python from R through Rcpp — Wush Wu
- Using iterators for sparse vectors and matrices — Soren Hojsgaard and Doug Bates
- A simple array class with specialized linear algebra routines — Fabian Scheipl