Project Euler 378: Count ordered triplets whose triangular numbers have decreasing number of factors

Question

Project Euler problem 378 says:

Let \$T(n)\$ be the \$n\$th triangle number, so \$T(n) = {n (n+1) \over 2}\$.

Let \$\mathit{dT}(n)\$ be the number of divisors of \$T(n)\$.
E.g.: \$T(7) = 28\$ and \$\mathit{dT}(7) = 6\$.

Let \$\mathit{Tr}(n)\$ be the number of triples \$(i, j, k)\$ such that \1ドル ≤ i < j < k ≤ n\$ and \$\mathit{dT}(i) > \mathit{dT}(j) > \mathit{dT}(k)\$.
\$\mathit{Tr}(20) = 14\$, \$\mathit{Tr}(100) = 5772\$ and \$\mathit{Tr}(1000) = 11174776\$.

Find \$\mathit{Tr}(60,000円,000円)\$.
Give the last 18 digits of your answer.

This is my attempt at solving it:

def t(n):
 tri_num = (n * (n+1))//2
 return tri_num
#finding nth triangle numbers
def dt(n):
 count = 0
 for i in range(1,t(n)+1):
 if t(n)%i == 0:
 count += 1
 return count
#finding nth triangle numbers' total number of dividers
def factors(n):
 factors = []
 for i in range(1,t(n)+1):
 if t(n)%i == 0:
 factors.append(i)
 return factors
#finding nth triangle number's dividers
def tr(n):
 triplesnum = 0
 triples = [(i, j, k) for i in range(n) for j in range(n) for k in range(n) if 1 <= i < j < k <= n and dt(i) > dt(j) > dt(k)]
 for i in triples:
 triplesnum += 1
 return triplesnum
while True:
 n = int(input("N =???"))
 print("triples number",tr(n))
#solution

I am new to python; I only know at most the loops,list comprehension, functions, and some class. I am sure this problem can be solved more simply.

Although I see this code as a success considering my knowledge, it works too slowly for me after 3-digit values of n.

Is there an improvement I can make with my current knowledge? If not, what should I be learning about next?

Note: I know some of the code is not necessary, but I used them as a testing material in the further functions.

• Peter Taylor: Have you solved Project Euler #12?
• Gareth Rees: There is an important clue in the question. It says, "Give the last 18 digits of your answer." This suggests that the answer is more than \10ドル^{18}\$. And that means that any program that tries to count the triples one at a time will have to execute triplesnum += 1 more than \10ドル^{18}\$ times, taking hundreds of years. So this can't be the way to do it.
• kerx1998: @PeterTaylor i tried to solve #12 but it also took way too much time to end so i searched for a similiar but different problem including triangular numbers therefore came across #378 idk what is the shortcut in both of my code.
• Martin R: PE #378 (solved by 635 people so far) is a bit harder than PE #12 (solved by 186711), so you might want to start with the easier one. Also you'll find several Q&A's about PE #12 here on Code Review, given essential hints about divisors of triangle numbers.
• Toby Speight: Welcome to Code Review! The original title, which stated your concerns about the code, applies to too many questions on this site to be useful. The site standard is for the title to simply state the task accomplished by the code. I've edited it to help you; please check that I haven't misrepresented your code!

Answer 1

We can improve your code in many areas.

First, a baseline on my laptop. Run time for N=70 is 16.4 seconds.

t(n) is only called inside dt(n), so we can avoid the extra function call and remove t(n) and move the calculation inside of dt(n). While doing so, we can see the calculation is actually being performed \$t(n)+1\$ times! We should save the result of the calculation in a local variable tn:

def dt(n):
 tn = n * (n+1) // 2
 count = 0
 for i in range(1, tn+1):
 if tn % i == 0:
 count += 1
 return count

Time for N=70 drops to 4.5 seconds.

Counting factors from 1 to tn seems like it does a lot of extra work. We can automatically include 1 and tn as factors, so we could start our count a 2 (assuming tn>1), and only count the factors between 2 and tn-1. But the largest possible factor below tn is tn/2, so we can cut that range in half.

def dt(n):
 tn = n * (n+1) // 2
 if tn == 1:
 return 1
 count = 2
 for i in range(2,tn//2+1):
 if tn % i == 0:
 count += 1
 return count

Computes N=70 in 2.1 seconds.

Our factors come in pairs. 1 matches tn, 2 would match tn//2 if tn is even, 3 would match tn//3, and so on. We could count by 2 for every pair, and only count up to \$sqrt(tn)\$. Of course, we need to take into account the possibility of t(n) being a perfect square; t(8) = 36 = 6*6, so 6 should only count for 1 factor, not 2.

def dt(n):
 tn = n * (n+1) // 2
 if tn == 1:
 return 1
 sqrt_tn = int(math.sqrt(tn))
 count = 2
 for i in range(2, sqrt_tn+1):
 if tn % i == 0:
 count += 2
 if sqrt_tn * sqrt_tn == tn:
 count -= 1
 return count

List comprehension can reduce the loop a little bit, too.

def dt(n):
 tn = n * (n+1) // 2
 if tn == 1:
 return 1
 sqrt_tn = int(math.sqrt(tn))
 count = 2 + sum(2 for i in range(2, sqrt_tn+1) if tn % i == 0)
 if sqrt_tn * sqrt_tn == tn:
 count -= 1
 return count

Computes N=70 in 0.3 seconds.

How many times are we calling dt(n) for any given value of n? Does the same value get returned each time? If so, why recalculate it? functools.lru_cache can do the caching for us:

@functools.lru_cache(maxsize=None)
def dt(n):
 tn = n * (n+1) // 2
 if tn == 1:
 return 1
 sqrt_tn = int(math.sqrt(tn))
 count = 2 + sum(2 for i in range(2, sqrt_tn+1) if tn % i == 0)
 if sqrt_tn * sqrt_tn == tn:
 count -= 1
 return count

Computes N=70 in 0.09 seconds.

---

Let's turn our attention to tr(n).

triples = [(i, j, k) for i in range(n)
 for j in range(n)
 for k in range(n)
 if 1 <= i < j < k <= n and dt(i) > dt(j) > dt(k)]

If 1 <= i < j < k <= n, why start i at 0, j at 0, and k at 0? j should start at i+1 and k at j+1. And k should actually reach n, not n-1 which is a bug.

def tr(n):
 triplesnum = 0
 triples = [(i, j, k) for i in range(1, n-1)
 for j in range(i+1, n)
 for k in range(j+1, n+1) # n+1: Bug fix
 if 1 <= i < j < k <= n and dt(i) > dt(j) > dt(k)]
 for i in triples:
 triplesnum += 1
 return triplesnum

Computes N=70 in 0.07 seconds.

Why generate the tuple(i, j, k) when we only count its existence? Why generate an array only to count every element in it?

def tr(n):
 return sum(1 for i in range(1, n-1)
 for j in range(i+1, n)
 for k in range(j+1, n+1)
 if dt(i) > dt(j) > dt(k))

Computes N=70 in 0.05 seconds.

Looks like a 328x speed-up, using just loops, list comprehensions and functions ... just by looking to avoid repeated work and overhead ... oh, and the @lru_cache decorator.

Unfortunately, N=700 now takes 16 seconds. So these incremental speedups aren't going to take you to \$Tr(60 000 000)\$ anytime in the near future.

---

Oh, one last optimization. Why loop over all the possible values of k when dt(i) > dt(j) is False?

def tr(n):
 return sum(1 for i in range(1, n-1)
 for j in range(i+1, n) if dt(i) > dt(j)
 for k in range(j+1, n+1) if dt(j) > dt(k))

That brings N=700 down to 5.9 seconds. We now reach 16 seconds at N=1000. 60 million is still in the far, far future.

Answer 2

I am new to python; I only know at most the loops,list comprehension, functions, and some class. I am sure this problem can be solved more simply.

...

Is there an improvement I can make with my current knowledge? If not, what should I be learning about next?

Loops and arrays are all you really need to solve almost any Project Euler challenge. The only exceptions are the ones which also require reading an input file. The point of PE is not to push your knowledge of language features, but to push your mathematical skills and knowledge of data structures and algorithms.

I would advise you to solve at least 40 of the first 50 challenges before you start looking at challenges beyond 100, because the early challenges introduce key techniques which you will use many times in later challenges. Trying to solve #378 before #12 is simply doing things in the wrong order.

The most important of the key techniques is dynamic programming. So many PE problems can only be solved in a reasonable time by exploiting shared structure. In this challenge, for example, you use the "number of divisors" function \$d(n)\$. The structure of this function can be summed up as $$\begin{align} d(p^a) &= a+1&\qquad \textrm{when }p\textrm{ is prime} \\ d(mn) &=d(m)d(n)&\qquad \textrm{when }\gcd(m,n)=1 \end{align}$$ Exploiting that structure by adapting the sieve of Eratosphenes will allow you to calculate an array containing dt(n) for every \1ドル \le n \le 60000000\$ in less than a minute.

Similarly, exploiting the structure of \1ドル ≤ i < j < k ≤ n\$ and \$\mathit{dT}(i) > \mathit{dT}(j) > \mathit{dT}(k)\$ will allow you to reduce the current \$O(n^3)\$ runtime of that accumulation phase to \$O(nw)\$ where \$w = \max \{dT(i) : 1 \le i \le n\}\$. That's enough to solve the problem in about 40 minutes, using only basic data structures.

To solve the problem in under two minutes you will need an advanced data structure, but that's more than enough spoiler. Solve the problem at all and you'll get access to the forum thread where people explain their solutions.