# How many squares in a three-dimensional $n \times n \times n$ cartesian grid?

This brings the classical question to three dimensions. Given a three-dimensional Cartesian grid of $$n \times n \times n$$ points (that is $$(n-1) \times (n-1) \times (n-1)$$ unit cubes), how many squares whose vertices are points of the grid are there? The post is quite long as I give some families of squares of the 3D grid, but the main question is:

Can we actually give a closed form expression of the number of squares whose vertices are points of a $$n \times n \times n$$ Cartesian grid?

I give below part of my attempt but it is highly inconclusive even though it generates some simple hidden squares.

A square can be represented using its sides as two perpendicular vectors of same positive length $$u,v$$ (that is $$u \cdot v = 0$$ and $$u \cdot u = v \cdot v$$), with corners $$0, u, v, u+v$$. This representation is not unique $$(u,v) = (v,u) = (-u,v) = \dots$$ for instance if we are looking at classes of squares identical by translation (essentially that is eight different ways to write the same class of squares).

In the plane given $$u = \begin{bmatrix} a & b \end{bmatrix}$$, $$v$$ must take the form $$v = k \begin{bmatrix} -b & a \end{bmatrix}$$ for some $$k$$ and yet have same length as $$u$$, that is $$k = \pm 1$$. If $$u,v$$ are two vectors forming a square, of the vectors $$u,v,-u,-v$$ exactly one of them must lie on the first orthant call it $$x = \begin{bmatrix} a & b \end{bmatrix}$$ with $$a >0$$ and $$b \geq 0$$. We can call $$y = \begin{bmatrix} -b & a \end{bmatrix}$$ the vector left such that $$(x,y)$$ is direct. This defines a unique representant for the class of squares $$(u,v)$$. We just need to add that such a square $$(x,y)$$ is tightly included in a straight square of dimensions $$(a+b) \times (a+b)$$, and then we can properly count and find that there are $$\frac{n^2(n^2-1)}{12}$$ different squares in a $$n\times n$$ planar Cartesian grid.

In three dimensions the question seems way trickier since we do not have this simple characterization. I was reasoning on the number of null coordinates to first get an idea of some simple hidden squares.

1. If $$u,v$$ have a total of six null coordinates we would not call that a square.
2. Surely $$u,v$$ cannot have five null coordinates in total.
3. The vectors $$u,v$$ could totalize four null coordinates, this brings the family of squares $$u = \begin{bmatrix} k & 0 & 0 \end{bmatrix}, v = \begin{bmatrix} 0 & k & 0 \end{bmatrix}$$ and so on, that is squares that are aligned with grid.
4. If $$u,v$$ have three null coordinates, it must be that say $$u = \begin{bmatrix} a & b & 0 \end{bmatrix}, v = \begin{bmatrix} 0 & 0 & c \end{bmatrix}$$ (not all $$0$$ could belong to $$u,v$$ but also if two $$0$$ were facing each other in $$u,v$$ we would have a problem with $$u \cdot v = 0$$), and $$a^2 + b^2 = c^2$$. These are essentially squares with one side sideways in a parallel plane to the grid but of integer length (Pythagorean triple) and another side perpendicular to that plane of the same integer length.
5. If $$u,v$$ have two null coordinates it must be that say $$u = \begin{bmatrix} a & b & 0 \end{bmatrix}, v = \begin{bmatrix} c & d & 0 \end{bmatrix}$$ since $$u \cdot v = 0$$, but again it must then be that $$v = \pm \begin{bmatrix} -b & a & 0 \end{bmatrix}$$, these squares are essentially the sideways square in the plane.
6. If there is exactly one null coordinate, we can write say $$u = \lambda \begin{bmatrix} pa & pb & c \end{bmatrix}, v = \pm \lambda \begin{bmatrix} -qb & qa & 0 \end{bmatrix}$$ where $$q > p > 0$$ are coprime, $$a,b$$ are coprime, $$c^2 = (q^2-p^2)(a^2+b^2)$$ and $$\lambda$$ is just a non-null integer factor. I am not sure what are all the solutions but when $$(x,y,z)$$ and $$(u,v,w)$$ are Pythagorean triples, then $$a=u, b=v, c=yw, p=x, q=z$$ is a solution.

I know that this does not encompass all solutions, for instance $$u = \begin{bmatrix} 2 & 2 & 1 \end{bmatrix}, v = \begin{bmatrix} 1 & -2 & 2 \end{bmatrix}$$ is another solution. So I was reasoning on the different numbers used. Here in $$u,v$$ above, the numbers used were $$1,2$$.

1. If only one number is used, $$u\cdot v = 0$$ is not be feasible.
2. If two distinct numbers are used say $$a,b$$, it must be that neither $$u$$ nor $$v$$ is constituted of the same one number (otherwise we end up on relations such as $$3a^2 = a^2 + 2b^2$$ for instance). Furthermore they can't contain different mixtures of $$a$$ and $$b$$, otherwise we would have something like $$a^2+2b^2 = b^2 + 2a^2$$ by comparing the lengths. Therefore $$u,v$$ contain the same mixture say two $$a$$ and one $$b$$ each, surely then the two $$b$$ can't be at the same position (we would have to solve $$0 = b^2$$ or $$a^2 = 2b^2$$ for the $$u \cdot v = 0$$ condition), so in the dot product we have the numbers $$ab,ab,a^2$$ with $$\pm$$ signs equalling $$0$$, surely $$a^2 = 2ab$$ that is $$a = 2b$$, we find the family of multiples of the above $$u,v$$ example and equivalents.
3. Now assume three distinct numbers $$a,b,c$$ are used. Working through all cases one can figure out that there are only three families $$u = \begin{bmatrix} a & b & c \end{bmatrix}, v = \begin{bmatrix} a & b & -c \end{bmatrix}$$ with $$a^2 + b^2 = c^2$$ and equivalents; the family of $$u = \begin{bmatrix} a & b & c \end{bmatrix}, v = \begin{bmatrix} a & -c & -b \end{bmatrix}$$ with $$a^2 = 2bc$$ and equivalents; and finally $$u = \begin{bmatrix} a & b & c \end{bmatrix}, v = \begin{bmatrix} -c & a & b \end{bmatrix}$$ where $$ab+bc = ac$$ that is $$\frac1b = \frac1a+\frac1c$$ (for example $$a = 4, b = 3, c = 12$$ or $$a = 10, b = 6, c = 15$$) and equivalents.
4. If there are four distinct numbers used, there are even more cases it only gets more complicated at this point.
• – Ethan Bolker Jul 27 at 21:55
• I published an OEIS sequence that might be of interest: oeis.org/A334881 – Peter Kagey Jul 27 at 22:38
• That is indeed exactly the answer, I assume it was computed? – Olivier Massicot Jul 27 at 22:50
• That's right—my algorithm was very brute-force. Pontus von Brömssen computed more terms, so you can contact him to ask if he did anything clever. – Peter Kagey Jul 27 at 23:46
• @Peter Kagey: With more brute-force ($5$ hours worth), here is the last term listed in the oeis sequence followed by $8$ more terms . . . $$91200,139338,206394,296832,417120,575556,779238,1037514,1359792$$ – quasi Jul 29 at 21:38

I computed the first 100 entries of the sequence in about half an hour.

The naive approach is to look at every triplet of 3 points, and see if they define 3 vertices of a square. This has a memory requirement of $$O(1)$$, but a time requirement of $$O(n^9)$$. My approach was instead to categorize every vector, and afterwards do a fairly simple computations to count the full number of squares. I implemented it with a memory requirement of $$O(n^3)$$ and a time requirement of $$O(n^5)$$.

Consider 2 vectors $$u=[a,b,c], v=[i,j,k]$$ . These vectors define the sides of a square if $$|u|=|v|$$ and $$u\cdot v=0$$. Iterated over every pair of vectors, each square has exactly 4 vector pairs which count it, corresponding to each of $$(\pm u, \pm v)$$.

With these vectors, we can then deduce how many 'copies' of this square are in the cube. The smallest 'box' aligned to the grid that contains the square has dimensions $$[|a|+|i|,|b|+|j|,|c|+|k|]$$; since the square is made of all 4 combinations of $$(\pm u, \pm v)$$ for each coordinate there is one where the components of the vectors have the same sign. Then, the number of these 'boxes' that fit in the first dimension of the cube is $$(n-(|a|+|i|))(n-(|b|+|j|))(n-(|c|+|k|))$$, unless one of the components is larger than the grid, i.e. if n=5 and |a|+|i|=6 then there are no copies of that square in the grid.

So, the ultimate computation that my code computes is $$\frac{1}{4} \sum_{\substack{u=[a,b,c],v=[i,j,k]\\ a,b,c,i,j,k\in[-n,…,n] \\ |u|=|v|\\ u\cdot v=0\\ |a|+|i|,|b|+|j|,|c|+|k|\leq n}} (n-(|a|+|i|))(n-(|b|+|j|))(n-(|c|+|k|))$$.

It should be possible to very carefully choose which vectors to check to make each square counted uniquely or exactly twice by taking symmetry into account, but it would need to be very careful.

My code does the calculation in a different order than the straight summation so that every vector is considered exactly once, not twice. First, it generates every vector, and adds it to a dictionary with vectors of the same length. For every vector in that dictionary that it is orthogonal too, it adds a count to the 'box' which it fits in, and keeps track of how many squares are counted in each box. Then, for every term up to $$n$$, it does the product and sum to determine how many copies of that box (and associated squares) are in the cube, and divides them by 4, and finally spits out all of the terms less than $$n$$ at once.

For the time calculation: Generating the vectors is $$O(n^3)$$.There are $$O(l^2)$$ vectors of length $$l$$, so each vector is compared to at most $$O(n^2)$$ vectors, and these are multiplied together resulting in a computation time of $$O(n^5)$$. The remaining computations are cheap.

My code is on github, here: https://github.com/fibbooo/OEIS/blob/master/A334881.py