Summation of primes

Find the sum of the first 2 million primes.

Answer: We use the sieve of Eratosthenes. Let us start small and work up. Find the sum of primes less than 50…

n.top <- 50
sieve <- 2:n.top 
p.current <- 2
while (p.current < sqrt(n.top)) {
  # eliminate multiples of p.current but not p.current:
  sieve <- sieve[ sieve %% p.current != 0 | sieve == p.current]
  # now we want to go to the next prime
  p.current <- sieve[which(sieve == p.current) + 1]
}
sieve
##  [1]  2  3  5  7 11 13 17 19 23 29 31 37 41 43 47
sum(sieve)
## [1] 328

So that is pretty good. Of course, we waste a lot of time by adding the even numbers and then eliminating them. We modify the above code slightly to avoid this inefficiency:

n.top <- 50
sieve <- c(2,seq(3,n.top,2)) 
p.current <- 3
while (p.current < sqrt(n.top)) {
  # eliminate multiples of p.current but not p.current:
  sieve <- sieve[ sieve %% p.current != 0 | sieve == p.current]
  # now we want to go to the next prime
  p.current <- sieve[which(sieve == p.current) + 1]
}
sieve
##  [1]  2  3  5  7 11 13 17 19 23 29 31 37 41 43 47
sum(sieve)
## [1] 328

OK, we are ready to sum the primes less than 2 million!

n.top <- 2000000
sieve <- c(2,seq(3,n.top,2)) 
p.current <- 3
while (p.current < sqrt(n.top)) {
  # eliminate multiples of p.current but not p.current:
  sieve <- sieve[ sieve %% p.current != 0 | sieve == p.current]
  # now we want to go to the next prime
  p.current <- sieve[which(sieve == p.current) + 1]
}
sum(sieve)
## [1] 142913828922