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$

Question

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
    memory.add(tuple_b)
    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$.

Answer 1

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.

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