Integer Optimization with Constraints

As your main objective is computer efficiency, I advise you to consider the following continuous setting (with correspondence: $x \leftrightarrow n, \ y \leftrightarrow m$):

How to maximize quantity:

$$q=yx-x^2 \tag{1}$$

under the following constraint:

$$0<x^2+y^2 \le d \ ? \tag{2}$$

(Please note that I have written (2) with an $\le$ instead of a $<$ for an easier handling).

(1) can be written

$$y=x+\frac{q}{x}\tag{1'}$$

Condition (2) means that the points must be in the disk centered at the origin with radius $\sqrt{d}$.

Condition (1') can be interpreted, for a given $q$, as the equation of a hyperbola $(H_q)$, all of them sharing the same vertical asymptote $x=0$ and slant asymptote $y=x$.

Fig. 1 : Some hyperbolas $(H_q)$, one of them being the optimal one with point of tangency $(x_0,y_0)$ (in green) given by formulas (4).

Therefore, among all those "russian dolls" hyperbolas, a single one is tangent to the disk, as can be seen on at two symmetrical points of tangency $\pm(x_0,y_0)$. An explicit formulation of $\pm(x_0,y_0)$ can be obtained with the technique of Lagrange multiplier. Do you know it ? See appendix below.

Now, coming back to the discrete setting, the optimal tuple $(n,m)$ is in the close vicinity of $(x_0,y_0)$. This is adapted for a computer approach, lowering drastically the number of points of $\mathbb{N^2}$ to be investigated (you should be in the "relatively quickly" range...).

Of course, a specific mathematical study using some kind of integer programming is surely possible, but it isn't in my "field of expertise"...

Appendix : Lagrange method amounts to express that the gradients (obtained by computing partial derivatives) of the two curves with equations:

$$yx-x^2=q \ \ \text{and} \ \ x^2+y^2=d$$

are proportional at the point(s) of tangency. This gives rise to the system of three equations with 3 unknowns :

$$\begin{cases}y-2x&=&\lambda 2x\\x&=&\lambda 2y\\x^2+y^2&=&d\end{cases}\tag{3}$$

whose positive solutions are:

$$\begin{cases} x_0&=&e\\y_0&=&(1+\sqrt{2})e\end{cases} \ \text{with} \ e=\tfrac12\sqrt{(2-\sqrt{2})d}\tag{4}$$

(the value of $\lambda$ is unimportant). Please note that all potential minimization points lie on a same line.

Figure 1 has been generated by the following Sage program ; see in particular the last lines implementing the solution of system (3).

 var('x y d L a e')
 a=5
 d=2*(1+sqrt(2))
 g=line(((-a,0),(a,0)),color='black')
 g+=line(((0,-a),(0,a)),color='black')
 g+=line(((-a,-a),(a,a)),color='black')
 g+=implicit_plot(x^2+y^2-d,(x,-a,a),(y,-a,a),color='red')
 for k in range(1,6) :
    g+=implicit_plot(x*y-x^2-k,(x,-a,a),(y,-a,a))
 e=sqrt((2-sqrt(2))*d)/2
 g+=point((e,(1+sqrt(2))*e),color='green',size=50)
 show(g)
 s=solve([y-2*x==L*2*x,x==L*2*y,x^2+y^2==d],x,y,L)
 show(s)