Function such that at each point the nearest point whose value is larger than the current point is known I'm trying to work out (if possible) what the function would be that would satisfy the following.
Given two integers x, y
Function f(x, y) is a real number. f(x, y) does not tend to infinity as x or y does. Nor does it have any asymptotes. (I'd like the output of f to seem random following a power law distribution. Clearly, if it was random the function g could not exist, at least not without an unbounded search)
g(x,y) = x', y' such that f(x', y')>f(x,y) minimising (x-x')^2+(y-y')^2
Is it possible to construct the function f such that the function g can be calculated without searching the area around the point?
What I've thought of so far:

*

*Some form of gradient-based construction so that the search is not just surround but can have direction, although anything with this I can think of doesn't allow for looks random to the eye. It seems to me the result would always be (x,y) + normalized (dv/dx, dv/dy) so that the point was always adjacent. I guess what I want is for the larger point to not often and certainly not always be the adjacent point.

*A surround search but caching which doesn't really eliminate the surrounding search but might make a series of inspections of g in an area slightly better.

 A: To get a function with the desired properties, I would recommend using the composition of three functions.
$$f(x,y) = f_3\circ f_2\circ f_1(x,y)$$
The middle function $f_2:\mathbb Z^2\to[-2,2]$ does the heavy lifting of the "randomness". This function is of the form
$$f_2(x,y) = \cos(p_1(x))+\cos(p_2(y))\text{.}$$
where $p_1$ and $p_2$ are polynomials. This function ensures that the output is bounded but not asymptotic anywhere on the domain. The discrete nature of the domain means that for appropriate polynomials, the underlying periodicity is obscured. With some testing in Matlab, I rather like the results of the polynomials
$$p_1(x) = 15x^2-8x,p_2(y)=14y^2+5y\text{.}$$
Now for the most part, this function is pretty random-looking. However, due to the cartesian way these are combined, when viewed along an axis it is clear that for any fixed value of $x$, varying $y$ always has the same effect as for any other value of $x$, just translated up or down.
This will be the basis of our $g$ function--to find a higher value, one only needs to look along the coordinate axes. At a minimum, this converts a 2d search into four 1d searches, which should be much faster. Additionally, there is likely some work in number theory and rational approximations for $\pi$ which could make this process faster, if not deterministic.
However, these cartesian dependencies are unsightly for a function which is supposed to appear random. Correcting this will be the purpose of the first function $f_1:\mathbb Z^2\to\mathbb Z^2$. This could be anything which keeps points fairly close to their starting locations, while disrupting the cartesian dependence. For this case, I would suggest tesselling the plane with rectangular maps, moving points based on their coordinate $\text{modulo }n$. For example,
$$\begin{array}{r|cccc}
_x\setminus^y& 0& 1& 2& 3\\\hline
0&(0,0) &(2,1) &(3,2) &(1,3) \\
1&(1,2) &(3,3) &(2,0) &(0,1) \\
2&(2,3) &(0,2) &(1,1) &(3,0) \\
3&(3,1) &(1,0) &(0,3) &(2,2) \\
\end{array}$$
scrambles the numbers on a grid of size $n=4$.
Now for $g$, it's simple enough to find the row/column the input point belongs to, and follow the points that turn into the straight lines after applying the $f_1$ transform. It does add some extra checks once a larger value has been found, since following the straight line no longer guarantees the closest match, but it should be within the same square or an immediate neighbor. Note that larger squares do mean more checks, so there's a trade-off there.
Finally, we need the function to have an output with the desired distribution. This is where the final function $f_3:[-2,2]\to(0,\infty)$ comes in. This function simply matches values from the output of $f_2$ to the desired output space. $f_3$ is also monotonic increasing, so that the "larger" criterion of $g$ can be calculated directly from $f_2$.
If the CDF's of the power law distribution and $f_2$ distribution are known (respectively $F_P(z)$ and $F_{f_2}(z)$), then $f_3(z) = F_P^{-1}(F_{f_2}(z))$. Furthermore, the $f_2$ CDF is equivalent to the CDF of $\cos(X)+\cos(Y)$ over $X,Y\in[0,2\pi)$.
