Chapter 3 Fibonacci Sequences Two

Example in Rcpp book.

R

mfibR <- local({
    memo <- c(1, 1, rep(NA, 1000))
    f <- function(x) {
    if (x == 0) return(0)
    if (x < 0) return(NA)
    if (x > length(memo))
    stop("x too big for implementation")
    if (!is.na(memo[x])) return(memo[x])
    ans <- f(x-2) + f(x-1)
    memo[x] <<- ans
    ans
}
})

C++

#include <Rcpp.h>
#include <algorithm>
#include <vector>
#include <stdexcept>
#include <cmath>
#include <iostream>

class Fib {
public:
    Fib(unsigned int n = 1000) {
        memo.resize(n);
        std::fill( memo.begin(), memo.end(), NAN );
        memo[0] = 0.0; 
        memo[1] = 1.0;
    }
    double fibonacci(int x) {
        if (x < 0) 
        return( (double) NAN );
        if (x >= (int) memo.size())
        throw std::range_error(\"x too large for implementation\");
        if (! std::isnan(memo[x]))
        return(memo[x]); 

        memo[x] = fibonacci(x-2) + fibonacci(x-1);
        return( memo[x] ); 
    }
private:
    std::vector< double > memo;
};

// [[Rcpp::export]]
double mfib_cpp(SEXP xs){
    int x = Rcpp::as<int>(xs);
    Fib f;
    return f.fibonacci(x);
}

Rust

pub use std::f64::*;

struct Fib {
    memo : Vec<f64>
}
impl Fib {
    fn new(cache:usize) ->Fib{
        let mut r = Vec::new();
        r.resize(cache, NAN);
        r[0] = 0.0;
        r[1] = 1.0;
        Fib{ memo : r}
    }
    fn eval(&mut self, x:usize)->RResult<f64>{
        if x < 0 { return Ok(NAN) }
        
        if x >= self.memo.len() {
            return rraise("x too large for implementation")
        }
        Ok(self.inner(x))
    }
    fn inner(&mut self, x:usize)->f64{
        if !self.memo[x].is_nan()  { return self.memo[x] }
        self.memo[x] = self.inner(x-2) as f64 + self.inner(x-1) as f64;
        self.memo[x]
    }
    
}

// #[rustr_export]
pub fn mfib_rs(x:usize)-> RResult<f64>{
    let mut fib = Fib::new(1000);
    fib.eval(x)
}

R script

library(Rcpp)
library(rustinr)
library(microbenchmark)

mfib_r <- local({
    memo <- c(1, 1, rep(NA, 1000))
    f <- function(x) {
    if (x == 0) return(0)
    if (x < 0) return(NA)
    if (x > length(memo))
    stop("x too big for implementation")
    if (!is.na(memo[x])) return(memo[x])
    ans <- f(x-2) + f(x-1)
    memo[x] <<- ans
    ans
}
})

sourceCpp(code = '
#include <Rcpp.h>
#include <algorithm>
#include <vector>
#include <stdexcept>
#include <cmath>
#include <iostream>

class Fib {
public:
    Fib(unsigned int n = 1000) {
        memo.resize(n);
        std::fill( memo.begin(), memo.end(), NAN );
        memo[0] = 0.0; 
        memo[1] = 1.0;
    }
    double fibonacci(int x) {
        if (x < 0) 
        return( (double) NAN );
        if (x >= (int) memo.size())
        throw std::range_error(\"x too large for implementation\");
        if (! std::isnan(memo[x]))
        return(memo[x]); 

        memo[x] = fibonacci(x-2) + fibonacci(x-1);
        return( memo[x] ); 
    }
private:
    std::vector< double > memo;
};

// [[Rcpp::export]]
double mfib_cpp(SEXP xs){
    int x = Rcpp::as<int>(xs);
    Fib f;
    return f.fibonacci(x);
}
')


rust(code = '
pub use std::f64::*;

struct Fib {
    memo : Vec<f64>
}
impl Fib {
    fn new(cache:usize) ->Fib{
        let mut r = Vec::new();
        r.resize(cache, NAN);
        r[0] = 0.0;
        r[1] = 1.0;
        Fib{ memo : r}
        
    }
    fn eval(&mut self, x:usize)->RResult<f64>{
        if x < 0 { return Ok(NAN) }
        
        if x >= self.memo.len() {
            return rraise("x too large for implementation")
        }
        Ok(self.inner(x))
    }
    fn inner(&mut self, x:usize)->f64{
        if !self.memo[x].is_nan()  { return self.memo[x] }
        self.memo[x] = self.inner(x-2) as f64 + self.inner(x-1) as f64;
        self.memo[x]
    }
    
}

// #[rustr_export]
pub fn mfib_rs(x:usize)-> RResult<f64>{
    let mut fib = Fib::new(1000);
    fib.eval(x)
}

')

microbenchmark(mfib_r(30L),mfib_cpp(30L),mfib_rs(30L))
#> Unit: microseconds
#>           expr  min   lq mean median   uq   max neval
#>    mfib_r(30L) 1.54 1.74 4.98   1.95 2.24 256.7   100
#>  mfib_cpp(30L) 2.40 2.99 3.87   3.22 3.59  39.0   100
#>   mfib_rs(30L) 1.54 1.83 2.82   2.06 2.68  53.2   100