For a given N, I am trying to find every positive integers a and b such that N = a*b.
I start decomposing into prime factors using sympy.ntheory.factorint, it gives me a dict factor -> exponent.
I have this code already, but I don't want to get duplicates (a and b play the same role):
import itertools
from sympy.ntheory import factorint
def find_decompositions(n):
prime_factors = factorint(n)
cut_points = {f: [i for i in range(1+e)] for f, e in prime_factors.items()}
cuts = itertools.product(*cut_points.values())
decompositions = [((a := np.prod([f**e for f, e in zip(prime_factors, cut)])), n//a) for cut in cuts]
return decompositions
Example:
In [235]: find_decompositions(12)
Out[235]: [(1, 12), (3, 4), (2, 6), (6, 2), (4, 3), (12, 1)]
What I would like to get instead:
Out[235]: [(1, 12), (3, 4), (2, 6)]
I tried to reduce halve the range in cut_points with range extends such as e//2, 1 + e//2, (1+e)//2, 1 + (1+e)//2. None of it ended up working.
A simple solution is obviously to compute the same and return:
decompositions[:(len(decompositions)+1)//2]
but I am looking for an eventual solution that reduces the number of computations instead.
You're using module sympy, which already has a divisors function:
from sympy import divisors
print(divisors(12))
# [1, 2, 3, 4, 6, 12]
def find_decompositions(n):
divs = divisors(n)
half = (len(divs) + 1) // 2 # +1 because perfect squares have odd number of divisors
return list(zip(divs[:half], divs[::-1]))
print(find_decompositions(12))
# [(1, 12), (2, 6), (3, 4)]
print(find_decompositions(100))
# [(1, 100), (2, 50), (4, 25), (5, 20), (10, 10)]