Detailed mathematical approach

Problem Summary

Let \(S(n)\) denote the sum, over every lattice cube that can be translated to lie inside \([0,n]^3\), of the number of lattice points contained in that cube. A lattice cube means that all eight vertices lie in \(\mathbb{Z}^3\).

A direct search over all vertex sets would be hopelessly redundant. The implementation instead enumerates primitive cube frames, extends them by every admissible integer scale, counts how many integer translations fit inside the box, and adds the exact lattice-point count of each scaled cube.

Mathematical Approach

Describe one cube by three edge vectors \(u,v,w\in\mathbb{Z}^3\) issuing from a single vertex. They must satisfy

$u\cdot v=u\cdot w=v\cdot w=0,\qquad u\cdot u=v\cdot v=w\cdot w.$

The key is that such integer triples can be generated systematically from integer quaternions, and once one primitive frame is known, all larger cubes in the same orientation follow by integer scaling.

Step 1: Generate primitive orthogonal edge triples from integer quaternions

Start from integers \(a,b,c,d\), not all zero, and define

$\begin{aligned} u&=\frac{1}{g}\left(a^2+b^2-c^2-d^2,\ 2(bc-ad),\ 2(bd+ac)\right),\\ v&=\frac{1}{g}\left(2(bc+ad),\ a^2-b^2+c^2-d^2,\ 2(cd-ab)\right),\\ w&=\frac{1}{g}\left(2(bd-ac),\ 2(cd+ab),\ a^2-b^2-c^2+d^2\right). \end{aligned}$

The divisor \(g\) removes the compulsory common factor forced by parity:

$g=\begin{cases} 4,& \text{if }a,b,c,d\text{ are all odd},\\ 2,& \text{if exactly two of them are odd},\\ 1,& \text{otherwise}. \end{cases}$

When \(\gcd(a,b,c,d)=1\), these vectors are primitive in the sense needed by the search, and they satisfy

$u\cdot v=u\cdot w=v\cdot w=0,\qquad u\cdot u=v\cdot v=w\cdot w=\left(\frac{a^2+b^2+c^2+d^2}{g}\right)^2.$

If we set

$N_0=a^2+b^2+c^2+d^2,\qquad L_0=\frac{N_0}{g},$

then \(L_0\) is the Euclidean edge length of the primitive cube frame.

Step 2: Pass from one primitive frame to every integer scale

For any integer scale \(s\ge 1\), the same orientation gives a larger cube with edges

$su,\qquad sv,\qquad sw,$

and edge length

$L=sL_0.$

For a lattice vector \(x=(x_1,x_2,x_3)\), write

$\delta(x)=\gcd(|x_1|,|x_2|,|x_3|).$

The segment from \(0\) to \(sx\) then contains \(s\,\delta(x)+1\) lattice points. For the three primitive edges define

$\Delta_0=\delta(u)+\delta(v)+\delta(w),\qquad \Delta=s\Delta_0.$

These quantities measure the lattice step size along the three edge directions and are exactly what the point-counting polynomial needs.

Step 3: Derive the exact lattice-point count of one scaled cube

A standard lattice-parallelepiped formula says that for edge vectors \(p,q,r\in\mathbb{Z}^3\),

$\#\bigl((\text{parallelepiped}(p,q,r))\cap\mathbb{Z}^3\bigr)=|\det(p,q,r)|+\delta(p\times q)+\delta(q\times r)+\delta(r\times p)+\delta(p)+\delta(q)+\delta(r)+1.$

Apply this with \(p=su\), \(q=sv\), \(r=sw\). Because the edges are mutually orthogonal and all have length \(L\), we get

$|\det(su,sv,sw)|=L^3.$

Also, orthogonality and equal edge lengths imply

$ (su)\times(sv)=\pm L(sw),\qquad (sv)\times(sw)=\pm L(su),\qquad (sw)\times(su)=\pm L(sv).$

Therefore

$\delta((su)\times(sv))=L\,\delta(sw),\qquad \delta((sv)\times(sw))=L\,\delta(su),\qquad \delta((sw)\times(su))=L\,\delta(sv).$

Summing the three face terms gives \(L\Delta\), and summing the edge terms gives \(\Delta\). Hence the lattice-point count of one scaled cube is

$\boxed{\operatorname{points}(s)=L^3+(L+1)\Delta+1.}$

As a sanity check, for the axis-aligned cube of side \(m\) we have \(L=m\) and \(\Delta=3m\), so the formula becomes \((m+1)^3\), exactly as expected.

Step 4: Count how many translations fit inside \([0,n]^3\)

For the primitive frame define the axis spans

$\sigma_x=|u_1|+|v_1|+|w_1|,\qquad \sigma_y=|u_2|+|v_2|+|w_2|,\qquad \sigma_z=|u_3|+|v_3|+|w_3|.$

These are the side lengths of the smallest axis-aligned box containing the primitive cube. The reason is that each vertex coordinate is a subset sum of the edge coordinates, so the total range along one axis is the sum of the absolute values on that axis.

After scaling by \(s\), the spans become \(s\sigma_x\), \(s\sigma_y\), \(s\sigma_z\). The number of admissible integer translations is therefore

$T(s)=(n+1-s\sigma_x)(n+1-s\sigma_y)(n+1-s\sigma_z),$

provided every factor is positive. Equivalently, the maximum allowed scale is

$s_{\max}=\min\left(\left\lfloor\frac{n}{\sigma_x}\right\rfloor,\left\lfloor\frac{n}{\sigma_y}\right\rfloor,\left\lfloor\frac{n}{\sigma_z}\right\rfloor\right)=\left\lfloor\frac{n}{\max(\sigma_x,\sigma_y,\sigma_z)}\right\rfloor.$

Step 5: Sum every frame, every scale, then divide by the 24-fold symmetry

For one primitive frame, the raw contribution is

$\sum_{s=1}^{s_{\max}} T(s)\left((sL_0)^3+(sL_0+1)s\Delta_0+1\right).$

The quaternion parameterization enumerates ordered orthogonal edge triples, not unlabeled geometric cubes. The same geometric cube appears in \(24\) equivalent orientations, corresponding to the rotational symmetry group of the cube. So if \(R(n)\) is the total raw sum over all primitive frames, then

$\boxed{S(n)=\frac{R(n)}{24}.}$

This is why the implementation checks that the raw total is divisible by \(24\) before reporting the final answer.

Step 6: Worked example for \(n=2\)

Take the primitive axis frame

$u=(1,0,0),\qquad v=(0,1,0),\qquad w=(0,0,1).$

Then \(L_0=1\), \(\Delta_0=3\), and \(\sigma_x=\sigma_y=\sigma_z=1\), so \(s_{\max}=2\).

For \(s=1\), we have

$T(1)=2^3=8,\qquad \operatorname{points}(1)=1^3+(1+1)\cdot 3+1=8,$

giving a contribution of \(8\cdot 8=64\).

For \(s=2\), we have

$T(2)=1,\qquad \operatorname{points}(2)=2^3+(2+1)\cdot 6+1=27,$

giving a contribution of \(27\).

No larger scale fits, so

$S(2)=64+27=91,$

which matches the checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations all use the same underlying search. The core solver splits the quaternion enumeration by parity pattern, because the parity determines the normalization divisor \(g\) and the norm bound \(N_0\le gn\). That lets it prebuild only the integer values that can actually occur in each case.

During enumeration, partial sums of \(a^2+b^2+c^2+d^2\) are checked early, so branches that already exceed the bound are abandoned immediately. The search also rejects non-primitive quadruples by gcd tests and uses canonical sign restrictions to avoid the trivial duplication coming from replacing the parameter quadruple by its negation.

For every surviving frame, the implementation computes the three edge vectors, rejects frames whose axis spans already exceed the box, derives \(L_0\), \(\Delta_0\), and \(s_{\max}\), and then loops over all valid scales \(s\). At each scale it multiplies the translation count by the cube point-count formula, adds the result to a thread-local accumulator, and only after all threads finish divides the raw total by \(24\).

The Python and Java versions are thin wrappers around the same compiled solver, so all three languages share the same mathematics, the same checkpoints, and the same final value.

Complexity Analysis

The dominant cost is the quaternion search. A loose geometric upper bound comes from the four-dimensional ball \(a^2+b^2+c^2+d^2\le 4n\), which contains \(O(n^2)\) integer tuples, so the worst-case search space is roughly quadratic in \(n\).

In practice the implementation is much faster than that crude bound suggests. Partial norm checks stop many branches early, gcd filters remove non-primitive tuples, span checks discard frames that can never fit in the box, and each surviving frame only iterates up to \(s_{\max}\). The method is therefore strongly prune-driven rather than a naive enumeration of all cubes.

Memory usage is modest. The precomputed even and odd value tables only reach \(O(\sqrt{n})\) in magnitude, and the parallel version stores only task metadata plus one accumulator per worker thread, so the overall memory footprint stays small compared with the arithmetic work.

Footnotes and References

1. Problem page: Project Euler 579
2. Quaternion background: Wikipedia - Quaternion
3. Euler-Rodrigues parameterization: Wikipedia - Euler-Rodrigues formula
4. Lattice polytopes: Wikipedia - Lattice polytope
5. Ehrhart theory: Wikipedia - Ehrhart polynomial