# Topological Tomfoolery in R: Plotting a Möbius Strip

Topology is, according to Clifford Pickover, the “silly putty of mathematics”. This branch of maths studies the transformation of shapes, knots and other complex geometry problems.  One of the most famous topics in topology is the Möbius strip. This shape has some unusual properties which have inspired many artists, inventors, mathematicians and magicians.

You can make a Möbius strip by taking a strip of paper, giving it one twist and glue the ends together to form a loop. If you now cut this strip lengthwise in half, you don’t end-up with two separate strips, but with one long one.

The Möbius strip can also be described with the following parametric equations (where $0 \leq u \leq 2\pi$, $-1 \leq v \leq 1$ and $R$ is the radius of the loop):

$x(u,v)= \left(R+\frac{v}{2} \cos \frac{u}{2}\right)\cos u$
$y(u,v)= \left(R+\frac{v}{2} \cos\frac{u}{2}\right)\sin u$
$z(u,v)= \frac{v}{2}\sin \frac{u}{2}$

The mathematics of this set of parametric equations is not as compex as it looks. $R$ is the radius of the ring, $u$ is the polar angle of each point and $v$ indicates the width of the strip. The polar angle $u/2$ indicates the number of half twists. To make the ring twist twice, change the anlge to $u$.

For my data science day job, I have to visualise some three-dimensional spaces so I thought I best learn how to do this by visualising a Möbis strip, using these three equations.

## Plotting a Möbius Strip

The RGL package provides the perfect functionality to play with virtual Möbius strips. This package produces interactive three-dimensional plots that you can zoom and rotate. This package has many options to change lighting, colours, shininess and so on. The code to create for plotting a Möbius strip is straightforward.

The first section defines the parameters and converts the $u$ and $v$ sequences to a mesh (from the plot3D package). This function creates two matrices with every possible combination of $u$ and $v$ which are used to calculate the $x, y, z$ points.

The last three lines define a 3D window with a white background and plot the 3D surface in blue. You can explore the figure with your mouse by zooming and rotating it. Parametric equations can be a bit of fun, play with the formula to change the shape and see what happens.

# Moebius strip
library(rgl)
library(plot3D)

# Define Parameters
R <- 5
u <- seq(0, 2 * pi, length.out = 100)
v <- seq(-1, 1, length.out = 100)
m <- mesh(u, v)
u <- m$x v <- m$y

# Móbius strip parametric equations
x <- (R + v/2 * cos(u /2)) * cos(u)
y <- (R + v/2 * cos(u /2)) * sin(u)
z <- v/2 * sin(u / 2)

# Visualise
bg3d(color = "white")
surface3d(x, y, z, color= "blue")


Plotting a Möbius Strip: RGL output.

We can take it to the next level by plotting a three-dimensional Möbius strip, or a Klein Bottle. The parametric equations for the bottle are mind boggling:

$x(u,v) = -\frac{2}{15} \cos u (3 \cos{v}-30 \sin{u}+90 \cos^4{u} \sin{u} -60 \cos^6{u} \sin{u} +5 \cos{u} \cos{v} \sin{u})$

$y(u,v) = -\frac{1}{15} \sin u (3 \cos{v}-3 \cos^2{u} \cos{v}-48 \cos^4{u} \cos{v} + 48 \cos^6{u} \cos{v} - 60 \sin{u}+5 \cos{u} \cos{v} \sin{u}-5 \cos^3{u} \cos{v} \sin{u}-80 \cos^5{u} \cos{v} \sin{u}+80 \cos^7{u} \cos{v} \sin{u})$

$z(u,v) = \frac{2}{15} (3+5 \cos{u} \sin{u}) \sin{v}$

Where: $0 \leq u \leq \pi$ and $0 \leq v \leq 2\leq$.

The code to visualise this bottle is essentially the same, just more complex equations.

u <- seq(0, pi, length.out = 100)
v <- seq(0, 2 * pi, length.out = 100)
m <- mesh(u, v)
u <- m$x v <- m$y
x <- (-2 / 15) * cos(u) * (3 * cos(v) - 30 * sin(u) + 90 * cos(u)^4 * sin(u) - 60 * cos(u)^6 * sin(u) + 5 * cos(u) * cos(v) * sin(u))
y <- (-1 / 15) * sin(u) * (3 * cos(v) - 3 * cos(u)^2 * cos(v) - 48 * cos(u)^4 * cos(v) + 48 * cos(u)^6 * cos(v) - 60 * sin(u) + 5 * cos(u) * cos(v) * sin(u) - 5 * cos(u)^3 * cos(v) * sin(u) - 80 * cos(u)^5 * cos(v) * sin(u) + 80 * cos(u)^7 * cos(v) * sin(u))
z <- (+2 / 15) * (3 + 5 * cos(u) * sin(u)) * sin(v)

bg3d(color = "white")
surface3d(x, y, z, color= "blue", alpha = 0.5)


Plotting a Klein Bottle in RGL. Click to view RGL widget.

The RGL package has some excellent facilities to visualise three-dimensional objects, far beyond simple strips. I am still learning and am working toward using it to visualise bathymetric surveys of water reservoirs. Möbius strips are, however, a lot more fun.

## Creating Real Möbius Strips

Even more fun than playing with virtual Möbius strips is to make some paper versions and start cutting, just like August Möbius did when he did his research. If you like to create a Möbius strip, you can recycle then purchase a large zipper from your local haberdashery shop, add some hook-and-loop fasteners to the ends and start playing. If you like to know more about the mathematics of the topological curiosity, then I can highly recommend Clifford Pickover’s book on the topic.

Möbius strip zipper.

## The Möbius Strip in Magic

In the first half of the twentieth century, many magicians used the Möbius strip as a magic trick. The great Harry Blackstone performed it regularly in his show.

If you are interested in magic tricks and Möbius strips, then you can read my ebook on the Afghan bands.

# Pandigital Products: Euler Problem 32

Euler Problem 32 returns to pandigital numbers, which are numbers that contain one of each digit. Like so many of the Euler Problems, these numbers serve no practical purpose whatsoever, other than some entertainment value. You can find all pandigital numbers in base-10 in the Online Encyclopedia of Interegers (A050278). The Numberhile video explains everything you ever wanted to

The Numberhile video explains everything you ever wanted to know about pandigital numbers but were afraid to ask.

## Euler Problem 32 Definition

We shall say that an n-digit number is pandigital if it makes use of all the digits 1 to n exactly once; for example, the 5-digit number, 15234, is 1 through 5 pandigital.

The product 7254 is unusual, as the identity, 39 × 186 = 7254, containing multiplicand, multiplier, and product is 1 through 9 pandigital.

Find the sum of all products whose multiplicand/multiplier/product identity can be written as a 1 through 9 pandigital.

HINT: Some products can be obtained in more than one way so be sure to only include it once in your sum.

## Proposed Solution

The pandigital.9 function tests whether a string classifies as a pandigital number. The pandigital.prod vector is used to store the multiplication.

The only way to solve this problem is brute force and try all multiplications but we can limit the solution space to a manageable number. The multiplication needs to have ten digits. For example, when the starting number has two digits, the second number should have three digits so that the total has four digits, e.g.: 39 × 186 = 7254. When the first number only has one digit, the second number needs to have four digits.

pandigital.9 <- function(x) # Test if string is 9-pandigital
(length(x)==9 & sum(duplicated(x))==0 & sum(x==0)==0)

t <- proc.time()
pandigital.prod <- vector()
i <- 1
for (m in 2:100) {
if (m < 10) n_start <- 1234 else n_start <- 123
for (n in n_start:round(10000 / m)) {
# List of digits
digs <- as.numeric(unlist(strsplit(paste0(m, n, m * n), "")))
# is Pandigital?
if (pandigital.9(digs)) {
pandigital.prod[i] <- m * n
i <- i + 1
print(paste(m, "*", n, "=", m * n))
}
}
}


Numbers can also be checked for pandigitality using mathematics instead of strings.

You can view the most recent version of this code on GitHub.

# Digit fifth powers: Euler Problem 30

Euler problem 30 is another number crunching problem that deals with numbers to the power of five. Two other Euler problems dealt with raising numbers to a power. The previous problem looked at permutations of powers and problem 16 asks for the sum of the digits of $2^{1000}$.

Numberphile has a nice video about a trick to quickly calculate the fifth root of a number that makes you look like a mathematical wizard.

## Euler Problem 30 Definition

Surprisingly there are only three numbers that can be written as the sum of fourth powers of their digits:

$1634 = 1^4 + 6^4 + 3^4 + 4^4$

$8208 = 8^4 + 2^4 + 0^4 + 8^4$

$9474 = 9^4 + 4^4 + 7^4 + 4^4$

As $1 = 1^4$ is not a sum, it is not included.

The sum of these numbers is $1634 + 8208 + 9474 = 19316$. Find the sum of all the numbers that can be written as the sum of fifth powers of their digits.

## Proposed Solution

The problem asks for a brute-force solution but we have a halting problem. How far do we need to go before we can be certain there are no sums of fifth power digits? The highest digit is $9$ and $9^5=59049$, which has five digits. If we then look at $5 \times 9^5=295245$, which has six digits and a good endpoint for the loop. The loop itself cycles through the digits of each number and tests whether the sum of the fifth powers equals the number.

largest <- 6 * 9^5
for (n in 2:largest) {
power.sum <-0
i <- n while (i > 0) {
d <- i %% 10
i <- floor(i / 10)
power.sum <- power.sum + d^5
}
if (power.sum == n) {
print(n)
}
}


# Euler Problem 29: Distinct Powers

Euler Problem 29 is another permutation problem that is quite easy to solve using brute force. The MathBlog site by Kristian Edlund has a nice solution using only pen and paper.

Raising number to a power can have interesting results. The video below explains why this pandigital formula approximates $e$ to billions of decimals:

$(1 + 9^{-4^{6 \times 7}})^{3^{2^{85}}} \approx e$

## Euler Problem 29 Definition

Consider all integer combinations of: $a^b$ for $2 \leq a \leq 5$ and $\leq b \leq 5$.

$2^2=4, \quad 2^3 = 8,\quad 2^4 = 16,\quad 2^5 = 32$

$3^2 = 9,\quad 3^3 = 27,\quad 3^4 = 81,\quad 3^5 = 243$

$4^2 = 16,\quad 4^3 = 64,\quad 4^4 = 256, \quad 4^5 = 1024$

$5^2 = 25,\quad 5^3 = 125,\quad 5^4 = 625,\quad 5^5 = 3125$

If they are then placed in numerical order, with any repeats removed, we get the following sequence of 15 distinct terms:

$4, \ 8, \ 9, \ 16, \ 25, \ 27, \ 32, \ 64, \ 81, \ 125, \ 243, \ 256,\ 625, \ 1024, \ 3125$

How many distinct terms are in the sequence generated by $a^b$ for $2 \leq a \leq 100$ and $2 \leq b \leq 100$?

## Brute Force Solution

This code simply calculates all powers from $2^2$ to $2^{1000}$ and determines the number of unique values. Since we are only interested in their uniqueness and not the precise value, there is no need to use Multiple Precision Arithmetic.

# Initialisation
target <- 100
terms <- vector()
i <- 1
# Loop through values of a and b and store powers in vector
for (a in 2:target) {
for (b in 2:target) {
terms[i] <- a^b
i <- i + 1
}
}
# Determine the number of distinct powers


# The Ulam Spiral: Euler Problem 28

Euler Problem 28 takes us to the world of the Ulam Spiral. This is a spiral that contains sequential positive integers in a square spiral, marking the prime numbers. Stanislaw Ulam discovered that a lot of primes are located along the diagonals. These diagonals can be described as polynomials. The Ulam Spiral is thus a way of generating quadratic primes (Euler Problem 27).

Ulam Spiral (WikiMedia).

## Euler Problem 28 Definition

Starting with the number 1 and moving to the right in a clockwise direction a 5 by 5 spiral is formed as follows:

21 22 23 24 25
20 07 08 09 10
19 06 01 02 11
18 05 04 03 12
17 16 15 14 13

It can be verified that the sum of the numbers on the diagonals is 101. What is the sum of the numbers on the diagonals in a 1001 by 1001 spiral formed in the same way?

## Proposed Solution

To solve this problem we do not need to create a matrix. This code calculates the values of the corners of a matrix with size $n$. The lowest number in the matrix with size $n$ is $n(n-3)+4$. The numbers increase by $n-1$.

The code steps through all matrices from size 3 to 1001. The solution uses only the uneven sized matrices because these have a centre. The answer to the problem is the sum of all numbers.

size <- 1001 # Size of matrix
answer <- 1 # Starting number
# Define corners of subsequent matrices
for (n in seq(from = 3, to = size, by = 2)) {
corners <- seq(from = n * (n - 3) + 3, by = n - 1, length.out = 4)
}


## Plotting the Ulam Spiral

We can go beyond Euler Problem 28 and play with the mathematics. This code snippet plots all the prime numbers in the Ulam Spiral. Watch the video for an explanation of the patterns that appear along the diagonals.

Ulam Spiral prime numbers.

The code creates a matrix of the required size and fills it with the Ulam Spiral. The code then identifies all primes using the is.prime function from Euler Problem 7. A heat map visualises the results.

# Ulam Spiral
size <- 201 # Size of matrix
ulam <- matrix(ncol = size, nrow = size)
mid <- floor(size / 2 + 1)
ulam[mid, mid] <- 1
for (n in seq(from = 3, to = size, by = 2)) {
numbers <- (n * (n - 4) + 5) : ((n + 2) * ((n + 2) - 4) + 4)
d <- mid - floor(n / 2)
l <- length(numbers)
ulam[d, d:(d + n - 1)] <- numbers[(l - n + 1):l]
ulam[d + n - 1, (d + n - 1):d] <- numbers[(n - 1):(n - 2 + n)]
ulam[(d + 1):(d + n - 2), d] <- numbers[(l - n):(l - 2 * n + 3)]
ulam[(d + 1):(d + n - 2), d + n - 1] <- numbers[1:(n - 2)]
}
ulam.primes &lt;- apply(ulam, c(1, 2), is.prime)

# Visualise
library(ggplot2)
library(reshape2)
ulam.primes <- melt(ulam.primes)
ggplot(ulam.primes, aes(x=Var1, y=Var2, fill=value)) + geom_tile() +
scale_fill_manual(values = c("white", "black")) +
guides(fill=FALSE) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank()) +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank()
)


# Generating Quadratic Primes: Euler Problem 27

Prime numbers are the engine of the Internet economy. One of the reasons prime numbers are so useful is that they cannot be generated through an algorithm. This impossibility has not stopped mathematicians from trying to find formulas to generate prime numbers.

Euler problem 27 deals with quadratic formulas that can be used to generate sets of prime numbers. We have already discussed this in the post about the Ulam Spiral. This Numerphile video discusses quadratic primes.

## Euler Problem 27 Definition

Euler discovered the remarkable quadratic formula:

$n^2+n+41$

It turns out that the formula will produce $40$ primes for the consecutive integer values $0 \leq n \leq 39$. However, when $n=40$, $40^2+40+41=40(40+1)+41$ is divisible by $41$, and certainly when $n=41$, $41^2+41+41$ is clearly divisible by $41$.

The incredible formula $n^2-79n+1601$ was discovered, which produces 80 primes for the consecutive values $0 \leq n \leq 79$. The product of the coefficients, $-79$ and $1601$, is $-126479$.

Considering quadratics of the form: $n^2+an+bn^2+an+b$,

where $|a| < 1000$ and $|b| \leq 1000$, where $|n|$ is the modulus/absolute value of $n$, e.g. $|11|=11$ and $|-4|=4$.

Find the product of the coefficients, $a$ and $b$, for the quadratic expression that produces the maximum number of primes for consecutive values of $n$, starting with $n=0$.

## Proposed Solution

The only way to solve this problem is through brute force and reduce the solution space to optimise it for speed (source: mathblog.dk). Because the outcome of the equation must be prime for $n = 0$, $b$ also has to be prime. We can use the prime sieve from Euler Problem 3, which reduces it from 2000 to 168 options. When we insert $n = 1$ it follows that a has to be an even number. If $1+a+b$ has to be prime and $b$ has to be a prime number, then a can only be an odd number. However, when $b = 2$, a has to be even.

## Euler Problem 27 code

seq.a <- seq(-999, 1001, 2) # a has to be odd
seq.b <- (esieve(1000)) # b has to be prime
max.count <- 0
for (a in seq.a) {
if (a == 2)
seq.a <- seq(-1000, 1000, 2) # a has to be even
for (b in seq.b) {
n <- 0 # Find sequence of primes for a and b
while (is.prime(n^2 + a * n + b)) {
n <- n + 1 } # Store maximum values if (n > max.count) {
max.count <- n
max.a <- a
max.b <- b
}
}
}


# The Viral Recurring Decimal: Euler Problem 26

A few years ago a fraction broke the internet. What happens when you divide 1 by 998001?

$\frac{1}{998001} = 0.000001002003004005006007008009010011012013014015 \ldots$

What is special about this fraction is that it lists every three-decimal number except for 998. Look carefully at the sequence to see that is 000, 001, 0002, 003, 004, 005 and so on. After it has reached 999, the sequence continues from the start. This fraction thus has 2997 recurring decimals. James Grime from Numberphile explains this mathematical oddity with his usual enthusiasm.

The decimal fraction of 1/998001 is a recurring decimal. These are decimal numbers with periodic digits (repeating its values at regular intervals). Euler problem 26 asks us to analyse recurring decimals (reciprocal cycles).

## Euler Problem 26 Definition

A unit fraction contains 1 in the numerator. The decimal representation of the unit fractions with denominators 2 to 10 are given:

$1/2 = 0.5$
$1/3 = 0.(3)$
$1/4 = 0.25$
$1/5 = 0.2$
$1/6 = 0.1(6)$
$1/7 = 0.(142857)$
$1/8 = 0.125$
$1/9 = 0.(1)$
$1/10 = 0.1$

Where 0.1(6) means 0.166666…, and has a 1-digit recurring cycle. It can be seen that 1/7 has a 6-digit recurring cycle. Find the value of d < 1000 for which 1/d contains the longest recurring cycle in its decimal fraction part.

## Solution

A051626 describes the length of the recurring numbers in 1/n in the On-Line Encyclopaedia of Integer Sequences. To solve Euler Problem 26, we need to generate the first 1000 numbers of this sequence and find out which number has the longest recurring cycle.

R can only display up to 22 decimals by using options(digits=22). The base R capability is unsuitable for solving this problem, so I wrote some code to perform long division the old-fashioned way.

The recur function divides 1 by any arbitrary integer. The code continues until the decimal terminates, for example 1/4 = 0.25, or when a recurring pattern emerges, e.g. 1/7 = 0.(142857).

The function has two arguments: n is the input number. The output argument determines the outcome of the function: “len” for the length of the recurring decimals. Any other value shows the result of the calculation. The output of the function is a string. Using the European notation, the recurring part of the decimals is shown between brackets, e.g. 1/14 = 0.0(714285).

recur <- function(x, output = "") {
# Prepare variable
if (x == 0) return(NaN)
if (x == 1) return(0)
x <- floor(abs(x))
# Initiate vectors to store decimals and remainders
dec <- vector()
rem <- vector()
# Initiate values
i <- 1
r <- 10
rem <- r
# Long division
repeat {
dec[i] <- floor(r / x)
r <- 10 * (r %% x)
# Test wether the number is terminating or repeating
if (r == 0 | r %in% rem) break
rem[i + 1] <- r
i <- i + 1
}
# Determine number of recurring digits
rep <- ifelse(r != 0, length(rem) - which(r == rem) + 1, 0)
# Output
if (output == "len")
return(rep)
else {
if (rep != 0) {
if (rep == length(dec))
l <- "("
else
l <- c(dec[1:(length(dec) - rep)], "(")
dec <- c(l, dec[(length(dec) - rep + 1):length(dec)], ")")
}
return(paste0("0.", paste0(dec, collapse = "", sep = "")))
}
}

A051626 <- sapply(1:1000, recur, "len")

recur(998001)


# Large integers in R: Fibonacci number with 1000 digits, Euler Problem 25

The Fibonacci Sequence occurs in nature: The nautilus shell.

Euler Problem 25 takes us back to the Fibonacci sequence and the problems related to working with very large integers.

The Fibonacci sequence follows a simple mathematical rule but it can create things of great beauty. This pattern occurs quite often in nature, like to nautilus shell shown in the image. The video by Arthur Benjamin at the end of this post illustrates some of the magic of this sequence.

## Large Integers in R

By default, numbers with more than 7 digits are shown in scientific notation in R, which reduces the accuracy of the calculation. You can change the precision of large integers with the options function but R struggles with integers with more than 22 digits. This example illustrates this issue.

factorial(24)
[1] 6.204484e+23
> options(digits=22)
> factorial(24)
[1] 620448401733239409999872


However, finding a number of 1000 digits is a problem with using special functions.

## Euler Problem 25 Definition

The Fibonacci sequence is defined by the recurrence relation:

$F_n = F_{n-1} + F_{n-2}$, where $F_1 = 1$ and $F_2 = 1$. The 12th term, $F_{12}$, is the first term to contain three digits.

What is the index of the first term in the Fibonacci sequence to contain 1000 digits?

## Proposed Solutions

The first solution uses the GMP library to manage very large integers. This library also contains a function to generate Fibonacci numbers. This solution cycles through the Fibonacci sequence until it finds a number with 1000 digits.

library(gmp) # GNU Multiple Precision Arithmetic Library
n <- 1
fib <- 1
while (nchar(as.character(fib)) < 1000) {
fib <- fibnum(n) # Determine next Fibonacci number
n <- n + 1
}


This is a very fast solution but my aim is to solve the first 100 Project Euler problems using only base-R code. The big.add function I developed to solve Euler Problem 13.

t <- proc.time()
fib <- 1 # First Fibonaci number
cur <- 1 # Current number in sequence
pre <- 1 # Previous number in sequence
index <- 2
while (nchar(fib) < 1000) {
fib <- big.add(cur, pre) # Determine next Fibonacci number
pre <- cur
cur <- fib
index <- index + 1
}


This code is much slower than the GMP library but it was fun to develop.

# Lexicographic Permutations: Euler Problem 24

Euler Problem 24 asks to develop lexicographic permutations which are ordered arrangements of objects in lexicographic order. Tushar Roy of Coding Made Simple has shared a great introduction on how to generate lexicographic permutations.

## Euler Problem 24 Definition

A permutation is an ordered arrangement of objects. For example, 3124 is one possible permutation of the digits 1, 2, 3 and 4. If all of the permutations are listed numerically or alphabetically, we call it lexicographic order. The lexicographic permutations of 0, 1 and 2 are:

012 021 102 120 201 210

What is the millionth lexicographic permutation of the digits 0, 1, 2, 3, 4, 5, 6, 7, 8 and 9?

## Brute Force Solution

The digits 0 to 9 have $10! = 3628800$ permutations (including combinations that start with 0). Most of these permutations are, however, not in lexicographic order. A brute-force way to solve the problem is to determine the next lexicographic permutation of a number string and repeat this one million times.

nextPerm <- function(a) {
# Find longest non-increasing suffix
i <- length(a) while (i > 1 && a[i - 1] >= a[i])
i <- i - 1
# i is the head index of the suffix
# Are we at the last permutation?
if (i <= 1) return (NA)
# a[i - 1] is the pivot
# Find rightmost element that exceeds the pivot
j <- length(a)
while (a[j] <= a[i - 1])
j <- j - 1
# Swap pivot with j
temp <- a[i - 1]
a[i - 1] <- a[j]
a[j] <- temp
# Reverse the suffix
a[i:length(a)] <- rev(a[i:length(a)])
return(a)
}

numbers <- 0:9
for (i in 1:(1E6 - 1)) numbers <- nextPerm(numbers)


This code takes the following steps:

1. Find largest index $i$ such that $a_{i-1} < a_i$.
1. If no such index exists, then this is already the last permutation.
2. Find largest index $j$ such that $j \geq i$ and $a_j > a_{i-1}$.
3. Swap $a_j$ and $a_{i-1}$.
4. Reverse the suffix starting at $a_i$.

## Combinatorics

A more efficient solution is to use combinatorics, thanks to MathBlog. The last nine digits can be ordered in $9! = 362880$ ways. So the first $9!$ permutations start with a 0. By extending this thought, it follows that the millionth permutation must start with a 2.

$\lfloor (1000000 - 1) / 9! \rfloor = 2$

From this rule, it follows that the 725761st permutation is 2013456789. We now need 274239 more lexicographic permutations:

$(1000000 - 1) - (2 \times 9!) = 274239$

We can repeat this logic to find the next digit. The last 8 digits can be ordered in 40320 ways. The second digit is the 6th digit in the remaining numbers, which is 7 (2013456789).

$\lfloor 274239 / 8! \rfloor = 6$

$274239 - (6 \times 7!) = 32319$

This process is repeated until all digits have been used.

numbers <- 0:9
n <- length(numbers)
remain <- 1E6 - 1
for (i in 1:n) {
j <- floor(remain / factorial(n - i))
remain <- remain %% factorial(n - i)
numbers <- numbers[-(j + 1)]
}


R blogger Tony’s Bubble Universe created a generalised function to solve this problem a few years ago.

# Euler Problem 23: Non-Abundant Sums

A demonstration of the abundance of the number 12 using Cuisenaire rods (Wikipedia).

Euler problem 23 asks us to solve a problem with abundant or excessive numbers.

These are numbers for which the sum of its proper divisors is greater than the number itself.

12 is an abundant number because the sum of its proper divisors (the aliquot sum) is larger than 12: (1 + 2 + 3 + 4 + 6 = 16).

All highly composite numbers or anti-primes greater than six are abundant numbers. These are numbers that have so many divisors that they are considered the opposite of primes, as explained in the Numberphile video below.

## Euler Problem 23 Definition

A perfect number is a number for which the sum of its proper divisors is exactly equal to the number. For example, the sum of the proper divisors of 28 would be 1 + 2 + 4 + 7 + 14 = 28, which means that 28 is a perfect number.

A number n is called deficient if the sum of its proper divisors is less than n and it is called abundant if this sum exceeds n.

As 12 is the smallest abundant number, 1 + 2 + 3 + 4 + 6 = 16, the smallest number that can be written as the sum of two abundant numbers is 24. By mathematical analysis, it can be shown that all integers greater than 28123 can be written as the sum of two abundant numbers. However, this upper limit cannot be reduced any further by analysis, even though it is known that the greatest number that cannot be expressed as the sum of two abundant numbers is less than this limit.

Find the sum of all the positive integers which cannot be written as the sum of two abundant numbers.

## Solution

This solution repurposes the divisors function that determines the proper divisors for a number, introduced for Euler Problem 21. The first code snippet creates the sequence of all abundant numbers up to 28123 (sequence A005101 in the OEIS). An abundant number is one where its aliquot sum is larger than n.

# Generate abundant numbers (OEIS A005101)
A005101 <- function(x){
abundant <- vector()
a <- 1
for (n in 1:x) {
aliquot.sum <- sum(proper.divisors(n)) - n
if (aliquot.sum > n) {
abundant[a] <- n
a <- a + 1
}
}
return(abundant)
}

abundant <- A005101(28123)


The solution to this problem is also a sequence in the Online Encyclopedia of Integer Sequences (OEIS A048242). This page states that the highest number in this sequence is 20161, not 28123 as stated in the problem definition.

The second section of code creates a list of all potential numbers not the sum of two abundant numbers. The next bit of code sieves any sum of two abundant numbers from the list. The answer is determined by adding remaining numbers in the sequence.

# Create a list of potential numbers that are not the sum of two abundant numbers
A048242 <- 1:20161

# Remove any number that is the sum of two abundant numbers
for (i in 1:length(abundant)) {
for (j in i:length(abundant)) {
if (abundant[i] + abundant[j] <= 20161) {
A048242[abundant[i] + abundant[j]] <- NA
}
}
}
A048242 <- A048242[!is.na(A048242)]