Based on the answer by <PERSON> 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. 614
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
```python 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] ``` 603