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.
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.
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.
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.
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.
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.$$
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.
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.
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.
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.
Sei \(S(n)\) die Summe der Gitterpunktzahlen aller Gitterwürfel, die sich durch eine Translation vollständig in \([0,n]^3\) legen lassen. Ein Gitterwürfel bedeutet hier, dass alle acht Ecken in \(\mathbb{Z}^3\) liegen.
Eine direkte Suche über alle möglichen Eckmengen wäre extrem redundant. Die Implementierung enumeriert stattdessen primitive Würfelrahmen, erweitert jeden Rahmen um alle zulässigen ganzzahligen Skalierungen, zählt die passenden ganzzahligen Translationen im Kasten und addiert jeweils die exakte Anzahl der Gitterpunkte des skalierten Würfels.
Ein Würfel wird durch drei Kantenvektoren \(u,v,w\in\mathbb{Z}^3\) beschrieben, die von derselben Ecke ausgehen. Sie müssen erfüllen
$$u\cdot v=u\cdot w=v\cdot w=0,\qquad u\cdot u=v\cdot v=w\cdot w.$$
Der zentrale Trick besteht darin, solche ganzzahligen orthogonalen und gleich langen Tripel systematisch aus ganzzahligen Quaternionen zu erzeugen. Kennt man einen primitiven Rahmen, dann liefern alle ganzzahligen Vielfachen davon die größeren Würfel derselben Orientierung.
Man startet mit ganzen Zahlen \(a,b,c,d\), nicht alle null, und definiert
$$\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}$$
Der Divisor \(g\) entfernt den durch die Parität erzwungenen gemeinsamen Faktor:
$$g=\begin{cases} 4,& \text{falls }a,b,c,d\text{ alle ungerade sind},\\ 2,& \text{falls genau zwei davon ungerade sind},\\ 1,& \text{sonst}. \end{cases}$$
Gilt zusätzlich \(\gcd(a,b,c,d)=1\), dann liefern diese Formeln die primitiven Rahmen, die in der Suche benötigt werden, und es gilt
$$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.$$
Setzt man
$$N_0=a^2+b^2+c^2+d^2,\qquad L_0=\frac{N_0}{g},$$
dann ist \(L_0\) genau die euklidische Kantenlänge des primitiven Würfelrahmens.
Für jede ganze Skalierung \(s\ge 1\) erhält man einen größeren Würfel mit Kanten
$$su,\qquad sv,\qquad sw,$$
und Kantenlänge
$$L=sL_0.$$
Für einen Gittervektor \(x=(x_1,x_2,x_3)\) schreiben wir
$$\delta(x)=\gcd(|x_1|,|x_2|,|x_3|).$$
Dann enthält das Segment von \(0\) nach \(sx\) genau \(s\,\delta(x)+1\) Gitterpunkte. Für die drei primitiven Kanten definieren wir
$$\Delta_0=\delta(u)+\delta(v)+\delta(w),\qquad \Delta=s\Delta_0.$$
Diese Größen messen die Gitter-Schrittweiten in den drei Kantenrichtungen und tauchen unmittelbar im Punktzählpolynom auf.
Für einen Gitter-Parallelepiped mit Kantenvektoren \(p,q,r\in\mathbb{Z}^3\) gilt die Standardformel
$$\#\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.$$
Hier setzen wir \(p=su\), \(q=sv\), \(r=sw\). Da die Kanten paarweise orthogonal sind und alle dieselbe Länge \(L\) haben, folgt
$$|\det(su,sv,sw)|=L^3.$$
Außerdem gilt
$$ (su)\times(sv)=\pm L(sw),\qquad (sv)\times(sw)=\pm L(su),\qquad (sw)\times(su)=\pm L(sv).$$
Somit erhält man
$$\delta((su)\times(sv))=L\,\delta(sw),\qquad \delta((sv)\times(sw))=L\,\delta(su),\qquad \delta((sw)\times(su))=L\,\delta(sv).$$
Die drei Flächenterme ergeben also zusammen \(L\Delta\), und die drei Kantenterme ergeben \(\Delta\). Damit lautet die exakte Gitterpunktzahl eines skalierten Würfels
$$\boxed{\operatorname{Punkte}(s)=L^3+(L+1)\Delta+1.}$$
Zur Kontrolle: Beim achsenparallelen Würfel der Seitenlänge \(m\) gilt \(L=m\) und \(\Delta=3m\), also wird daraus \((m+1)^3\), genau die bekannte Anzahl.
Für den primitiven Rahmen definieren wir die Achsenspannen
$$\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|.$$
Das sind die Kantenlängen der kleinsten achsenparallelen Box, die den primitiven Würfel enthält. Der Grund ist, dass jede Eckkoordinate eine Teilsumme der Kantenkoordinaten ist und der gesamte Wertebereich entlang einer Achse daher die Summe der Absolutbeträge ist.
Nach Skalierung mit \(s\) werden daraus \(s\sigma_x\), \(s\sigma_y\), \(s\sigma_z\). Die Zahl der zulässigen ganzzahligen Translationen ist daher
$$T(s)=(n+1-s\sigma_x)(n+1-s\sigma_y)(n+1-s\sigma_z),$$
sofern alle Faktoren positiv bleiben. Äquivalent ist die maximal erlaubte Skalierung
$$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.$$
Der rohe Beitrag eines primitiven Rahmens ist
$$\sum_{s=1}^{s_{\max}} T(s)\left((sL_0)^3+(sL_0+1)s\Delta_0+1\right).$$
Die Quaternion-Parametrisierung zählt geordnete orthogonale Kantentripel und nicht unlabeled geometrische Würfel. Derselbe geometrische Würfel erscheint in \(24\) äquivalenten Orientierungen, entsprechend der Rotationssymmetrie des Würfels. Bezeichnet \(R(n)\) die rohe Gesamtsumme über alle primitiven Rahmen, so gilt daher
$$\boxed{S(n)=\frac{R(n)}{24}.}$$
Genau deshalb prüft die Implementierung vor der Ausgabe, dass die rohe Summe durch \(24\) teilbar ist.
Betrachte den primitiven Achsenrahmen
$$u=(1,0,0),\qquad v=(0,1,0),\qquad w=(0,0,1).$$
Dann sind \(L_0=1\), \(\Delta_0=3\) und \(\sigma_x=\sigma_y=\sigma_z=1\), also \(s_{\max}=2\).
Für \(s=1\) erhält man
$$T(1)=2^3=8,\qquad \operatorname{Punkte}(1)=1^3+(1+1)\cdot 3+1=8,$$
also den Beitrag \(8\cdot 8=64\).
Für \(s=2\) erhält man
$$T(2)=1,\qquad \operatorname{Punkte}(2)=2^3+(2+1)\cdot 6+1=27,$$
also den Beitrag \(27\).
Größere Skalen passen nicht mehr, also
$$S(2)=64+27=91,$$
genau der Kontrollwert der Implementierungen.
Die Implementierungen in C++, Python und Java verwenden alle dieselbe zugrunde liegende Suche. Der Kernsolver zerlegt die Quaternion-Enumeration in Paritätsmuster, weil die Parität sowohl den Normierungsdivisor \(g\) als auch die Schranke \(N_0\le gn\) bestimmt. Dadurch werden nur die tatsächlich relevanten ganzzahligen Werte vorbereitet.
Während der Enumeration werden partielle Summen von \(a^2+b^2+c^2+d^2\) früh geprüft, sodass Zweige oberhalb der Schranke sofort wegfallen. Zusätzlich werden nicht primitive Vierer durch ggT-Tests verworfen, und kanonische Vorzeichenregeln vermeiden die triviale Doppelzählung durch \((a,b,c,d)\) und \(-(a,b,c,d)\).
Für jeden verbleibenden Rahmen berechnet die Implementierung die drei Kantenvektoren, verwirft Rahmen mit zu großen Achsenspannen, bestimmt \(L_0\), \(\Delta_0\) und \(s_{\max}\) und durchläuft dann alle zulässigen Skalen \(s\). Auf jeder Skala wird die Zahl der Translationen mit der Punktzahl des Würfels multipliziert, zu einem threadlokalen Akkumulator addiert und erst nach Abschluss aller Threads die rohe Gesamtsumme durch \(24\) geteilt.
Die Python- und Java-Versionen sind nur schlanke Hüllen um denselben kompilierten Solver. Deshalb teilen alle drei Sprachen dieselbe Mathematik, dieselben Kontrollwerte und denselben Endwert.
Der dominierende Aufwand ist die Quaternion-Suche. Eine grobe geometrische Obergrenze ergibt sich aus der vierdimensionalen Kugel \(a^2+b^2+c^2+d^2\le 4n\), die \(O(n^2)\) ganzzahlige Tupel enthält. Als lockere Worst-Case-Abschätzung ist der Suchraum also ungefähr quadratisch in \(n\).
In der Praxis ist die Laufzeit deutlich kleiner als diese grobe Schranke. Partielle Normprüfungen schneiden viele Zweige früh ab, ggT-Filter entfernen nicht primitive Tupel, Spannweitenprüfungen verwerfen unpassende Rahmen, und jeder überlebende Rahmen läuft nur bis \(s_{\max}\). Die Methode ist daher stark von Pruning bestimmt und keineswegs eine naive Enumeration aller Würfel.
Der Speicherbedarf bleibt überschaubar. Die vorberechneten geraden und ungeraden Wertelisten reichen nur bis Größenordnung \(O(\sqrt{n})\), und die parallele Version speichert im Wesentlichen nur Task-Metadaten sowie einen Akkumulator pro Worker-Thread. Im Vergleich zur arithmetischen Arbeit ist der Speicherverbrauch daher klein.
\(S(n)\), \([0,n]^3\) kutusu içine bir öteleme ile yerleştirilebilen bütün örgü küplerin içerdiği örgü noktalarının toplamıdır. Burada örgü küp, sekiz köşesinin de \(\mathbb{Z}^3\) içinde olduğu küp anlamına gelir.
Bütün köşe kümelerini doğrudan taramak çok fazla tekrar üretir. Bunun yerine uygulama, ilkel küp çerçevelerini üretir, her çerçeve için izin verilen tüm tamsayı ölçekleri dener, kutuya sığan tamsayı ötelemeleri sayar ve ölçeklenmiş küpün tam örgü nokta sayısını toplar.
Bir küpü, tek bir köşeden çıkan üç kenar vektörüyle yazalım: \(u,v,w\in\mathbb{Z}^3\). Bunların sağlaması gereken koşullar
$$u\cdot v=u\cdot w=v\cdot w=0,\qquad u\cdot u=v\cdot v=w\cdot w.$$
Temel fikir, böyle tamsayı ortogonal ve eş uzunluklu üçlüleri tamsayı kuaterniyonlardan sistematik olarak üretmek ve sonra aynı yönelimin bütün daha büyük küplerini sadece tamsayı ölçeklerle elde etmektir.
Sıfır olmayan bir \((a,b,c,d)\in\mathbb{Z}^4\) dörtlüsü için
$$\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}$$
tanımlarını yapalım. Buradaki \(g\), paritenin zorunlu kıldığı ortak çarpanı temizler:
$$g=\begin{cases} 4,& \text{eğer }a,b,c,d\text{ dördü de tekse},\\ 2,& \text{eğer bunlardan tam ikisi tekse},\\ 1,& \text{diğer durumlarda}. \end{cases}$$
\(\gcd(a,b,c,d)=1\) olduğunda bu formüller aramanın ihtiyaç duyduğu ilkel çerçeveleri verir ve
$$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$$
eşitliği sağlanır. Yani
$$N_0=a^2+b^2+c^2+d^2,\qquad L_0=\frac{N_0}{g}$$
yazarsak, \(L_0\) ilkel küpün Öklidyen kenar uzunluğudur.
Her \(s\ge 1\) tamsayısı için aynı yönelimde daha büyük bir küp elde ederiz:
$$su,\qquad sv,\qquad sw,$$
ve kenar uzunluğu
$$L=sL_0$$
olur. Bir örgü vektörü \(x=(x_1,x_2,x_3)\) için
$$\delta(x)=\gcd(|x_1|,|x_2|,|x_3|)$$
tanımını kullanalım. O zaman \(0\) ile \(sx\) arasındaki doğru parçası üzerinde tam olarak \(s\,\delta(x)+1\) örgü noktası bulunur. Üç ilkel kenar için
$$\Delta_0=\delta(u)+\delta(v)+\delta(w),\qquad \Delta=s\Delta_0$$
tanımlarını yaparız. Bu büyüklükler, üç kenar yönündeki örgü adım yoğunluğunu ölçer ve nokta sayma formülünde doğrudan görünür.
Standart bir örgü paralelyüz formülü şunu söyler: \(p,q,r\in\mathbb{Z}^3\) kenar vektörleri için
$$\#\bigl((\text{paralelepiped}(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.$$
Burada \(p=su\), \(q=sv\), \(r=sw\) alıyoruz. Kenarlar birbirine dik ve hepsinin uzunluğu \(L\) olduğundan
$$|\det(su,sv,sw)|=L^3$$
olur. Ayrıca
$$ (su)\times(sv)=\pm L(sw),\qquad (sv)\times(sw)=\pm L(su),\qquad (sw)\times(su)=\pm L(sv)$$
olduğu için
$$\delta((su)\times(sv))=L\,\delta(sw),\qquad \delta((sv)\times(sw))=L\,\delta(su),\qquad \delta((sw)\times(su))=L\,\delta(sv)$$
elde edilir. Böylece üç yüz terimi toplamda \(L\Delta\), üç kenar terimi de \(\Delta\) verir. Sonuç olarak ölçeklenmiş küpün tam örgü nokta sayısı
$$\boxed{\operatorname{nokta}(s)=L^3+(L+1)\Delta+1.}$$
Kontrol için eksenlere paralel, kenarı \(m\) olan küpte \(L=m\) ve \(\Delta=3m\) olur; formül \((m+1)^3\) verir. Bu da beklenen klasik sonuçtur.
İlkel çerçeve için eksen yayılımlarını
$$\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|$$
şeklinde tanımlayalım. Bunlar, ilkel küpü saran en küçük eksen-paralel kutunun kenar uzunluklarıdır. Sebebi, her köşe koordinatının kenar koordinatlarının bir alt toplamı olması ve dolayısıyla bir eksendeki toplam aralığın mutlak değerler toplamına eşit olmasıdır.
\(s\) ile ölçeklendikten sonra yayılımlar \(s\sigma_x\), \(s\sigma_y\), \(s\sigma_z\) olur. Bu yüzden geçerli tamsayı öteleme sayısı
$$T(s)=(n+1-s\sigma_x)(n+1-s\sigma_y)(n+1-s\sigma_z)$$
olur; tabii bütün çarpanlar pozitif kalmalıdır. Buna eşdeğer olarak en büyük ölçek
$$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$$
şeklindedir.
Tek bir ilkel çerçevenin ham katkısı
$$\sum_{s=1}^{s_{\max}} T(s)\left((sL_0)^3+(sL_0+1)s\Delta_0+1\right)$$
olur. Kuaterniyon parametrizasyonu, etiketsiz geometrik küpleri değil, sıralı ortogonal kenar üçlülerini sayar. Aynı geometrik küp, küpün dönme simetri grubu nedeniyle \(24\) eşdeğer yönelimle görünür. Bu nedenle bütün ilkel çerçeveler üzerindeki ham toplam \(R(n)\) ise
$$\boxed{S(n)=\frac{R(n)}{24}.}$$
Uygulamanın son aşamada ham toplamın \(24\)'e bölünebildiğini kontrol etmesinin nedeni budur.
Eksenlere paralel ilkel çerçeveyi alalım:
$$u=(1,0,0),\qquad v=(0,1,0),\qquad w=(0,0,1).$$
Burada \(L_0=1\), \(\Delta_0=3\) ve \(\sigma_x=\sigma_y=\sigma_z=1\), dolayısıyla \(s_{\max}=2\).
\(s=1\) için
$$T(1)=2^3=8,\qquad \operatorname{nokta}(1)=1^3+(1+1)\cdot 3+1=8$$
ve katkı \(8\cdot 8=64\) olur.
\(s=2\) için
$$T(2)=1,\qquad \operatorname{nokta}(2)=2^3+(2+1)\cdot 6+1=27$$
ve katkı \(27\)'dir.
Daha büyük bir ölçek sığmadığı için
$$S(2)=64+27=91$$
elde edilir; bu da uygulamadaki kontrol değeriyle birebir aynıdır.
C++, Python ve Java uygulamalarının hepsi aynı temel aramayı kullanır. Çekirdek çözücü, kuaterniyon taramasını parite desenlerine ayırır; çünkü hem normlaştırma böleni \(g\) hem de \(N_0\le gn\) üst sınırı pariteye bağlıdır. Böylece her durumda yalnızca gerçekten gerekebilecek tamsayı değerleri hazırlanır.
Tarama sırasında \(a^2+b^2+c^2+d^2\) için kısmi toplamlar erkenden kontrol edilir; sınırı aşan dallar anında budanır. Ayrıca ebob testleriyle ilkel olmayan dörtlüler atılır ve \((a,b,c,d)\) ile \(-(a,b,c,d)\) arasındaki anlamsız tekrarlar kanonik işaret kurallarıyla engellenir.
Ayakta kalan her çerçeve için uygulama üç kenar vektörünü hesaplar, kutuya zaten sığamayacak eksen yayılımlı çerçeveleri eler, \(L_0\), \(\Delta_0\) ve \(s_{\max}\) değerlerini çıkarır, sonra da bütün geçerli \(s\) ölçeklerini dolaşır. Her ölçekte öteleme sayısı ile küpün örgü nokta sayısı çarpılır, thread-yerel toplama eklenir ve bütün iş parçacıkları bittikten sonra ham toplam \(24\)'e bölünür.
Python ve Java sürümleri, aynı derlenmiş çözücünün ince sarmalayıcılarıdır. Bu yüzden üç dilde de aynı matematik, aynı kontrol noktaları ve aynı nihai sonuç kullanılır.
Baskın maliyet kuaterniyon aramasıdır. Kaba bir geometrik üst sınır, \(a^2+b^2+c^2+d^2\le 4n\) dört boyutlu küresi üzerinden gelir; bu bölge \(O(n^2)\) kadar tamsayı dörtlüsü içerir. Dolayısıyla gevşek bir en kötü durum bakışında arama uzayı yaklaşık olarak \(n\)'de ikinci derecedir.
Pratikte süre bu kaba sınırdan belirgin biçimde düşüktür. Kısmi norm kontrolleri çok sayıda dalı erken keser, ebob filtreleri ilkel olmayan dörtlüleri temizler, yayılım kontrolleri kutuya asla sığmayacak çerçeveleri atar ve her yaşayan çerçeve yalnızca \(s_{\max}\)'e kadar ilerler. Bu yüzden yöntem, bütün küpleri körlemesine saymaktan çok yoğun budamaya dayalıdır.
Bellek kullanımı düşüktür. Önceden oluşturulan çift ve tek değer tabloları yalnızca \(O(\sqrt{n})\) büyüklüğüne kadar gider; paralel sürüm de esas olarak görev meta verilerini ve her iş parçacığı için bir toplama değişkenini saklar. Dolayısıyla bellek maliyeti aritmetik işe kıyasla küçüktür.
Sea \(S(n)\) la suma, sobre todos los cubos de red que pueden trasladarse para quedar contenidos en \([0,n]^3\), del número de puntos de red contenidos en cada cubo. Un cubo de red significa aquí que sus ocho vértices pertenecen a \(\mathbb{Z}^3\).
Una búsqueda directa sobre conjuntos de vértices sería extremadamente redundante. La implementación enumera en cambio marcos primitivos de cubo, los amplía por todas las escalas enteras admisibles, cuenta cuántas traslaciones enteras caben dentro de la caja y suma el número exacto de puntos de red del cubo escalado.
Describimos un cubo mediante tres vectores de arista \(u,v,w\in\mathbb{Z}^3\) que salen de un mismo vértice. Deben cumplir
$$u\cdot v=u\cdot w=v\cdot w=0,\qquad u\cdot u=v\cdot v=w\cdot w.$$
La idea clave es que tales ternas enteras, ortogonales y de igual longitud, pueden generarse sistemáticamente a partir de cuaterniones enteros. Una vez conocido un marco primitivo, todos los cubos mayores con la misma orientación aparecen simplemente multiplicando por escalas enteras.
Partimos de enteros \(a,b,c,d\), no todos nulos, y definimos
$$\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}$$
El divisor \(g\) elimina el factor común obligado por la paridad:
$$g=\begin{cases} 4,& \text{si }a,b,c,d\text{ son todos impares},\\ 2,& \text{si exactamente dos de ellos son impares},\\ 1,& \text{en los demás casos}. \end{cases}$$
Si además \(\gcd(a,b,c,d)=1\), estas fórmulas producen los marcos primitivos necesarios en la búsqueda y se cumple
$$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.$$
Si escribimos
$$N_0=a^2+b^2+c^2+d^2,\qquad L_0=\frac{N_0}{g},$$
entonces \(L_0\) es la longitud euclídea de la arista del cubo primitivo.
Para cada escala entera \(s\ge 1\), obtenemos un cubo mayor con aristas
$$su,\qquad sv,\qquad sw,$$
y longitud de arista
$$L=sL_0.$$
Para un vector de red \(x=(x_1,x_2,x_3)\), definimos
$$\delta(x)=\gcd(|x_1|,|x_2|,|x_3|).$$
Entonces el segmento entre \(0\) y \(sx\) contiene exactamente \(s\,\delta(x)+1\) puntos de red. Para las tres aristas primitivas definimos
$$\Delta_0=\delta(u)+\delta(v)+\delta(w),\qquad \Delta=s\Delta_0.$$
Estas cantidades miden el paso de red a lo largo de las tres direcciones de arista y son precisamente las que aparecen en el polinomio de conteo.
Una fórmula estándar para un paralelepípedo de red con aristas \(p,q,r\in\mathbb{Z}^3\) dice
$$\#\bigl((\text{paralelepipedo}(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.$$
Aquí tomamos \(p=su\), \(q=sv\), \(r=sw\). Como las aristas son mutuamente ortogonales y todas tienen longitud \(L\), resulta
$$|\det(su,sv,sw)|=L^3.$$
Además,
$$ (su)\times(sv)=\pm L(sw),\qquad (sv)\times(sw)=\pm L(su),\qquad (sw)\times(su)=\pm L(sv).$$
Por tanto,
$$\delta((su)\times(sv))=L\,\delta(sw),\qquad \delta((sv)\times(sw))=L\,\delta(su),\qquad \delta((sw)\times(su))=L\,\delta(sv).$$
La suma de los tres términos de cara es \(L\Delta\), y la suma de los tres términos de arista es \(\Delta\). Así, el número exacto de puntos de red del cubo escalado es
$$\boxed{\operatorname{puntos}(s)=L^3+(L+1)\Delta+1.}$$
Como verificación, en el cubo alineado con los ejes de lado \(m\) se tiene \(L=m\) y \(\Delta=3m\), con lo cual la fórmula se convierte en \((m+1)^3\), exactamente el valor esperado.
Para el marco primitivo definimos las extensiones por eje
$$\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|.$$
Estas son las longitudes del menor contenedor alineado con los ejes que encierra al cubo primitivo. La razón es que cada coordenada de un vértice es una suma parcial de las coordenadas de las aristas, de modo que el rango total sobre un eje es la suma de los valores absolutos en ese eje.
Tras escalar por \(s\), las extensiones pasan a ser \(s\sigma_x\), \(s\sigma_y\), \(s\sigma_z\). El número de traslaciones enteras válidas es entonces
$$T(s)=(n+1-s\sigma_x)(n+1-s\sigma_y)(n+1-s\sigma_z),$$
siempre que todos los factores sean positivos. De forma equivalente, la escala máxima permitida es
$$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.$$
La contribución bruta de un marco primitivo es
$$\sum_{s=1}^{s_{\max}} T(s)\left((sL_0)^3+(sL_0+1)s\Delta_0+1\right).$$
La parametrización por cuaterniones enumera ternas ordenadas de aristas ortogonales, no cubos geométricos sin etiquetas. El mismo cubo geométrico aparece en \(24\) orientaciones equivalentes, correspondientes al grupo de simetrías rotacionales del cubo. Si \(R(n)\) es la suma bruta total sobre todos los marcos primitivos, entonces
$$\boxed{S(n)=\frac{R(n)}{24}.}$$
Por eso la implementación comprueba antes de terminar que la suma bruta sea divisible entre \(24\).
Tomemos el marco axial primitivo
$$u=(1,0,0),\qquad v=(0,1,0),\qquad w=(0,0,1).$$
Entonces \(L_0=1\), \(\Delta_0=3\) y \(\sigma_x=\sigma_y=\sigma_z=1\), así que \(s_{\max}=2\).
Para \(s=1\), tenemos
$$T(1)=2^3=8,\qquad \operatorname{puntos}(1)=1^3+(1+1)\cdot 3+1=8,$$
y la contribución es \(8\cdot 8=64\).
Para \(s=2\), tenemos
$$T(2)=1,\qquad \operatorname{puntos}(2)=2^3+(2+1)\cdot 6+1=27,$$
y la contribución es \(27\).
No cabe ninguna escala mayor, de modo que
$$S(2)=64+27=91,$$
en acuerdo exacto con el punto de control usado por las implementaciones.
Las implementaciones en C++, Python y Java utilizan la misma búsqueda subyacente. El solucionador principal separa la enumeración de cuaterniones por patrones de paridad, porque la paridad determina tanto el divisor de normalización \(g\) como la cota \(N_0\le gn\). Así solo se preparan los valores enteros que realmente pueden intervenir en cada caso.
Durante la enumeración se revisan pronto las sumas parciales de \(a^2+b^2+c^2+d^2\), de forma que las ramas que ya exceden la cota se descartan inmediatamente. Además, se eliminan las cuádruplas no primitivas con pruebas de mcd y se aplican restricciones canónicas de signo para evitar la duplicación trivial entre \((a,b,c,d)\) y \(-(a,b,c,d)\).
Para cada marco que sobrevive, la implementación calcula los tres vectores de arista, descarta los marcos cuya extensión axial ya supera la caja, obtiene \(L_0\), \(\Delta_0\) y \(s_{\max}\), y luego recorre todas las escalas válidas \(s\). En cada escala multiplica el número de traslaciones por la fórmula exacta de puntos del cubo, acumula el resultado en un total local por hilo y, una vez terminan todos los hilos, divide la suma bruta entre \(24\).
Las versiones de Python y Java son envoltorios ligeros sobre el mismo solucionador compilado. Por eso las tres comparten la misma matemática, los mismos puntos de control y el mismo valor final.
El coste dominante es la búsqueda de cuaterniones. Una cota geométrica grosera proviene de la bola de cuatro dimensiones \(a^2+b^2+c^2+d^2\le 4n\), que contiene \(O(n^2)\) tuplas enteras. En una visión conservadora, el espacio de búsqueda es por tanto aproximadamente cuadrático en \(n\).
En la práctica el tiempo real es sensiblemente menor que esa cota grosera. Las comprobaciones parciales de norma cortan muchas ramas pronto, los filtros por mcd eliminan tuplas no primitivas, las comprobaciones de extensión descartan marcos que nunca cabrán en la caja y cada marco superviviente solo itera hasta \(s_{\max}\). El método está por tanto guiado por una poda fuerte, no por un barrido ingenuo de todos los cubos.
El uso de memoria es reducido. Las tablas precomputadas de valores pares e impares solo llegan hasta magnitud \(O(\sqrt{n})\), y la versión paralela almacena esencialmente metadatos de tareas y un acumulador por hilo de trabajo. El coste de memoria queda pequeño frente al trabajo aritmético.
记 \(S(n)\) 为这样一个总和：把所有能够通过平移完整放进 \([0,n]^3\) 的晶格立方体都取出来，对每个立方体统计其包含的晶格点个数，再把这些数全部相加。这里的“晶格立方体”是指八个顶点都落在 \(\mathbb{Z}^3\) 上的立方体。
如果直接按顶点集合去枚举，重复会非常严重。实现采用的做法是先枚举“原始立方体框架”，再对每个框架枚举所有允许的整数放大倍数，计算每个放大后立方体在盒子中的整数平移个数，并乘上该立方体本身精确的晶格点数。
一个立方体可以由同一顶点发出的三条边向量 \(u,v,w\in\mathbb{Z}^3\) 描述。它们必须满足
$$u\cdot v=u\cdot w=v\cdot w=0,\qquad u\cdot u=v\cdot v=w\cdot w.$$
也就是说，这三条边必须两两正交且长度相同。关键点在于：这样的整数向量三元组可以由整数四元数系统地生成出来；而一旦找到一个原始框架，同方向的更大立方体都只需做整数倍缩放即可得到。
从整数 \(a,b,c,d\) 出发，要求它们不全为零，并定义
$$\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}$$
这里的 \(g\) 不是随意加入的，而是用来消去由奇偶性强制产生的公共因子：
$$g=\begin{cases} 4,& \text{如果 }a,b,c,d\text{ 全都是奇数},\\ 2,& \text{如果四个数中恰有两个是奇数},\\ 1,& \text{其他情况}. \end{cases}$$
再加上 \(\gcd(a,b,c,d)=1\) 这一原始性条件，上述公式就会产生程序真正需要枚举的“原始立方体框架”。并且可以直接验证
$$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.$$
因此若记
$$N_0=a^2+b^2+c^2+d^2,\qquad L_0=\frac{N_0}{g},$$
那么 \(L_0\) 就是这个原始立方体框架的欧几里得边长。
对任意整数 \(s\ge 1\)，同一方向上更大的立方体边向量就是
$$su,\qquad sv,\qquad sw,$$
对应的边长为
$$L=sL_0.$$
对于一个晶格向量 \(x=(x_1,x_2,x_3)\)，定义
$$\delta(x)=\gcd(|x_1|,|x_2|,|x_3|).$$
这样一来，从 \(0\) 到 \(sx\) 的线段上恰好有 \(s\,\delta(x)+1\) 个晶格点。对三条原始边，记
$$\Delta_0=\delta(u)+\delta(v)+\delta(w),\qquad \Delta=s\Delta_0.$$
\(\Delta_0\) 和 \(\Delta\) 描述的是三条边方向上的晶格步长总量，而这正是后面点数公式里会出现的量。
对任意晶格平行六面体，只要三条边向量是 \(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.$$
现在把它应用到 \(p=su\)、\(q=sv\)、\(r=sw\)。由于三条边两两垂直且长度都等于 \(L\)，体积项立刻变成
$$|\det(su,sv,sw)|=L^3.$$
另一方面，正交且等长意味着
$$ (su)\times(sv)=\pm L(sw),\qquad (sv)\times(sw)=\pm L(su),\qquad (sw)\times(su)=\pm L(sv).$$
因此三个面项分别满足
$$\delta((su)\times(sv))=L\,\delta(sw),\qquad \delta((sv)\times(sw))=L\,\delta(su),\qquad \delta((sw)\times(su))=L\,\delta(sv).$$
把三项加起来就是 \(L\Delta\)。而三条边本身的晶格项之和正好是 \(\Delta\)。所以缩放后一个立方体的精确晶格点个数为
$$\boxed{\operatorname{points}(s)=L^3+(L+1)\Delta+1.}$$
这个公式可以用轴对齐立方体做快速自检。若边长为 \(m\)，则 \(L=m\)，\(\Delta=3m\)，上式就化成 \((m+1)^3\)，和最熟悉的轴对齐结果完全一致。
对原始框架，定义三个坐标轴方向的跨度
$$\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|.$$
它们就是包含该原始立方体的最小轴对齐包围盒在三个方向上的边长。原因是：立方体八个顶点在某一坐标轴上的坐标，都是三条边该坐标分量的若干子和，因此该轴上的总范围恰好等于绝对值之和。
把立方体放大到比例 \(s\) 后，三个跨度变成 \(s\sigma_x\)、\(s\sigma_y\)、\(s\sigma_z\)。于是可行的整数平移数就是
$$T(s)=(n+1-s\sigma_x)(n+1-s\sigma_y)(n+1-s\sigma_z),$$
前提是这三个因子都仍然为正。等价地，允许的最大缩放倍数为
$$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.$$
对于一个固定的原始框架，其“原始贡献”是
$$\sum_{s=1}^{s_{\max}} T(s)\left((sL_0)^3+(sL_0+1)s\Delta_0+1\right).$$
但这里还有一个重要的重计数问题。四元数参数化枚举的是“有序的正交边三元组”，而不是“不带标签的几何立方体”。同一个几何立方体会以 \(24\) 种等价的方向出现，这正对应于立方体的旋转对称群。所以如果把所有原始框架的原始贡献都加起来得到 \(R(n)\)，那么真正需要的答案是
$$\boxed{S(n)=\frac{R(n)}{24}.}$$
这也解释了为什么实现会在最后先检查原始总和是否能被 \(24\) 整除，再输出最终结果。
取最简单的轴对齐原始框架
$$u=(1,0,0),\qquad v=(0,1,0),\qquad w=(0,0,1).$$
此时 \(L_0=1\)，\(\Delta_0=3\)，且 \(\sigma_x=\sigma_y=\sigma_z=1\)，因此 \(s_{\max}=2\)。
当 \(s=1\) 时，
$$T(1)=2^3=8,\qquad \operatorname{points}(1)=1^3+(1+1)\cdot 3+1=8,$$
所以这一层的贡献是 \(8\cdot 8=64\)。
当 \(s=2\) 时，
$$T(2)=1,\qquad \operatorname{points}(2)=2^3+(2+1)\cdot 6+1=27,$$
所以贡献是 \(27\)。
更大的尺度已经放不进盒子，于是
$$S(2)=64+27=91,$$
这正好与实现中的校验值一致。
C++、Python 和 Java 三个实现共享同一套底层搜索逻辑。核心求解器会先按照奇偶模式拆分四元数参数空间，因为奇偶性同时决定了归一化因子 \(g\) 和范数上界 \(N_0\le gn\)。这样每一种情况只需要准备真正可能出现的整数取值范围。
在枚举过程中，\(a^2+b^2+c^2+d^2\) 的部分和会被尽早检查；一旦某个分支已经超过上界，就立刻停止向下展开。接着还会用 gcd 条件剔除非原始四元组，并施加规范化的符号限制，避免 \((a,b,c,d)\) 与 \(-(a,b,c,d)\) 这种显然等价的重复。
对每一个通过筛选的框架，实现都会计算三条边向量，先排除那些轴向跨度已经超过盒子的候选，再求出 \(L_0\)、\(\Delta_0\) 和 \(s_{\max}\)。随后程序遍历所有可行的缩放 \(s\)，在每个尺度上把平移个数与立方体点数公式相乘，加到线程局部累加器中。等所有线程结束后，才把原始总和统一除以 \(24\)。
Python 与 Java 版本本质上只是同一个已编译求解器的轻量包装，因此三种语言使用完全相同的数学推导、相同的校验点以及相同的最终答案。
主导成本来自四元数搜索。一个保守但比较粗的上界，是考虑四维球体 \(a^2+b^2+c^2+d^2\le 4n\) 中的整数点个数；该数量级为 \(O(n^2)\)。因此从几何体积的角度看，候选参数空间可以视为大约二次规模。
不过真实运行时间通常明显小于这个粗上界。部分范数判断会很早砍掉大量分支，gcd 过滤会剔除非原始参数，跨度判断会直接排除根本放不进盒子的框架，而每个保留下来的框架也只会循环到 \(s_{\max}\)。因此这套方法的效率主要来自强剪枝，而不是机械地枚举所有立方体。
内存占用相对温和。预生成的奇数表和偶数表只需要覆盖到 \(O(\sqrt{n})\) 的规模；并行部分也主要保存任务信息以及每个工作线程自己的累加器。所以与整体算术工作相比，内存开销并不大。
Пусть \(S(n)\) обозначает сумму числа узлов решетки по всем решетчатым кубам, которые можно сдвигом целочисленного вектора целиком поместить внутрь \([0,n]^3\). Под решетчатым кубом понимается куб, все восемь вершин которого лежат в \(\mathbb{Z}^3\).
Прямой перебор по наборам вершин привел бы к огромному числу повторов. Поэтому реализация перечисляет примитивные каркасы куба, затем для каждого каркаса перебирает все допустимые целые масштабы, считает число допустимых целочисленных сдвигов внутри коробки и добавляет точное число решетчатых точек в соответствующем масштабированном кубе.
Один куб можно описать тремя векторами ребер \(u,v,w\in\mathbb{Z}^3\), выходящими из одной вершины. Они должны удовлетворять условиям
$$u\cdot v=u\cdot w=v\cdot w=0,\qquad u\cdot u=v\cdot v=w\cdot w.$$
Иными словами, ребра должны быть попарно ортогональны и иметь одинаковую длину. Ключевой факт состоит в том, что такие целочисленные тройки можно систематически генерировать с помощью целочисленных кватернионов. Как только найден один примитивный каркас, все большие кубы той же ориентации появляются простым умножением на целый масштаб.
Возьмем целые числа \(a,b,c,d\), не все нулевые, и определим
$$\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}$$
Делитель \(g\) убирает общий множитель, который неизбежно возникает из-за четности:
$$g=\begin{cases} 4,& \text{если }a,b,c,d\text{ все нечетные},\\ 2,& \text{если ровно два из них нечетные},\\ 1,& \text{в остальных случаях}. \end{cases}$$
При условии \(\gcd(a,b,c,d)=1\) эти формулы дают именно те примитивные каркасы, которые нужны поиску, и выполняется
$$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.$$
Если ввести обозначения
$$N_0=a^2+b^2+c^2+d^2,\qquad L_0=\frac{N_0}{g},$$
то \(L_0\) есть евклидова длина ребра примитивного каркаса куба.
Для любого целого масштаба \(s\ge 1\) получаем больший куб с ребрами
$$su,\qquad sv,\qquad sw,$$
и длиной ребра
$$L=sL_0.$$
Для решетчатого вектора \(x=(x_1,x_2,x_3)\) введем
$$\delta(x)=\gcd(|x_1|,|x_2|,|x_3|).$$
Тогда отрезок от \(0\) до \(sx\) содержит ровно \(s\,\delta(x)+1\) решетчатых точек. Для трех примитивных ребер положим
$$\Delta_0=\delta(u)+\delta(v)+\delta(w),\qquad \Delta=s\Delta_0.$$
Эти величины измеряют шаг решетки вдоль трех реберных направлений и именно они появляются в формуле для числа точек.
Для решетчатого параллелепипеда с ребрами \(p,q,r\in\mathbb{Z}^3\) существует стандартная формула
$$\#\bigl((\text{параллелепипед}(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.$$
Применим ее к \(p=su\), \(q=sv\), \(r=sw\). Так как ребра попарно ортогональны и все имеют длину \(L\), объемный член равен
$$|\det(su,sv,sw)|=L^3.$$
Кроме того,
$$ (su)\times(sv)=\pm L(sw),\qquad (sv)\times(sw)=\pm L(su),\qquad (sw)\times(su)=\pm L(sv).$$
Следовательно,
$$\delta((su)\times(sv))=L\,\delta(sw),\qquad \delta((sv)\times(sw))=L\,\delta(su),\qquad \delta((sw)\times(su))=L\,\delta(sv).$$
Сумма трех гранных членов равна \(L\Delta\), а сумма трех реберных членов равна \(\Delta\). Поэтому точное число решетчатых точек в масштабированном кубе равно
$$\boxed{\operatorname{points}(s)=L^3+(L+1)\Delta+1.}$$
Проверка на осевом кубе проста: если сторона равна \(m\), то \(L=m\) и \(\Delta=3m\), а формула превращается в \((m+1)^3\), что полностью совпадает с очевидным результатом.
Для примитивного каркаса введем осевые размахи
$$\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|.$$
Это длины ребер наименьшего осевого контейнера, содержащего примитивный куб. Причина в том, что координата любой вершины по одной оси есть сумма некоторых координат ребер, и полный диапазон значений вдоль оси равен сумме модулей.
После масштабирования на \(s\) размахи становятся равны \(s\sigma_x\), \(s\sigma_y\), \(s\sigma_z\). Поэтому число допустимых целочисленных сдвигов равно
$$T(s)=(n+1-s\sigma_x)(n+1-s\sigma_y)(n+1-s\sigma_z),$$
если все множители положительны. Эквивалентно максимальный допустимый масштаб равен
$$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.$$
Грубый вклад одного примитивного каркаса имеет вид
$$\sum_{s=1}^{s_{\max}} T(s)\left((sL_0)^3+(sL_0+1)s\Delta_0+1\right).$$
Однако кватернионная параметризация перечисляет упорядоченные тройки ортогональных ребер, а не геометрические кубы без меток. Один и тот же геометрический куб возникает в \(24\) эквивалентных ориентациях, что соответствует группе вращательных симметрий куба. Поэтому, если \(R(n)\) обозначает суммарный грубый вклад по всем примитивным каркасам, то окончательный ответ равен
$$\boxed{S(n)=\frac{R(n)}{24}.}$$
Именно поэтому реализация перед выводом убеждается, что сырая сумма делится на \(24\).
Возьмем примитивный осевой каркас
$$u=(1,0,0),\qquad v=(0,1,0),\qquad w=(0,0,1).$$
Тогда \(L_0=1\), \(\Delta_0=3\), а \(\sigma_x=\sigma_y=\sigma_z=1\), откуда \(s_{\max}=2\).
Для \(s=1\) получаем
$$T(1)=2^3=8,\qquad \operatorname{points}(1)=1^3+(1+1)\cdot 3+1=8,$$
то есть вклад \(8\cdot 8=64\).
Для \(s=2\) получаем
$$T(2)=1,\qquad \operatorname{points}(2)=2^3+(2+1)\cdot 6+1=27,$$
то есть вклад \(27\).
Более крупный масштаб уже не помещается, поэтому
$$S(2)=64+27=91,$$
и это в точности совпадает с контрольным значением, используемым реализациями.
Реализации на C++, Python и Java используют один и тот же базовый поиск. Основной решатель разбивает перечисление кватернионов по шаблонам четности, потому что четность определяет и нормализующий делитель \(g\), и ограничение \(N_0\le gn\). Благодаря этому для каждого случая заранее подготавливаются только те целые значения, которые действительно могут понадобиться.
Во время перебора частичные суммы \(a^2+b^2+c^2+d^2\) проверяются как можно раньше, поэтому ветви, уже вышедшие за предел, немедленно отсекаются. Далее с помощью условий на gcd удаляются непримитивные четверки, а канонические ограничения на знаки устраняют тривиальное дублирование между \((a,b,c,d)\) и \(-(a,b,c,d)\).
Для каждого уцелевшего каркаса реализация вычисляет три реберных вектора, сразу отбрасывает каркасы со слишком большим осевым размахом, находит \(L_0\), \(\Delta_0\) и \(s_{\max}\), а затем проходит по всем допустимым масштабам \(s\). На каждом масштабе число сдвигов умножается на точную формулу числа точек в кубе, результат добавляется в потоковый локальный аккумулятор, и только после завершения всех потоков грубая сумма делится на \(24\).
Версии на Python и Java являются тонкими оболочками над тем же скомпилированным решателем. Поэтому все три языка используют одну и ту же математику, одни и те же контрольные точки и одно и то же окончательное значение.
Главная стоимость связана с перебором кватернионов. Грубая геометрическая верхняя оценка получается из четырехмерного шара \(a^2+b^2+c^2+d^2\le 4n\), содержащего \(O(n^2)\) целочисленных четверок. В таком консервативном приближении пространство поиска имеет примерно квадратичный размер по \(n\).
На практике время работы заметно меньше этой грубой оценки. Проверки частичной нормы очень рано обрубают множество ветвей, фильтры gcd убирают непримитивные параметры, проверки осевых размахов отбрасывают каркасы, которые никогда не поместятся в коробке, а каждый оставшийся каркас итерируется только до \(s_{\max}\). Поэтому производительность определяется в первую очередь агрессивным отсечением, а не наивным перебором всех кубов.
Память расходуется умеренно. Предварительно подготовленные таблицы четных и нечетных значений доходят только до масштаба \(O(\sqrt{n})\), а параллельная часть хранит главным образом метаданные задач и по одному аккумулятору на рабочий поток. По сравнению с объемом арифметической работы это небольшой расход памяти.
لنرمز بـ \(S(n)\) إلى مجموع أعداد نقاط الشبكة في جميع المكعبات الشبكية التي يمكن نقلها بعملية إزاحة كاملة لتقع داخل الصندوق \([0,n]^3\). والمقصود بالمكعب الشبكي هنا أن تكون رؤوسه الثمانية كلها ضمن \(\mathbb{Z}^3\).
العد المباشر على مجموعات الرؤوس يكرر الحالات نفسها مرارًا. لذلك تعتمد الخوارزمية على تعداد الأطر الأولية للمكعب، ثم توسيع كل إطار بكل معاملات القياس الصحيحة الممكنة، ثم حساب عدد الإزاحات الصحيحة التي تبقي المكعب داخل الصندوق، وأخيرًا ضرب ذلك في العدد الدقيق لنقاط الشبكة داخل المكعب بعد القياس.
يمكن وصف المكعب بثلاثة متجهات حافة \(u,v,w\in\mathbb{Z}^3\) تخرج من رأس واحد. ويجب أن تحقق
$$u\cdot v=u\cdot w=v\cdot w=0,\qquad u\cdot u=v\cdot v=w\cdot w.$$
أي إن الحواف متعامدة زوجيًا ومتساوية الطول. والفكرة الجوهرية هي أن هذه الثلاثيات الصحيحة يمكن توليدها نظاميًا من كواتيرنيونات صحيحة، ثم إن معرفة إطار أولي واحد تكفي لتوليد جميع المكعبات الأكبر ذات الاتجاه نفسه بواسطة القياس الصحيح.
نبدأ بأعداد صحيحة \(a,b,c,d\) ليست كلها صفرًا، ثم نعرّف
$$\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}$$
القاسم \(g\) يزيل العامل المشترك الذي تفرضه الفردية والزوجية:
$$g=\begin{cases} 4,& \text{all odd},\\ 2,& \text{two odd},\\ 1,& \text{otherwise}. \end{cases}$$
ومع الشرط \(\gcd(a,b,c,d)=1\)، تعطينا هذه الصيغ الأطر الأولية التي يحتاجها البحث، كما تحقق
$$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.$$
وإذا كتبنا
$$N_0=a^2+b^2+c^2+d^2,\qquad L_0=\frac{N_0}{g},$$
فإن \(L_0\) هو طول حافة الإطار الأولي في المعنى الإقليدي.
لكل عدد صحيح \(s\ge 1\)، نحصل على مكعب أكبر بالحواف
$$su,\qquad sv,\qquad sw,$$
وطول الحافة
$$L=sL_0.$$
ولأي متجه شبكي \(x=(x_1,x_2,x_3)\)، نعرّف
$$\delta(x)=\gcd(|x_1|,|x_2|,|x_3|).$$
وعندئذ يحتوي المقطع من \(0\) إلى \(sx\) على \(s\,\delta(x)+1\) نقطة شبكة بالضبط. وللحواف الأولية الثلاث نضع
$$\Delta_0=\delta(u)+\delta(v)+\delta(w),\qquad \Delta=s\Delta_0.$$
هذه الكمية تقيس كثافة الخطوات الشبكية على اتجاهات الحواف الثلاثة، وهي التي تظهر مباشرة في صيغة عد النقاط.
توجد صيغة قياسية لمتوازي المستطيلات الشبكي ذي الحواف \(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.$$
نطبّقها هنا على \(p=su\)، \(q=sv\)، \(r=sw\). وبما أن الحواف متعامدة ومتساوية الطول \(L\)، فإن حد الحجم يساوي
$$|\det(su,sv,sw)|=L^3.$$
كما أن
$$ (su)\times(sv)=\pm L(sw),\qquad (sv)\times(sw)=\pm L(su),\qquad (sw)\times(su)=\pm L(sv).$$
ومن ثم
$$\delta((su)\times(sv))=L\,\delta(sw),\qquad \delta((sv)\times(sw))=L\,\delta(su),\qquad \delta((sw)\times(su))=L\,\delta(sv).$$
وبجمع حدود الوجوه نحصل على \(L\Delta\)، وبجمع حدود الحواف نحصل على \(\Delta\). إذًا العدد الدقيق لنقاط الشبكة داخل المكعب المقاس هو
$$\boxed{\operatorname{points}(s)=L^3+(L+1)\Delta+1.}$$
ويمكن اختبار الصيغة سريعًا في حالة المكعب الموازي للمحاور ذي الضلع \(m\): عندها \(L=m\) و\(\Delta=3m\)، فتتحول الصيغة إلى \((m+1)^3\)، وهي النتيجة المعتادة تمامًا.
للإطار الأولي نعرّف امتدادات المحاور
$$\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|.$$
وهذه هي أطوال أضلاع أصغر صندوق موازٍ للمحاور يحتوي المكعب الأولي. والسبب أن إحداثي أي رأس على محور معين هو مجموع جزئي لإحداثيات الحواف، وبالتالي يكون المدى الكامل على ذلك المحور مساويًا لمجموع القيم المطلقة.
بعد القياس بمعامل \(s\)، تصبح الامتدادات \(s\sigma_x\)، \(s\sigma_y\)، \(s\sigma_z\). وعندئذ يكون عدد الإزاحات الصحيحة المسموح بها هو
$$T(s)=(n+1-s\sigma_x)(n+1-s\sigma_y)(n+1-s\sigma_z),$$
بشرط أن تبقى العوامل كلها موجبة. وبصورة مكافئة فإن أكبر قياس ممكن هو
$$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.$$
المساهمة الخام لإطار أولي واحد هي
$$\sum_{s=1}^{s_{\max}} T(s)\left((sL_0)^3+(sL_0+1)s\Delta_0+1\right).$$
لكن ترميز الكواتيرنيونات يعد ثلاثيات الحواف المتعامدة المرتبة، لا المكعبات الهندسية غير المعلَّمة. والمكعب نفسه يظهر في \(24\) اتجاهًا مكافئًا تمثل مجموعة دورانات المكعب. لذلك إذا كانت \(R(n)\) هي المحصلة الخام على جميع الأطر الأولية، فإن الجواب النهائي هو
$$\boxed{S(n)=\frac{R(n)}{24}.}$$
ولهذا تتحقق الخوارزمية في النهاية من أن المجموع الخام يقبل القسمة على \(24\) قبل إخراج النتيجة.
لنأخذ الإطار الأولي الموازي للمحاور
$$u=(1,0,0),\qquad v=(0,1,0),\qquad w=(0,0,1).$$
عندئذٍ \(L_0=1\)، و\(\Delta_0=3\)، و\(\sigma_x=\sigma_y=\sigma_z=1\)، ومن ثم \(s_{\max}=2\).
عندما \(s=1\)، نحصل على
$$T(1)=2^3=8,\qquad \operatorname{points}(1)=1^3+(1+1)\cdot 3+1=8,$$
فتكون المساهمة \(8\cdot 8=64\).
وعندما \(s=2\)، نحصل على
$$T(2)=1,\qquad \operatorname{points}(2)=2^3+(2+1)\cdot 6+1=27,$$
فتكون المساهمة \(27\).
ولا يوجد قياس أكبر يمكن أن يلائم الصندوق، لذلك
$$S(2)=64+27=91,$$
وهو تمامًا مقدار نقطة التحقق المستخدمة في التنفيذ.
تنفيذات C++ وPython وJava تعتمد جميعها على البحث الرياضي نفسه. يقسم الحل الأساسي تعداد الكواتيرنيونات بحسب نمط الفردية والزوجية، لأن هذا النمط يحدد كلًا من عامل التطبيع \(g\) والحد \(N_0\le gn\). وبهذا لا تُجهز إلا القيم الصحيحة التي يمكن أن تظهر فعلًا في كل حالة.
أثناء التعداد تُفحص المجاميع الجزئية لـ \(a^2+b^2+c^2+d^2\) مبكرًا، ولذلك تُقطع الفروع التي تجاوزت الحد مباشرة. ثم تُزال الرباعيات غير الأولية بواسطة اختبارات القاسم المشترك الأكبر، وتُستخدم قيود معيارية على الإشارات لتفادي التكرار البديهي بين \((a,b,c,d)\) و\(-(a,b,c,d)\).
وبالنسبة لكل إطار ينجو من التصفية، يحسب التنفيذ متجهات الحواف الثلاثة، ويستبعد الأطر التي تتجاوز امتداداتها المحورية أبعاد الصندوق، ثم يستخرج \(L_0\) و\(\Delta_0\) و\(s_{\max}\)، وبعدها يمر على جميع قيم \(s\) الصالحة. في كل قياس يضرب عدد الإزاحات في صيغة عدد نقاط المكعب، ويضيف الناتج إلى مجموع محلي خاص بالخيط، ثم بعد انتهاء جميع الخيوط يقسم المجموع الخام على \(24\).
نسختا Python وJava ليستا سوى أغلفة خفيفة فوق الحل المترجم نفسه، ولذلك فالثلاثة تشترك في الرياضيات ذاتها وفي نقاط التحقق نفسها وفي القيمة النهائية نفسها.
الكلفة المسيطرة تأتي من بحث الكواتيرنيونات. ويمكن الحصول على حد أعلى هندسي تقريبي بالنظر إلى الكرة رباعية الأبعاد \(a^2+b^2+c^2+d^2\le 4n\)، والتي تحتوي \(O(n^2)\) رباعيات صحيحة. لذا فمن منظور محافظ يمكن اعتبار فضاء البحث تقريبًا تربيعيًا في \(n\).
لكن الزمن الفعلي يكون أصغر بكثير من هذا الحد الخام. فاختبارات المعيار الجزئي تقطع عددًا كبيرًا من الفروع مبكرًا، ومرشحات القاسم المشترك الأكبر تزيل الرباعيات غير الأولية، واختبارات الامتداد تستبعد الأطر التي لا يمكن أن تدخل الصندوق أصلًا، كما أن كل إطار باقٍ لا يكرر إلا حتى \(s_{\max}\). لذلك فهذه الطريقة معتمدة على التقليم القوي، لا على تعداد ساذج لجميع المكعبات.
أما الذاكرة فتبقى محدودة. فالجداول المسبقة للقيم الزوجية والفردية لا تمتد إلا حتى رتبة \(O(\sqrt{n})\)، والنسخة المتوازية تخزن في الأساس بيانات المهام ومجمّعًا واحدًا لكل خيط عامل. ولهذا يظل استهلاك الذاكرة صغيرًا مقارنة بالعمل الحسابي.