Project Euler #12: Highly divisible triangular number

Question

I want to receive advice on my code. It is too slow... The problem I am trying to solve is:

What is the value of the first triangle number to have over five hundred divisors?

import math
import numpy as np
import copy
from operator import add
memoi_table = {}
#check whether the x is prime or not
#to do a prime factorization
def prime_finder(x):
 cnt = 0
 a = math.sqrt(x)
 sqrt_number = int(math.floor(a))
 for i in range(2,sqrt_number+1):
 if x % i == 0:
 cnt = cnt +1
 break
 if cnt == 0:
 return True
 else:
 return False
#Do a prime_factorization
#Do a dynamic programing with memoization table
#Do a recursive method to cover the triangular number
#ex)Triangular number = 1,3,6,10,15,21,28
#When we treat 10, we couldn't utilize dynamic programming => Recursive method
def prime_factorization(x):
 global memoi_table
 memoi_table[x] = {}
 if prime_finder(x) == True:
 memoi_table[x][x] = 1
 for i in range(2,math.floor(math.sqrt(x)+1)):
 if (prime_finder(i) == True) and (x % i == 0):
 new_x = x / i
 if not new_x in memoi_table.keys():
 prime_factorization(new_x)
 #Recursive way
 if new_x in memoi_table.keys():
 memoi_table[x] = memoi_table[new_x].copy()
 #to shallow copy
 if i in memoi_table[x].keys():
 memoi_table[x][i] = memoi_table[x][i] + 1
 else:
 memoi_table[x][i] = 1
 return (x,memoi_table[x])
#ex) 20 => return (20,(2:2,5:1))
# To calculate the number of divisor based on prime factorization
def find_number_of_divisor(x):
 global memoi_table
 number_of_divisor = 1
 for key in memoi_table[x].keys():
 number_of_divisor = number_of_divisor * (memoi_table[x][key]+1)
 return number_of_divisor
#Let's do it
def main():
 global memoi_table
 n = 0
 while True:
 n = n + 1
 Triangular_number = n * (n + 1) / 2
 prime_factorization(Triangular_number)
 num_divisor = find_number_of_divisor(Triangular_number)
 if num_divisor > 500:
 return (Triangular_number,num_divisor)
 break
print(main())

• ErikR: codereview.stackexchange.com/search?q=euler+12+python

Answer 1

The Algorithm

Let's just look at your top level algorithm:

n = 0
while True:
 n = n + 1
 Triangular_number = n * (n + 1) / 2
 prime_factorization(Triangular_number)
 num_divisor = find_number_of_divisor(Triangular_number)

Factoring is hard. It's a difficult problem in the best of cases. You are taking a difficult problem and making it harder. You are starting with n * (n + 1) / 2 - that is already a partial factorization. You are then throwing out that information and factoring the whole number. It is much easier to factor 7 * 4 than it is to factor 28, but you are choosing to do the latter.

It would be much more efficient to do something like:

factor1 = prime_factorization(n/2 if n % 2 == 0 else n)
factor2 = prime_factorization(n+1 if n % 2 == 0 else (n+1)/2)
full_factorization = merge(factor1, factor2)

If you memoize prime_factorization, then you're only doing a little bit of work each iteration. Let's consider just going from the 49th triangular number to the 50th and assume we have a very good memoized prime factorization (more on this later). Factoring the full number:

(1275)
--> 3 * (425)
--> 3 * 5 * (85)
--> 3 * 5 * 5 * (17)

None of 1275, 425, or 85 were in the table. 17 was the first number we would have found. Consider the alternative, we are instead finding the prime factorization of:

(25)
--> memoized as 5*5 (already done when finding the 24th triangular number) 
(51)
--> 3 * (17)

And 17 is memoized already (from the 16th triangular number). So we did one division. That time difference adds up. The full solution I present later in this answer takes just under 1.0s (when run 10x) if I factor each factor independently, but 4.46s when I factor the whole number. That is a lot slower.

And when you find yourself doing a counting loop, you should really use itertools.count() so that the main loop becomes:

for n in itertools.count(start=1):
 factors1 = prime_factorization(n/2 if n % 2 == 0 else n)
 factors2 = prime_factorization(n+1 if n % 2 == 0 else (n+1)/2)
 full_factorization = merge_factors(factors1, factors2)
 if num_factors(full_factorization) >= 500:
 return n * (n+1) / 2

Memoization

Now, you're kind of memoizing, but then not really, since the first thing you do is wipe memoi_table[x]. If should definitely start with:

def prime_factorization(x):
 if x in memoi_table:
 return memoi_table[x]
 ...

In any event, it's simpler to just write a decorator:

def memoize(f):
 cache = {}
 def wrapper(*args):
 if not args in cache:
 cache[args] = f(*args)
 return cache[args]
 return wrapper
@memoize
def prime_factorization(x):
 ...

Prime Finder

If you do your loop correctly, you don't actually need prime_finder(i). What does that function do? It looks over all the numbers smaller than i... but we're already doing that in prime_factorization! You're just doubling work for no reason. The whole body should just be something like:

for i in range(2,math.floor(math.sqrt(x)+1)):
 if x%i == 0:
 old_factorization = prime_factorization(x / i).copy()
 old_factorization[i] += 1
 return old_factorization
# still here? we didn't find a factor, must be prime
res = collections.defaultdict(int)
res[x] = 1
return res

That's it.

Putting it all together

All we need is merging two dictionaries, which is an easy loop:

def merge_factors(f1, f2):
 res = f1.copy()
 for k, v in f2.items():
 res[k] += v
 return res 

And finding the number of factors, which is just a product:

def num_factors(factorization):
 product = 1 
 for v in factorization.values():
 product *= (v+1)
 return product

The times based on finding the first triangular number with at least X factors:

 mine OP
 5 0.0005s 0.0007s
 10 0.0009s 0.0041s
 25 0.0039s 0.0860s
100 0.0256s 6.0304s
500 0.9849s ???????