**Tony's bubble universe » R**, and kindly contributed to R-bloggers)

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.

^{?}View Code RSPLUS

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.

**leave a comment**for the author, please follow the link and comment on his blog:

**Tony's bubble universe » R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...