Let \(M = 10^9\). Clearly \(D_{min} \leq M\).

We do a binary search on \(D\), and look for a solution with penalty not exceeding \(D\). Initially set \(x_i\) to smallest possible values, and \(y_i\) to values certainly larger than any solution, e.g. \(x_i = i\), \(y_j = D + M + j\).

Define \(L_{i, j} = \max(1, a_{i, j} - D)\), \(R_{i, j} = a_{i, j} + D\) — bounds on \(x_i y_j\). We will increase \(x_i\) and decrease \(y_j\) when we have to, that is:

• if \(x_i y_j > R_{i, j}\), we can put \(y_j = \lfloor R_{i, j} / x_i \rfloor\),
• if \(x_i y_j < L_{i, j}\), we can put \(x_i = \lceil L_{i, j} / y_j \rceil\),
• if \(x_i \geq x_{i + 1}\), put \(x_{i + 1} = x_i + 1\),
• if \(y_j \geq y_{j + 1}\), put \(y_j = y_{j + 1} - 1\).

Note that at any point \(x_i\) are unconditional lower bounds to any solution with penalty at most \(D\), similarly \(y_j\) are unconditional upper bounds. Therefore, with this process we will always find a solution if there exists one.

Assume additionally that in the solution \(x_n \leq y_m\) (the opposite case is symmetrical). We will terminate the process once \(x_n > y_m\), or \(y_j \leq 0\) for any \(j\). Since \(x_n y_m \leq D + M\), we have \(x_n \leq \sqrt{D + M}\). Thus throughout the process we will make at most \(n \sqrt{D + M}\) increases to \(x_i\), thus at most \(n m \sqrt{D + M}\) decreases to \(y_j\), and \(n^2m \sqrt{D + M}\) operations in total (note that for each decrease to \(y_j\), we need to check \(O(n)\) conditions).

The total complexity is \(O(nm(n + m) \sqrt{M} \log{M})\). There is also an \(O(nm(n + m) \sqrt{M})\) variant of the above solution without the binary search.