Project Euler Problem 12 in Python

Question

I'm trying to do Project Euler Problem 12, which reads as:

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?

My code is as follows:

import math
from time import time
t = time()
def divisors(n):
    number_of_factors = 0
    for i in range(1, int(math.ceil(math.sqrt(n)))):
        if n % i == 0:
            number_of_factors +=2
        else:
            continue
    return number_of_factors

x=1
for y in range(2,1000000):
    x += y
    if divisors(x) >= 500:
        print x
        break
tt = time()-t
print tt

My code takes 3.8 seconds to run, and I know this can be done in around 0.2 seconds in Java. Any suggestions on how to make this run faster?

Answer 1

Improving project Euler solutions is usually done by improving the algorithms itself rather than just optimising the code.

First things first, your divisors() function is wrong and does not work for perfect squares. To fix this, you need to take care of the sqrt(n) case separately:

def divisors(n):
    number_of_factors = 0
    for i in xrange(1, int(math.ceil(math.sqrt(n)))+1):
        if n % i == 0:
            number_of_factors +=2
        if i*i==n:
            number_of_factors -=1
    return number_of_factors

Now we come to the main part of the code. The important observation to improve your code is that the nth triangular number is given by T(n) = n(n+1)/2. Now, n and n+1 are coprime. If a and b are coprime numbers, the number of divisors of a*b is just the product of number of divisors of a and b.

This suggests the following improvement for the main part of the code:

for n in xrange(1,1000000):
    Tn=(n*(n+1))/2
    if n%2==0:
        cnt=divisors(n/2)*divisors(n+1)
    else:
        cnt=divisors(n)*divisors((n+1)/2)
    if cnt >= 500:
        print Tn
        break

The improvement this provides is very significant. You can improve the performance further by modifying the divisor function to use the same technique:

def divisors(n,start):
    if n==1:
        return 1
    for i in xrange(st, int(math.ceil(math.sqrt(n)))+1):
        if n % i == 0:
            cnt=1
            while(n%i==0):
                n/=i
                cnt+=1
            return divisors(n,i+1)*cnt
    return 2

Essentially, we find p, the first prime factor of n. If p^k is the maximum power of p that divides n, (k+1)*divisors(n/p^k) is the number of divisors of n. start is just a starting point for checking prime divisors.

Overall, these modifications seem to reduce running time from ~5s to 0.15s on my laptop.

Answer 2

You can get a 28% speed-up if you use xrange instead of range, in Python 2 range used to create a full list consuming time and memory. (It has been fixed in Python 3).

---

Also I would like to suggest longer names for readibility sake.

triangle = 1
for natural in range(2,1000000):
    triangle += natural
    if divisors(triangle) >= 500:
        print triangle
        break

---

The following is useless and should be removed to reduce clutter (it also speeds the programme up a tiny bit

else:
    continue

---

Avoid magic numbers: DIVISORS_WANTED = 500 would be easier to change than a number buried inside the code.

Answer 3

Your divisors function is wrong. As noted by Raziman T V, it fails on perfect squares, but this is not the only problem.

You iterate over range(1, int(math.ceil(math.sqrt(n)))).

Let's see what it does on sa few values. To make it clear, I took Raziman T V's correction, and added a print (don't mind the range, I'm using Python 3):

def divisors(n):
    number_of_factors = 0
    for i in range(1, int(math.ceil(math.sqrt(n)))+1):
        if n % i == 0:
            number_of_factors +=2
            print(i, n//i)
        if i*i==n:
            number_of_factors -=1
    return number_of_factors

So here we go.

>>> divisors(9)
1 9
3 3
# Alright

>>> divisors(10)
1 10
2 5
#Alright

>>> divisors(11)
1 11
# Alright    

>>> divisors(12)
1 12
2 6
3 4
4 3
# Oops

You took 3 and 4 twice.

---

So what happened here?

As I said, you iterate until int(math.ceil(math.sqrt(n))). As a consequence, if n is divisible by both floor(sqrt(n)) and ceil(sqrt(n)), the iteration will continue, and (ceil(sqrt(n)), floor(sqrt(n))) will be found as another couple of divisors.

So if n is a product of two consecutive numbers p and p+1, it's likely that p and p+1 will be each counted twice.

Demonstration:

>>> divisors(25*26)
1 650
2 325
5 130
10 65
13 50
25 26
26 25

>>> divisors(45*46)
1 2070
2 1035
3 690
5 414
6 345
9 230
10 207
15 138
18 115
23 90
30 69
45 46
46 45

So your count will exceed by two the actual number of divisors. Fortunately for you, the result is the same with 500... But the logic is still off.

Therefore, you must iterate over floor(sqrt(n)) instead. I suggest that you change your xrange to xrange(1, int(math.sqrt(n))).