How to enumerate 2D integer coordinates ordered by Euclidean distance?

The square of Euclidean distance between $(x, y)\in\mathbb{Z}^2$ and origin is $d = x^2+y^2$. How to enumerate the coordinates $(x, y)$ in ascending order of $d$?

For example, the first 14 sets of coordinates are:

d=0: { (0,0) }
d=1: { (1,0), (0,1), (0,-1), (-1,0) }
d=2: { (1,-1), (-1,-1), (-1,1), (1,1) }
d=4: { (2,0), (0,2), (-2,0), (0,-2) }
d=5: { (1,2), (-1,2), (1,-2), (-2,1), (-2,-1), (2,-1), (2,1), (-1,-2) }
d=8: { (2,2), (-2,2), (-2,-2), (2,-2) }
d=9: (0,3), (-3,0), (0,-3), (3,0)
d=10: (1,3), (-1,3), (3,1), (-3,1), (-3,-1), (1,-3), (-1,-3), (3,-1)
d=13: (2,-3), (-3,-2), (3,-2), (-2,-3), (-3,2), (3,2), (-2,3), (2,3)
d=16: (0,-4), (-4,0), (4,0), (0,4)
d=17: (-4,1), (-4,-1), (4,1), (1,-4), (4,-1), (-1,-4), (-1,4), (1,4)
d=18: (3,-3), (-3,-3), (-3,3), (3,3)
d=20: (4,2), (4,-2), (-4,-2), (2,-4), (-4,2), (-2,-4), (-2,4), (2,4)
d=25: (-3,-4), (-5,0), (5,0), (4,3), (-3,4), (-4,3), (0,-5), (4,-3), (-4,-3), (3,-4), (3,4), (0,5)


The first 14 iterations are depicted as:

          13
131210 9101213
1311 8 7 6 7 81113
12 8 5 4 3 4 5 812
10 7 4 2 1 2 4 710
13 9 6 3 1 0 1 3 6 913
10 7 4 2 1 2 4 710
12 8 5 4 3 4 5 812
1311 8 7 6 7 81113
131210 9101213
13


For a finite range of $(x, y)$, a trivial algorithm is to store all coordinates within the range into a list, and then sort the coordinates with $d$. However, this will need $O(n\log n)$ time and $O(n)$ space.

Another possible approach is to solve the diophantine equation $x^2 + y^2 = m$ for $m = 0, 1, \ldots, ... d_\mathrm{max}$. But it seems also a hard problem.

Is there any simpler way with lower time/space complexity?

Here is a C++ code of the trivial solution for reference.

• As far as storing goes, you can certainly take advantage of symmetries and perhaps more generally use the Sum of Squares function $r_2(n)$ (or $r'_2(n)$ if you're not interested in order/signs). I will also link this mathoverflow question as it seems to answer precisely what you desire. – Fimpellizieri Jun 28 '16 at 22:54
• Why the complexity of the trivial algorithm is not $O(n^2)$? – user84976 Jun 29 '16 at 10:35
• Your $n$ seems to be the number of points in your grid, not the edge length of the grid, right? So in terms of the edge length $a$, you get $O(a^2\log a)$ time complexity and $O(a^2)$ space complexity. But you can't go asymptotically lower than $O(a^2)$ time complexity since you have to enumerate this many points. Are you interested in tweaking constants here, or only in getting rid of the $\log a$ term or reducing the space complexity? – MvG Jul 5 '16 at 12:29

You can employ symmetry to only generating one eighth of a circle, e.g. $0\le x\le y$. It should be obvious how to do this.

You can reduce the space complexity to $O(r)=O(\sqrt n)$, i.e. linear in the radius of your disk, not in the number of grid points it covers. The way to do this is by considering parallel lines of grid points. Then on each line you know that the points will be enumerated starting closest to the center and moving outwards, so it is enough to keep track of the next candidate on each line.

The most appropriate data structure for this keeping track is most likely a priority queue, which will give you quick access to the element with minimal $d$. But even then, elementary operations in a priority queue have $O(\log n)$ time complexity. Which in this case means you are still at $O(r^2\log r)=O(n\log n)$ overall time complexity since you enumerate $O(r^2)$ elements and maintain a data structure where you have to perform two $O(\log r)$ operations for each of them.

Here is an implementation of these ideas. The core portions are these:

struct Pt { // essentially your struct Point, just shorter names
int x, y, d;
Pt(int x, int y) : x(x), y(y), d(x*x + y*y) {}
bool operator<(const Pt& p) const { return d > p.d; } // flipped!
};

void enumerateGridPoints(int r) {
priority_queue<Pt> q;
q.push(Pt(0, 0));
int d, level = -1;
while ((d = q.top().d) <= r*r) {
reportLevel(++level, d);
do {
Pt p = q.top();
q.pop();
reportPoint(p);
if (p.y < p.x) q.push(Pt(p.x, p.y + 1));
if (p.y == 0) q.push(Pt(p.x + 1, 0));
} while (q.top().d == d);
}
}


Since the C++ std::priority_queue is a max queue by default, and since I'm lazy, I just flipped the sign in the comparison operation so that the element with smallest d is the maximal element and thus extracted next.

One benefit of the above algorithm compared to yours is that it's easy to modify it so that enumerates points without any preset limit. That might be useful if you'd like to e.g. inspect sites until you find something you are looking for, but you don't know up front how long you will have to search until you find it. Both time and space complexity scale with the generated output, so small initial portions can be obtained quickly while still allowing you to proceed to larger results.

Edsger Dijkstra gives an elegant recursive solution to the problem of finding solutions to $x^2 + y^2 = m$ in his 1976 book The Discipline of Programming. The solution is outlined here, and involves sweeping potential $x$ values down from $\sqrt m$ and $y$ values up from 0.

Consider the function $B(x, y)$ that returns solutions between $x$ and $y$, defined recursively as follows:

$$B(x, y) = \begin{cases} \varnothing, & \text{if x < y (base case)} \\ (x, y) \cup B(x-1, y+1), & \text{if x^2 + y^2 = m (satisfying case)} \\ B(x, y+1), & \text{if x^2 + y^2 < m} \\ B(x-1, y), & \text{if x^2 + y^2 > m} \\ \end{cases}$$

Using this definition, you can calculate $B(\sqrt m, 0)$ for each $m=0,1,…,d_{max}$. Written as a generator or iterator, this algorithm will have constant space complexity.

• Would you agree that the time complexeity is $O(r^3)$ here (with $r=\sqrt{d_{\text{max}}}$), since for each of the $O(r^2)$ possible values for $m$ you'd have $O(r)$ steps of evaluating $B$ till you reach the base case? – MvG Jul 6 '16 at 11:25