Little useless-useful R functions – Goldbach’s conjecture and Sieve of Sundaram

This is fun It is also O(MAX) complexity. But first some background. Since the problem is super old, we are not intending to solve it, merely to play with it. In the number theory of mathematics, the Goldbach’s conjecture states that for every even integer (greater than 2) can be expressed with the sum of two prime numbers. There are also far cries from this theory. For example, prove that every even number can be written as the sum of not more than 300.000 primes (by Schnirelman (1939)).

To put this to test, we created a simple Prime function, to determine if integer is a prime or not.

# sieve of sundaram
sieve_of_sundaram <- function(limit) {
  n <- (limit - 1) %/% 2
  sieve <- rep(TRUE, n + 1)
  
  for (i in 1:n) {
    j <- 1
    while (i+j+2*i*j <= n) {
      sieve[i+j+2*i*j] <- FALSE
      j <- j + 1
    }
  }
  primes <- c(2,(2*(1:n)+1)[sieve])
  return(primes)
}

# is prime
is_prime <- function(n) {
  if (n <= 1)  return(FALSE)
  if (n <= 3)  return(TRUE)
  if (n %% 2 == 0 || n %% 3 == 0) return(FALSE)
  i <- 5
  while (i*i <= n) {
    if (n %% i == 0 || n %% (i + 2) == 0) return(FALSE)
    i <- i + 6
  }
  return(TRUE)
}
  

For fun, I am adding the famous Sieve of Sundaram algorithm for finding multiple primes in an array of integers.

For the next step, we find sum of all two primes for every even number as input; with other words, the Goldbach’s conjecture.

goldbach_conjecture <- function(even_num) {
  if (even_num <= 2 || even_num %% 2 != 0) {
    return("Number must be even and greater than 2.")
  }
  c <- NULL
  for (i in 2:(even_num / 2)) {
    if (is_prime(i) && is_prime(even_num - i)) {
      #cat("Goldbach's pairs for", even_num, "are:", i, "+", even_num - i, "\n")
      c <- cbind(c,i) # nof solutions
    } 
  }
  return(length(c))
}

But not to end here, let’s put the Goldbach’s conjecture to test and see how the solutions are generated:

#make some 1000 solutions
sol <- NULL
for (i in seq(4,1000, by=2)){
  nof_solutions <- goldbach_conjecture(i)
  sol <- rbind(sol, data.frame(n=i, nof=nof_solutions))
}

# plot solutions; alternating solutions
plot(sol$n, sol$nof, type = "p", xlab = "Even number", ylab = "Number of Solutions", main = "Goldbach's Conjecture") 
reg<-lm(nof ~ n, data = sol)
abline(reg, col="red")

We see the “alternating” pattern of solutions for every even number between 4 and 1000 as sum of two primes.

And the distribution of prime frequencies is equally “interesting”

As always, code is available on Github in  Useless_R_function repository. The sample file in this repository is here (filename: Goldbach_conjecture.R) Check the repository for future updates.

Happy R-coding and stay healthy!