Optimize simple r code for Project Euler 12

Question

The idea of Project Euler question 12 is to find the smallest triangular number with a specified number of divisors(https://projecteuler.net/problem=12). As an attempt to solve this problem, I wrote the following code:

# This function finds the number of divisors of a number and returns it. 
FUN <- function(x) {
    i = 1
    lst = integer(0)
    while(i<=x)
    {
        if(x %% i ==0)
        {
          lst = c(lst, i)
        }
      i = i +1
    }

  return(lst)
}

and

n = 1 
i=1
while (length(FUN(n))<500)
{ 
   i = i + 1
   n = n + i
}

This code is producing the correct answer for few smaller test cases: length(FUN(n))<4 will produce 6, and length(FUN(n))<6 will produce 28.

However, this simple looking code is taking over 24 hours to run (and still running) for length(FUN(n))<500. I understand that for a number to have 500 divisors, the number is probably very big, but I am wondering why is it taking so long to run.

You're growing vectors element by element, which is inherently slow in R because of memory allocation. You get rid of the `while` loop in `FUN` by vectorizing the operations to avoid the issue in part, though that still leaves the outer loop. So faster, but not fast (and tidied a bit) for the first 10 triangle numbers: `triangles <- data.frame(number = cumsum(1:10)); triangles$n_factors <- sapply(triangles$number, function(x){sum(x %% seq.int(x) == 0)}); triangles` — alistaire, Aug 31 '17 at 05:51
Have a read of [this post](https://stackoverflow.com/questions/6424856/r-function-for-returning-all-factors) — Mist, Aug 31 '17 at 06:05
@alistaire your function could be `function(x){sum(x %% seq.int(sqrt(x)) == 0)*2})`. Would be a lot quicker but numbers that are squared get one more factor than they actually should have. Only an issue if you find a number with exactly 500 divisors. — Mist, Aug 31 '17 at 06:40
@Mist Yeah, I remember when I first solved this question I used an approach which was significantly mathematically optimized. There are ways to make a solution really fast, but it takes some clever code. — alistaire, Aug 31 '17 at 15:06

score 2 · Accepted Answer · answered Aug 31 '17 at 09:14

You FUN is much too inefficient for this task. As the first triangular number is above the 12,000th with a value of 75,000,000 and FUN runs through all these numbers ... the number of iterations to perform is almost

12000 * 75000000 / 2 = 450 * 10^9

This is clearly more than R's relatively slow for-loop can do in a reasonable time frame.

Instead, you could apply the divisors function from the numbers package that employs a prime factor decomposition. The following code need about 5-6 seconds (on my machine) to find the triangular number.

library(numbers)

t <- 0
system.time(
    for (i in 1:100000) {
        t <- t + i
        d <- length( divisors(t) )
        if (d > 500) {
            cat(i, t, d, '\n')
            break
        }
    }
)
## 12375 76576500 576 
##    user  system elapsed 
##   5.660   0.000   5.658

Instead of calculating the i-th triangular number, here i is added to the last triangular number. The time saving is minimal.

This is amazing. I need to look into `divisors()` and `primeFactors()` and `gmp` package and learn the details. — user101998, Sep 01 '17 at 00:32

score 1 · Answer 2 · answered Aug 31 '17 at 06:54

Here's my attempt:

library(gmp)
library(plyr)

get_all_factors <- function(n)
{
  prime_factor_tables <- lapply(
    setNames(n, n), 
    function(i)
    {
      if(i == 1) return(data.frame(x = 1L, freq = 1L))
      plyr::count(as.integer(gmp::factorize(i)))
    }
  )
  lapply(
    prime_factor_tables, 
    function(pft)
    {
      powers <- plyr::alply(pft, 1, function(row) row$x ^ seq.int(0L, row$freq))
      power_grid <- do.call(expand.grid, powers)
      sort(unique(apply(power_grid, 1, prod)))
    }
  )
}

for (i in 99691200:100000) {
  if (length(get_all_factors(i)[[1]])>500) print(paste(i, length(get_all_factors(i)[[1]])))
  if (i %% 100000 == 0) print(paste("-",i,"-"))
}

Let it run as long as you can be bothered...

Optimize simple r code for Project Euler 12

2 Answers2