How does Mathematica solve this 2D piecewise function from constraints?

I am interested in learning the technique(s) used to solve the following 2D piecewise equation for $$x, y$$ when given $$z$$, where $$z = f(x, y)$$.

$$\text{Given that x \in \mathbb{Z}, y \in \mathbb{Z}}$$ $$f(x, y) = \begin{cases} x(4x - 3) + y & \text{y < x and y > -x} \\ y(4y - 3) + x & \text{y < x and y \leq -x} \\ y(4y - 1) - x & \text{y \geq x and y > -x} \\ x(4x - 1) - y & \text{y \geq x and y \leq -x} \\ \end{cases} \\$$

Mathematica is able to solve it easily using the following code:

f[x_, y_] = ConditionalExpression[Piecewise[{
{x (4 x - 3) + y, y < x && y > -x},
{y (4 y - 3) + x, y < x && y <= -x},
{y (4 y - 1) - x, y >= x && y > -x},
{x (4 x - 1) - y, y >= x && y <= -x}
}, Undefined],
x \[Element] Integers && y \[Element] Integers];

Solve[f[x, y] == 1000000000000000000000000000000000000, {x, y}] // Timing
(*Output: {0.03125, {{x -> -500000000000000000, y -> 500000000000000000}}}*)


My main question is how Mathematica solves this problem, as well as general mathematical techniques one could use to solve these types of problems.

I was not able to find any hints from Google, as the search results mostly yielded calculus related problems. Rather than having a second equation to solve for $$x$$ and $$y$$, Mathematica solves the equation from using the inequalities and the restriction of $$x, y$$ to $$\mathbb{Z}$$ instead.

The large $$z$$ value in Solve was chosen to see if Mathematica was simply iterating through all possible solutions, which was the approach I had attempted. It seems that Mathematica is employing a smarter technique, as the timing remains fairly consistent regardless of input size. Iterating through that many values would also probably take longer than 30–50 ms. The ideal speed, however, should be faster by a factor of 10.

I would like to learn the techniques for solving this equation, in order to apply the solution in other programs for visualization purposes.

If I had do guess, I'd say that Mathematica grabs each of the 4 cases and applies Solve to each individual one.

rhs = 1000000000000000000000000000000000000;
cond = {
{x (4 x - 3) + y == rhs, y < x, y > -x},
{y (4 y - 3) + x == rhs, y < x, y <= -x},
{y (4 y - 1) - x == rhs, y >= x, y > -x},
{x (4 x - 1) - y == rhs, y >= x, y <= -x}
};
Join @@ (Solve[#, {x, y}, Integers] & /@ cond) // Timing


{0.009538, {{x -> -500000000000000000, y -> 500000000000000000}}}

Being not an expert in diophantine equations, I cannot tell exactly how Mathematica solves each of these cases. But I guess, the first equation is used to eliminate one of the variables. Since only one unknown variable is left afterwards, quick numerical methods can be used to provide coarse bounds for the regions to so search for. If done clever enough, these regions are so small that they can be search exhaustively.