Search code examples
pythonalgorithmpython-3.xprimesprime-factoring

How can I generate all possible divisor products for a number?


I'm struggling to implement an algorithm which would give me possible products for a number. For example, for N=24 these are:

24*1, 12*2, 8*3, 6*4, 4*3*2, 3*2*2*2

I've implemented a function which calculates the prime factors for a given number with their powers (e.g. 2^3 and 3^1 for N=24). But I can't figure out how to get the divisor combinations from the prime factors.

EDIT: Here's what I've tried:

def divisors(factors): # prime factors, e.g. [2,2,2,3] for 24
    yield list(factors)

    d = factors.pop()

    for i in range(len(factors)):
        m = [d*factors[i]] + factors[:i] + factors[i+1:]
        yield from divisors(m)

Solution

  • You don't say what size numbers you need this to work for, or whether speed is a concern, but here's a very simple unoptimized solution that should work just fine for small inputs (n smaller than 10**7, say).

    def products(n, min_divisor=2):
        """Generate expressions of n as a product of ints >= min_divisor."""
        if n == 1:
            yield []
        for divisor in range(min_divisor, n+1):
            if n % divisor == 0:
                for product in products(n // divisor, divisor):
                    yield product + [divisor]
    

    To explain the code, think about how you might do this methodically by hand. You could start with the products that contain a 2. If n is odd, there aren't any. If n is even, then we can do a simple recursion: find all possible decompositions of n // 2, and then output decomposition * 2 for each one. Once we've exhausted all the products containing a 2, we move on to products involving a 3. But here there's an additional complication: in the first step, we've already found all products involving a 2, so to avoid repeated solutions, we want to restrict to products in which every divisor is at least 3. So our recursive call needs to keep track of the minimum permitted divisor, min_divisor above. Finally, we need a base case: 1 is expressible as the empty product.

    And here's the output for n=24, including the 6*2*2 case that you missed:

    >>> for product in products(24):
    ...     print('*'.join(map(str, product)))
    ... 
    3*2*2*2
    6*2*2
    4*3*2
    12*2
    8*3
    6*4
    24
    

    This is less than satisfying, though: as other commenters have pointed out, it should be possible to compute the multiplicative partitions of n from the prime factorisation, and even just from the list of exponents in the prime factorisation, using the primes only when it's necessary to reconstruct the factors. Here's a variant of the above that works with an existing prime factorisation. There's still plenty of scope for efficiency improvements: in particular, the itertools.product call and subsequent filtering to ignore everything lexicographically smaller than min_exponents should be replaced with a custom iterator that starts from min_exponents. But this should serve as a starting point.

    import itertools
    
    def exponent_partitions(exponents, min_exponents):
        """Generate all vector partitions of 'exponents', each of whose
        entries is lexicographically at least 'min_exponents'."""
        if all(exponent == 0 for exponent in exponents):
            yield []
        else:
            for vector in itertools.product(*(range(v+1) for v in exponents)):
                if vector >= min_exponents:
                    remainder = tuple(x - y for x, y in zip(exponents, vector))
                    for partition in exponent_partitions(remainder, vector):
                        yield partition + [vector]
    
    
    def divisor_from_exponents(primes, exponent_vector):
        """Reconstruct divisor from the list of exponents."""
        divisor = 1
        for p, e in zip(primes, exponent_vector):
            divisor *= p**e
        return divisor
    
    
    def multiplicative_partitions(primes, exponents):
        """Generate all multiplication partitions of
        product(p**e for p, e in zip(primes, exponents))"""
        if len(exponents) == 0:
            # Corner case for partitions of 1.
            yield []
        else:
            initial_vector = (0,) * (len(exponents) - 1) + (1,)
            for partition in exponent_partitions(exponents, initial_vector):
                yield [divisor_from_exponents(primes, vector) for vector in partition]
    

    And the output for 24, again: we're writing 24 as 2**3 * 3**1, so the tuple of primes is (2, 3) and the corresponding tuple of exponents is (3, 1).

    >>> for product in multiplicative_partitions((2, 3), (3, 1)):
    ...     print('*'.join(map(str, product)))
    ... 
    2*2*2*3
    4*2*3
    8*3
    6*2*2
    12*2
    4*6
    24
    

    There's plenty of literature on generating and counting the multiplicative partitions of an integer. See for example the links from OEIS A001055, and the SAGE functions for computing vector partitions.