## Introduction

I’ll begin by stating that I am by no means an expert in program optimisation. If you are and you cringe at the things I write in this post then feel free to let me know. I won’t be offended. Also, if speed is absolutely of the essence then R is probably not the language for you. Rather, the purpose of this post is to show, using a case study, how one might speed up R code to get a slower task done much quicker using some of the more advanced features in R. This post was part inspired by reading Hadley Whickham’s excellent Advanced R site and partly inspired by attempting some of the challenges on the Project Euler website. Of course, any bugs in the code that follows are entirely my own doing.

The basic goal is this: write an R function that enables the user to quickly determine whether each of a potentially large number of potentially large numbers are prime or composite numbers. We’ll begin with the obvious then gradually move on to more advanced code (edge-case bugs may come and go along the way). Hopefully, however, things should remain relatively intuitive. While illustrating by example, some of the improvements will be applicable to other problems one might face using R.

## The base: divide and conquer

The most obvious method to test if a number, *N*, is a prime is to test for divisors – an approach known, for obvious reasons, as trial division. Naively one may test for any divisors between 2 and *N*-1; slightly less naively we could test numbers up to *N*/2. But for every factor a number has above its square root, it has a compatriot below the square root. For example, the factors of 24 (square root ~4.9) are 1 and 24 (obviously), 2 and 12, 3 and 8 and 4 and 6. So to show if a number, *N*, is a prime we only have to test if any integer between 2 and the root of *N* divide exactly in to it. This makes things a lot quicker. Of course if *N* is even and bigger than 2 we need not go to any effort to determine it isn’t a prime. Similarly, there’s little point looking for integers by dividing an odd number by an even one. So, if *N* is odd, we *just* need to do trial division with the other odd numbers above 1 and up to the root of *N*.

It’s simple enough to write some R code that implements this:

1 2 3 4 5 6 7 8 9 10 11 12 |
isPrime0 <- function(val){ #Special cases with sqrt < 3 if(val%in%c(2,3,5,7)){return(TRUE)} #Return FALSE immediately if val not a whole number if(val%%2==0 | val<2 | floor(val)!=val){return(FALSE)} top <- floor(sqrt(val)) return(all(val%%seq(from=3,by=2,to=top)>0)) } |

This code works but is it fast? For repeated testing purposes here are four datasets of random data:

1 2 3 4 5 |
set.seed(2014) da <- round(runif(10000,0,1000000)) db <- round(runif(1000000,0,1000000)) dc <- round(runif(10000,0,1000000000)) dd <- round(runif(1000000,0,1000000000)) |

These sets contain, in order, what could loosely be describes as a small set of small numbers, a large set of small numbers, a small set of large numbers and a large set of large numbers (though the last two sets also contain small numbers).

We can use the microbenchmark package (though it’s really meant for small chunks of code) to repeatedly test how quickly the function can check these arrays, for example:

1 2 |
library(microbenchmark) microbenchmark(sapply(da,isPrime0),times=30) |

The use of sapply (or similar) is required because the code isn’t vectorised (more on that shortly). For the da dataset the median total time (from 30 trials) to compute whether each value in the vector is a prime is just over 350 ms on my 2011 MacBook Pro. This is bearable, but for dataset db we’re looking at a time of 35-40 s while the smaller dc set takes about 3.4 s. Clearly we want to improve on this.

## Step 1: use a bit of maths and cache results

The fundamental theorem of arithmetic tells us that every integer greater than 1 is either a prime or can be represented as a product of primes. Because of this it is sufficient to establish that a number is not divisible by any prime number equal to or less than its square root to determine it is a prime.

To show how far we’ve come from our earliest considerations above, let’s consider the reasonably large number 999,983 – it’s the largest prime below 1 million. To check this by taking the most naive of approaches to trial division we’d try all 999,981 integers from 2 to 999,982. A more reasonable approach, all odd integers from three up to the (floor of the) square root of 999,982 – 999 – would reduce this to 499 trial divisions. But if we do trial division using only prime numbers below 1000 this reduces further, to just 168.

Of course this would require knowledge of primes below 1000. We can get them by using trial division. But then have we really progressed? Well, actually, yes. If we calculate and store a list of small primes once using trial division with odd numbers we can then use those primes to do a much more efficient trial division using only primes of larger numbers. It doesn’t take much effort to take the code for isPrime0 and create a different function that calculates all primes up to some specified number n. It should look something like this:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
doTD <- function(n){ #Primes with square root less than 3 hard-coded to make loop easier primes <- c(2,3,5,7) if(n<11){return(primes[primes<=n])} counter <- 5 num <- 11 while(num<=n){ top <- floor(sqrt(num)) if(all(num%%seq(from=3,to=top,by=2)>0)){ primes[counter] <- num counter <- counter + 1 } num <- num + 2 } return(primes) } |

If this is actually going to be used effectively and efficiently then we need some “safe” place to store the results of a call to this function.

### Closure time

Wikipedia has a detailed article describing exactly what a closure is in computer programming. It’s a bit verbose and unnecessarily complicated for what we need here. For our purposes, all we need to understand is that R can return a new function from inside a function and that all the variables created inside the outer function (including, potentially, more functions) can be seen by the inner, returned, function. Those additional variables can be used again and again and even changed (using the double arrow assignment operator) by the returned function.

We can utilise a closure to both store a list (well, actually a vector in R-speak) of prime numbers that our function for checking prime numbers can consult and to hold a function, such as the doTD function above, for finding new prime numbers to add to the list. This new environment is self-contained – we don’t have to worry about some external function deleting the list of primes because it won’t have direct access to it.

### Putting it together

The function below utilises a closure to store the doTD function above, a set of prime numbers and a single number that marks the maximum value, used as the argument to the doTD function, to create the vector of prime numbers. For convenience all primes below 1,000 are found as part of the function initiation.

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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 |
isPrime1 <- (function(){ doTD <- function(n){ #Primes with square root less than 3 hard-coded to make loop simpler primes <- c(2,3,5,7) num <- 11 if(n<num){return(primes[primes<=n])} counter <- 5 while(num<=n){ top <- floor(sqrt(num)) if(all(num%%seq(from=3,to=top,by=2) > 0)){ primes[counter] <- num counter <- counter + 1 } num <- num + 2 } maxPrime <<- n primes <<- primes } maxPrime <- 1000 primes <- c() doTD(maxPrime) return (function(val){ #Return FALSE immediately if val not a whole number if(floor(val) != val){return(FALSE)} #If val <= maxPrime just need to check whether val is in list of primes if(val <= maxPrime){ return(val%in%primes) } #If composite, val must have prime factor <= sqrt(val) top <- floor(sqrt(val)) #Two cases... top<=maxPrimes and top>maxPrimes #The first of these if(top <= maxPrime){ #Just check whether or not there are any prime factors of top return(all(val%%primes[primes<=top] > 0)) } #We can still check whether top has any know prime factors if(!all((val%%primes) > 0)){ return(FALSE) } #But if it doesn't then we don't know whether val is a prime or it just has a larger prime factor #We have to extend or knowledge of primes... oldLength <- length(primes) doTD(top) #Then check again #Repeating test for final number in prior vector of primes ensures no warnings if new primes is same as old primes return(all(val%%(primes[oldLength:length(primes)]) > 0)) }) #End of returned function })() |

The code in the returned function is more complex than that of isPrime0. Now than we have a list of primes up to a certain number, maxPrime, we no longer need to perform any trial division on numbers less than or equal to it. We can just return TRUE or FALSE depending on whether the value is or is not in the vector of primes. For remaining numbers whose floor of the square root is below maxPrime we can return the result of trial division with those primes. Finally, for larger numbers we can try trial division with these primes. However, a result of TRUE does not necessarily mean the number is a prime since it could have a prime factor outside the vector of primes. In that case we have to recompute our list of primes and try again. The new vector of primes is stored for future reference.

We can test this new function out with datasets da and db and it won’t ever have to calculate a new set of primes. Unfortunately it’s still slow: around 190 ms seconds for the former case and a little over 20 seconds for the latter case. So only a modest improvement. For the dc set, with larger numbers, and a freshly sourced isPrime1 function we’re looking at a time of around 4.5 s, over a second slower than the original function. However, once the list of primes has been added to in the initial call to sapply, repeating the calculations is quicker, only taking around 1.7 s.

There are two obvious ways to move forward: replace the doTD function with something better and call that replacement less often. The next section tackles the first of these. After that we’ll move on to the latter.

## Improvement 2: use the work of an ancient Greek

### Enter Erastothenes

The Sieve of Erastothenes is an algorithm for finding all primes up to a certain number by systematically marking off the products of smaller primes. In essence, start at 2 then mark – as not prime – all the products of 2. Then move on to the next unmarked number – 3 – and repeat. Then the next – 5 – and so on. The numbers unmarked after this process has been completed are the primes.

Here’s some R code that implements this:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
doSieve <- function(n){ #Mark all odd numbers as prime and even numbers as not prime isAPrime <- c(rep(c(TRUE,FALSE), length=n)) #Swap special cases of 1 and 2 isAPrime[c(1,2)] <- c(FALSE,TRUE) p <- 3 #Now cross-off all the odd numbers that are products of known primes while((pp<-p^2)<=n){ if(isAPrime[p]){ #First value that needs crossing off is p^2, not 2*p because... #earier values will already have been crossed off... #eg for 5: 10 and 20 are products of 2 and 15 is product of 3 #Step through by 2*p since p^2 + odd p will be even isAPrime[seq(from=pp,to=n,by=(2*p))] <- FALSE } p <- p + 2 } return((1:n)[isAPrime]) } |

Using doSieve to find prime numbers up to a specified number is much faster as this chart shows (note log-log scale):

We can swap the doTD function for the doSieve function in our prime-checking function (with a couple of minor changes):

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 41 42 43 44 45 46 47 48 49 50 51 |
isPrime2 <- (function(){ doSieve <- function(n){ isAPrime <- c(rep(c(TRUE,FALSE), length=n)) isAPrime[c(1,2)] <- c(FALSE,TRUE) p <- 3 while((pp<-p^2)<=n){ if(isAPrime[p]){ isAPrime[seq(from=pp,to=n,by=(2*p))] <- FALSE } p <- p + 2 } maxPrime <<- n primes <<- ((1:n)[isAPrime]) } maxPrime <- 1000 primes <- c() doSieve(maxPrime) return (function(val){ if(floor(val) != val){return(FALSE)} if(val <= maxPrime){ return(val%in%primes) } top <- floor(sqrt(val)) if(top <= maxPrime){ return(all(val%%primes[primes<=top] > 0)) } if(!all((val%%primes) > 0)){ return(FALSE) } oldLength <- length(primes) doSieve(top) return(all(val%%(primes[oldLength:length(primes)]) > 0)) }) #End of returned function })() |

Testing the freshly sourced isPrime2 code with the dc dataset drops the calculation time to about 1.9 s.

If you never need to check large numbers then isPrime1 and isPrime2 should be indistinguishable in speed (since this inner functions are essentially identical). However, if the vector of primes needs frequent updates then clearly isPrime2 has the edge.

## Step 3: vectorisation

At the minute, if you want to determine if many numbers are prime you’d have to use isPrime2 inside a call to sapply (or a for loop). This is unnecessarily slow because:

- Loop and function calls are slow.
- The function has no way of knowing if a bigger value will shortly be checked. This leads to unnecessary calls to doSieve and is especially problematic if the numbers you want to check are large and in ascending order.

Let’s back statement 2 up with some data (again, your results will vary). From above we know the median time to find all primes in the dc data set is about 1.9 s when freshly sourced. If we source the function again and sort the data into descending order this drops to about 1.8 s (little different from when the function is executed a second time with the same data). But if we sort it into ascending order instead then running the code takes about 2.4 s as we have to call the doSieve function repeatedly.

We can do better by vectorising the code so that multiple values can be passed in at once. This is also as good a time as any to include an explicit isWhole function to test whether a number is or is not whole to within an acceptable tolerance. (As is often the case, this came from a stackoverflow question which is well worth a read if you’re confused about why this might be deemed necessary.) My initial solution looked like this:

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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 |
isPrime3 <- (function(){ doSieve <- function(n){ isAPrime <- c(rep(c(TRUE,FALSE), length=n)) isAPrime[c(1,2)] <- c(FALSE,TRUE) p <- 3 while((pp<-p^2)<=n){ if(isAPrime[p]){ isAPrime[seq(from=pp,to=n,by=(2*p))] <- FALSE } p <- p + 2 } maxPrime <<- n primes <<- ((1:n)[isAPrime]) } maxPrime <- 1000 primes <- c() doSieve(maxPrime) isWhole <- function(x, tol = .Machine$double.eps^0.5){return(abs(x - round(x)) < tol)} return (function(vals){ vals <- ifelse(isWhole(vals),vals,1) #Even values > 2 and all values >2 must be FALSE out <- ifelse((vals%%2==0 & vals>2)|vals<2,FALSE,NA) indices <- which(is.na(out)) nRemain <- length(indices) remain <- vals[indices] #Check whether remaining values are found in set of stored primes out[indices] <- ifelse(remain%in%primes,TRUE,ifelse(remain%%2==0 | remain<=maxPrime, FALSE,NA)) #Do we still have NAs? if(any(is.na(out))){ indices <- which(is.na(out)) nRemain <- length(indices) remain <- vals[indices] tops <- floor(sqrt(remain)) results <- rep(NA,nRemain) #Loop and do trial division by known primes for(i in 1:nRemain){ top <- tops[i] val <- remain[i] test <- all(val%%primes[primes<=top] > 0) results[i] <- ifelse(test,ifelse(top<=maxPrime,TRUE,NA),FALSE) } out[indices] <- results #Odd numbers that might have primes > maxPrimes will still be NA if(any(is.na(out))){ indices <- which(is.na(out)) nRemain <- length(indices) remain <- vals[indices] tops <- floor(sqrt(remain)) #Our new maximum: theMax <- max(tops) oldLength <- length(primes) doSieve(theMax) newPrimes <- primes[oldLength:length(primes)] results <- rep(NA,nRemain) #Loop and do trial division by new primes for(i in 1:nRemain){ val <- remain[i] top <- tops[i] primeSet <- newPrimes[which(newPrimes<=top)] results[i] <- all(val%%primeSet > 0) } out[indices] <- results } } return(out) }) #End of returned function })() |

We can now pass a single vector with – say – a million numbers to check for primality and the doSieve function will be called at most once. (Of course if you pass a vector with even bigger numbers in a subsequent call you’ll have to re-sieve.)

To my disappointment this code was still pretty slow. This isn’t helped by the presence of for loops. I was unable to vectorise these out of the code (I’m not saying that’s not possible, I just couldn’t work out how to do it).

And another result I didn’t expect: when freshly sourced the dc set would be computed in about 200 ms, but on repeated calculations this went up to about a second! This must be because the function is carrying around and subsetting a larger set of primes after the first use. For most tested numbers that isn’t really necessary.

We could of course not save the larger set of primes. But that is making a big assumption about how the function will be utilised. If, for example, it’s called frequently with small vectors of large numbers then things could get very slow. Before I try to tackle that problem, however, we’ll dispose of the for loops – or rather delegate them to another programming language altogether.

## Step 4: “cheat”

As just mentioned, I was unable to get rid of the slow for loops using R code. So at this point, you may say, I cheated.

### Rcpp

The Rcpp package lets you write and compile functions in C++ and then use them as if they were R functions in R. Using loops in R is notoriously slow. This isn’t the case in C++.

I first learnt to program in C++, though things may have got a bit rusty in the last year or two. Thankfully, Rcpp makes the whole experience quite pleasant. I’ll defer to Hadley Wickham for a more detailed explanation, but a few little points for those unfamiliar:

- Each variable needs to have an explicitly declared type (int, double, bool etc plus some really useful types defined by the Rcpp package that mimic R types);
- Every statement ends with a semi-colon (90% of the time this is why my C++ code fails to compile);
- Arrays start at 0;
- You must explicitly return a value at the end of a function. (I nearly always do this in R anyway since I think it makes the code clearer.)

After that, most of the other bits I needed I got by modifying examples from Whickham’s web page and answers on stackoverflow.

Anyway, this C++ function does the work of the loops in isPrime4 if you pass it your values and your list of prime numbers:

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 |
#include <Rcpp.h> using namespace Rcpp; // [[Rcpp::export]] LogicalVector checkPrimality(IntegerVector vals, IntegerVector primes) { const int32_t nVals = vals.size(); const int32_t nPrimes = primes.size(); LogicalVector out(nVals); for(int32_t i=0; i<nVals; i++){ int32_t val = vals[i]; int32_t top = sqrt(val); bool isAPrime = true; if(val < 2){ //If val < 2 no need to loop isAPrime = false; } else{ for(int32_t j=0; j<nPrimes; j++){ int32_t prime = primes[j]; if(prime > top){ //We have a prime break; } if(val%prime == 0){ //We have a composite isAPrime = false; break; } } } out[i] = isAPrime; } return(out); } |

Implementing it in a new addition of isPrime looks like this:

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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 |
isPrime4 <- (function(){ doSieve <- function(n){ isAPrime <- c(rep(c(TRUE,FALSE), length=n)) isAPrime[c(1,2)] <- c(FALSE,TRUE) p <- 3 while((pp<-p^2)<=n){ if(isAPrime[p]){ isAPrime[seq(from=pp,to=n,by=(2*p))] <- FALSE } p <- p + 2 } maxPrime <<- n primes <<- ((1:n)[isAPrime]) } maxPrime <- 1000 primes <- c() doSieve(maxPrime) require("Rcpp") sourceCpp("checkPrimality.cpp",env=environment()) isWhole <- function(x, tol = .Machine$double.eps^0.5){return(abs(x - round(x)) < tol)} return (function(vals){ vals <- ifelse(isWhole(vals),vals,1) out <- ifelse((vals%%2==0 & vals>2)|vals<2,FALSE,NA) indices <- which(is.na(out)) nRemain <- length(indices) remain <- vals[indices] out[indices] <- ifelse(remain%in%primes,TRUE,ifelse(remain%%2==0 | remain<=maxPrime, FALSE,NA)) if(any(is.na(out))){ indices <- which(is.na(out)) remain <- vals[indices] tops <- floor(sqrt(remain)) out[indices] <- ifelse(checkPrimality(remain,primes),NA,FALSE) if(any(is.na(out))){ indices <- which(is.na(out)) remain <- vals[indices] tops <- floor(sqrt(remain)) theMax <- max(floor(sqrt(max(remain)))) oldLength <- length(primes) doSieve(theMax) newPrimes <- primes[oldLength:length(primes)] out[indices] <- checkPrimality(remain,primes) } } return(out) }) #End of returned function })() |

Now the factory-fresh isPrime4 will calculate the results for the dc dataset in about 22 ms, while a repeated calculation will take about 26 ms. For the bigger dd data set we’re looking at around 2.2 s and 2.7 s, respectively (nicely in-line with the factor of 100 increase in vector size). Generally, things are much improved: we’ve roughly gained an order of magnitude in speed from isPrime3.

Profiling the execution of isPrime4 with the profr package and the dd dataset showed that more than 50% of the execution time was spent outside the checkPrimality function. The looping is no longer dominating execution time now that we’ve introduced Rcpp.

## Step 5: cheat some more

Using C++ brought back to mind one thing that I never really worry about using R: numbers need to be stored in some format or other. Most pertinently, integer types have a limit to the maximum number they can store. R actually has a specific integer type, a 32-bit signed integer, whose maximum value can be found via .Machine$integer.max. That’s 2,147,483,647. That’s also the limit to the size of the integer in the Rcpp defined IntegerVector. So if we test isPrime4 with an odd number larger than that (which R stores as a floating-point number) we get a warning and an NA value.

We can swap the “IntegerVector vals” entry in the argument list of checkPrimality with “std::vector<int64_t> vals” and then swap the “int32_t” type declarations for the val variable with an “int64_t” type declaration. Much to my surprise this actually seems to work, though I haven’t rigorously tested it. It is, however, 50% slower for the dc and dd datasets.

If we don’t care about testing really big numbers there’s an alternative that makes the code much simpler and much faster. The floor of the square root of .Machine$integer.max is “just” 46,340 and doSieve can calculate all primes up to that on my computer in just ~5 ms. It also creates a vector of just 4,792 prime numbers – not problematic for storing in memory.

This means that we can sieve just once and be guaranteed to have all the primes we need to test all numbers below our 32-bit integer ceiling. And to make things really quick we can do pre-flight checks using C++ too:

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 |
#include <Rcpp.h> using namespace Rcpp; // [[Rcpp::export]] IntegerVector convertFloats(NumericVector vals, double tol, int32_t lim) { const int32_t nVals = vals.size(); IntegerVector out(nVals); for(int32_t i=0; i<nVals; i++){ double val = vals[i]; int32_t rVal = round(val); //Check val. In checkPrimality function a no. less that 2 quickly returns false if(fabs(val) > lim){ //Out of range of int32_t if(val>lim){ //Too big: error stop("One or more numbers is out of range"); } else{ //Too small: but can but set to 0 to get false out[i] = 0; } } else if((val == rVal) || fabs(val - rVal) < tol){ //Just right out[i] = rVal; } else{ //Not integer or within integer tolerance so set to 0 out[i] = 0; } } return(out); } |

This is then used in the final (for now at least) isPrime function, isPrime5. For this I’ve added another couple of arguments. The first additional argument checks what the user wants to do with NA values (preserve them in the final output or return FALSE in their place). The second argument allows the user to specify their own tolerance level when checking whether a number is an integer. The end result is this:

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 41 42 43 44 45 46 47 48 49 50 51 52 |
isPrime5 <- (function(){ lim <- .Machine$integer.max primes <- (function(n = floor(sqrt(lim))){ isAPrime <- c(rep(c(TRUE,FALSE), length=n)) isAPrime[c(1,2)] <- c(FALSE,TRUE) p <- 3 while((pp<-p^2) <= n){ if(isAPrime[p]){ isAPrime[seq(from=pp,to=n,by=(2*p))] <- FALSE } p <- p + 2 } return((1:n)[isAPrime]) })() require("Rcpp") sourceCpp("convertFloats.cpp",env=environment()) sourceCpp("checkPrimality.cpp",env=environment()) return (function(vals, preserveNAs=TRUE, tol=.Machine$double.eps^0.5){ indicesNA <- c() if(preserveNAs & any(is.na(vals))){ #Ugly hack to make sure NA values are kept rather than returning FALSE indicesNA <- which(is.na(vals)) #store indices of NA values } if(!is.integer(vals)){ #If we have integers already then we're good to go if(is.numeric(vals)){ #If we have doubles though do custom conversion vals <- convertFloats(vals, tol, lim) } else{ #Otherwise we have something we can't use: throw an error stop("vals must be numeric") } } out <- checkPrimality(vals, primes) if(length(indicesNA) > 0){ #Second part of ugly hack out[indicesNA] <- NA #add NA values back in } return(out) }) #End of returned function })() |

Testing with dataset db produces a median calculation time of 105 milliseconds. Compared with the original isPrime0 function we’ve gained a factor of about ~350. Testing with the dd dataset gives a median time of ~0.8 s, roughly a factor of three improvement on isPrime4. As importantly, the code is somewhat simpler (I think) and there are added checks and options.

At this point I ran out of ideas for making any more noticeable changes, so now we’ll move on to comparing with pre-existing functions from several R packages I found.

## Comparisons with existing packages

### Sieving

The sfsmisc package provides a primes function and the numbers package provides a similarly named Primes function. Both essentially implement the Sieve of Erastothenes. These and doSieve are comparable in speed (see the chart below) but (to my slight disappointment) primes is marginally quicker. Because I don’t like losing I created doSieve2: like doSieve but with the while loop written as a C++ function. For larger input arguments, doSieve2 is about a factor of 2.5 quicker than the other functions:

### Prime-checking

I found a couple of R packages with functions that allow for prime-number checking in a similar manner to isPrime5. The numbers package (see above) provides the isPrime function while the gmp package provides the isprime function. The former calls the Primes function just discussed and then uses trial division with the resulting set of primes. It’s about a factor of 2 slower than isPrime3 (but the code is somewhat neater).

The isprime function is a more complex prime-number checker. I tried to follow through the source code to understand exactly how it worked but got lost in a maze of C. The function help file, however, tells us that trial division is tried first and then followed by some number of Miller-Rabin probabilistic tests. Because the tests are only probabilistic it’s not always certain whether a candidate is or isn’t prime. Hence the function returns 0 for not prime, 1 for probably prime and 2 for definitely prime.

Rather than compare isPrime5 and is prime using the usual data sets I instead created six new ones. Each new data set contained 1 million integers between 0 and some upper limit. The values for the upper limit chosen went in powers of 10 from 100,000 (with repetitions of course) to 1 billion and then a final set that ranged up to .Machine$integer.max (ie around 2 billion). Using these data I compared isPrime5 with isprime (using three different values for the number of Miller-Rabin tests to perform). The results are plotted below (note the vertical axis is *not* logarithmic here):

Below ~1 billion, isPrime5 is faster than isprime for any number of Miller-Rabin tests (including 0), this is generally the case with smaller vectors of numbers too. isPrime5 also gives definite answers. Having said that, with both 5 and 40 Miller-Rabin tests, all primes identified as probably prime were actually prime. The definite advantage that isprime has over isPrime is that it can deal with 64-bit integers.

## Summary

Admittedly we were starting from a fairly low base to begin with, but by using some “advanced” features of R (including the ability to use another programming language entirely) and a bit of maths, we’ve managed to speed up a (I guess) fairly common task by about two orders of magnitude. The main resulting function is now, for vectors of numbers below ~ a billion, faster than the equivalent functions I found in existent packages. For larger numbers you may be better off with the isprime function in the gmp package.

There are still some open questions: Have I missed any other notable functions in existent packages? Are there any edge-cases not properly dealt with? Are there any obvious (or obscure) means to significantly speed up the code while keeping the mathematics relatively intuitive? Can isPrime5 be extended in some way to cope with 64-bit integers without suffering a huge speed penalty?