x <- rnorm(1000000,10,3) y <- 4.1*x+2.2 + rnorm(1000000,0,.5) beta1 <- lm(y~x)$coef beta1 library(Rcpp) sourceCpp("lm.cpp") X <- as.matrix(cbind(1,x)) Y <- matrix(y) beta2 <- lmA(X,Y) beta2 sourceCpp("lmgood.cpp") beta3 <- lmB(X,Y) beta3 library(microbenchmark) microbenchmark(lm(y~x),lmA(X,Y),lmB(X,Y),time=10L)