Find all integer points that lay in a 3-ball with a given radius

Question

How can I efficiently find all lattice points in the cubic lattice $Z^3$ (that is to say, all integer points in a 3-space) that lay in a closed ball of radius $R$ centred at the origin?

Essentially,

Let $dist(p)$ be a function denoting the euclidean distance between a point in n-space and the origin point of that space, so $dist(p)=\sqrt{p_12+p_22+p_32\ldots p_n2}$.
How might I efficiently iterate over $\{p \in \mathrm{Z}^3\ | \ dist(p) \le R\}$?

I'm aware that this is trivial to do in $\mathbb{O}(R^3)$ time by iterating over all lattice points that lay inside the minimum bounding box of the ball and filtering out every point $p$ where $dist(p) > R$, and I'm also aware that this can be optimised by squaring both sides of the distance function, but this algorithm is still too slow for my needs.

• D.W.: What's $n$? Do you mean $R$? What counts as "efficiently"? Do you care about asymptotic running time, or about constant factors?
• bluelhf: Sorry, I mean $R,ドル and I care about constant factors (otherwise I'd settle for the naive algorithm, even an ideal algorithm would need to iterate over approximately $\frac{4}{3}\pi R^3$ points, which is $\frac{\pi}{6}$ times the iterations needed by the naive algorithm)

Answer 1

One straightforward approach is to iterate over $x,y$ in the bounding box, then find $z_\max$ so $(x,y,z)$ is in the sphere iff $|z| \le z_\max$. You can find $z_\max$ via the formula

$$z_\max = \sqrt{R - x^2 - y^2}.$$

If you can compute squares and square roots in constant time, then the running time is proportional to the number of pixels in the ball, i.e., you iterate over $\frac43 \pi R^3$ points.

A similar but slightly better algorithm is to iterate over $x,y$ in the bounding box, then iterate over increasing values of $z$ (i.e., $z=0,1,2,\dots$) until $z^2 > R - x^2 - y^2$. This is a little more efficient, because it doesn't require computing square roots, and because you can use the identity

$$(z+1)^2 = z^2 + 2z + 1$$

to avoid the need to compute squares in the inner loop either, instead updating the value of $z^2$ in each iteration. (Make sure to add in the negative values of $z$ too.)