# Find the number of points of distance n away from origin as function of n

I came across a seemingly simple problem the other day and I thought I'd share it with anyone interested.

Say you have a point in 3 dimensions.

The number of points that are of distance $0$ away is $1$ (itself).

$1 \rightarrow 6$ (up/down/left/right/in/out).

$\sqrt{2} \rightarrow 12$ (across all 2D diagonals).

$\sqrt{3} \rightarrow 8$ (across all 3D diagonals).

$\sqrt{4}=2 \rightarrow 6$

...

$\sqrt{7} \rightarrow 0$

etc...

Does anyone know how to create an expression or recurrence relation to find the number of points distance $\sqrt{n}$ away for any given $n$?

• OEIS A005875 – Henry Apr 15 '16 at 15:50
• That ended up being much more complicated than I thought. Thank you for the link! – Brian Apr 15 '16 at 16:07
• Complicated? We have $n = \sqrt{ x^2 + y^2 + z^2 }$ where $n \in \mathbb{N}$, and $x, y, z \in \mathbb{Z}$. So, you just need to count the number of integer triplets $x, y, z$ whose squares summed equals $n$ squared, $x^2 + y^2 + z^2 = n^2$. So, actually straightforward -- only hard if you wanted to do it algebraically. (It is also surprisingly easy to construct the sequence programmatically, using a triple loop and a histogram for $n^2$. Use symmetries to make it faster. Works in other dimensions (2D, 4D, etc.), too.) – Nominal Animal Apr 16 '16 at 19:27
• @NominalAnimal I realize that it's straightforward, but I was asking if there was an algebraic way to do it. I originally did it manually (part of some homework for a class, find the number of points distance n away from the center of a cubic lattice). Then I tried doing it programmatically (to see if I could find a pattern) and I was curious if there was an algebraic way to solve it. – Brian Apr 18 '16 at 12:55
• I think if you want more than that link, you need to be more specific about what you want. The OEIS entry is highly likely to contain (links to) everything that's known about this sequence. It contains expressions for computing elements of the sequence. If you want more than that, or you want explanations of things in the entry, you should write exactly what you're missing. – joriki May 6 '16 at 22:16

Consider this naïve C implementation of counting the cubic lattice points as a function of squared distance, the sequence $s(n)$ for $n = 0 .. max$, where the largest distance considered is $\sqrt{max-1}$:

#include <stdlib.h>
#include <math.h>

void histogram_3d(unsigned long sqr[], unsigned long max)
{
const unsigned long n = 1UL + (unsigned long)ceil(sqrt((double)max));
unsigned long       z, y;

for (z = 0UL; z <= n; z++) {
for (y = 0UL; y <= n; y++) {
const unsigned long zzyy = z*z + y*y;
if (zzyy < max) {
const unsigned long c = 2UL * (y == 0UL ? 1UL : 2UL) * (z == 0UL ? 1UL : 2UL);
unsigned long       i = zzyy + 1UL;
unsigned long       d = 3UL;
sqr[zzyy] += c / 2UL;
while (i < max) {
sqr[i] += c;
i += d;
d += 2UL;
}
}
}
}
}


The innermost loop is over the squares of $x$. It applies the fact that the difference between consecutive squares ($1,4,9,16,25,\dots$) increments by 2 every step -- the differences are (starting at 1) $3,5,7,9,\dots$.

Given $max \ge 10$, the ten first elements in sqr array are incremented by 1,6,12,8,6,24,24,0,12,30. On a HP laptop with a Core i5-4200U CPU (x86-64, unsigned long is 64-bit), this calculates the first million terms of the sequence in about three seconds; the first ten million terms in about 180 seconds (three minutes). Using 32-bit unsigned int type instead, on the same machine and compiler, calculating the first million terms takes about two seconds, and first ten million terms about 160 seconds.

It is also possible to modify the above to consider only some window $n = min..max$. That would allow better cache usage (if the histogram for the entire window fits in cache), but I have not tested in practice whether this yields significant speedups or not for very large $n$.

Simply put, the naïve method may yield quite acceptable performance in practice, even though the algorithm technically has $O(n^3)$ time complexity.

The OEIS A005875 page mentions that the generating function for this sequence is $$\left(\sum_{i=0} x^{i\cdot i}\right)^3 = 1 + 6x + 12x^2 + 8x^3 + 6x^4 + 24x^5 + 24x^6 + 0x^7 + 12x^8 + 30x^9 + O(x^{10})$$ In other words, expanding the above sum, the coefficients of the increasing powers of $x$ yield the (nonzero) terms in the sequence.

The hard part here is the fact that the sum is cubed. (The powers of $x$ in the sum are $0, 1, 4, 9, 16, 25, 36, \dots$, so they are not a problem.) Using a triple loop over this sequence (instead of the lattice coordinates within the $\sqrt{max}$-radius ball as in the example code above) might yield a constant factor speedup, but optimized code would be much more complex. This, too, would have $O(n^3)$ time complexity.

The OEIS page also mentions that the sequence is an Euler transform of period 4 sequence $[6, -9, 6, -3, \dots]$. This means that $$s(n) = \sum_{k=0}^n -1^k \binom{n}{k} a_k = \sum_{k=0}^n -1^k \frac{n!}{k!(n-k)!} a_k$$ where $a_0 = a_{4i} = 6$, $\;a_1 = a_{4i+1} = -9$, $\; a_2 = a_{4i+2} = 6$, $\; a_3 = a_{4i+3} = -3$, $\;0 \le i \in \mathbb{Z}$. The negations occur in pairs, so actually $$s(n) = \sum_{k=0}^n \binom{n}{k} b_k$$ where $b_i = \lvert a_i \rvert$.

(In C, $b_i = \lvert a_k\rvert$ can be evaluated using 3*( ((i+1)&2) | ((i+2)&3)) ), which may be faster than a lookup on some architectures; a lookup may have unwanted cache effects. Another option is (13974 >> ((i&3)<<2)) & 15, which uses only bit shifts and two binary and operations, essentially packing the sequence into a 16-bit hexadecimal number, 13974 = 0x3696.)

This is a very promising approach in general, as its time complexity looks to be $O(n^2)$, if we disregard the complexity of computing the binomial coefficients. It may very well be possible to calculate the coefficients using Pascal's triangle, although extra memory may be needed to store the binomial coefficients for each innermost loop.

Essentially, this approach is mathematically fully defined, but the efficiency of the approach depends on the implementation details -- and as such, is really a programming question, not a mathematical one.