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) == [1, 1, 1]
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$.