This is a notebook for Project Euler #12.

The sequence of triangle numbers is generated by adding the natural numbers. So the 7th triangle number would be \[1 + 2 + 3 + 4 + 5 + 6 + 7 = 28.\] The first ten terms would be:

\[1, 3, 6, 10, 15, 21, 28, 36, 45, 55,\ldots\]

Let us list the factors of the first seven triangle numbers:

1: 1 3: 1,3 6: 1,2,3,6 10: 1,2,5,10 15: 1,3,5,15 21: 1,3,7,21 28: 1,2,4,7,14,28

We can see that 28 is the first triangle number to have over five divisors.

What is the value of the first triangle number to have over five hundred divisors?

Solution: We first note the following important remark

Result: Given a number \(n\), if we write it in terms of prime factorization: \[n=p_1^{k_1}\cdots p_m^{k_m},\] then the number divisors is simply \((k_1 +1)\times \cdots \times (k_m +1)\)

Proof: Any divisor of \(n\) is of the form \(p_1^{j_1}\cdots p_m^{j_m}\) where \[0\le j_1 \le k_1,\ldots,0 \le j_m \le k_m. \] Thus, we need to define a function that gives us the vector of exponents of the prime divisors. Here is the function that for a given integer, it outputs the factorization:

# The following will give the necessarily prime first divisor 
# of a number:

first.divisor <- function(n){
  # because we start out testing at 2 and go to 
  # square root of n, we need to manually give the answers 
  # for 2 and 3, both of which have square roots < 2:
  if (n==2) return(2)
  if (n==3) return(3)  
  top.range <- as.integer(sqrt(n))
  for (i in 2L:top.range){
    if (n %% i == 0) return(i)
  }
  return(n)
}

# We use the previous function to define the following:

prime.factorization <- function(n) {
  if (n==2) return(2)
  pd <- c()
  while (n > 1){
    fd <- first.divisor(n)
    pd <- c(pd,fd)
    n <- n %/% fd    # integer division
  }
  return(pd)
}

Let’s check this

print(prime.factorization(2205))
## [1] 3 3 5 7 7

Now we want, for a given integer, the multiplicities of the prime factors. Here is a function that will give the prime divisors without repetition:

members <- function(a){          # a is a vector of values
  m <- c()    # start with the empty set
  for (i in a){
    # if i is a member that hasn't been listed, then add it...
    if (!is.element(i,m)) m <- c(m,i)
  }
  return(m)
}

Let’s try it:

print(members(c(3,3,5,7,7,13)))
## [1]  3  5  7 13

Next we need a function that takes a vector of number and gives the vector of multiplicities of the elements, thus for the vector, [3, 3, 5, 7, 7, 13], it will give [2, 1, 2, 1], which is the number of times that 3, 5, 7, and 13 appear, respectively.

# a will be a vector of entries with repetitions
# like [3,3,5,7,7,13]

multip <- function(a) {
  b <- members(a)
  l.b <- length(b)
  m <- rep(1,l.b)     # initialize multip with 1's
  for (i in 1:l.b) {
    # the function, which(5==a), gives us the indices of 
    # where 5 is the element of the vector a. We can use
    # the length of this to give the multiplicity of 5. 
    m[i] <- length(which(b[i]==a))
  }
  return(m)
}

Let’s test this out…

multip(c(3,3,5,7,7,13))
## [1] 2 1 2 1

Now the number of divisors is given by the product of the multiplicities of the prime factorization plus 1. For example, the number 12 has prime factorization \(2\times 2\times 3\). the multiplicities are 2 and 1. The product of multiplicities plus 1 is \(3\times 2 = 6\) and indeed the divisors are \(1,2,3,4,6,12\).

Therefore finally we can define a function that gives the number of divisors (whew!):

n.of.divisors <- function(n){
  p.f <- prime.factorization(n)
  # m <- members(p.f)
  # print(m)
  mul <- multip(p.f)
  # print(mul)
  return(prod(mul + 1))
  #
  # note we could have written this in terms of 
  # one long function. we spread it out and used 
  # unnecessary, intermediate variables for clarity.
}

Test this. We consider the number of divisors of 10 and 12. The divisors are 1, 2, 5, 10 and 1, 2, 3, 4, 6, 12, respectively. Hence the number of divisors should be 4 and 6 respectively:

n.of.divisors(10)
## [1] 4
n.of.divisors(12)
## [1] 6

Let’s look at number of divisors for the first 20 triangular numbers:

for (i in 1:20) {
  n.triangle <- ( i*(i+1) )%/% 2
  print(c(n.triangle,n.of.divisors(n.triangle)))
}
## [1] 1 1
## [1] 3 2
## [1] 6 4
## [1] 10  4
## [1] 15  4
## [1] 21  4
## [1] 28  6
## [1] 36  9
## [1] 45  6
## [1] 55  4
## [1] 66  8
## [1] 78  8
## [1] 91  4
## [1] 105   8
## [1] 120  16
## [1] 136   8
## [1] 153   6
## [1] 171   6
## [1] 190   8
## [1] 210  16

Now let’s find the answer to the question! Want the first triangular number with the number of divisors to be more than 500.

nd <- 1
i <- 1
while (nd<500){
  tri <- (i*(i+1)) %/% 2
  nd <- n.of.divisors(tri)
  if (nd>500) {
    print(c(tri,nd))
    break
  }
  i <- i +1
  
}
## [1] 76576500      576

And the number, 76,576,500 (the 12375th triangular number) is indeed the first triangular number with more than 500 divisors.