Problem 10 - Summation of Primes

Train tracks leading to Nine Arches Bridge & and surrounding village
Demodara, Ella, Sri Lanka - April 30th, 2025

The sum of the primes below \(10\) is \(2 + 3 + 5 + 7 = 17\).
Find the sum of all the primes below two million.

https://projecteuler.net/problem=10

Development

R

Answer 1 - Trial Division

sieve <- function(n) {
  x <- 2:n
  p <- c()
  while(length(x) > 0) {
    p <- c(p, x[1])
    # incremental elimination
    x <- x[x %% x[1] != 0]
  }
  return(p)
}

tic()
x <- sieve(2e6)
sum(x)
toc()

[1] 142913828922
78.201 sec elapsed

Answer 2 - {primes} package

?primes::generate_primes tells us that it uses a:

“a fast implementation of the Sieve of Eratosthenes.”

tic()
x <- primes::generate_primes(min = 2, max = 2e6)
sum(x)
toc()

[1] 142913828922
0.039 sec elapsed

This is much faster than Answer 1. How can we make our implementation faster?

Answer 3 - Sieve of Eratosthenes

This idea is very similar to our implementation, but an order of magnitude more efficient. The Sieve of Eratosthenes does appear to implement an elimination logic, but it is not done by trial division.

The x <- x[x %% x[1] != 0] line is technically running trial division by the \(j^{th}\) prime on every number in the set. However, this is inefficient.

Instead, multiples of the \(j^{th}\) prime starting at \(j^2\) can be marked as composite (i.e. *not prime). Logically, this gives us the same result at each increment of the loop, but progressively runs an order of magnitude less calculations.

1st Attempt

sieve <- function(n) {
  x <- 2:n
  i <- 1
  while(x[x != 0][i]^2 <= n) {
    x[seq(x[x != 0][i]^2-1, n-1, x[x != 0][i])] <- 0
    i <- i + 1
  }
  return(x[x != 0])
}

tic()
x <- sieve(2e6)
sum(x)
toc()

[1] 142913828922
3.934 sec elapsed

2nd Optimization

This solution can be optimized a bit as well with a for loop & then writing a single assignmnet to the non-eliminated prime index (x != 0), making less calculations:

sieve <- function(n) {
  x <- 2:n
  for (i in 1:(n-1)) {
    j <- x != 0
    if (x[j][i]^2 > n) {
      return(x[j])
    }
    x[seq(x[j][i]^2-1, n-1, x[j][i])] <- 0
  }
}

tic()
x <- sieve(2e6)
sum(x)
toc()

[1] 142913828922
2.582 sec elapsed

However, this is stil not optimal. We can simply work with primitives (logicals & integers) & loops, not making consecutive index operation calls (e.g. j <- x!= 0), and nested calls (e.g. x[seq(x[j][i]^2-1, n-1, x[j][i])] <- 0).

Optimal

Below is optimized.

sieve <- function(n) {
  # prime candidates
  x <- rep(TRUE, n)
  # number 1 is not prime
  x[1] <- FALSE
  i = 2
  n0 <- floor(sqrt(n))
  for (i in 2:n0) {
    if (x[i]) {
      x[seq(i^2, n, i)] <- FALSE
    }
  }
  return((1:n)[x])
}

tic()
x <- sieve(2e6)
sum(x)
toc()

[1] 142913828922
0.058 sec elapsed

Python

Answer 1 - Trial Division

This is extremely slow.

import time

def sieve(n: float) -> list:
  x = list(range(2, int(n)+1))
  p = []
  while len(x) > 0:
    p.append(x[0])
    # incremental elimination
    x = [i for i in x if i % x[0] != 0]
  return(p)

t0 = time.perf_counter()
x = sieve(2e6)
print(sum(x))
t1 = time.perf_counter()
t = t1 - t0
print(f'{t:.3f} sec elapsed')

142913828922
398.397 sec elapsed

Answer 2 - Sieve of Eratosthenes - Optimized

import time
from math import floor

def sieve(n: float) -> list:
  n = int(n)
  # prime candidates
  x = [True] * n
  # number 1 is not prime
  x[0] = False
  n0 = floor(n**0.5)
  for i in range(2, n0+1):
    if x[i-1]:
      for j in range(i**2, n+1, i):
        x[j-1] = False
  return([i for i, j in zip(list(range(1,n+1)), x) if j])

t0 = time.perf_counter()
x = sieve(2e6)
print(sum(x))
t1 = time.perf_counter()
t = t1 - t0
print(f'{t:.3f} sec elapsed')

142913828922
0.168 sec elapsed