# Project Euler #12: Highly divisible triangular number

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
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())


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   ???????