For integers \(R \gt r \gt 0\), the hypocycloid is defined by
$$x(t)=(R-r)\cos t+r\cos\!\left(\frac{R-r}{r}t\right), \qquad y(t)=(R-r)\sin t-r\sin\!\left(\frac{R-r}{r}t\right).$$
We only keep parameters \(t\) for which \(\sin t\) and \(\cos t\) are rational, and among those points we keep only the ones with integer coordinates. For a fixed pair \((R,r)\), let \(S(R,r)\) be the sum of \(|x|+|y|\) over the distinct lattice points obtained this way. The global target is
$$T(N)=\sum_{3\le R\le N}\sum_{1\le r\lt R/2} S(R,r).$$
The implementations verify the checkpoints \(S(3,1)=10\), \(T(10)=524\), \(T(100)=580442\), and \(T(1000)=583108600\).
Write
$$R=gp,\qquad r=gq,\qquad \gcd(p,q)=1,\qquad e=p-q.$$
The condition \(r\lt R/2\) becomes \(q \lt p/2\), so \(e \gt q\). Now set \(t=q\theta\) and define \(z=e^{i\theta}\). Then \(z^q=e^{it}\), and the hypocycloid coordinates combine into one complex expression:
$$x+iy=g\left(ez^q+q\overline{z}^{\,e}\right).$$
This identity is the key simplification. For each reduced pair \((p,q)\), the geometry is fixed, and changing \(g\) only scales the final lattice point linearly.
Every rational point on the unit circle can be written as
$$z=\frac{A+iB}{D},\qquad A^2+B^2=D^2,\qquad \gcd(A,B)=1.$$
The trivial case \(D=1\) gives the four roots of unity \((\pm1,0)\) and \((0,\pm1)\). Every nontrivial rational point comes from a primitive Pythagorean triple, together with sign changes and swapping the two legs.
Substituting this parametrization into the complex form gives
$$x+iy=\frac{g}{D^e}\left(e(A+iB)^qD^{e-q}+q(A-iB)^e\right).$$
If we write
$$U_x+iU_y=e(A+iB)^qD^{e-q}+q(A-iB)^e,$$
then the entire arithmetic question becomes: after cancelling the common factors of \(D^e\), when does the remaining denominator divide \(g\)?
Let
$$g_0=\gcd(D^e,U_x,U_y),\qquad \delta=\frac{D^e}{g_0}.$$
After reducing the fraction, the point becomes
$$\left(x,y\right)=\frac{g}{\delta}\left(\frac{U_x}{g_0},\frac{U_y}{g_0}\right).$$
Therefore the coordinates are integers if and only if
$$\delta \mid g.$$
If \(g=m\delta\), then this one geometric shape contributes
$$m\left(\left|\frac{U_x}{g_0}\right|+\left|\frac{U_y}{g_0}\right|\right).$$
Summing over all admissible multiples of \(\delta\) produces a triangular number:
$$\sum_{m=1}^{m_{\max}} m=\frac{m_{\max}(m_{\max}+1)}{2},\qquad m_{\max}=\left\lfloor\frac{N/p}{\delta}\right\rfloor.$$
This is why the global computation can be organized by reduced shapes instead of looping over all \((R,r)\) one by one.
When \(D=1\), the only rational unit-circle points are the four roots of unity. For a fixed reduced pair \((p,q)\) and a scale \(g\), their total Manhattan contribution is
$$C_0(p,q,g)=g\begin{cases} 4p-2q,& p \text{ odd},\\ 4p,& 4\mid p,\\ 4p-4q,& p\equiv 2 \pmod 4. \end{cases}$$
The parity of \(p\) determines whether one axis pair has norm \(p\) or \(p-2q\). A small example is \((R,r)=(3,1)\), where \((g,p,q,e)=(1,3,1,2)\). The four points are
$$ (3,0),\quad (-1,0),\quad (-1,2),\quad (-1,-2), $$
so
$$S(3,1)=3+1+3+3=10.$$
For a fixed \(p\), define
$$A(p)=\#\{1\le q \lt p/2:\gcd(p,q)=1\}=\frac{\varphi(p)}2,$$
and
$$Q(p)=\sum_{\substack{1\le q \lt p/2\\ \gcd(p,q)=1}} q.$$
If \(M_p=\lfloor N/p\rfloor\), then the full \(D=1\) contribution is
$$T_{\mathrm{main}}(N)=\sum_{p=3}^{N} C_{\mathrm{main}}(p)\frac{M_p(M_p+1)}{2},$$
with
$$C_{\mathrm{main}}(p)=\begin{cases} 4pA(p)-2Q(p),& p \text{ odd},\\ 4pA(p),& 4\mid p,\\ 4pA(p)-4Q(p),& p\equiv 2 \pmod 4. \end{cases}$$
The nontrivial part is the coprime sum \(Q(p)\). Möbius inversion removes the gcd condition:
$$Q(p)=\sum_{d\mid p}\mu(d)\,d\,\frac{m_d(m_d+1)}{2},\qquad m_d=\left\lfloor\frac{p-1}{2d}\right\rfloor.$$
This is exactly the divisor sum accumulated by the implementation before the final outer summation over \(p\).
All remaining lattice points come from nontrivial rational points on the unit circle, so \(D \gt 1\). For each reduced pair \((p,q)\), each primitive Pythagorean triple, and each sign or coordinate-swap symmetry, Step 3 gives a required denominator \(\delta\). Whenever \(\delta\le M_p\), that one family contributes
$$\left(\left|\frac{U_x}{g_0}\right|+\left|\frac{U_y}{g_0}\right|\right)\frac{m_{\max}(m_{\max}+1)}{2},\qquad m_{\max}=\left\lfloor\frac{M_p}{\delta}\right\rfloor.$$
The direct single-pair routine removes coincident points explicitly before summing \(|x|+|y|\). For the bulk computation of \(T(N)\), the implementations combine the closed form for \(D=1\) with this explicit \(D \gt 1\) correction. At \(N=10^6\), the correction is validated to vanish for \(p \gt 12\), so the expensive branch can be truncated safely.
The C++, Python, and Java implementations all follow the same plan. They first build tables for the Möbius function and Euler's totient function up to \(N\). Those tables are then used to accumulate the divisor sums needed for \(Q(p)\), which yields the dominant arithmetic contribution from the \(D=1\) case.
After that, the implementation enumerates primitive Pythagorean triples only for the nontrivial \(D \gt 1\) branch. For each admissible reduced pair and each rational unit-circle point generated from the triple symmetries, it computes the reduced denominator \(\delta\) and adds the corresponding triangular-scale contribution. A separate direct evaluator for one pair \((R,r)\) is used for checkpoint verification and deduplicates geometric coincidences explicitly.
The sieve for \(\mu\) and \(\varphi\) is \(O(N)\) time and \(O(N)\) memory. The divisor accumulation that evaluates the Möbius formula for every \(p\le N\) is near-linear, with the usual \(O(N\log\log N)\) sieve-like behavior. The \(D \gt 1\) correction is much smaller in practice because only small reduced \(p\) values survive the validation cutoff and most primitive triples are eliminated immediately by the denominator test. Overall the method is dominated by arithmetic precomputation with \(O(N)\) storage.
Für ganze Zahlen \(R \gt r \gt 0\) ist die Hypozykloide durch
$$x(t)=(R-r)\cos t+r\cos\!\left(\frac{R-r}{r}t\right), \qquad y(t)=(R-r)\sin t-r\sin\!\left(\frac{R-r}{r}t\right)$$
gegeben. Es werden nur solche Parameter \(t\) betrachtet, für die \(\sin t\) und \(\cos t\) rational sind, und darunter nur die Punkte mit ganzzahligen Koordinaten. Für ein festes Paar \((R,r)\) sei \(S(R,r)\) die Summe von \(|x|+|y|\) über alle verschiedenen Gitterpunkte. Gesucht ist
$$T(N)=\sum_{3\le R\le N}\sum_{1\le r\lt R/2} S(R,r).$$
Die Implementierungen bestätigen unter anderem \(S(3,1)=10\), \(T(10)=524\), \(T(100)=580442\) und \(T(1000)=583108600\).
Schreibe
$$R=gp,\qquad r=gq,\qquad \gcd(p,q)=1,\qquad e=p-q.$$
Aus \(r\lt R/2\) folgt \(q \lt p/2\), also \(e \gt q\). Setze nun \(t=q\theta\) und \(z=e^{i\theta}\). Dann gilt \(z^q=e^{it}\), und die beiden Koordinaten lassen sich zu einem komplexen Ausdruck zusammenfassen:
$$x+iy=g\left(ez^q+q\overline{z}^{\,e}\right).$$
Damit ist die Geometrie für jedes reduzierte Paar \((p,q)\) fest vorgegeben; der Faktor \(g\) skaliert anschließend nur noch linear.
Jeder rationale Punkt auf dem Einheitskreis besitzt die Form
$$z=\frac{A+iB}{D},\qquad A^2+B^2=D^2,\qquad \gcd(A,B)=1.$$
Der triviale Fall \(D=1\) liefert die vier vierten Einheitswurzeln. Jeder nichttriviale rationale Punkt entsteht aus einem primitiven pythagoreischen Tripel, ergänzt um Vorzeichenwechsel und Vertauschung der Katheten.
Einsetzen in die komplexe Form ergibt
$$x+iy=\frac{g}{D^e}\left(e(A+iB)^qD^{e-q}+q(A-iB)^e\right).$$
Mit
$$U_x+iU_y=e(A+iB)^qD^{e-q}+q(A-iB)^e$$
reduziert sich die Frage auf das verbleibende Nennerproblem nach dem Kürzen.
Definiere
$$g_0=\gcd(D^e,U_x,U_y),\qquad \delta=\frac{D^e}{g_0}.$$
Nach dem Kürzen lautet der Punkt
$$\left(x,y\right)=\frac{g}{\delta}\left(\frac{U_x}{g_0},\frac{U_y}{g_0}\right).$$
Die Koordinaten sind also genau dann ganzzahlig, wenn
$$\delta \mid g$$
gilt. Für \(g=m\delta\) liefert dieser geometrische Punkt den Beitrag
$$m\left(\left|\frac{U_x}{g_0}\right|+\left|\frac{U_y}{g_0}\right|\right).$$
Über alle zulässigen Vielfachen summiert entsteht wieder eine Dreieckszahl:
$$\sum_{m=1}^{m_{\max}}m=\frac{m_{\max}(m_{\max}+1)}2,\qquad m_{\max}=\left\lfloor\frac{N/p}{\delta}\right\rfloor.$$
Darum kann die Gesamtsumme nach reduzierten Formen organisiert werden, statt jedes Paar \((R,r)\) einzeln zu behandeln.
Für \(D=1\) bleiben nur die vier Einheitswurzeln übrig. Für ein festes reduziertes Paar \((p,q)\) und einen Maßstab \(g\) ist die gesamte Manhattan-Summe
$$C_0(p,q,g)=g\begin{cases} 4p-2q,& p \text{ odd},\\ 4p,& 4\mid p,\\ 4p-4q,& p\equiv 2 \pmod 4. \end{cases}$$
Der Grund ist eine reine Paritätsbetrachtung: je nach Restklasse von \(p\) hat ein Achsenpaar Norm \(p\) oder \(p-2q\). Für \((R,r)=(3,1)\) mit \((g,p,q,e)=(1,3,1,2)\) erhält man die vier Punkte
$$ (3,0),\quad (-1,0),\quad (-1,2),\quad (-1,-2), $$
also
$$S(3,1)=3+1+3+3=10.$$
Für festes \(p\) setzen wir
$$A(p)=\#\{1\le q \lt p/2:\gcd(p,q)=1\}=\frac{\varphi(p)}2,$$
sowie
$$Q(p)=\sum_{\substack{1\le q \lt p/2\\ \gcd(p,q)=1}} q.$$
Mit \(M_p=\lfloor N/p\rfloor\) lautet der gesamte Beitrag aus dem Fall \(D=1\)
$$T_{\mathrm{main}}(N)=\sum_{p=3}^{N} C_{\mathrm{main}}(p)\frac{M_p(M_p+1)}{2},$$
wobei
$$C_{\mathrm{main}}(p)=\begin{cases} 4pA(p)-2Q(p),& p \text{ odd},\\ 4pA(p),& 4\mid p,\\ 4pA(p)-4Q(p),& p\equiv 2 \pmod 4. \end{cases}$$
Die schwierige Größe ist \(Q(p)\). Per Möbius-Inversion verschwindet die Teilerfremdheitsbedingung:
$$Q(p)=\sum_{d\mid p}\mu(d)\,d\,\frac{m_d(m_d+1)}{2},\qquad m_d=\left\lfloor\frac{p-1}{2d}\right\rfloor.$$
Genau diese Divisorensumme wird in der Implementierung effizient für alle \(p\le N\) aufgebaut.
Alle übrigen Gitterpunkte stammen von nichttrivialen rationalen Punkten des Einheitskreises, also von \(D \gt 1\). Für jedes reduzierte Paar \((p,q)\), jedes primitive Tripel und jede Vorzeichen- oder Vertauschungssymmetrie liefert Schritt 3 einen erforderlichen Nenner \(\delta\). Falls \(\delta\le M_p\), dann trägt diese Familie
$$\left(\left|\frac{U_x}{g_0}\right|+\left|\frac{U_y}{g_0}\right|\right)\frac{m_{\max}(m_{\max}+1)}{2},\qquad m_{\max}=\left\lfloor\frac{M_p}{\delta}\right\rfloor$$
zur Summe bei. Die direkte Auswertung eines einzelnen Paares entfernt geometrische Kollisionen ausdrücklich. Für die Gesamtfunktion \(T(N)\) kombinieren die Implementierungen den geschlossenen \(D=1\)-Term mit dieser expliziten \(D \gt 1\)-Korrektur. Bei \(N=10^6\) ist verifiziert, dass für \(p \gt 12\) kein solcher Zusatzterm mehr überlebt.
Die C++-, Python- und Java-Implementierungen verwenden denselben Aufbau. Zuerst werden Tabellen für die Möbius-Funktion und die Eulersche Phi-Funktion bis \(N\) erzeugt. Aus diesen Tabellen werden anschließend die Divisorensummen für \(Q(p)\) akkumuliert; damit erhält man den dominanten arithmetischen Anteil aus dem Fall \(D=1\).
Danach werden primitive pythagoreische Tripel nur noch für den nichttrivialen Zweig \(D \gt 1\) durchlaufen. Für jedes zulässige reduzierte Paar und jeden daraus erzeugten rationalen Punkt auf dem Einheitskreis berechnet die Implementierung den reduzierten Nenner \(\delta\) und addiert den passenden Dreieckszahl-Beitrag. Für Kontrollwerte einzelner Paare \((R,r)\) gibt es zusätzlich eine direkte Auswertung, die doppelte Punkte explizit entfernt.
Das Sieb für \(\mu\) und \(\varphi\) benötigt \(O(N)\) Zeit und \(O(N)\) Speicher. Die Divisorakkumulation für die Möbius-Formel über alle \(p\le N\) ist nahezu linear und verhält sich wie ein klassisches Sieb mit Aufwand \(O(N\log\log N)\). Die Korrektur für \(D \gt 1\) ist in der Praxis deutlich kleiner, weil nur kleine reduzierte \(p\)-Werte den geprüften Grenzfall überstehen und die meisten Tripel sofort am Nennerkriterium scheitern. Insgesamt dominiert also eine arithmetische Vorverarbeitung mit \(O(N)\) Speicherbedarf.
\(R \gt r \gt 0\) tam sayıları için hiposikloid
$$x(t)=(R-r)\cos t+r\cos\!\left(\frac{R-r}{r}t\right), \qquad y(t)=(R-r)\sin t-r\sin\!\left(\frac{R-r}{r}t\right)$$
ile tanımlanır. Yalnızca \(\sin t\) ve \(\cos t\) rasyonel olan parametreler alınır; bunlar arasından da sadece tam sayı koordinat veren noktalar tutulur. Sabit bir \((R,r)\) çifti için \(S(R,r)\), elde edilen farklı kafes noktalarının \(|x|+|y|\) toplamıdır. İstenen toplam
$$T(N)=\sum_{3\le R\le N}\sum_{1\le r\lt R/2} S(R,r)$$
şeklindedir. Uygulamalar \(S(3,1)=10\), \(T(10)=524\), \(T(100)=580442\) ve \(T(1000)=583108600\) denetimlerini doğrular.
Şöyle yazalım:
$$R=gp,\qquad r=gq,\qquad \gcd(p,q)=1,\qquad e=p-q.$$
\(r\lt R/2\) koşulu \(q \lt p/2\) demektir; dolayısıyla \(e \gt q\). Şimdi \(t=q\theta\) ve \(z=e^{i\theta}\) alalım. Böylece \(z^q=e^{it}\) olur ve iki koordinat tek bir kompleks ifadede birleşir:
$$x+iy=g\left(ez^q+q\overline{z}^{\,e}\right).$$
Bu dönüşüm temel adımdır. Çünkü indirgenmiş \((p,q)\) sabit olduğunda geometri sabitlenir; \(g\) yalnızca sonucu doğrusal olarak büyütür.
Birim çember üzerindeki her rasyonel nokta
$$z=\frac{A+iB}{D},\qquad A^2+B^2=D^2,\qquad \gcd(A,B)=1$$
şeklinde yazılabilir. \(D=1\) durumu dört kök birliği verir. \(D \gt 1\) olan tüm durumlar, ilkel Pisagor üçlülerinden ve bunların işaret ile bacak değiştirme simetrilerinden gelir.
Bu parametrizasyonu kompleks forma yerleştirince
$$x+iy=\frac{g}{D^e}\left(e(A+iB)^qD^{e-q}+q(A-iB)^e\right)$$
elde edilir. Eğer
$$U_x+iU_y=e(A+iB)^qD^{e-q}+q(A-iB)^e$$
yazarsak, bütün problem payda sadeleştikten sonra geriye hangi bölenin kaldığını anlamaya dönüşür.
Tanımlayalım:
$$g_0=\gcd(D^e,U_x,U_y),\qquad \delta=\frac{D^e}{g_0}.$$
Sadeleştirme sonrası nokta
$$\left(x,y\right)=\frac{g}{\delta}\left(\frac{U_x}{g_0},\frac{U_y}{g_0}\right)$$
olur. Demek ki koordinatların tam sayı olması için gerekli ve yeterli şart
$$\delta \mid g$$
eşitliğidir. Eğer \(g=m\delta\) ise, bu geometrik şeklin katkısı
$$m\left(\left|\frac{U_x}{g_0}\right|+\left|\frac{U_y}{g_0}\right|\right)$$
olur. Tüm uygun katlar toplanınca üçgensel toplam çıkar:
$$\sum_{m=1}^{m_{\max}}m=\frac{m_{\max}(m_{\max}+1)}2,\qquad m_{\max}=\left\lfloor\frac{N/p}{\delta}\right\rfloor.$$
Böylece toplamı tek tek tüm \((R,r)\) çiftleri üzerinden değil, indirgenmiş şekiller üzerinden toplamak mümkün olur.
\(D=1\) olduğunda birim çemberde sadece dört kök birliği kalır. Sabit bir indirgenmiş \((p,q)\) ve ölçek \(g\) için Manhattan toplamı
$$C_0(p,q,g)=g\begin{cases} 4p-2q,& p \text{ odd},\\ 4p,& 4\mid p,\\ 4p-4q,& p\equiv 2 \pmod 4 \end{cases}$$
olur. Bunun nedeni paritedir: bazı kalıntı sınıflarında bir eksen çifti \(p\) normuna, bazılarında ise \(p-2q\) normuna iner. Küçük örnek olarak \((R,r)=(3,1)\) için \((g,p,q,e)=(1,3,1,2)\) alınır ve dört nokta
$$ (3,0),\quad (-1,0),\quad (-1,2),\quad (-1,-2) $$
şeklinde gelir. Dolayısıyla
$$S(3,1)=3+1+3+3=10.$$
Sabit \(p\) için
$$A(p)=\#\{1\le q \lt p/2:\gcd(p,q)=1\}=\frac{\varphi(p)}2,$$
ve
$$Q(p)=\sum_{\substack{1\le q \lt p/2\\ \gcd(p,q)=1}} q$$
tanımlarını yapalım. \(M_p=\lfloor N/p\rfloor\) ile, \(D=1\) dalının tamamı
$$T_{\mathrm{main}}(N)=\sum_{p=3}^{N} C_{\mathrm{main}}(p)\frac{M_p(M_p+1)}{2}$$
olur; burada
$$C_{\mathrm{main}}(p)=\begin{cases} 4pA(p)-2Q(p),& p \text{ odd},\\ 4pA(p),& 4\mid p,\\ 4pA(p)-4Q(p),& p\equiv 2 \pmod 4. \end{cases}$$
Zor kısım \(Q(p)\)'dir. Möbius terslemesi ile aralarında asal olma koşulu kaldırılır:
$$Q(p)=\sum_{d\mid p}\mu(d)\,d\,\frac{m_d(m_d+1)}{2},\qquad m_d=\left\lfloor\frac{p-1}{2d}\right\rfloor.$$
Uygulama, dış toplama başlamadan önce tam olarak bu bölen toplamlarını tüm \(p\le N\) için biriktirir.
Kalan bütün kafes noktaları, birim çemberin \(D \gt 1\) olan rasyonel noktalarından gelir. Her indirgenmiş \((p,q)\), her ilkel Pisagor üçlüsü ve her işaret ya da koordinat değişim simetrisi için, Adım 3 gerekli paydayı \(\delta\) olarak verir. Eğer \(\delta\le M_p\) ise, o aile
$$\left(\left|\frac{U_x}{g_0}\right|+\left|\frac{U_y}{g_0}\right|\right)\frac{m_{\max}(m_{\max}+1)}{2},\qquad m_{\max}=\left\lfloor\frac{M_p}{\delta}\right\rfloor$$
katkısını yapar. Tek bir \((R,r)\) çifti doğrudan hesaplanırken çakışan noktalar açıkça elenir. \(T(N)\)'nin toplu hesabında ise uygulamalar, \(D=1\) kapalı formunu bu açık \(D \gt 1\) düzeltmesiyle birleştirir. \(N=10^6\) için doğrulanan sonuç, \(p \gt 12\) olduğunda artık hiçbir ek terimin kalmadığıdır.
C++, Python ve Java uygulamaları aynı iskeleti izler. Önce \(N\)'ye kadar Möbius fonksiyonu ile Euler phi fonksiyonunun tabloları kuruluyor. Ardından \(Q(p)\) için gereken bölen toplamları hesaplanıyor; böylece \(D=1\) durumundan gelen baskın aritmetik katkı elde ediliyor.
Daha sonra yalnızca \(D \gt 1\) dalı için ilkel Pisagor üçlüleri taranıyor. Her uygun indirgenmiş çift ve bu çiftin ürettiği rasyonel birim çember noktası için, uygulama indirgenmiş paydayı \(\delta\) hesaplayıp ilgili üçgensel ölçek katkısını ekliyor. Tek bir \((R,r)\) çifti için kullanılan ayrı doğrulama yordamı ise geometrideki çakışmaları kaldırarak \(|x|+|y|\) toplamını doğrudan alıyor.
\(\mu\) ve \(\varphi\) süzgeci \(O(N)\) zaman ve \(O(N)\) bellek ister. Möbius formülünün tüm \(p\le N\) için hesaplandığı bölen biriktirme adımı, klasik süzgeç mantığında yaklaşık \(O(N\log\log N)\) davranır. \(D \gt 1\) düzeltmesi ise pratikte çok daha küçüktür; çünkü doğrulanmış kesim nedeniyle yalnızca küçük indirgenmiş \(p\) değerleri kalır ve çoğu üçlü payda testiyle hemen elenir. Sonuç olarak yöntem esasen aritmetik ön-hesap tarafından belirlenir ve \(O(N)\) bellek kullanır.
Para enteros \(R \gt r \gt 0\), la hipocicloide viene dada por
$$x(t)=(R-r)\cos t+r\cos\!\left(\frac{R-r}{r}t\right), \qquad y(t)=(R-r)\sin t-r\sin\!\left(\frac{R-r}{r}t\right).$$
Solo se consideran los parámetros \(t\) para los que \(\sin t\) y \(\cos t\) son racionales, y dentro de ellos solo se conservan los puntos con coordenadas enteras. Para un par fijo \((R,r)\), definimos \(S(R,r)\) como la suma de \(|x|+|y|\) sobre los puntos de la red distintos que aparecen así. El objetivo global es
$$T(N)=\sum_{3\le R\le N}\sum_{1\le r\lt R/2} S(R,r).$$
Las implementaciones verifican los valores de control \(S(3,1)=10\), \(T(10)=524\), \(T(100)=580442\) y \(T(1000)=583108600\).
Escribimos
$$R=gp,\qquad r=gq,\qquad \gcd(p,q)=1,\qquad e=p-q.$$
La condición \(r\lt R/2\) equivale a \(q \lt p/2\), luego \(e \gt q\). Ahora ponemos \(t=q\theta\) y \(z=e^{i\theta}\). Entonces \(z^q=e^{it}\), y las dos coordenadas se combinan en una sola expresión compleja:
$$x+iy=g\left(ez^q+q\overline{z}^{\,e}\right).$$
Esta identidad es la simplificación central: para cada par reducido \((p,q)\) la forma geométrica queda fijada, y cambiar \(g\) solo multiplica linealmente el punto final.
Todo punto racional del círculo unitario puede escribirse como
$$z=\frac{A+iB}{D},\qquad A^2+B^2=D^2,\qquad \gcd(A,B)=1.$$
El caso trivial \(D=1\) produce las cuatro raíces de la unidad. Todo punto racional no trivial proviene de una terna pitagórica primitiva, junto con cambios de signo y el intercambio de los dos catetos.
Al sustituir esta parametrización en la fórmula compleja obtenemos
$$x+iy=\frac{g}{D^e}\left(e(A+iB)^qD^{e-q}+q(A-iB)^e\right).$$
Si escribimos
$$U_x+iU_y=e(A+iB)^qD^{e-q}+q(A-iB)^e,$$
la cuestión aritmética se reduce a determinar qué denominador queda después de cancelar factores comunes con \(D^e\).
Definimos
$$g_0=\gcd(D^e,U_x,U_y),\qquad \delta=\frac{D^e}{g_0}.$$
Después de simplificar, el punto queda
$$\left(x,y\right)=\frac{g}{\delta}\left(\frac{U_x}{g_0},\frac{U_y}{g_0}\right).$$
Por tanto, las coordenadas son enteras si y solo si
$$\delta \mid g.$$
Si \(g=m\delta\), la contribución de esa forma geométrica es
$$m\left(\left|\frac{U_x}{g_0}\right|+\left|\frac{U_y}{g_0}\right|\right).$$
Al sumar todos los múltiplos admisibles de \(\delta\) aparece un número triangular:
$$\sum_{m=1}^{m_{\max}}m=\frac{m_{\max}(m_{\max}+1)}2,\qquad m_{\max}=\left\lfloor\frac{N/p}{\delta}\right\rfloor.$$
Ese es el motivo por el que el cálculo global puede agruparse por formas reducidas y no por cada par \((R,r)\) por separado.
Cuando \(D=1\), los únicos puntos racionales del círculo unitario son las cuatro raíces de la unidad. Para un par reducido \((p,q)\) y una escala \(g\), la contribución total en norma Manhattan es
$$C_0(p,q,g)=g\begin{cases} 4p-2q,& p \text{ odd},\\ 4p,& 4\mid p,\\ 4p-4q,& p\equiv 2 \pmod 4. \end{cases}$$
La razón es paritaria: según la clase de \(p\), uno de los pares sobre los ejes tiene norma \(p\) o \(p-2q\). Un ejemplo pequeño es \((R,r)=(3,1)\), donde \((g,p,q,e)=(1,3,1,2)\). Los cuatro puntos son
$$ (3,0),\quad (-1,0),\quad (-1,2),\quad (-1,-2), $$
y por eso
$$S(3,1)=3+1+3+3=10.$$
Para un \(p\) fijo definimos
$$A(p)=\#\{1\le q \lt p/2:\gcd(p,q)=1\}=\frac{\varphi(p)}2,$$
y también
$$Q(p)=\sum_{\substack{1\le q \lt p/2\\ \gcd(p,q)=1}} q.$$
Si \(M_p=\lfloor N/p\rfloor\), el aporte completo del caso \(D=1\) es
$$T_{\mathrm{main}}(N)=\sum_{p=3}^{N} C_{\mathrm{main}}(p)\frac{M_p(M_p+1)}{2},$$
donde
$$C_{\mathrm{main}}(p)=\begin{cases} 4pA(p)-2Q(p),& p \text{ odd},\\ 4pA(p),& 4\mid p,\\ 4pA(p)-4Q(p),& p\equiv 2 \pmod 4. \end{cases}$$
La cantidad delicada es \(Q(p)\). La inversión de Möbius elimina la condición de coprimalidad:
$$Q(p)=\sum_{d\mid p}\mu(d)\,d\,\frac{m_d(m_d+1)}{2},\qquad m_d=\left\lfloor\frac{p-1}{2d}\right\rfloor.$$
Esa es exactamente la suma sobre divisores que las implementaciones acumulan antes de la suma exterior sobre \(p\).
Todos los demás puntos de la red proceden de puntos racionales no triviales del círculo unitario, es decir, de \(D \gt 1\). Para cada par reducido \((p,q)\), cada terna pitagórica primitiva y cada simetría de signos o intercambio de coordenadas, el Paso 3 produce un denominador requerido \(\delta\). Siempre que \(\delta\le M_p\), esa familia aporta
$$\left(\left|\frac{U_x}{g_0}\right|+\left|\frac{U_y}{g_0}\right|\right)\frac{m_{\max}(m_{\max}+1)}{2},\qquad m_{\max}=\left\lfloor\frac{M_p}{\delta}\right\rfloor.$$
La rutina directa para un solo par elimina coincidencias geométricas de manera explícita. Para el cálculo masivo de \(T(N)\), las implementaciones combinan la forma cerrada del caso \(D=1\) con esta corrección explícita para \(D \gt 1\). En \(N=10^6\) se valida además que no sobrevive ningún término no trivial cuando \(p \gt 12\).
Las implementaciones en C++, Python y Java siguen la misma estructura. Primero construyen tablas para la función de Möbius y la función phi de Euler hasta \(N\). Con esas tablas acumulan las sumas sobre divisores que permiten obtener \(Q(p)\), y con ello evalúan el término aritmético dominante procedente del caso \(D=1\).
Después enumeran ternas pitagóricas primitivas solo para la rama no trivial \(D \gt 1\). Para cada par reducido admisible y para cada punto racional del círculo unitario generado por las simetrías de la terna, la implementación calcula el denominador reducido \(\delta\) y suma la contribución triangular correspondiente. Un evaluador separado para un par \((R,r)\) sirve como verificación de puntos de control y elimina duplicados geométricos antes de sumar \(|x|+|y|\).
El cribado de \(\mu\) y \(\varphi\) cuesta \(O(N)\) tiempo y \(O(N)\) memoria. La acumulación de divisores que evalúa la fórmula de Möbius para todos los \(p\le N\) es casi lineal, con el comportamiento típico \(O(N\log\log N)\) de los cribados armónicos. La corrección para \(D \gt 1\) es mucho menor en la práctica porque solo sobreviven valores reducidos pequeños de \(p\) y la prueba del denominador elimina la mayoría de las ternas de inmediato. En conjunto, el coste está dominado por la precomputación aritmética y el uso de memoria es \(O(N)\).
对整数 \(R \gt r \gt 0\)，题目的 hypocycloid 参数方程为
$$x(t)=(R-r)\cos t+r\cos\!\left(\frac{R-r}{r}t\right), \qquad y(t)=(R-r)\sin t-r\sin\!\left(\frac{R-r}{r}t\right).$$
我们只保留那些满足 \(\sin t,\cos t\in\mathbb{Q}\) 的参数值，并且只统计其中坐标同时为整数的点。对固定的 \((R,r)\)，令 \(S(R,r)\) 表示所有不同整点的 \(|x|+|y|\) 之和。全局要求的是
$$T(N)=\sum_{3\le R\le N}\sum_{1\le r\lt R/2} S(R,r).$$
程序中的小规模校验值包括 \(S(3,1)=10\)、\(T(10)=524\)、\(T(100)=580442\) 和 \(T(1000)=583108600\)。这些值说明推导出来的计数方式与实现完全一致。
写成
$$R=gp,\qquad r=gq,\qquad \gcd(p,q)=1,\qquad e=p-q.$$
条件 \(r\lt R/2\) 等价于 \(q \lt p/2\)，因此 \(e \gt q\)。再令 \(t=q\theta\)，并记 \(z=e^{i\theta}\)。于是 \(z^q=e^{it}\)，原来的两个实坐标可以合并成一个复数表达式：
$$x+iy=g\left(ez^q+q\overline{z}^{\,e}\right).$$
这一步非常关键。因为对固定的既约参数对 \((p,q)\) 来说，曲线的“形状”已经固定，而 \(g\) 只负责把最终点按整数比例放大。也就是说，先分析既约形状，再对所有允许的缩放统一求和，会比直接枚举所有 \((R,r)\) 高效得多。
单位圆上的任意有理点都可写成
$$z=\frac{A+iB}{D},\qquad A^2+B^2=D^2,\qquad \gcd(A,B)=1.$$
\(D=1\) 时就是四个平凡点 \((\pm1,0)\)、\((0,\pm1)\)。所有 \(D \gt 1\) 的非平凡有理点，都来自原始勾股三元组，再配合符号翻转与两条直角边互换的对称性。
把这个表示代入上面的复数公式，可以得到
$$x+iy=\frac{g}{D^e}\left(e(A+iB)^qD^{e-q}+q(A-iB)^e\right).$$
若记
$$U_x+iU_y=e(A+iB)^qD^{e-q}+q(A-iB)^e,$$
那么问题就被压缩成一个纯算术判定：在 \(D^e\) 与分子中的公共因子约掉之后，还剩下什么分母，它是否整除 \(g\)。
定义
$$g_0=\gcd(D^e,U_x,U_y),\qquad \delta=\frac{D^e}{g_0}.$$
约分之后的点坐标就是
$$\left(x,y\right)=\frac{g}{\delta}\left(\frac{U_x}{g_0},\frac{U_y}{g_0}\right).$$
因此坐标同时为整数，当且仅当
$$\delta \mid g$$
成立。如果 \(g=m\delta\)，那么这个几何形状的贡献就是
$$m\left(\left|\frac{U_x}{g_0}\right|+\left|\frac{U_y}{g_0}\right|\right).$$
对所有允许的 \(m\) 求和，就会出现三角数
$$\sum_{m=1}^{m_{\max}}m=\frac{m_{\max}(m_{\max}+1)}2,\qquad m_{\max}=\left\lfloor\frac{N/p}{\delta}\right\rfloor.$$
这正是整体算法能把“缩放因子求和”提前压成封闭公式的原因。
当 \(D=1\) 时，单位圆上只有四个四次单位根，因此一个既约参数对 \((p,q)\) 在尺度 \(g\) 下的总曼哈顿贡献可以直接写成
$$C_0(p,q,g)=g\begin{cases} 4p-2q,& p \text{ odd},\\ 4p,& 4\mid p,\\ 4p-4q,& p\equiv 2 \pmod 4. \end{cases}$$
这里的分段来自纯粹的奇偶分析：某些余数类下，轴上的一个点对长度保持为 \(p\)；另一些余数类下，它会缩成 \(p-2q\)。最简单的例子是 \((R,r)=(3,1)\)。此时 \((g,p,q,e)=(1,3,1,2)\)，四个点恰好是
$$ (3,0),\quad (-1,0),\quad (-1,2),\quad (-1,-2), $$
所以
$$S(3,1)=3+1+3+3=10.$$
固定 \(p\) 后，记
$$A(p)=\#\{1\le q \lt p/2:\gcd(p,q)=1\}=\frac{\varphi(p)}2,$$
再记
$$Q(p)=\sum_{\substack{1\le q \lt p/2\\ \gcd(p,q)=1}} q.$$
设 \(M_p=\lfloor N/p\rfloor\)，则全部 \(D=1\) 贡献为
$$T_{\mathrm{main}}(N)=\sum_{p=3}^{N} C_{\mathrm{main}}(p)\frac{M_p(M_p+1)}{2},$$
其中
$$C_{\mathrm{main}}(p)=\begin{cases} 4pA(p)-2Q(p),& p \text{ odd},\\ 4pA(p),& 4\mid p,\\ 4pA(p)-4Q(p),& p\equiv 2 \pmod 4. \end{cases}$$
真正需要技巧的是 \(Q(p)\)。用 Möbius 反演去掉互素条件，可得
$$Q(p)=\sum_{d\mid p}\mu(d)\,d\,\frac{m_d(m_d+1)}{2},\qquad m_d=\left\lfloor\frac{p-1}{2d}\right\rfloor.$$
实现正是先把这个除数和对所有 \(p\le N\) 累积出来，再做最外层的求和。因此主项完全由数论预处理支撑，而不需要对每个 \(q\) 单独扫描。
剩余的整点全部来自非平凡有理点，也就是 \(D \gt 1\) 的情况。对每个既约对 \((p,q)\)、每个原始勾股三元组以及每个符号或坐标交换对称，步骤 3 都会给出一个所需分母 \(\delta\)。只要 \(\delta\le M_p\)，这一族点就贡献
$$\left(\left|\frac{U_x}{g_0}\right|+\left|\frac{U_y}{g_0}\right|\right)\frac{m_{\max}(m_{\max}+1)}{2},\qquad m_{\max}=\left\lfloor\frac{M_p}{\delta}\right\rfloor.$$
单独计算某个 \((R,r)\) 时，程序会显式去重，以防不同参数化得到同一个几何点。计算总和 \(T(N)\) 时，则把 \(D=1\) 的封闭主项与这部分 \(D \gt 1\) 修正拼接起来。对于 \(N=10^6\)，实现还验证了当 \(p \gt 12\) 时不会再出现任何 \(D \gt 1\) 的有效项，因此昂贵分支可以安全截断。
C++、Python 和 Java 实现遵循同一条主线。首先建立到 \(N\) 为止的 Möbius 函数与 Euler phi 函数表。然后利用这些表批量累积 \(Q(p)\) 所需的除数和，从而一次性得到主导性的 \(D=1\) 算术部分。
接下来，程序只在 \(D \gt 1\) 的分支上枚举原始勾股三元组。对每个允许的既约参数对，以及由三元组对称性生成的每个单位圆有理点，程序计算约分后的分母 \(\delta\)，再把对应的三角数缩放贡献加入总和。另有一个针对单个 \((R,r)\) 的直接求值器，用于校验小样例，并在求 \(|x|+|y|\) 之前显式去除重复点。
\(\mu\) 与 \(\varphi\) 的筛法是 \(O(N)\) 时间、\(O(N)\) 空间。对所有 \(p\le N\) 累积 Möbius 反演除数和的过程是近线性的，具有典型的 \(O(N\log\log N)\) 筛法行为。\(D \gt 1\) 修正项在实际运行中小得多，因为经过验证后只剩下很小的既约 \(p\) 范围，而且大多数原始三元组会被分母判定立即剪枝。整体来看，主要成本来自数论预处理，空间复杂度保持在 \(O(N)\)。
Для целых \(R \gt r \gt 0\) гипоциклоида задается формулами
$$x(t)=(R-r)\cos t+r\cos\!\left(\frac{R-r}{r}t\right), \qquad y(t)=(R-r)\sin t-r\sin\!\left(\frac{R-r}{r}t\right).$$
Рассматриваются только такие параметры \(t\), для которых \(\sin t\) и \(\cos t\) рациональны, и среди них выбираются лишь точки с целыми координатами. Для фиксированной пары \((R,r)\) обозначим через \(S(R,r)\) сумму \(|x|+|y|\) по всем различным решетчатым точкам. Требуется вычислить
$$T(N)=\sum_{3\le R\le N}\sum_{1\le r\lt R/2} S(R,r).$$
Реализации подтверждают контрольные значения \(S(3,1)=10\), \(T(10)=524\), \(T(100)=580442\) и \(T(1000)=583108600\).
Пишем
$$R=gp,\qquad r=gq,\qquad \gcd(p,q)=1,\qquad e=p-q.$$
Условие \(r\lt R/2\) превращается в \(q \lt p/2\), поэтому \(e \gt q\). Теперь положим \(t=q\theta\) и \(z=e^{i\theta}\). Тогда \(z^q=e^{it}\), и обе координаты удобно собрать в одну комплексную формулу:
$$x+iy=g\left(ez^q+q\overline{z}^{\,e}\right).$$
В этом и состоит главная редукция. Для фиксированной приведенной пары \((p,q)\) геометрическая форма уже задана, а параметр \(g\) только линейно растягивает конечную точку.
Любая рациональная точка единичной окружности имеет вид
$$z=\frac{A+iB}{D},\qquad A^2+B^2=D^2,\qquad \gcd(A,B)=1.$$
Случай \(D=1\) дает четыре тривиальные корня единицы. Все нетривиальные рациональные точки с \(D \gt 1\) получаются из примитивных пифагоровых троек вместе со сменой знаков и перестановкой катетов.
После подстановки в комплексную форму получаем
$$x+iy=\frac{g}{D^e}\left(e(A+iB)^qD^{e-q}+q(A-iB)^e\right).$$
Если обозначить
$$U_x+iU_y=e(A+iB)^qD^{e-q}+q(A-iB)^e,$$
то вся арифметическая часть задачи сводится к вопросу: какой знаменатель остается после сокращения общих множителей с \(D^e\).
Введем
$$g_0=\gcd(D^e,U_x,U_y),\qquad \delta=\frac{D^e}{g_0}.$$
После сокращения координаты равны
$$\left(x,y\right)=\frac{g}{\delta}\left(\frac{U_x}{g_0},\frac{U_y}{g_0}\right).$$
Следовательно, точка лежит в решетке тогда и только тогда, когда
$$\delta \mid g.$$
Если \(g=m\delta\), то вклад этой геометрической конфигурации равен
$$m\left(\left|\frac{U_x}{g_0}\right|+\left|\frac{U_y}{g_0}\right|\right).$$
Суммирование по всем допустимым масштабам снова дает треугольное число:
$$\sum_{m=1}^{m_{\max}}m=\frac{m_{\max}(m_{\max}+1)}2,\qquad m_{\max}=\left\lfloor\frac{N/p}{\delta}\right\rfloor.$$
Именно поэтому глобальную сумму удобно перестраивать по приведенным формам, а не по всем парам \((R,r)\) отдельно.
При \(D=1\) остаются только четыре корня единицы. Для фиксированной приведенной пары \((p,q)\) и масштаба \(g\) суммарный вклад в манхэттенской норме равен
$$C_0(p,q,g)=g\begin{cases} 4p-2q,& p \text{ odd},\\ 4p,& 4\mid p,\\ 4p-4q,& p\equiv 2 \pmod 4. \end{cases}$$
Это простой эффект четности: в зависимости от класса \(p\) одна осевая пара имеет норму либо \(p\), либо \(p-2q\). Для примера \((R,r)=(3,1)\) имеем \((g,p,q,e)=(1,3,1,2)\), и четыре точки равны
$$ (3,0),\quad (-1,0),\quad (-1,2),\quad (-1,-2), $$
поэтому
$$S(3,1)=3+1+3+3=10.$$
Для фиксированного \(p\) положим
$$A(p)=\#\{1\le q \lt p/2:\gcd(p,q)=1\}=\frac{\varphi(p)}2,$$
а также
$$Q(p)=\sum_{\substack{1\le q \lt p/2\\ \gcd(p,q)=1}} q.$$
Если \(M_p=\lfloor N/p\rfloor\), то весь вклад случая \(D=1\) записывается как
$$T_{\mathrm{main}}(N)=\sum_{p=3}^{N} C_{\mathrm{main}}(p)\frac{M_p(M_p+1)}{2},$$
где
$$C_{\mathrm{main}}(p)=\begin{cases} 4pA(p)-2Q(p),& p \text{ odd},\\ 4pA(p),& 4\mid p,\\ 4pA(p)-4Q(p),& p\equiv 2 \pmod 4. \end{cases}$$
Тонкое место здесь именно \(Q(p)\). Условие взаимной простоты снимается обращением Мёбиуса:
$$Q(p)=\sum_{d\mid p}\mu(d)\,d\,\frac{m_d(m_d+1)}{2},\qquad m_d=\left\lfloor\frac{p-1}{2d}\right\rfloor.$$
Именно эту сумму по делителям реализации аккумулируют заранее для всех \(p\le N\), а затем выполняют внешнее суммирование.
Все остальные решетчатые точки происходят из нетривиальных рациональных точек единичной окружности, то есть из случаев \(D \gt 1\). Для каждой приведенной пары \((p,q)\), каждой примитивной пифагоровой тройки и каждой симметрии по знакам или перестановке координат шаг 3 дает требуемый знаменатель \(\delta\). Если \(\delta\le M_p\), то эта семья дает вклад
$$\left(\left|\frac{U_x}{g_0}\right|+\left|\frac{U_y}{g_0}\right|\right)\frac{m_{\max}(m_{\max}+1)}{2},\qquad m_{\max}=\left\lfloor\frac{M_p}{\delta}\right\rfloor.$$
При прямом вычислении для одной пары \((R,r)\) совпадающие геометрические точки удаляются явно. Для массового вычисления \(T(N)\) реализации соединяют замкнутую формулу для \(D=1\) с этой явной поправкой для \(D \gt 1\). Для \(N=10^6\) дополнительно проверено, что при \(p \gt 12\) нет ни одного выжившего нетривиального вклада.
Реализации на C++, Python и Java построены одинаково. Сначала они строят таблицы функции Мёбиуса и функции Эйлера \(\varphi\) до \(N\). Затем эти таблицы используются для накопления сумм по делителям, которые дают \(Q(p)\), то есть главный арифметический вклад из случая \(D=1\).
После этого перечисляются примитивные пифагоровы тройки только для нетривиальной ветви \(D \gt 1\). Для каждой допустимой приведенной пары и для каждой рациональной точки единичной окружности, возникающей из симметрий тройки, вычисляется сокращенный знаменатель \(\delta\), после чего добавляется соответствующий вклад с треугольным масштабированием. Отдельный прямой вычислитель для одной пары \((R,r)\) используется для контрольных проверок и явно удаляет дубликаты перед суммированием \(|x|+|y|\).
Решето для \(\mu\) и \(\varphi\) требует \(O(N)\) времени и \(O(N)\) памяти. Накопление сумм по делителям для формулы Мёбиуса по всем \(p\le N\) близко к линейному и ведет себя как классическое решето с трудоемкостью \(O(N\log\log N)\). Поправка для \(D \gt 1\) на практике значительно меньше, потому что после проверенного ограничения остаются только малые приведенные значения \(p\), а большинство троек сразу отбрасывается проверкой знаменателя. В итоге основную стоимость задает арифметическая предобработка при памяти \(O(N)\).
لأعداد صحيحة \(R \gt r \gt 0\)، يُعطى منحنى hypocycloid بالمعادلتين
$$x(t)=(R-r)\cos t+r\cos\!\left(\frac{R-r}{r}t\right), \qquad y(t)=(R-r)\sin t-r\sin\!\left(\frac{R-r}{r}t\right).$$
نحتفظ فقط بقيم \(t\) التي تجعل \(\sin t\) و\(\cos t\) عددين نسبيين، ثم نختار من بينها النقاط التي تكون إحداثياتها صحيحة. لزوج ثابت \((R,r)\)، نرمز بـ \(S(R,r)\) إلى مجموع \(|x|+|y|\) على جميع نقاط الشبكة المختلفة الناتجة بهذه الطريقة. والمطلوب هو
$$T(N)=\sum_{3\le R\le N}\sum_{1\le r\lt R/2} S(R,r).$$
تؤكد التطبيقات قيم الفحص \(S(3,1)=10\) و\(T(10)=524\) و\(T(100)=580442\) و\(T(1000)=583108600\).
نكتب
$$R=gp,\qquad r=gq,\qquad \gcd(p,q)=1,\qquad e=p-q.$$
الشرط \(r\lt R/2\) يكافئ \(q \lt p/2\)، ومن ثم \(e \gt q\). بعد ذلك نضع \(t=q\theta\) ونعرف \(z=e^{i\theta}\). عندئذٍ يكون \(z^q=e^{it}\)، ويمكن جمع الإحداثيين في صيغة مركبة واحدة:
$$x+iy=g\left(ez^q+q\overline{z}^{\,e}\right).$$
هذه هي الفكرة الأساسية. فلكل زوج مختزل \((p,q)\) تصبح البنية الهندسية ثابتة، بينما يقوم \(g\) فقط بتكبير النقطة النهائية تكبيرًا خطيًا.
كل نقطة نسبية على دائرة الوحدة يمكن كتابتها بالصورة
$$z=\frac{A+iB}{D},\qquad A^2+B^2=D^2,\qquad \gcd(A,B)=1.$$
الحالة البسيطة \(D=1\) تعطي الجذور الرابعة للوحدة. أما كل نقطة نسبية غير تافهة بحيث \(D \gt 1\)، فتأتي من ثلاثية فيثاغورية أولية مع تبديل الإشارات وتبادل الضلعين القائمين.
بالتعويض في الصيغة المركبة نحصل على
$$x+iy=\frac{g}{D^e}\left(e(A+iB)^qD^{e-q}+q(A-iB)^e\right).$$
وإذا كتبنا
$$U_x+iU_y=e(A+iB)^qD^{e-q}+q(A-iB)^e,$$
فإن السؤال الحسابي كله يتحول إلى معرفة المقام المتبقي بعد اختزال العوامل المشتركة مع \(D^e\).
نعرف
$$g_0=\gcd(D^e,U_x,U_y),\qquad \delta=\frac{D^e}{g_0}.$$
بعد الاختزال تصبح النقطة
$$\left(x,y\right)=\frac{g}{\delta}\left(\frac{U_x}{g_0},\frac{U_y}{g_0}\right).$$
إذن تكون الإحداثيات صحيحة إذا وفقط إذا تحقق
$$\delta \mid g.$$
وإذا كان \(g=m\delta\)، فإن مساهمة هذا الشكل الهندسي هي
$$m\left(\left|\frac{U_x}{g_0}\right|+\left|\frac{U_y}{g_0}\right|\right).$$
وعند الجمع على جميع المضاعفات المسموح بها لـ \(\delta\) يظهر عدد مثلثي:
$$\sum_{m=1}^{m_{\max}}m=\frac{m_{\max}(m_{\max}+1)}2,\qquad m_{\max}=\left\lfloor\frac{N/p}{\delta}\right\rfloor.$$
ولهذا يمكن إعادة تنظيم المجموع الكلي بحسب الأشكال المختزلة بدل المرور على جميع الأزواج \((R,r)\) مباشرة.
عندما يكون \(D=1\)، لا تبقى إلا الجذور الأربع للوحدة. ولزوج مختزل ثابت \((p,q)\) ومعامل قياس \(g\)، يكون مجموع مسافات Manhattan هو
$$C_0(p,q,g)=g\begin{cases} 4p-2q,& p \text{ odd},\\ 4p,& 4\mid p,\\ 4p-4q,& p\equiv 2 \pmod 4. \end{cases}$$
السبب هنا مرتبط بتكافؤ الأعداد: فحسب باقي \(p\) قد يكون أحد الزوجين الواقعين على المحورين بطول \(p\) أو بطول \(p-2q\). المثال الأبسط هو \((R,r)=(3,1)\)، حيث \((g,p,q,e)=(1,3,1,2)\)، فنحصل على النقاط الأربع
$$ (3,0),\quad (-1,0),\quad (-1,2),\quad (-1,-2), $$
ومن ثم
$$S(3,1)=3+1+3+3=10.$$
لـ \(p\) ثابت نعرّف
$$A(p)=\#\{1\le q \lt p/2:\gcd(p,q)=1\}=\frac{\varphi(p)}2,$$
وكذلك
$$Q(p)=\sum_{\substack{1\le q \lt p/2\\ \gcd(p,q)=1}} q.$$
إذا كان \(M_p=\lfloor N/p\rfloor\)، فإن كل مساهمة الحالة \(D=1\) تساوي
$$T_{\mathrm{main}}(N)=\sum_{p=3}^{N} C_{\mathrm{main}}(p)\frac{M_p(M_p+1)}{2},$$
حيث
$$C_{\mathrm{main}}(p)=\begin{cases} 4pA(p)-2Q(p),& p \text{ odd},\\ 4pA(p),& 4\mid p,\\ 4pA(p)-4Q(p),& p\equiv 2 \pmod 4. \end{cases}$$
الكمية غير المباشرة هنا هي \(Q(p)\). ويزيل عكس Möbius شرط التباين المشترك:
$$Q(p)=\sum_{d\mid p}\mu(d)\,d\,\frac{m_d(m_d+1)}{2},\qquad m_d=\left\lfloor\frac{p-1}{2d}\right\rfloor.$$
وهذا بالضبط هو مجموع القواسم الذي تبنيه التطبيقات مسبقًا قبل تنفيذ الجمع الخارجي على جميع قيم \(p\).
جميع نقاط الشبكة المتبقية تأتي من النقاط النسبية غير التافهة على دائرة الوحدة، أي من الحالات التي فيها \(D \gt 1\). لكل زوج مختزل \((p,q)\)، ولكل ثلاثية فيثاغورية أولية، ولكل تناظر ناتج عن تغيير الإشارة أو تبديل الإحداثيين، تعطي الخطوة 3 مقامًا مطلوبًا \(\delta\). وإذا تحقق \(\delta\le M_p\)، فإن تلك العائلة تضيف
$$\left(\left|\frac{U_x}{g_0}\right|+\left|\frac{U_y}{g_0}\right|\right)\frac{m_{\max}(m_{\max}+1)}{2},\qquad m_{\max}=\left\lfloor\frac{M_p}{\delta}\right\rfloor.$$
عند الحساب المباشر لزوج واحد \((R,r)\)، يتم حذف النقاط المتطابقة صراحةً قبل جمع \(|x|+|y|\). أما في الحساب الكلي لـ \(T(N)\)، فتجمع التطبيقات بين الصيغة المغلقة للحالة \(D=1\) وبين هذا التصحيح الصريح للحالة \(D \gt 1\). ولدى \(N=10^6\) تم التحقق من أن أي مساهمة غير تافهة تختفي عندما \(p \gt 12\)، ولذلك يمكن قطع هذا الفرع بأمان.
تتبع تطبيقات C++ وPython وJava البنية نفسها. أولًا تُبنى جداول دالتي Möbius وEuler phi حتى \(N\). ثم تُستخدم هذه الجداول لتجميع مجاميع القواسم التي تعطي \(Q(p)\)، ومن ثم يُحسب الحد الحسابي الغالب القادم من الحالة \(D=1\).
بعد ذلك تُعدَّد الثلاثيات الفيثاغورية الأولية فقط في الفرع غير التافه \(D \gt 1\). ولكل زوج مختزل مسموح به، ولكل نقطة نسبية على دائرة الوحدة يولدها تناظر الثلاثية، تحسب الشيفرة المقام المختزل \(\delta\) وتضيف مساهمة القياس المثلثية الموافقة. كما يوجد مُقيِّم مباشر لزوج واحد \((R,r)\) للتحقق من قيم الاختبار، ويزيل التكرارات الهندسية قبل جمع \(|x|+|y|\).
الغربال الخاص بـ \(\mu\) و\(\varphi\) يعمل في \(O(N)\) زمنًا و\(O(N)\) ذاكرةً. أما خطوة تجميع القواسم اللازمة لتطبيق صيغة Möbius على جميع \(p\le N\)، فهي شبه خطية وتتصرف كسلوك غرابيل الأعداد المعتاد من رتبة \(O(N\log\log N)\). وتصحيح \(D \gt 1\) أصغر بكثير عمليًا، لأن حد التحقق يبقي فقط قيم \(p\) المختزلة الصغيرة، كما أن اختبار المقام يستبعد معظم الثلاثيات مباشرة. لذلك تهيمن مرحلة المعالجة الحسابية المسبقة، بينما تبقى الذاكرة من رتبة \(O(N)\).