Euler Problem 12

Back to overview.

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, \dots$

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?

Let's try brute forcing.

In [5]:
def get_divisors(n):
    return [i for i in range(1, int(n/2) + 1) if n % i == 0] + [n]

assert(get_divisors(28) == [1,2,4,7,14,28])

def triangle_number_generator_function():
    c = 0
    for i in range(1, 100000000000):
        c += i
        yield c

#ts = triangle_number_generator_function()
#for t in ts:
#    if len(get_divisors(t)) > 500:
#        print(t)
#        break

That failed miserably. We have to come up with something more efficient. I looked at my old Haskell solution which looks like this:

divisor_count' x = product . map (succ . length) . group $ prim_factors x

Add first I did not understand what it does at all. But the algorithm is as follows:

# get prime factors, for example for 28
2 * 2 * 7
# group the primes
[[2, 2], [7]]
# get the length of each group and increment it by one
[3, 2]
# get the product of all values
6

Honestly, I have no idea why this gives as the number of divisors. I will implement the solution and then try to come up with an explanation.

In [6]:
def is_prime(n, smaller_primes):
    for s in smaller_primes:
        if n % s == 0:
            return False
        if s * s > n:
            return True
    return True

def prime_generator_function():
    primes = [2, 3, 5, 7]
    for p in primes:
        yield p
    while True:
        p += 2
        if is_prime(p, primes):
            primes.append(p)
            yield p
            
def get_prime_factors(number):
    prime_generator = prime_generator_function()
    remainder = number
    factors = []
    for p in prime_generator:
        while remainder % p == 0:
            remainder /= p
            factors.append(p)
        if remainder == 1 or p * p > number:
            break
    if remainder != 1:
        factors.append(remainder)
    return factors

This are the prime numbers related functions we already now. Now we implement a group function, the product function we already know, and based on that the algorithm to get the number of divisors mentioned above.

In [7]:
def group(xs):
    from functools import reduce
    def f(xss, x):
        if xss and x in xss[-1]:
            xss[-1].append(x)
        else:
            xss.append([x])
        return xss
    return reduce(f, xs, [])

def product(xs):
    from functools import reduce
    from operator import mul
    return reduce(mul, xs, 1)

assert(group([1, 1, 3, 2, 2]) == [[1, 1], [3], [2, 2]])

def get_number_of_divisors(n):
    # get prime factors, for example for 28
    ps = get_prime_factors(n)
    # group the primes
    ps = group(ps)
    # get the length of each group and increment it by one
    ps = map(lambda x: len(x) + 1, ps)
    # get the product of all values
    return product(ps)
            
assert(get_number_of_divisors(28) == 6)

Now we are ready to do another brute force attempt.

In [8]:
ts = triangle_number_generator_function()
for t in ts:
    if get_number_of_divisors(t) > 500:
        print(t)
        break

76576500

Now the only question is why this crazy algorithm works. Okay, I got it with the help of this page. The problem is actually an instance of the multiplication principle for counting things. Each prime can be used to calculate (n + 1) other divisors. The incrementation by one is required because we can also choose to not use a certain prime. Of course, to get the potential combinations we have to get the product of the potential divisors for all primes. The web page explains it better. The only question is whether I came up with this algorithm myself back then.

In [ ]: