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

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_1²+p_2²+p_3²\ldots p_n²}$$.
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.

• What's $n$? Do you mean $R$? What counts as "efficiently"? Do you care about asymptotic running time, or about constant factors?
– D.W.
Jan 6 at 21:51
• 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) Jan 6 at 22:11

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.)