The sequence of triangle numbers is generated by adding the natural numbers. So the \(7^\text{th}\) 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, \dots\]
Let us list the factors of the first seven triangle numbers:
\[\begin{aligned} \mathbf 1 &\colon 1\\ \mathbf 3 &\colon 1,3\\ \mathbf 6 &\colon 1,2,3,6\\ \mathbf{10} &\colon 1,2,5,10\\ \mathbf{15} &\colon 1,3,5,15\\ \mathbf{21} &\colon 1,3,7,21\\ \mathbf{28} &\colon 1,2,4,7,14,28 \end{aligned}\]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?
Let \(d(n)\) be the number of divisors of \(n\). The first thing to note is that if \(n\) can be written as \(p_1^{a_1} p_2^{a_2} \cdots p_k^{a_k}\) for distinct primes \(p_1,p_2,\ldots,p_k\), then \(d(n) = (a_1+1)(a_2+1)\cdots(a_k+1)\). This is because each divisor of \(n\) can be written as \(p_1^{b_1}p_2^{b_2}\cdots p_k^{b_k}\), where \(0 \leq b_i \leq c_i\) for all \(i\); thus, there are \(a_1+1\) choices for \(b_1\), \(a_2+1\) choices for \(b_2\), and so forth.
The \(n\)-th triangular number is equal to \(\frac{n(n+1)}2\). Because \(n\) and \(n+1\) have no prime factors in common, we have
\[d \left( \frac{n(n+1)}2 \right) = \begin{cases} d(n) d(\frac{n+1}2) & \text{if } n \text{ is odd} \\ d(\frac n2) d(n+1) & \text{if } n \text{ is even} \end{cases} \]My code iterates along the positive integers, finds \(d(\frac n2)\) for even \(n\) and finds \(d(n)\) for odd \(n\) as it goes, and keeps track of the previous value of \(d\) as it goes. If we're examining \(n\) and the two most recent values of \(d\) multiply to more than five hundred, we know the answer is \(\frac{(n-1)n}2\).
from sympy import primerange
from itertools import count
def p12(a=500):
primes = list(primerange(20000)) # Should be enough
prev_d = 1
for n in count(2):
m = n # copy used for dividing
d = 1
if m % 2 == 0:
m //= 2
for p in primes:
t = 1
while m % p == 0:
t += 1
m //= p
d *= t
if m == 1:
break
else:
raise ValueError("Prime list not long enough")
if prev_d * d > a:
return (n-1) * n // 2
prev_d = d