Finding antipodes using the globe and ggmap packages

The antipodes of each point on the Earth's surface

The antipode of each point on the Earth’s surface—the points where the blue and yellow overlap, are the land antipodes.

When I was a kid, I was was fascinated by the conundrum of what happens when you drill a hole straight through the centre of the earth. I always believed that I would turn up in Australia. But is this really the case?

The antipodes of any place on earth is the place that is diametrically opposite to it. A pair of antipodes are connected by a straight line running through the centre of the Earth. These points are as far away from each other as is possible on this planet. Two people are antipodes when they live on opposite sides of the globe. Their feet (πούς/pous in Ancient Greek) are directly opposite each other.

How can we use coding in R to solve this conundrum?

Using R to find antipodes

We can roughly recreate the antipodean map on Wikipedia with the globe package. This package, written by Adrian Baddeley, plots 2D and 3D views of the earth. The package contains a data file with major coastlines that can be used to create a flipped map of the world.

The package contains a data file with major coastlines that can be used to create a flipped map of the world. To turn a spatial location into its antipode you subtract 180 degrees from the longitude and reverse the sign of the latitude, shown below.

# Create antipodean map
flip <- earth
flip$coords[,1] <- flip$coords[,1]-180
flip$coords[,2] <- -flip$coords[,2]
par(mar=rep(0,4)) # Remove plotting margins
globeearth(eye=c(90,0), col="blue")
globepoints(loc=flip$coords, eye=c(90,0), col="red", pch=".")

We can also use the ggmap package to visualise antipodes. This package, developed by David Kahle antipodean R-guru Hadley Wickham, has a neat geocoding function to obtain a spatial location. The antipode function takes the description of a location and a zoom level to plot a dot on the antipode location. The gridExtra package is used to created a faceted map, which is not otherwise possible in ggmap.


antipode <- function(location, zm=6) {
    # Map location
    lonlat <- geocode(location)
    loc1 <- get_map(lonlat, zoom=zm)
    map1 <- ggmap(loc1) + geom_point(data=lonlat, aes(lon, lat, col="red", size=10)) + 
    # Define antipode
    lonlat$lon <- lonlat$lon-180
    if (lonlat$lon < -180) lonlat$lon <- 360 + lonlat$lon
    lonlat$lat <- -lonlat$lat
    loc2 <- get_map(lonlat, zoom=zm)
    map2 <- ggmap(loc2) + geom_point(data=lonlat, aes(lon, lat, col="red", size=10)) + 
    grid.arrange(map1, map2, nrow=1)

antipode("Rector Nelissenstraat 47 Hoensbroek", 4)

This code solves the problem I was thinking about as a child. Running the code shows that the antipodes location of the home I grew up in is not in Australia, but quite a way south of New Zealand. Another childhood fantasy shattered …

Antipodes using ggmap and gridExtra.

Euler Problem 8: Largest Product in a Series

Euler Problem 8 is a combination of mathematics and text analysis. The problem is defined as follows:

Euler Problem 8 Definition

The four adjacent digits in the 1,000-digit number below that have the greatest product are 9 \times 9 \times 8 \times 9 = 5832.


Find the thirteen adjacent digits in the 1,000-digit number that have the greatest product. What is the value of this product?


The first step is to define the digits as a character string. The answer is found by cycling through all 13-character n-grams to find the highest product using the prod function.

# Define digits
digits <- "7316717653133062491922511967442657474235534919493496983520312774506326239578318016984801869478851843858615607891129494954595017379583319528532088055111254069874715852386305071569329096329522744304355766896648950445244523161731856403098711121722383113622298934233803081353362766142828064444866452387493035890729629049156044077239071381051585930796086670172427121883998797908792274921901699720888093776657273330010533678812202354218097512545405947522435258490771167055601360483958644670632441572215539753697817977846174064955149290862569321978468622482839722413756570560574902614079729686524145351004748216637048440319989000889524345065854122758866688116427171479924442928230863465674813919123162824586178664583591245665294765456828489128831426076900422421902267105562632111110937054421750694165896040807198403850962455444362981230987879927244284909188845801561660979191338754992005240636899125607176060588611646710940507754100225698315520005593572972571636269561882670428252483600823257530420752963450"

ngram <- 13  # Define length
answer <- 0
# Clycle through digits
for (i in 1:(nchar(digits)-ngram+1)) {
    # Pick 13 consecutive digits
    adjecent <- substr(digits, i, i + ngram - 1)
    # Define product
    mult <- prod(as.numeric(unlist(strsplit(adjecent, "")))) # Largest? if (mult > answer) 
        answer <- mult

Finding n-grams is a basic function of text analysis called tokenization. The Google Ngram Viewer provides the ability to search for word strings (n-grams) in the massive Google books collection. This type of data can be used to define writing styles and analyse the evolution of language.

Euler Problem 7: 10,001st Prime

Euler Problem 7 Definition

By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13. What is the 1,0001st prime number?


The function determines whether a number is a prime number by checking that it is not divisible by any prime number up to the square root of the number.

The Sieve of used in Euler Problem 3 can be reused to generate prime numbers.

This problem can only be solved using brute force because prime gaps (sequence A001223 in the OEIS) do not follow a predictable pattern. <- function(n) {
    primes <- esieve(ceiling(sqrt(n)))

i <- 2 # First Prime
n <- 1 # Start counter
while (n<10001) { # Find 10001 prime numbers
    i <- i + 1 # Next number
    if( { # Test next number
        n <- n + 1 # Increment counter
        i <- i + 1 # Next prime is at least two away

answer <- i-1

The largest prime gap for the first 10,001 primes is 72. Sexy primes with a gap of 6 are the most common and there are 1270 twin primes.

Euler Problem 7: Prime gap frequency distribution for the first 10001 primes.

Prime gap frequency distribution for the first 10001 primes.

Euler Problem 6: Sum Square Difference

Euler Problem 6 Definition

The sum of the squares of the first ten natural numbers is:

1^2 + 2^2 + \ldots + 10^2 = 385

The square of the sum of the first ten natural numbers is:

(1 + 2 + \ldots + 10)^2 = 552 = 3025

The difference between the sum of the squares of the first ten natural numbers and the square of the sum is 3025 - 385 = 2640. Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.


This is a straightforward problem for vector processing capabilities in R.

answer <- sum(1:100)^2-sum((1:100)^2)

This problem can also be solved arithmetically. When Carl Friedrich Gauss (1777-1855) when he was a child his teacher challenged his students to add all numbers from 1 to 100. All of his friends struggled to add all the 100 numbers one by one but Carl completed the task in a few seconds.

The same principle applies to computing. The first solution is like Carl’s classmates who slavishly add all numbers. This solution is based on arithmetic progressions.

The sum of natural numbers can be expressed as:

\frac{n(n + 1)}{2}
The sum of the squares of the first n natural numbers is:

With these two formulas, a fast solution without having to use loops can be written.

n <- 100
answer <- ((n*(n+1))/2)^2 - (n*(n+1)*(2*n+1))/6

Euler Problem 5: Smallest Multiple

Euler Problem 5 relates to the divisibility of numbers.

Euler Problem 5

2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder. What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?


The solution will also be divisible by the number 1 to 10 so we can start at 2520 and increment by 2520. The loop checks whether the number is divisible by the numbers 1 to 20.

# Start as high as possible
i <- 2520
# Check consecutive numbers for divisibility by 1:20
while (sum(i%%(1:20)) != 0) {
    i <- i + 2520 # Increase by smallest number divisible by 1:10
answer <- i

Euler Problem 4: Largest Palindromic Product

Euler Problem 4 Definition

A palindromic number reads the same both ways. The largest palindrome made from the product of two 2-digit numbers is 9009 = 91 × 99. Find the largest palindrome made from the product of two 3-digit numbers.


This code searches fo palindromic numbers, starting at the highest values. The palindromes are tested by converting the number to a character string. When the first palindromic number is found, the loop is broken.

# Cycles through all number combinations, starting at 999
for (i in 999:900) {
    for (j in 990:900) {
        word <- as.character(i * j)
        # Create reverse
        reverse <- paste(rev(unlist(strsplit(word, ""))), collapse = "")
        # Check whether palindrome
        palindrome <- word == reverse
        if (palindrome) 
    if (palindrome) {
answer <- i * j

Searching a bit further results in twelve palindromic numbers. Sequence A002113 in the OEIS lists palindromes in base 10. Mathematicians used these numbers only for fun, without any purpose in reality. The graph below shows the number of palindromic numbers between 0 and 1 million.

Euler Problem 4 : Palindromic numbers

Palindromic number count

Euler Problem 3: Largest Prime Factor

Euler Problem 3 Definition

The prime factors of 13195 are 5, 7, 13 and 29. What is the largest prime factor of the number 600851475143?

Generating Prime Numbers

This solution relies on two functions that can be used for multiples problems. The Sieve of Eratosthenes generates prime numbers from 2 to n. The code is commented to explain the sieve and the image shows how numbers from 1 to 100 are sieved to find the primes.

Sieve of Eratosthenes.

Sieve of Eratosthenes. Green are multiples of 2, Blue are multiples of 3, the orange coloured numbers are multiples of 5 and purple the multiples of 7. The remaining numbers (except 1) are the prime numbers (black).

The prime.factors function generates the list of unique prime divisors and then generates the factors. The factors are identified by dividing the number by the candidate prime factors until the result is 1.


# Sieve of Eratosthenes for generating primes 2:n
esieve <- function(n) {
    if (n==1) return(NULL)
    if (n==2) return(n)
    # Create a list of consecutive integers {2,3,…,N}.
    l <- 2:n
    # Start counter
    i <- 1
    # Select p as the first prime number in the list, p=2.
    p <- 2
    while (p^2<=n) {
        # Remove all multiples of p from the l.
        l <- l[l==p | l%%p!=0]
        # set p equal to the next integer in l which has not been removed.
        i <- i+1 # Repeat steps 3 and 4 until p2 > n, all the remaining numbers in the list are primes
        p <- l[i]

# Prime Factors
prime.factors <- function (n) {
    factors <- c() # Define list of factors
    primes <- esieve(floor(sqrt(n))) # Define primes to be tested
    d <- which(n%%primes == 0) # Idenitfy prime divisors
    if (length(d) == 0) # No prime divisors
    for (q in primes[d]) { # Test candidate primes
        while (n%%q == 0) { # Generate list of factors
            factors <- c(factors, q)
            n <- n/q } } if (n > 1) factors <- c(factors, n)


The solution can also be found by using the primeFactors function in the numbers package. This package provides a range of functions related to prime numbers and is much faster than the basic code provided above.


Euler Problem 2: Even Fibonacci Numbers

Euler Problem 2 Definition

Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first 10 terms will be:

1, 2, 3, 5, 8, 13, 21, 34, 55, 89, \ldots

By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms.


The code generates Fibonacci numbers until it reaches the value of four million. The code then sums the even numbers in the sequence.

fib <- c(1, 2)  #Define first two numbers
while (max(fib) < 4e+06) {
    # Generate Fibonacci numbers until limit is reached
    len <- length(fib)
    fib <- c(fib, fib[len - 1] + fib[len])
answer <- sum(fib[fib%%2 == 0])

A range of packages exists that generate Fibonacci numbers. The gmp package for Multiple Precision Arithmetic and the numbers package supplies functions to calculate the nth Fibonacci number. This package is also able to work with very large numbers.

i <- 1
answer <- 0
fib <- fibnum(1)
while (fibnum(i) <= 4e6) {
    fib <- fibnum(i)
    if (fib%%2==0) answer <- answer + fib
    i <- i + 1

Fibonacci Numbers as a magic trick

Fibonacci numbers describe many natural processes and can also be used to create magic tricks. The Missing Square Puzzle is based on this principle.

By Trekky0623 at English Wikipedia (Transferred from en.wikipedia to Commons.) [Public domain], via Wikimedia Commons

By Trekky0623 at English Wikipedia (Transferred from en.wikipedia to Commons.) [Public domain], via Wikimedia Commons

Euler Problem 1: Multiples of 3 or 5

Euler Problem 1 Definition

If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23. Find the sum of all the multiples of 3 or 5 below 1000.


There are four ways to solve this problem in R.

  1. Brute force loop through all numbers from 1 to 999 and test whether they are divisible by 3 or by 5 using the modulus function.
  2. Vector arithmetics.
  3. Sequences of 3 and 5, excluding duplicates (numbers divisible by 15).
  4. Using an arithmetic approach.

By the way, the problem definition on the Project Euler website is not consistent: the title mentions multiples of 3 AND 5, while the description asks for multiples of 3 OR 5.

# Solution 1
answer <- 0
for (i in 1:999) {
    if (i%%3 == 0 | i%%5 == 0) 
        answer <- answer + i

# Solution 2
sum((1:999)[((1:999)%%3 == 0) | ((1:999)%%5 == 0)])

# Solution 3
sum(unique(c(seq(3, 999, 3), seq(5, 999, 5))))

These three brute-force solutions all take time to process, especially when using targets much igher than 999. An analytical solution will significantly reduce the processing time.

Analytical Solution

The sum of an arithmetic progression , where n is the number of elements and a1 and an are the lowest and highest value, is:

\mathrm{sum}= \frac{n(a_{1} + a_n)}{2}

The numbers divisible by n=3 can be expressed as:

\mathrm{sum}_3(999)=3+6+9+12+ \ldots +999 = 3(1+2+3+4+ \ldots +333)

We can now calculate the sum of all divisors by combining the above progression with the formula for arithmetic progressions as expressed in the above code, where m is the divisor and n the extent of the sequence.

p is the highest number less than n divisible by m. In the case of 5, this number is 995.

p = n \lfloor (m/n) \rfloor

Substitution gives:

\mathrm{sum}_m(n) =  p \frac{1+(p/m)}{2}

# Solution 4
SumDivBy <- function(m, n) {
    p <- floor(n/m)*m # Round to multiple of n
    return (p*(1+(p/m))/2)

answer <- SumDivBy(3, 999) + SumDivBy(5, 999) - SumDivBy(15, 999)