Project Euler Problem 45

Question

As a self-teaching Python beginner for almost 4 months, I have mostly been doing online challenges including Project Euler problems.

Problem 45 asks:

Triangle, pentagonal, and hexagonal numbers are generated by the following formulae:

$$ \begin{array}{lll} \textrm{Triangle} & T_n=n(n+1)/2 & 1, 3, 6, 10, 15, \ldots \\ \textrm{Pentagonal} & P_n=n(3n−1)/2 & 1, 5, 12, 22, 35, \ldots \\ \textrm{Hexagonal} & H_n=n(2n−1) & 1, 6, 15, 28, 45, \ldots \\ \end{array} $$

It can be verified that \$T_{285} = P_{165} = H_{143} = 40755\$.

Find the next triangle number that is also pentagonal and hexagonal.

I encountered this problem and was able to write the code without much problem although most of the solutions I write are mostly brute force. However I am not satisfied with the time the coding took to find the answer as it took 530.7 seconds to find the solution with this code:

import time

startTime = time.time()

limit = 1000000
triangle = []
pentagonal = []
hexagonal = []
triangle_number = []

class Shape:
    def __init__(self, term):
        self.term = term

    def triangle(self):
        return int(self.term * (self.term + 1) / 2)

    def pentagonal(self):
        return int(self.term * (3 * self.term -1) / 2)

    def hexagonal(self):
        return int(self.term * (2 * self.term - 1))

def generate():
    for terms in range(1, limit + 1):
        product = Shape(terms)

        triangle.append(product.triangle())
        pentagonal.append(product.pentagonal())
        hexagonal.append(product.hexagonal())

def main():
    generate()

    for _, terms in enumerate(hexagonal):
        if len(triangle_number) == 3:
            break
        elif terms in triangle and terms in pentagonal:
            triangle_number.append(terms)
            print(terms)

    print(triangle_number)
    print(time.time() - startTime, "seconds")

main()

I am aware that using the class is unecessary for this problem, but I want to include it as a practice of getting used to writing solutions that require classes.

Alternatively I have re-written another solution of this problem without any classes included but keeping the same style, however it still took 328.2 seconds.

So what suggestions are there to improve the code to run faster? I have tried to review other solutions from the solution page however I do not understand how it was simplified so much to make it run more efficient.

Answer 1

Arithmetic

Project Euler questions are meant to educate you about both mathematics and programming. It would be a good idea to understand what these triangular, pentagonal, and hexagonal numbers actually are, rather than blindly applying the given formulas.

One performance improvement would be to find a way to generate successive elements of each sequence without plugging \$n\$ into the formulas, which involve division. (Division tends to be slow. Floating-point division, using the / instead of the // operator, is even slower, and it also causes you to have to cast the result back to an int.)

If you take the formula for pentagonal numbers \$P_n\$, you can figure out another formula for the difference between successive elements in the sequence.

$$ \begin{align} P_n &= \frac{n(3n-1)}{2} = \frac{3n^2-n}{2} \\ P_{n+1} &= \frac{(n+1)(3(n+1)-1)}{2} = \frac{(n+1)(3n+2)}{2} = \frac{3n^2+5n+2}{2} \\ P_{n+1} - P_n &= \frac{(3n^2+5n+2)-(3n^2-n)}{2} = \frac{6n+2}{2} = 3n + 1 \end{align} $$

If you do the same for triangular, square, and hexagonal numbers, you'll find:

$$ \begin{align} T_{n+1} - T_n &= n+1 \\ P_{n+1} - P_n &= 3n+1 \\ H_{n+1} - H_n &= 4n+1 \end{align} $$

Considering that square numbers would be \$S_{n+1} - S_n = 2n+1\$, you can see a pattern to produce polygonal numbers in general.

Algorithm

Your strategy is to generate a million elements of each sequence and find the elements that exist in common.

First of all, one million is an arbitrary limit. You might need fewer than a million to find the next element in common (in which case you've wasted execution time), or you might need more than a million (in which case you would have to raise the limit and run your code again). It would be nice if your algorithm did not have to rely on a guess.

Secondly, the millionth hexagonal number is certainly going to be much larger than the millionth triangular number. There is no way that the millionth hexagonal number is going to coincide with anything, so that's wasted work.

Thirdly, you store the sequences as lists. Searching a list (e.g. terms in triangle) involves inspecting every element in that list (a so-called O(n) operation). Searching a set takes only O(1) time. Therefore, simply changing

triangle = []
pentagonal = []
hexagonal = []

and

        triangle.append(product.triangle())
        pentagonal.append(product.pentagonal())
        hexagonal.append(product.hexagonal())

to

triangle = set()
pentagonal = set()
hexagonal = set()

and

        triangle.add(product.triangle())
        pentagonal.add(product.pentagonal())
        hexagonal.add(product.hexagonal())

brings the execution time down from hundreds of seconds down to about 2 seconds. Better yet, your main() function could be simplified using the set intersection operator &:

def main():
    generate()
    triangle_number = triangle & pentagonal & hexagonal
    print(sorted(triangle_number))
    print(time.time() - startTime, "seconds")

Pythonicity

Be consistent with your naming. If you write "pentagonal" and "hexagonal", then use "triangular" rather than "triangle".

The Shape class shouldn't exist at all. It's just a very weird and cryptic way to call three functions that take a numerical parameter.

for _, terms in enumerate(hexagonal) is a nonsensical use of enumerate. If you're going to throw away the index anyway, why not just write for terms in hexagonal? And why is your iteration variable pluralized (terms rather than term)?

Your code would be much more expressive if you could say "give me the next pentagonal number". A good way to do that in Python is to define a generator, so that you can write next(pentagonal_numbers).

Suggested solution

from itertools import count

def polygonal_numbers(sides):
    result = 0
    for n in count():
        yield result
        result += (sides - 2) * n + 1

tt, pp, hh = polygonal_numbers(3), polygonal_numbers(5), polygonal_numbers(6)
t = p = 0 
for h in hh:
    while p < h: p = next(pp)
    while t < h: t = next(tt)
    if t == p == h > 40755:
        print(h)
        break

If you take into account that every hexagonal number is also a triangular number, you can ignore the triangular numbers altogether:

from itertools import count

def polygonal_numbers(sides):
    result = 0
    for n in count():
        yield result
        result += (sides - 2) * n + 1

pentagonal_numbers = polygonal_numbers(5)
p = 0
for h in polygonal_numbers(6):
    while p < h: p = next(pentagonal_numbers)
    if p == h > 40755:
        print(h)
        break

My last solution takes about 50 milliseconds to run on my machine.

Answer 2

Code

limit = 1000000
triangle = []
pentagonal = []
hexagonal = []
triangle_number = []

Global variables do not help readability.

What's the difference between triangle and triangle_number? Those names don't help me understand what they represent.

---

class Shape:
    def __init__(self, term):
        self.term = term

    def triangle(self):
        return int(self.term * (self.term + 1) / 2)

    def pentagonal(self):
        return int(self.term * (3 * self.term -1) / 2)

    def hexagonal(self):
        return int(self.term * (2 * self.term - 1))

A shape doesn't have a term: it has sides. Concretely, since we're talking regular shapes, it has two properties: the number of sides and the length of each side. This class doesn't make sense to me.

If you really want to practise structuring code with classes, the class should probably be a Solver.

---

    for _, terms in enumerate(hexagonal):
        if len(triangle_number) == 3:
            break
        elif terms in triangle and terms in pentagonal:
            triangle_number.append(terms)
            print(terms)

If you're testing x in ys then ys had better be a set, not a list, or you have to do a linear search.

---

Algorithm

The current algorithm can be summed up as such:

fix a large limit
generate `limit` terms in each of the sequences
for term in first_sequence
    if term in second_sequence and term in third_sequence:
        term is a candidate solution

The limit is guesswork, so it might be too small and not find the solution, or be too large and waste lots of time generating terms.

If you note that all of the sequences are strictly increasing, you can instead do a kind of merge:

while problem not solved:
    initialise each of the sequences at the first term
    if all sequences have the same current term:
        term is a candidate solution
    advance one of the sequences which has the smallest current term

---

Project Euler is more about maths than programming. Let's look again at those term formulae: $$T_n = \frac{n(n+1)}{2} \\ H_n = n(2n−1)$$ We can rewrite the latter as $$H_n = \frac{(2n−1)(2n)}{2}$$ Can you spot a major simplification which you can make to the search?

There are more sophisticated mathematical improvements, but this isn't the place. Check out the Project Euler discussion thread which you gain access to after solving the problem, and if you can distill out questions from that then ask them on our sister site math.stackexchange.com.

Answer 3

I want to include it as a practice of getting used to writing solutions that require classes.

No solution "requires" classes, although some situations are better represented with classes than with other techniques.

In this particular case, Shape doesn't really need to exist - as you've already identified. Since each of those methods only depends on term, you can simply have three functions all accepting one integer.

Some other things that can improve:

Int division

This:

return int(self.term * (self.term + 1) / 2)

can be

return self.term * (self.term + 1) // 2

Enumerate-to-/dev/null

You don't need to call enumerate - you aren't using the index. Just use for terms in hexagonal.

Answer 4

The main reason why your code is so slow is because your for loop in main spends most of the time checking for things that are logically impossible to be true. You have one million elements in every number group, every iteration of the for loop you are comparing one value to 2 million other values, when most of them can not be true.

for _, terms in enumerate(hexagonal):
    if len(triangle_number) == 3:
        break
    elif terms in triangle and terms in pentagonal:
        triangle_number.append(terms)
        print(terms)

With changing as little as possible and taking into account that as @Peter Taylor said all of the sequences are increasing. I was able to get the runtime of the program from 452.8s to 3.2s

t_iterator = 0
p_iterator = 0
for _, terms in enumerate(hexagonal):
    if len(triangle_number) == 3:
        break
    while (pentagonal[p_iterator] <= terms):
        if pentagonal[p_iterator] == terms:
            while (triangle[t_iterator] <= terms):
                if triangle[t_iterator] == terms:
                    print(terms)
                    triangle_number.append(terms)
                t_iterator += 1
        p_iterator += 1

This code I gave is still far from nice, it can be optimised further. You really should consider what the other answers have highlighted. It's not just that it does not look good and would be hard to understand. If you add another arbitrary 0 to the arbitrary "limit" you will probably run out of memory very quickly. You do not actually need to pre-calculate any of the values to solve the problem.

Answer 5

Python is a great choice of language for a challenge like this, mainly because of how easy it is to use sets. Basically, any challenge which states "find a number that matches these criteria" can be thought of as a set intersection problem. We want to find \$T \cap P \cap H\$, the intersection of triangular, pentagonal, and hexagonal numbers.

Depending on your Python version, you might not have access to the "walrus operator" :=. But for this case, it's quite handy.

upper_limit = 10000
hexs = set()
while len(hexs) < 2:
    tris = {n*(n+1)//2 for n in range(2, upper_limit)}
    pents = {v for n in range(upper_limit) if (v := n*(3*n-1)//2) in tris}
    hexs = {v for n in range(upper_limit) if (v:= n*(2*n-1)) in pents}
    upper_limit *= 10

print(hexs)

To find the number, we create larger and larger sets of triangular numbers, pentagonal numbers, and hexagonal numbers. However, we can filter out the pentagonal numbers which are not triangular, and the filter out the hexagonal numbers which are neither pentagonal nor triangular.

By using the optimizations presented in other answers, this could also be written as:

upper_limit = 10000
pents = set()
while len(pents) < 2:
    hexs = {n*(2*n-1) for n in range(2, upper_limit) if n*(2*n-1)}
    pents = {n*(3*n-1)//2 for n in range(upper_limit) if n*(3*n-1)//2 in hexs}
    upper_limit *= 10

print(pents)

The advantage of this approach is that it can easily be adapted for multiple different problems, and provides a lot of abstraction along with its performance.

Answer 6

Laziness is a programmer virtue.

These folks are spending more time on theory than doing it the lazy way would take in practice.

The big problem here is you're doing a bunch of unnecessary work, particularly around storage.

You're constantly appending to arrays for no apparent reason. You don't need the history of old numbers. Throw them away. Appending to arrays is expensive - it often has to reallocate the whole range into a new array then copy from the old one.

The question is "find the first 3 that's also a 5 and a 6."

Notice that you can rearrange those with no impact. That's the same question as "find the first 5 that's also a 3 and a 6," because they all go straight upwards, meaning the first of any will be the first of all.

Fine. You don't need to track anything. Iterate one of the columns. At each three, ask "is this one of the others?" If it is, ask "is this also the other other?" If it is, return success; otherwise keep marching.

So really, what you need is an efficient way to be able to ask "is this a foo-angular number?"

Right now you're doing that by building and searching a table. That's silly. Just reverse the math on the two "is it also" columns.

Almost all of this work can go.

Triangular and hexagonal are easy to reverse, so pentagonal is the one I'll keep in the original direction.

if triangular is "triangular x is (x * (x+1)) / 2,"

then in math, you have "n = x(x+1)/2".

solve for n, you get "x^2 + x - 2n = 0", or "x = (sqrt(1 + 8n) - 1)/2"

Therefore,

const triangle_from = x => (Math.sqrt(1+(x*8))-1)/2;

function is_triangular(x) {
  const pip = triangle_from(x);
  return pip == Math.floor(pip);
}

Now, you can actually throw that away; Euler is playing a trick on you, and I'm abusing that to show you how to do the work without actually doing it for you.

Why?

Because every hexagonal number is also a triangular number. By the time you've tested if it's hexagonal, that it's triangular is already in the bag. You can just skip that wholesale, as such.

You can retain the pentagonal in its existing notation, since we're using it to drive the bus. Also, TEST YOUR FUCKING CODE. I got this wrong because of order of operations differences between math and JS the first time. Just run it on 165 and see if it matches the problem description.

const to_pentagon = x => (x * ((3 * x)-1)) / 2;

Then it's just

function find_triple(cursor = 165, bail = 1000000) { // question says "after 40755," which is pentagonal #165
  while (true) {
    if (from_start >= bail) { throw new RangeException("exceeded safety cap"); }
    const current = to_pentagon(cursor);
    if (is_hexagonal(current)) { return cursor; }  
    ++cursor;
  }
}

Which, if you're feeling tricky, you could write as

function find_triple(cursor = 165, bail = 1000000) { // question says "after 40755," which is pentagonal #165
  while (true) {
    if (from_start >= bail) { throw new RangeException("exceeded safety cap"); }
    if (is_hexagonal(to_pentagon(cursor++))) { return --cursor; }  
  }
}

Answer 7

You can easily check that H(n) = T(2n - 1). So all hexagonal numbers are triangular numbers, meaning we can ignore the triangular numbers altogether.

To compute pentagonal numbers: Start with p = 1, dp = 4. To get the next pentagonal number, let p = p + dp, dp = dp + 3.

To compute hexagonal numbers: Start with h = 1, dh = 5. To get the next hexagonal number, let h = h + dh, dh = dh + 4.

So here is simple pseudo-code that finds all numbers that are both pentagonal and hexagonal (and therefore triangular):

Let p = 1, dp = 4
Let h = 1, dh = 5

Forever:
    If p = h then output p. 
    Let h = h + dh, dh = dh + 4
    Repeat
        Let p = p + dp, dp = dp + 3
    While p < h

Implemented in a low-level language like C or C++, this should find all 64 bit numbers in a few seconds. If you want to go further, change the code slightly: "diff" equals h - p from the previous code, but it will be a lot smaller, so you get much further with 64 bit integers.

Let dp = 4, np = 1
Let dh = 5, nh = 1
Let diff = 0

Forever:
    If diff = then output "np'th pentagonal number = nh'th pentagonal number". 
    Let diff = diff + dh, dh = dh + 4, nh = nh + 1
    Repeat
        Let diff = diff - dp, dp = dp + 3, np = np + 1
    While p < h

This outputs the indexes of pentagonal and hexagonal numbers that are equal. On an eight year old MacBook it takes less than six nanoseconds to examine each index, more than 10 billion indexes per minute, or about 150 trillion indexes per day. Hn with n = 1,042,188,953 is also pentagonal and triangular. There is another such Hn with n a bit over 201 billions; Hn is about 8.175 x 10^22. Finding another solution with this method will likely take a few days or weeks.

If you want to go further, solve p(m) = h(n), for integer n, calculating m as a real number. As n gets large, m as a function of n will be getting closer and closer to a linear function. Then you can use the GCD algorithm to quickly find values where m is close to an integer. You will need multi-precision arithmetic, but it should get you arbitrary large solutions very quickly. (If P(m) = H(n), then m ≈ n * sqrt (4/3) - (sqrt (1/12) - 1/6), with an error less than 4 / n^2, so starting with some n, you use the GCD algorithm to find the next n where m = n * sqrt (4/3) - (sqrt (1/12) - 1/6) is within 3.5 / n^2 of an integer).

Answer 8

As 200_success says, you can look at how the numbers are derived to generate the numbers one by one.

So rather than creating all the numbers and checking for an intersection, a simple algorithm is to look at a pentagon and hexagon number. If they are equal, you're done. If the pentagon number is larger than the hexagon number, then check whether the next hexagon number is equal. If the hexagon number is larger, then check the next pentagon number.

pentagon_index = 165
pentagon_number = 40755+3*165+1

hexagon_index = 143
hexagon_number = 40755+4*143+1

max_tries = 10**6

for i in range(max_tries):
    if pentagon_number < hexagon_number:
        pentagon_index +=1
        pentagon_number += 3*pentagon_index+1
    if pentagon_number > hexagon_number:
        hexagon_index +=1
        hexagon_number += 4*hexagon_index+1
    if pentagon_number == hexagon_number:
        break