Problem: What is the 10001st prime number?

Commentary: So it’s finally time to discuss my prime-searching algorithm, as promised back in Problem 5. We will find primes in the following way:

1. 2 and 3 are prime
2. We examine the odd integers (having already discovered the only even prime) starting with $i=5$.
3. First, check if $i$ is of the form $6k+3$ for some integer $k>0$. If so, then $i$ can’t be prime, because then $i=6k+3=3(2k+1)$, and since $k>0$, $2k+1 > 1$ (so $i eq 3$). In this case, don’t even bother proceeding; replace $i$ with $i+2$ and start over.
4. Calculate $s(i)=\sqrt{i}$. Among all known primes below $s(i)$ (going no further for the same reason as that outlined in Problem 3), check if these primes divide into $i$. If we find one that does, then there is no reason to proceed; $i$ is not prime. Go back to 3. and try again with $i$ replced with $i+2$. If none of the known primes below $s(i)$ divide $i$, then stop; $i$ must be prime. Add it to the list of known primes.
5. Repeat this process by replacing $i$ with $i+2$ until we have 10001 primes.

It’s nothing fancy, but it works. As for the implementation, believe it or not, I actually don’t construct a list of length 10001 and insert into it as I gather my primes like a good little R user. See the post-code commentary for an explanation.

R Code:

ElapsedTime <- system.time({
##########################
primes <- c(2, 3)
i <- 3
while (length(primes) < 10001){
flag <- 0
i <- i+2
if (i%%6 == 3){
flag <- 1
}
if (flag == 0){
s <- sqrt(i)+1
for (prime in primes){
if ((i%%prime == 0)){
flag <- 1
break
}
if (prime > s){
break
}
}
if (flag == 0){
primes <- c(primes, i)
}
}
}

##########################
})[3]
ElapsedMins <- floor(ElapsedTime/60)
ElapsedSecs <- (ElapsedTime-ElapsedMins*60)
cat(sprintf("\nThe answer is:  %d\nTotal elapsed time:  %d minutes and %f seconds\n",


Output: The answer is: 104743 Total elapsed time: 0 minutes and 1.191000 seconds

Post-Code Commentary: You may notice that I’m not pre-allocating the prime list and then storing new prime values as I go, like a good little R user. The reason is that it’s actually slower to do it this way.** Now maybe I’m an idiot and I don’t really understand how R does things, but if I take this algorithm and only modify the prime storage part of it to include preallocating the list and then injecting into that list (replacing a 0) as I go, then the time is nearly a full order of magnitude slower. Here’s the R code that you would produce if you were to do what I just described:

primes <- numeric(10001)
primes[1:2] <- c(2, 3)
i <- 3
prime.position <- 3

while (primes[10001] == 0){
flag <- 0
i <- i+2
if (i%%6 == 3){
flag <- 1
}
if (flag == 0){
s <- sqrt(i)+1
for (prime in primes[primes > 0]){
if ((i%%prime == 0)){
flag <- 1
break
}
if (prime > s){
break
}
}
if (flag == 0){
primes[prime.position] <- i
prime.position<-prime.position+1
}
}
}