Return to Question

Update

Based on the answer by D.W. I was able to implement an alternative algorithm. This effectively iterates over all permutations of $D$-length factorizations of $c$. Instead of utilizing the prime factorization to explicitly enumerate all the factorizations, I increment one element of $b$ to the next multiple of a divisor of $c$, then divide out that divisor from $c$ and repeat recursively. This effectively explores all $D$-length factorizations of $c$, but has the benefit of being able to ignore a large subset of the search space by abandoning a branch when the objective product becomes larger or equal to the smallest feasible objective found so far. Combined with some simple memoization, the resulting algorithm appears to be much faster on average than my original algorithm.

I've tested this on lots of random cases and tried to construct some pathological cases as well. I think there might be some pathological cases where this is slower than the original brute force, probably when there are tons of prime factors in $c$, but it seems to be much better on average.

My Python 3 implementation of the updated algorithm is

from functools import reduce
from operator import mul
from copy import deepcopy

def prod(x):
    return reduce(mul, x, 1)

def argsort(x, reverse=False):
    return sorted(range(len(x)), key=lambda idx: x[idx], reverse=reverse)

def divisors(v):
    """ does not include 1 """
    d = {v} if v > 1 else set()
    for n in range(2, int(v**0.5) + 1):
        if v % n == 0:
            d.add(n)
            d.add(v // n)
    return d

def update_memory(b, c_rem, memory):
    tuple_m = tuple(b + [c_rem])
    if tuple_m in memory:
        return False
    memory.add(tuple_m)
    return True

def dfs(a, b, c, c_rem, memory, p_best=float('inf'), b_best=None):
    ab = [ai + bi for ai, bi in zip(a, b)]
    p = prod(ab)
    if p >= p_best:
        return p_best, b_best
    elif p % c == 0:
        return p, deepcopy(b)

    dc = divisors(c_rem)
    for i in argsort(ab):
        for d in dc:
            db = (d - ab[i]) % d
            b[i] += db
            if update_memory(b, c_rem // d, memory):
                p_best, b_best = dfs(a, b, c, c_rem // d, memory, p_best, b_best)
            b[i] -= db

    return p_best, b_best

def solve(a, c):
    b = [0 for _ in range(len(a))]
    result = dfs(a, b, c, c, set())
    p_best, b_best = result
    return b_best


# a few test cases
assert solve([2, 3], 9) == [1, 0]
assert solve([2, 8], 6) == [0, 1]
assert solve([13, 17, 25], 8) in [[1, 1, 1], [0, 1, 3]]
assert solve([5, 13, 19], 6) == [0, 1, 2]

Based on the answer by D.W. I was able to implement an alternative algorithm. This effectively iterates over all permutations of $D$-length factorizations of $c$. Instead of utilizing the prime factorization to explicitly enumerate all the factorizations, I increment one element of $b$ to the next multiple of a divisor of $c$, then divide out that divisor from $c$ and repeat recursively. This effectively explores all $D$-length factorizations of $c$, but has the benefit of being able to be ignore a large subset of the search space by abandoning a branch when the objective product becomes larger or equal to the smallest feasible objective found so far. Combined with some simple memoization, the resulting algorithm appears to be much faster on average than my original algorithm.

Based on the answer by D.W. I was able to implement an alternative algorithm. This effectively iterates over all permutations of $D$-length factorizations of $c$. Instead of utilizing the prime factorization to explicitly enumerate all the factorizations, I increment one element of $b$ to the next multiple of a divisor of $c$, then divide out that divisor from $c$ and repeat recursively. This effectively explores all $D$-length factorizations of $c$, but has the benefit of being able to ignore a large subset of the search space by abandoning a branch when the objective product becomes larger or equal to the smallest feasible objective found so far. Combined with some simple memoization, the resulting algorithm appears to be much faster on average than my original algorithm.

from functools import reduce
from operator import mul
from copy import deepcopy

def prod(x):
    return reduce(mul, x, 1)

def argsort(x, reverse=False):
    return sorted(range(len(x)), key=lambda idx: x[idx], reverse=reverse)

def divisors(v):
    """ does not include 1 """
    d = {v}
    for n in range(2, int(v**0.5) + 1):
        if v % n == 0:
            d.add(n)
            d.add(v // n)
    return d

def update_memory(b, c_rem, memory):
    tuple_m = tuple(b + [c_rem])
    if tuple_m in memory:
        return False
    memory.add(tuple_m)
    return True

def dfs(a, b, c, c_rem, memory, p_best=float('inf'), b_best=None):
    ab = [ai + bi for ai, bi in zip(a, b)]
    p = prod(ab)
    if p >= p_best:
        return p_best, b_best
    elif p % c == 0:
        return p, deepcopy(b)

    dc = divisors(c_rem)
    for i in argsort(ab):
        for d in dc:
            db = (d - ab[i]) % d
            b[i] += db
            if update_memory(b, c_rem // d, memory):
                p_best, b_best = dfs(a, b, c, c_rem // d, memory, p_best, b_best)
            b[i] -= db

    return p_best, b_best

def solve(a, c):
    b = [0 for _ in range(len(a))]
    result = dfs(a, b, c, c, set())
    p_best, b_best = result
    return b_best


# a few test cases
assert solve([2, 3], 9) == [1, 0]
assert solve([2, 8], 6) == [0, 1]
assert solve([13, 17, 25], 8) in [[1, 1, 1], [0, 1, 3]]
assert solve([5, 13, 19], 6) == [0, 1, 2]

from functools import reduce
from operator import mul
from copy import deepcopy

def prod(x):
    return reduce(mul, x, 1)

def argsort(x, reverse=False):
    return sorted(range(len(x)), key=lambda idx: x[idx], reverse=reverse)

def divisors(v):
    """ does not include 1 """
    d = {v} if v > 1 else set()
    for n in range(2, int(v**0.5) + 1):
        if v % n == 0:
            d.add(n)
            d.add(v // n)
    return d

def update_memory(b, c_rem, memory):
    tuple_m = tuple(b + [c_rem])
    if tuple_m in memory:
        return False
    memory.add(tuple_m)
    return True

def dfs(a, b, c, c_rem, memory, p_best=float('inf'), b_best=None):
    ab = [ai + bi for ai, bi in zip(a, b)]
    p = prod(ab)
    if p >= p_best:
        return p_best, b_best
    elif p % c == 0:
        return p, deepcopy(b)

    dc = divisors(c_rem)
    for i in argsort(ab):
        for d in dc:
            db = (d - ab[i]) % d
            b[i] += db
            if update_memory(b, c_rem // d, memory):
                p_best, b_best = dfs(a, b, c, c_rem // d, memory, p_best, b_best)
            b[i] -= db

    return p_best, b_best

def solve(a, c):
    b = [0 for _ in range(len(a))]
    result = dfs(a, b, c, c, set())
    p_best, b_best = result
    return b_best


# a few test cases
assert solve([2, 3], 9) == [1, 0]
assert solve([2, 8], 6) == [0, 1]
assert solve([13, 17, 25], 8) in [[1, 1, 1], [0, 1, 3]]
assert solve([5, 13, 19], 6) == [0, 1, 2]