2
\$\begingroup\$

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())
\$\endgroup\$
2
\$\begingroup\$

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   ???????
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.