[Manual] Rcpp使用手册 – 南有嘉鱼

如果你是第一次阅读这篇文章,以下资源贴可能是你感兴趣的:

Rtools下载安装
Rcpp官网
Dirk Eddelbuettel大神的神坛
Rcpp-introduction
高级R语言编程指南Rcpp章节
Rcpp Quick Reference Guide


示例一 斐波那契数列

R解决方案之一

1
2
3
4
5
fibR=function(n) {
if(n==0) return(0)
if(n==1) return(1)
return(fibR(n-1)+fibR(n-2))
}

C++解决方案之一

装好Rtools最好重启一下RStudio,运行以下例子,如果编译成功,没有报错那就可以放心使用了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
library(inline)
library(Rcpp)
incltxt='
int fibonacci(const int x) {
if(x==0) return(0);
if(x==1) return(1);
return fibonacci(x-1)+fibonacci(x-2);
}
'
fibRcpp=cxxfunction(signature(xs='int'),
plugin = "Rcpp",
includes = incltxt,
body = '
int x=Rcpp::as<int>(xs);
return Rcpp::wrap(fibonacci(x));')

signature变量确认函数签名,定义变量进入body的接口;includes允许传递额外的指令,或者类似以下例子中的函数和类定义

R解决方案之二

前面的递归计算中,存在重复的计算,可以通过“存储”的方式更高效地计算。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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("你输入的值太大啦")
if(!is.na(memo[x])) return(memo[x])
ans<-f(x-2)+f(x-1)
memo[x]<<-ans
ans
}
}
)
mfibR(900)

可以看到这个办法中,会在递归之前查找是否已经计算过某一个值

C++解决方案之二

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
mincltxt='
#include <algorithm>
#include <vector>
#include <stdexcept>
#include <cmath>
#include <iostream>
class Fib {
public:
Fib(unsigned int n=1000) {
memo.resize(n); //预留n个元素
std::fill(memo.begin(),memo.end(),NAN); //全部设为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; //存储之前的计算结果
};
'
mfibRcpp=cxxfunction(signature(xs='int'),
plugin = "Rcpp",
includes = mincltxt,
body='
int x= Rcpp::as<int>(xs);
Fib f;
return Rcpp::wrap( f.fibonacci(x-1) );
')
mfibRcpp(600)

R解决方案之三

1
2
3
4
5
6
7
8
9
10
11
12
fibRiter=function(n) {
first=0
second=1
third=0
for (i in seq_len(n)) {
third=first+second
first=second
second=third
}
return(first)
}
fibRiter(1000)

由于迭代算法既不需要存储状态,也不需要递归,所以其可以在存储方法的基础上进一步提高。

C++解决方案之三

1
2
3
4
5
6
7
8
9
10
11
12
13
14
fibRcppIter=cxxfunction(signature(xs="int"),
plugin = 'Rcpp',
body='
int n=Rcpp::as<int>(xs);
double first=0;
double second=1;
double third=0;
for(int i=0;i<n;i++) {
third=first+second;
first=second;
second=third;
}
return Rcpp::wrap(first);
')

这是目前最快的版本。

数据结构

IntegerVector 类

返回完美数

完美数是指等于其因子之和的正整数,本例的函数签名为空,选择“Rcpp”作为plugin会引导cxxfunction()寻找Rcpp中的头文件和库:

1
2
3
4
5
6
7
8
9
10
src='
Rcpp::IntegerVector epn(4);
epn[0]=6;
epn[1]=14;
epn[2]=496;
epn[3]=8182;
return epn; //这里没有显式声明从C++向量转换为R向量
'
fun=cxxfunction(signature(),src,plugin = 'Rcpp')
fun()
  • 由于隐式调用wrap,返回向量不需要额外的代码,直接使用return epn;
  • 创建一个新向量和定义其大小可以调用Rcpp::

使用R向量输入

1
2
3
4
5
6
7
8
9
10
src='
Rcpp::IntegerVector vec(vx);
int prod=1;
for (int i=0;i<vec.size();i++) {
prod*=vec[i];
}
return Rcpp::wrap(prod); //这里的返回值用warp包裹起来了
'
fun=cxxfunction(signature(vx='integer'),src,plugin = 'Rcpp')
fun(1:10)

NumericVector类

使用两个输入

1
2
3
4
5
6
7
8
9
10
11
12
src='
Rcpp::NumericVector vec(vx);
double p=Rcpp::as<double>(dd);
double sum=0.0;
for(int i=0;i<vec.size();i++) {
sum+=pow(vec[i],p);
}
return Rcpp::wrap(sum);
'
fun=cxxfunction(signature(vx='numeric',dd='numeric'),src,plugin = 'Rcpp')
fun(1:4,2.2)