# R/C functions in R code

< R

This is a quick way to write a loop as a C function in R. This description assumes you are using Windows. A few small changes are required for other OSs.

First, install Rtools

The following R code illustrates a C function, its compilation and its inclusion in an R function.

'benchmark' is used to compare the speed of of the raw R function, the compiled R function, the C function using a 'double' accumulator, and the C function using a 'long double' accumulator.

The speed advantage of using the C function compared with the raw R function is in the order of 100. The compiled R function runs roughly twice as fast as the raw function.

```install.packages('rbenchmark')
library(rbenchmark)
library(compiler)

# Speeding up a function with a loop

Rfun <- function(t, N = 100000) {
tmp1 <- 0
tmp2 <- 0
tmp3 <- 0
for (j in N:1) {
tmp1 <- tmp1 + 0.5*log(j*(j+1)/(j*(j+1)-2*t))
tmp2 <- tmp2 + 1/(j*(j+1)-2*t)
tmp3 <- tmp3 + 2/((j*(j+1)-2*t))^2
}
c(tmp1, tmp2, tmp3)
}

# Using the R compiler

Rfun.comp <- cmpfun(Rfun)

# The function in C

cat('
#include <R.h>
#include <math.h>
void fun(double *t, double *N, double *ret)
{
double tmp1, tmp2, tmp3, j;
tmp1 = tmp2 = tmp3 = 0;
for (j = *N; j > 0.9 ; j--) {
tmp1 = tmp1 + 0.5*log(j*(j+1)/(j*(j+1)-2.0*t[0]));
tmp2 = tmp2 + 1.0/(j*(j+1)-2.0*t[0]);
tmp3 = tmp3 + 2.0/pow(((j*(j+1)-2.0*t[0])),2.0);
};
ret[0] = tmp1;
ret[1] = tmp2;
ret[2] = tmp3;
ret[3] = j;
}
', file = 'fun.c')
system( 'rm fun.dll')
system( "R CMD SHLIB fun.c" )

Cfun <- function(t, N = 100000){
ret <- double(4)
ret2 <- .C('fun', t = as.double(t), N = as.double(N),
ret = ret)
ret2\$ret
}

# Using long double as an accumulator to see what effect that has

cat('
#include <R.h>
#include <math.h>
void funl(double *t, double *N, double *ret)
{
long double tmp1, tmp2, tmp3, j, Nl;
tmp1 = tmp2 = tmp3 = 0;
Nl = *N;
for (j = Nl; j > 0.9 ; j--) {
tmp1 = tmp1 + 0.5*log(j*(j+1)/(j*(j+1)-2.0*t[0]));
tmp2 = tmp2 + 1.0/(j*(j+1)-2.0*t[0]);
tmp3 = tmp3 + 2.0/pow(((j*(j+1)-2.0*t[0])),2.0);
};
ret[0] = tmp1;
ret[1] = tmp2;
ret[2] = tmp3;
ret[3] = j;
}
', file = 'funl.c')
system( 'rm funl.dll')
system( "R CMD SHLIB funl.c" )

Cfunl <- function(t, N = 100000){
ret <- double(4)
ret2 <- .C('funl', t = as.double(t), N = as.double(N),
ret = ret)
ret2\$ret
}

# -475 is a plausible argument for this function
Rfun( -475) - Rfun.comp( -475)
Cfun( -475) - Cfunl( -475)
Rfun( -475) - Cfunl(-475)[1:3]

benchmark( Rfun( -475), Rfun.comp(-475),Cfun(-475),Cfunl(-475), replications=10)```