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