Project Euler Problem 12: Highly divisible triangular number

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, ...

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?

Version 1: Brute force

In [1]:
from six.moves import map, range, reduce
In [2]:
def triangle_num(n):
    "The nth triangle number"
    return n*(n+1)>>1
In [3]:
triangle_num(7)
Out[3]:
28
In [4]:
list(map(triangle_num, range(1, 11)))
Out[4]:
[1, 3, 6, 10, 15, 21, 28, 36, 45, 55]

We already implemented a function that finds all prime factors of a number which may be useful if we are required to optimize. For now, let's use something simple like this. After all,

Premature optimization is the root of all evil

  • Donald Knuth
In [5]:
factors = lambda n: list(filter(lambda x: n % x == 0, range(1, n+1)))
In [6]:
{n:factors(n) for n in map(triangle_num, range(1, 8))}
Out[6]:
{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]}

Function composition is going to be useful for this problem. Of course, Python is not a functional programming language, but we can nonetheless hack together something that closely resembles it.

In [7]:
compose = lambda *fns: reduce(lambda f, g: lambda *args, **kwargs: f(g(*args, **kwargs)), fns)

Here is a "semi-quine" to test our function composition.

In [8]:
compose(lambda x: x**2, len)('should be 169')
Out[8]:
169
In [9]:
map(compose(factors, triangle_num), range(1, 11))
Out[9]:
<map at 0x11046ab38>

Now we're ready to go

In [10]:
from itertools import count, islice
In [11]:
def first(pred, iterable):
    for x in iterable:
        if pred(x):
            return x

First, a quick composition with the generator version of map and count/islice instead of the usual range, since we're going to apply the higher-order function first later, instead of islice.

In [12]:
list(islice(map(compose(len, factors, triangle_num), count(1)), 10))
Out[12]:
[1, 2, 4, 4, 4, 4, 6, 9, 6, 4]

We can see that the first triangle number to have over 5 divisors has 6 divisors, but we don't know what the triangle number is. We need to do a bit more work to retain that information in our function composition pipeline.

In [13]:
first(lambda n: n >= 5, map(compose(len, factors, triangle_num), count(1)))
Out[13]:
6

Excellent. Now we have a pair of the triangle number, and its number of factors.

In [14]:
list(islice(map(compose(lambda n: (n, compose(len, factors)(n)), triangle_num), count(1)), 10))
Out[14]:
[(1, 1),
 (3, 2),
 (6, 4),
 (10, 4),
 (15, 4),
 (21, 4),
 (28, 6),
 (36, 9),
 (45, 6),
 (55, 4)]

Instead of using islice, we just need to use our higher-order function first with a predicate. What does this function do? The name speaks for itself, return the first item in an iterable that satisfies the predicate.

In [15]:
first(lambda l: l[1] >= 5, map(compose(lambda n: (n, compose(len, factors)(n)), triangle_num), count(1)))
Out[15]:
(28, 6)

Hold on to your hats...

In [16]:
# (lambda l: l[1] >= 5*10**2, imap(compose(lambda n: (n, compose(len, factors)(n)), triangle_num), count(1)))

Version 2: The Fundamental Theorem of Arithmetic

By the fundamental theorem of arithmetic, every integer $n$ has a unique prime factorization, say, $n = \prod_{i} p_i^{e_i}$ and it follows that $n$ has $\prod_{i} (e_i+1)$ factors. To see this, observe that each prime factor $p_i$ of $n$ can be used 0 to $e_i$ times to form a factor of $n$. This is where our prime_factors function comes in handy.

In [17]:
%load_ext autoreload
%autoreload 2
In [18]:
from common.utils import prime_factors
In [19]:
prime_factors(28)
Out[19]:
Counter({2: 2, 7.0: 1})

First, let's implement prod, which is the multiplication analog to sum, with the help of the built-in reduce function.

In [20]:
prod = lambda iterable: reduce(lambda x, y: x*y, iterable, 1)
In [21]:
prod([])
Out[21]:
1

Now, we need a way to add 1 to each item of the Counter which represents an integer's prime factors. We could do this for example

In [22]:
from collections import Counter
c = Counter({2: 2, 7: 1})
for i in c:
    c[i] += 1
c
Out[22]:
Counter({2: 3, 7: 2})

But that's not very Pythonic, nor does it utilize the power of Counter very well. Let's do it right - by making use of the fact that Counter objects supports addition.

In [23]:
c = Counter({2: 2, 7: 1})
d = Counter(c.keys())
d
Out[23]:
Counter({2: 1, 7: 1})
In [24]:
c + d
Out[24]:
Counter({2: 3, 7: 2})

Now we're ready to implement count_factors.

In [25]:
def count_factors(n):
    c = prime_factors(n)
    d = Counter(c.keys())
    return prod((c+d).values()) 
In [26]:
count_factors(10000)
Out[26]:
25
In [27]:
first(lambda l: l[1] >= 5, map(compose(lambda n: (n, count_factors(n)), triangle_num), count(2)))
Out[27]:
(28, 6)
In [28]:
first(lambda l: l[1] >= 50, map(compose(lambda n: (n, count_factors(n)), triangle_num), count(2)))
Out[28]:
(25200, 90)
In [29]:
first(lambda l: l[1] >= 500, map(compose(lambda n: (n, count_factors(n)), triangle_num), count(2)))
Out[29]:
(76576500, 576)

Version 3: Pretty basic number theory

Denote the $n$th triangle by $T_n$ and recall that $T_n = \frac{n(n+1)}{2}$.

Observe that $n$ and $n+1$ are coprime. That is, $\gcd{(n, n+1)} = 1$. Note that one of $n, n+1$ must be even and that in general $\gcd{(a, 2b)} = 1 \Rightarrow \gcd{(a, b)} = 1$ for integers $a, b$.

Say $n = 2k$. Then $T_n = \frac{n(n+1)}{2} = k(2k+1)$. Then, $\gcd{(n, n+1)} = 1 \Rightarrow \gcd{(\frac{n}{2}, n+1)} = \gcd{(k, 2k+1)} = 1$ - the only common positive factor of $k$ and $2k+1$ is $1$, so we need only find their respective factors and combine them together.

The same goes for the case where $n = 2k+1$. We have $T_n = \frac{n(n+1)}{2} =(k+1)(2k+1)$ and $\gcd{(n, n+1)} = 1\Rightarrow \gcd{(n, \frac{n+1}{2})} = \gcd{(2k+1, k+1)}=1$, so we factorize $2k+1$ and $k+1$.

Now that we have decomposed the problem to factorizing 2 integers, we essentially have the overlapping subproblem of factorizing $2k+1$. So we can even use dynamic programming to further optimize this solution.

In [30]:
for n in count(2):
    k = n // 2
    T_n = triangle_num(n)
    if n % 2:
        a, b = 2*k+1, k+1
    else:
        a, b = k, 2*k+1
    T_n_factors = count_factors(a)*count_factors(b)
    if T_n_factors >= 500:
        print(T_n, T_n_factors)
        break
76576500 576
In [31]:
for k in count(1):
    a = count_factors(k)
    b = count_factors(2*k+1)
    c = count_factors(k+1)
    n_even = triangle_num(2*k)
    n_odd = triangle_num(2*k+1)
    if a*b >= 500:
        print(n_even)
        break
    elif b*c >= 500:
        print(n_odd)
        break
76576500

Comments

Comments powered by Disqus