This is a solution for problem 21 on the Project Euler website. It consists of finding the sum of all the amicable numbers under 10000. This was pretty easy to solve, but the solution could probably be improved quite a bit.
Solution #1 in R is as follows (it calculates the proper divisors of each number using prop.divs, and then adds up the sequence of amicable numbers in the main function).
[sourcecode language=”matlab”]
prop.divs <- function(x) {
if (x == 1) return (1)
divs <- integer(30)
j <- 1
divs[j] <- 1
j <- j + 1
for (i in 2:(floor(x/2))) {
if ((x %% i) == 0) {
divs[j] <- i
j <- j + 1
}
}
sum(divs[1:(j-1)])
}
problem.21 <- function(N) {
s <- 0
for (i in 2:N) {
da <- prop.divs(i)
if (da == i) next
db <- prop.divs(da)
if ( db==i ) {
s <- s + da + db
}
}
s/2
}
[/sourcecode]
The s/2 is needed as each factor is added twice during the calculation.
This gives the correct answer, but the implementation is a bit naive. I remembered coming across an article about prime factors and proper divisors on PlanetMath a while ago, and this seemed like potentially a more efficient way to calculate the factors involved. Specifically, the sum of proper divisors of a number n can be given by:
\[ \prod_{i=1}^k\frac{p_i^{m_i+1}-1}{p_i – 1}-\prod_{i=1}^kp_i^{m_i} \]
The second attempt at this problem looked like the following:
[sourcecode language=”r”]
prime.sieve <- function(n) {
a <- seq.int(1,n)
p <- 1
M <- as.integer(sqrt(n))
while ((p <- p + 1) <= M) {
if (a[p] != 0)
a[seq.int(p*p, n, p)] <- 0
}
a[a>1 & a>0]
}
sum.proper.divisors <- function(x) {
primes <- prime.sieve( x )
primes <- primes[ x %% primes == 0]
geo.sum <- numeric(length(primes))
i <- 1
for (prime in primes) {
n <- x
curr <- 0
while (n %% prime == 0) {
curr <- curr + 1
n <- n %/% prime
}
geo.sum[i] <- ( (prime^(curr+1) – 1)/(prime – 1) )
i <- i + 1
}
prod(geo.sum)-x
}
problem.21_2 <- function(N) {
s <- 0
for (i in 2:N) {
da <- sum.proper.divisors(i)
if (da == i) next
db <- sum.proper.divisors(da)
if (db==i) s <- s + da +db
}
s/2
}
[/sourcecode]
This also gives the correct answer, but with much reduced runtime overhead:
> system.time(problem.21(10000))
user system elapsed
103.943 0.511 106.978
> system.time(problem.21_2(10000))
user system elapsed
24.834 0.160 26.565