It’s been over one month since my last post on Euler problem 20, when  I was planning to post at least one on either Euler project or visualization. So I am four posts behind; I’ll try to catch up. Tonight, I’ll solve the 21st Euler problem. Let’s take a look.

Let d(n) be defined as the sum of proper divisors of n (numbers less than n which divide evenly into n). If d(a) = b and d(b) = a, where a!=b, then a and b are an amicable pair and each of a and b are called amicable numbers. For example, the proper divisors of 220 are 1, 2, 4, 5, 10, 11, 20, 22, 44, 55 and 110; therefore d(220) = 284. The proper divisors of 284 are 1, 2, 4, 71 and 142; so d(284) = 220. Evaluate the sum of all the amicable numbers under 10000.

The new term ”amicable” above is just a definition. What is important is the how to find the amicable numbers. I’m doing it in two steps: first, to find out the divisors for a number from 1:9999 and sum them up; second, to see whether the sum of the divisors of the previous sum is equal to the given number. Sounds a little twisted. The code might explain itself well.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

DivisorsSearch <- function(n) {
  # helper function to find divisors
  if (n <= 3) {
    divisors <- 1
  }
  else {
    candidate <- 2:floor(sqrt(n))
    flag <- n %% candidate == 0
    divisors.half <- candidate[flag]
    divisors <- unique(c(1, divisors.half, n / divisors.half))
  }
  return(divisors)
}
 
# calculate the sum of divisors
num.lst <- 1:9999
divisors.sum <- numeric(length(num.lst))
for (i in 1:length(num.lst)) {
  divisors.sum[i] <- sum(DivisorsSearch(i))
}
 
# discard those with sum bigger than 10000
mt <- divisors.sum < 10000
num.lst.sub <- num.lst[mt]
divisors.sum.sub <- divisors.sum[mt]
flags <- logical(length(num.lst.sub))
amicable <- list()
j <- 1
for (i in 1:length(num.lst.sub)) {
  num1 <- num.lst.sub[i]
  num.divisors.sum <- divisors.sum.sub[i]
  num2 <- divisors.sum[num.divisors.sum]
  if (num1 == num2 & num1 !=num.divisors.sum) {
    flags[i] <- TRUE
    amicable[[j]] <- c(num1, num.divisors.sum)
	# this will give the actual amicable pair, although there is duplication
    j <- j + 1
  }
}
result <- sum(num.lst.sub[flags])
cat("The result is:", result, "\n")

I also draw the scatterplot between the number below 10000 and the sum of divisors. You would notice some divisors sum is much bigger than the number itself. So the amicable number pair should lie in the green area as shown in the following figure.