Functions in C++

2026-03-30

Outline

  • Good or bad code?
  • timing code
  • profiling tools

Connecting R and C++

  • lots of support for C++ programming from R community

  • on your machine: C++ compiler (check with which cpp in the Terminal)

A first program

library(Rcpp)

cppFunction('int add(int x, int y, int z) {
  int sum = x + y + z;
  return sum;
}')


add(1,2,3)
[1] 6

Vector input, scalar output

sumR <- function(x) {
  total <- 0
  for (i in seq_along(x)) {
    total <- total + x[i]
  }
  total
}

Vector input, scalar output in C++

cppFunction('double sumC(NumericVector x) {
  int n = x.size();
  double total = 0;
  for(int i = 0; i < n; ++i) {
    total += x[i];
  }
  return total;
}')

Package bench

x <- runif(1e3)
bench::mark(
  sum(x),
  sumC(x),
  sumR(x)
)[1:6]
# A tibble: 3 × 6
  expression      min   median `itr/sec` mem_alloc `gc/sec`
  <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
1 sum(x)       1.23µs   1.37µs   716014.        0B        0
2 sumC(x)      3.04µs   3.79µs   251875.        0B        0
3 sumR(x)     27.33µs  27.93µs    34059.      18KB        0

Let’s try this

Buffon’s needle is a Monte Carlo estimation of \(pi\) by counting the number of times that a needle of length \(\ell\) tossed on a piece of lined paper crosses a line, when the lines are \(\ell\) apart. The theoretical probability that any tossed needle crosses a line is \(2/\pi\).

Write a function rbuffon(n) that produces a binary vector with \(n\) 1s and 0s (1s for crossing needles and 0s otherwise).

How many digits of \(\pi\) can be reliably estimated with 50,000 tosses?

Double check the digits of \(\pi\) against Eve Andersson’s one million digits of pi

Let’s try this

Use Advanced R, Chapter 25 to implement a C version of rbuffon.

How do the functions compare in speed?

What about number of crossings in n throws?

Buffon’s Needle in R and C

rbuffon <- function(n) {
  # Simulate the center of the needle's distance from the nearest line
  # Uniform on [0, 0.5] since distance between lines = 1
  center_dist <- runif(n, min = 0, max = 0.5)
  
  # Simulate the acute angle the needle makes with the lines
  # Uniform on [0, pi/2]
  angle <- runif(n, min = 0, max = pi / 2)
  
  # A needle crosses a line if the horizontal projection of its half-length
  # is >= the distance from its center to the nearest line.
  # Half-length projection = 0.5 * sin(angle), needle length = 1
  crosses <- as.integer(0.5 * sin(angle) >= center_dist)
  
  return(crosses)
}

2/mean(rbuffon(50000))
[1] 3.138929
2/mean(rbuffon(50000))
[1] 3.159957
2/mean(rbuffon(50000))
[1] 3.1346

C implementation

download code file

library(Rcpp)

# Compile and load the C++ function
sourceCpp("../code/buffons-needle.cpp")

2/mean(rbuffon_cpp(50000))
[1] 3.149408
2/mean(rbuffon_cpp(50000))
[1] 3.14416
2/mean(rbuffon_cpp(50000))
[1] 3.140901

bench::mark(
  2/mean(rbuffon(100000))
)[1:6]
# A tibble: 1 × 6
  expression                  min   median `itr/sec` mem_alloc `gc/sec`
  <bch:expr>             <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
1 2/mean(rbuffon(1e+05))   3.46ms   4.24ms      226.    3.05MB     14.3
bench::mark(
  2/mean(rbuffon_cpp(100000))
)[1:6]
# A tibble: 1 × 6
  expression                      min   median `itr/sec` mem_alloc `gc/sec`
  <bch:expr>                 <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
1 2/mean(rbuffon_cpp(1e+05))   3.02ms   3.61ms      283.    1.91MB     13.9

bench::mark(
  sum(rbuffon(1000000))
)[1:6]
# A tibble: 1 × 6
  expression               min   median `itr/sec` mem_alloc `gc/sec`
  <bch:expr>          <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
1 sum(rbuffon(1e+06))   45.7ms   52.7ms      15.5    30.5MB     20.6
bench::mark(
  rbuffon_count_cpp(1000000)
)[1:6]
# A tibble: 1 × 6
  expression                    min   median `itr/sec` mem_alloc `gc/sec`
  <bch:expr>               <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
1 rbuffon_count_cpp(1e+06)     28ms   33.1ms      31.6    15.3MB     15.8

Example: Spigot algorithm for \(\pi\)

# install.packages("Rcpp")
library(Rcpp)

cppFunction('
#include <Rcpp.h>
#include <boost/multiprecision/cpp_int.hpp>
using namespace Rcpp;
using boost::multiprecision::cpp_int;

// [[Rcpp::export]]
std::string pi_spigot_boost(int digits) {
    cpp_int q = 1, r = 0, t = 1;
    cpp_int k = 1, n = 3, l = 3;

    int count = 0;
    std::string out = "";
    int newline = 0;

    while (count < digits) {
        if ((4*q + r - t) < n * t) {
            // emit digit
            out += n.convert_to<std::string>();
            count++;
            newline++;

            if (count == 1) {
              out += ".\\n";
              newline = 0;
            }
            if (newline==50) {
              out += "\\n";
              newline = 0;
            }

            cpp_int q_old = q;
            cpp_int r_old = r;

            q = 10 * q_old;
            r = 10 * (r_old - n * t);
            n = (10 * (3*q_old + r_old) / t) - 10 * n;
        } else {
            cpp_int q_old = q;
            cpp_int r_old = r;
            cpp_int t_old = t;
            cpp_int k_old = k;
            cpp_int l_old = l;

            q = q_old * k_old;
            r = (2*q_old + r_old) * l_old;
            t = t_old * l_old;
            k = k_old + 1;
            n = (q_old * (7*k_old) + 2 + r_old * l_old) / (t_old * l_old);
            l = l_old + 2;
        }
    }

    return out;
}
')

# test
cat(pi_spigot_boost(101))
3.
14159265358979323846264338327950288419716939937510
58209749445923078164062862089986280348253421170679