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?

Mathematics Behind a Solution

If $T(n)$ is the n-th triangular number, then we know that

$T(n)=1+2+\cdots+n=\sum_{i=1}^n=\frac{n(n+1)}{2}.$

Let’s assume that we know the prime factorizations of both $n$ and $n+1$, and we’ll write

$n = p_1^{e_1}p_2^{e_2}\cdots p_s^{e^s}\quad\text{and}\quad n+1 = q_1^{f_1}q_2^{f_2}\cdots q_t^{f_t}$

Notice that $n$ and $n+1$ cannot share any prime factors (as they are consecutive numbers), and so we know that all the $p_i$ and $q_i$ are distinct primes. Also, one and only one of $n$ and $n+1$ are divisible by 2. The exponents $e_i$ and $f_i$ are therefore the only things we really need to consider in order to determine the number of divisors of $T(n)$. The fact that $T(n)=(n(n+1))/2$ means that we’ll need to neglect a single power of two in the factorization of $n$ or $n+1$ (remember, only one is even). Let’s assume, without loss of generality, that $n$ is even and that $p_1=2$. Then, some quick combinatorics tell us that the total number of factors of $T(n)$ will be

$(e_1)(e_2+1)\cdots(e_s+1)(f_1+1)(f_2+1)\cdots(f_t+1)=e_1\prod_{i=2}^s(e_i+1)\prod_{i=1}^t(f_i+1)$

Even better, as we’re increasing our triangular numbers looking for the first one to satisfy the property in question, we only need to calculate the factorization of $n+1$ (as we already know that factorization of $n$, as it was previously $n+1$). This should decrease the runtime substantially.

Python Solution

The mathematics described above can be coded in Python as follows.

import time   def num_divisors(n): if n % 2 == 0: n = n/2 divisors = 1 count = 0 while n % 2 == 0: count += 1 n = n/2 divisors = divisors * (count + 1) p = 3 while n != 1: count = 0 while n % p == 0: count += 1 n = n/p divisors = divisors * (count + 1) p += 2 return divisors   def find_triangular_index(factor_limit): n = 1 lnum, rnum = num_divisors(n), num_divisors(n+1) while lnum * rnum < 500: n += 1 lnum, rnum = rnum, num_divisors(n+1) return n   start = time.time() index = find_triangular_index(500) triangle = (index * (index + 1)) / 2 elapsed = (time.time() - start)   print "result %s returned in %s seconds." % (triangle,elapsed)

When executed, we get the following result.

result 76576500 returned in 3.48622322083 seconds.

Cython Solution

We can recode the Python to Cython and achieve a significant speedup.

%cython   import time   cdef num_divisors(unsigned long int n): if n % 2 == 0: n = n/2 cdef unsigned long divisors = 1 cdef unsigned int count = 0 while n % 2 == 0: count += 1 n = n/2 divisors = divisors * (count + 1) cdef unsigned int p = 3 while n != 1: count = 0 while n % p == 0: count += 1 n = n/p divisors = divisors * (count + 1) p += 2 return divisors   cdef find_triangular_index(unsigned int factor_limit): cdef unsigned long n = 1 cdef unsigned long lnum, rnum lnum, rnum = num_divisors(n), num_divisors(n+1) while lnum * rnum < 500: n += 1 lnum, rnum = rnum, num_divisors(n+1) return n   start = time.time() index = find_triangular_index(500) triangle = (index * (index + 1)) / 2 elapsed = (time.time() - start)   print "result %s returned in %s seconds." % (triangle,elapsed)

When executed, we obtain the following.

result 76576500 returned in 0.122308969498 seconds.

Thus, the Cython executes roughly 28.5 times faster than the Python version.

Sage Solution

We could also solve the problem using Sage’s built-in functions.

import time   start = time.time()   i = 1 triangle = 1 while True: if len(triangle.divisors()) > 500: break i += 1 triangle += i   elapsed = (time.time() - start) print triangle print "found %s in %s seconds" % (i,elapsed)

This yields the following.

found 76576500 in 0.597285985947 seconds

Comments are closed