Find a vector of non-negative integers $b$ that minimizes $\prod_{i = 1}^{D}\left(a_i + b_i\right)$ such that the product is a multiple of $c$

I'm trying to come up an efficient algorithm that, given a list of positive integers $$a = \left(a_1, \ldots, a_D\right)$$ and positive integer $$c$$, finds a list of non-negative integers $$b = (b_1, \ldots, b_D)$$ that minimizes $$\prod_{i = 1}^{D}\left(a_i + b_i\right)$$ such that $$\prod_{i = 1}^{D}\left(a_i + b_i\right)$$ is a multiple of $$c$$.

The brute force search I came up with is

1. Let $$t$$ be the smallest multiple of $$c$$ that is $$\ge \prod_{i = 1}^{D} a_i$$.
2. Use depth-first-search to search for values $$\hat{b}$$, starting at $$\mathbf{0}$$ and incrementing one element of $$\hat{b}$$ at a time until $$\prod_{i = 1}^{D} \left(a_i + \hat{b}_i\right) \geq t$$.
3. If we found $$\hat{b}$$ such that $$\prod_{i = 1}^{D} \left(a_i + \hat{b}_i\right) = t$$ then $$\hat{b}$$ is the optimal solution. Otherwise increment $$t$$ by $$c$$ and go back to step 2.

The above algorithm does work, but if $$D$$ or $$c$$ are too large then it will potentially take a very very long time. I'm wondering if this maps to any well known algorithm or if there's a more efficient solution. I'm considering that the prime factors of $$c$$ and $$a$$ could play a large role in reducing the search space but I can't quite figure it out.

In case anyone wants to play with this, a Python 3 implementation of the brute force algorithm described above is

from functools import reduce
from operator import mul

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

def ceil_divide(num, denom):
return -(-num // denom)

def update_memory(b, memory):
tuple_b = tuple(b)
if tuple_b in memory:
return False
return True

def dfs(a, b, t, memory):
if not update_memory(b, memory):
return False

p = prod([ai + bi for ai, bi in zip(a, b)])
if p == t:
return True
elif p > t:
return False

for i in range(len(a)):
b[i] += 1
if dfs(a, b, t, memory):
return True
b[i] -= 1

def solve(a, c):
b = [0 for _ in range(len(a))]
t = c * ceil_divide(prod(a), c)
while not dfs(a, b, t, set()):
t = t + c
return b

# 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]


Note: One potential use-case of such a thing could be, for example, to find the minimum padding of a $$D$$-dimensional tensor such that the number of elements of the tensor is divisible by $$c$$.

• How about finding an optimal $B$ such that $(\prod_{n=1}^{D} a_i) + B = kc$ ?? Feb 5 '21 at 21:13
• Can you tell us anything about typical sizes for $c$, $D$, and the $a_i$'s?
– D.W.
Feb 5 '21 at 22:22
• @D.W. So I don't have exact bounds on $c$, $D$ or $a_i$ but i think if your method is feasible for $c \leq 2^{20}$ and $D \leq 10$ then that's probably good enough for any use case I can think of. Even if it doesn't work for such large numbers then I'm just mostly interested in better alternatives to what I've already got. Feb 6 '21 at 0:12
• @droptop This is the initial value $t$ but I'm not sure how this solves the problem. Feb 6 '21 at 0:22
• I think you can speed up the search by only considering the smallest $t$ for which $B = t - ( \prod_{n=1}^{10} a_i)$, with $B$ having at least one factor in $a$. So if no $a_i$ divides $B$, move on to the next multiple of $c$. Feb 6 '21 at 12:51

If $$c$$ doesn't have too many factors and is not too large, the following should work reasonably well for typical numbers.

Factor $$c$$ into its prime factorization, say $$c=p_1^{e_1} \cdots p_k^{e_k}$$. Enumerate all tuples $$(f_1,\dots,f_D)$$ of non-negative integers such that $$f_1 \cdots f_D = c$$. For each such tuple $$(f_1,\dots,f_D)$$, and each $$i$$, find the smallest $$b_i$$ such that $$a_i+b_i$$ is a multiple of $$f_i$$ (namely, $$b_i = f_i - a_i \bmod f_i$$); then compute $$\prod_i (a_i+b_i)$$ and keep the smallest found so far. By construction, this will explore all possible solutions, and output the best one.

You can use the prime factorization of $$c$$ to help you enumerate all tuples $$(f_1,\dots,f_D)$$. In particular, the prime factorization makes it easy to enumerate all factors of $$c$$, so we can enumerate all possibilities for $$f_1$$ (namely, all factors of $$c$$), then enumerate all possibilities for $$f_2$$ (namely, all factors of $$c/f_1$$), etc.

The running time will be proportional to the number of such tuples you have to enumerate. This number is $$\prod_i {e_i + D - 1 \choose D - 1},$$ which is not too bad if the number of prime factors of $$c$$ is small, but is horrible if $$c$$ has many prime factors.

• This is very interesting! It is along the lines of what I was thinking but I was missing the idea of enumerating all permutations of factorizations of $c$ of length $D$. I do still wonder if this problem can be solved more efficiently, perhaps by by considering only the prime factors of $c$ not already included in the $a_i$s. I would like to leave this question open for a few days before accepting to see if more efficient algorithms might be found. Feb 6 '21 at 0:18
• @jodag That sounds like a good idea! I too wonder if there might be a better way. I like your idea: as an optimization: if $a_i$ divides $c$, then only enumerate tuples $(f_1,\dots,f_D)$ where $a_i$ divides $f_i$. It should be easy to adjust the method above to include that optimization. I'm not 100% certain that this optimization is correct in all cases (it will never cause us to miss a better solution) but it seems quite intuitive and reasonable to me.
– D.W.
Feb 6 '21 at 0:24
• The approach in this answer worked really well for me. I implemented a variant of it that can stop early and updated the original question with the algorithm and code. I am definitely still interested in further opimizations. One thing I would point out is that adding parenthesis to $b_i = (f_i - a_i) \bmod f_i$ may improve clarity, as it's not clear to me what the OOO for the mod operator should be. Feb 7 '21 at 0:29

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:
return d

def update_memory(b, c_rem, memory):
tuple_m = tuple(b + [c_rem])
if tuple_m in memory:
return False
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]