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?


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.

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

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.