factorial_benchmarks

facbench.R

library(microbenchmark)
library(Rcpp)
library(ggplot2); theme_set(theme_bw())
library(colorspace)
library(tidyverse)

## the most naive: a for loop, in pure R
f_forloop <- function(n) {
  val = 0
  for (i in 2:n) {
     val = val + log(i)
  }
  return(val)
}

## vectorized summation in R
f_sum <- function(n) {
  sum(log(2:n))
}

## for loop in C
cppFunction('double f_forloop_cpp(int n) {
  double val = 0.0;
  for (int i = 2; i <= n; i++) {
     val += log((double) i );
  }
  return val;
}')

## sum in C (see https://dirk.eddelbuettel.com/code/rcpp/Rcpp-sugar.pdf)
cppFunction('double f_sum_cpp(int n) {
  IntegerVector vec =  seq(2, n);
  return sum(log(vec));
}')

## test
all.equal(f_forloop(40), lfactorial(40))
all.equal(f_sum(40), lfactorial(40))
all.equal(f_forloop_cpp(40), lfactorial(40))
all.equal(f_sum_cpp(40), lfactorial(40))

mm <- function(n) {
  microbenchmark(f_forloop(n),
                 f_sum(n),
                 f_forloop_cpp(n),
                 f_sum_cpp(n),
                 lfactorial(n))
}

dd <- (lst(3,10,20,40,100,500)
  %>% purrr::map_dfr(~as.data.frame(mm(.)), .id="n")
  %>% mutate(across(expr, ~forcats::fct_reorder(., time)),
             across(n, forcats::fct_inorder),
             across(time, ~ . /1000))
)

ggplot(dd, aes(x = time, y = expr, fill = n)) +
  geom_violin() +
  scale_x_log10(name="time (microseconds)") + ## , limits=c(NA, 10), oob=scales::squish) +
  scale_fill_discrete_qualitative(guide=guide_legend(reverse=TRUE)) +
  labs(y="")