Let \(F(n)\) be the number of nonzero squarefree Gaussian integers \(z=a+bi\) with norm
$$N(z)=a^2+b^2\le n,$$
where two values that differ only by multiplication by a unit in \(\{\pm 1,\pm i\}\) are regarded as the same object. The goal is to evaluate \(F(10^{14})\).
A direct scan over all lattice points in the disk \(a^2+b^2\le n\) is far too slow. The solution instead turns squarefreeness into a Möbius-style inclusion-exclusion problem in \(\mathbb{Z}[i]\), collapses that Gaussian sum to an ordinary integer sum indexed by norms, and then accelerates the remaining geometric counts by grouping equal quotients.
Define the lattice-point counting function
$$Q(x)=\#\{(a,b)\in\mathbb{Z}^2:a^2+b^2\le x\}.$$
This counts all Gaussian integers of norm at most \(x\), including \(0\). The final formula will use \(Q(x)-1\) so that the origin does not contribute.
In the unique factorization domain \(\mathbb{Z}[i]\), the squarefree indicator has the same shape as over the ordinary integers:
$$\mathbf{1}_{\mathrm{sqfree}}(z)=\sum_{d^2\mid z}\mu_{\mathbb{Z}[i]}(d),$$
where \(\mu_{\mathbb{Z}[i]}(d)\) is the Gaussian Möbius function. It is \(0\) if \(d\) is divisible by the square of a Gaussian prime, and otherwise \((-1)^k\) when \(d\) is a product of \(k\) distinct Gaussian primes up to units.
Summing over all nonzero Gaussian integers with norm at most \(n\) gives
$$4F(n)=\sum_{\substack{z\in\mathbb{Z}[i]\\0<N(z)\le n}}\mathbf{1}_{\mathrm{sqfree}}(z)=\sum_{d\ne 0}\mu_{\mathbb{Z}[i]}(d)\left(Q\!\left(\left\lfloor\frac{n}{N(d)^2}\right\rfloor\right)-1\right).$$
The factor \(4\) appears because every nonzero Gaussian integer has exactly four associates obtained by multiplying by the units.
Only the norm \(N(d)\) matters in the geometric term, so we group together all Gaussian divisors with the same norm. Define
$$W(m)=\sum_{N(d)=m}\mu_{\mathbb{Z}[i]}(d).$$
Then the previous identity becomes
$$4F(n)=\sum_{m\le \sqrt{n}} W(m)\left(Q\!\left(\left\lfloor\frac{n}{m^2}\right\rfloor\right)-1\right).$$
This already explains the overall structure of the programs: build the arithmetic weight \(W(m)\), evaluate the circle count \(Q(\cdot)\), and combine them.
The function \(W\) is multiplicative, so it is enough to know its values on prime powers. These come directly from the three possible behaviors of a rational prime in \(\mathbb{Z}[i]\).
For \(p=2\), the prime is ramified:
$$2=-i(1+i)^2.$$
There is only one Gaussian prime above \(2\), with norm \(2\). Therefore the local contribution is
$$1-t.$$
For \(p\equiv 1\pmod{4}\), the prime splits into two nonassociate conjugate Gaussian primes, both of norm \(p\). The local squarefree possibilities are: choose neither, choose one of the two, or choose both. Hence the local polynomial is
$$\left(1-t\right)^2=1-2t+t^2.$$
For \(p\equiv 3\pmod{4}\), the prime stays Gaussian-prime, but its norm is \(p^2\). Thus the local polynomial is
$$1-t^2.$$
Reading off the coefficients gives
$$W(p^e)= \begin{cases} -1, & p=2,\ e=1,\\ -2, & p\equiv 1\pmod{4},\ e=1,\\ 1, & p\equiv 1\pmod{4},\ e=2,\\ -1, & p\equiv 3\pmod{4},\ e=2,\\ 0, & \text{otherwise}. \end{cases}$$
Multiplicativity then determines \(W(m)\) for every \(m\).
The geometric part is the circle-counting function \(Q(x)\). The implementations use the exact identity
$$Q(x)=1+4\sum_{a=0}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\sqrt{x-a^2}\right\rfloor.$$
For each nonnegative \(a\), the inner floor gives the largest admissible \(b>0\). As \(a\) increases, that maximal \(b\) never increases, so a monotone two-pointer scan evaluates the whole sum in \(O(\sqrt{x})\) time rather than recomputing a square root search from scratch for every \(a\).
The quantity
$$k=\left\lfloor\frac{n}{m^2}\right\rfloor$$
is constant on an interval of consecutive \(m\)-values. If the current block starts at \(\ell\), its right endpoint is
$$r=\left\lfloor\sqrt{\frac{n}{k}}\right\rfloor.$$
If we also precompute the prefix sums
$$P(t)=\sum_{m\le t}W(m),$$
then the whole block contributes
$$\left(P(r)-P(\ell-1)\right)\left(Q(k)-1\right).$$
This is the key speedup: the expensive circle count is evaluated once per distinct quotient instead of once per \(m\).
Here \(\lfloor\sqrt{10}\rfloor=3\), so only \(m=1,2,3\) can appear. From the prime-power rule we get
$$W(1)=1,\qquad W(2)=-1,\qquad W(3)=0.$$
Now compute the two needed circle counts:
$$Q(10)=37,\qquad Q(2)=9.$$
Therefore
$$4F(10)=1\cdot(37-1)+(-1)\cdot(9-1)+0\cdot(Q(1)-1)=36-8=28,$$
so
$$F(10)=7.$$
This matches the small checkpoint used by the implementations.
The C++, Python, and Java implementations follow the same plan. First they set \(M=\lfloor\sqrt{n}\rfloor\) and build a smallest-prime-factor table up to \(M\). This makes it easy to factor every integer \(m\le M\) and to evaluate the multiplicative weight \(W(m)\) from the prime-power rules above.
Next they build prefix sums of \(W\), so any block sum \(\sum_{\ell\le m\le r}W(m)\) is available in constant time. They then partition the interval \(1\le m\le M\) into maximal blocks on which \(\lfloor n/m^2\rfloor\) is constant. For each block they evaluate the lattice count \(Q(k)\), multiply by the corresponding weight sum, add everything together, and divide the final total by \(4\).
The C++ version computes the independent circle counts for different quotient blocks in parallel, while the Python and Java versions perform the same arithmetic sequentially. The mathematical formula is identical in all three languages.
Let \(M=\lfloor\sqrt{n}\rfloor\). Building the smallest-prime-factor data and the multiplicative weights needs \(O(M)\) memory and near-linear time in \(M\). The prefix sums are also \(O(M)\).
The number of distinct values of \(\lfloor n/m^2\rfloor\) is much smaller than \(M\); it is \(O(n^{1/3})\). Each such value requires one circle-count evaluation \(Q(k)\), and that evaluation costs \(O(\sqrt{k})\) time with constant extra memory. In practice the geometric counts dominate the running time, which is why quotient grouping and, in C++, parallel evaluation are the important optimizations.
Sei \(F(n)\) die Anzahl der von \(0\) verschiedenen quadratfreien gaußschen ganzen Zahlen \(z=a+bi\) mit Norm
$$N(z)=a^2+b^2\le n,$$
wobei zwei Werte, die sich nur um eine Einheit aus \(\{\pm 1,\pm i\}\) unterscheiden, als dasselbe Objekt gezählt werden. Gesucht ist also \(F(10^{14})\).
Ein direktes Durchmustern aller Gitterpunkte in der Kreisscheibe \(a^2+b^2\le n\) ist viel zu teuer. Die Lösung verwandelt die Quadratfreiheit daher zuerst in eine Möbius-artige Inklusions-Exklusions-Formel in \(\mathbb{Z}[i]\), fasst dann alle gaußschen Teiler mit gleicher Norm zusammen und beschleunigt die verbleibenden Kreiszählungen durch Intervallgruppierung.
Definiere die Gitterpunktfunktion
$$Q(x)=\#\{(a,b)\in\mathbb{Z}^2:a^2+b^2\le x\}.$$
Sie zählt alle gaußschen ganzen Zahlen mit Norm höchstens \(x\), einschließlich \(0\). In der Endformel verwendet man deshalb \(Q(x)-1\), um den Ursprung auszuschließen.
Im faktoriellen Ring \(\mathbb{Z}[i]\) hat die Indikatorfunktion für Quadratfreiheit dieselbe Form wie bei den gewöhnlichen ganzen Zahlen:
$$\mathbf{1}_{\mathrm{sqfree}}(z)=\sum_{d^2\mid z}\mu_{\mathbb{Z}[i]}(d).$$
Dabei ist \(\mu_{\mathbb{Z}[i]}(d)\) die gaußsche Möbius-Funktion. Sie ist \(0\), falls \(d\) durch das Quadrat einer gaußschen Primzahl teilbar ist, und sonst \((-1)^k\), wenn \(d\) das Produkt von \(k\) verschiedenen gaußschen Primzahlen bis auf Einheiten ist.
Summiert man über alle von \(0\) verschiedenen gaußschen ganzen Zahlen mit Norm höchstens \(n\), so erhält man
$$4F(n)=\sum_{\substack{z\in\mathbb{Z}[i]\\0<N(z)\le n}}\mathbf{1}_{\mathrm{sqfree}}(z)=\sum_{d\ne 0}\mu_{\mathbb{Z}[i]}(d)\left(Q\!\left(\left\lfloor\frac{n}{N(d)^2}\right\rfloor\right)-1\right).$$
Der Faktor \(4\) kommt daher, dass jede von \(0\) verschiedene gaußsche ganze Zahl genau vier Assoziierte besitzt.
Im geometrischen Teil hängt alles nur von der Norm \(N(d)\) ab. Deshalb bündeln wir alle gaußschen Teiler mit derselben Norm und definieren
$$W(m)=\sum_{N(d)=m}\mu_{\mathbb{Z}[i]}(d).$$
Dann wird die vorige Identität zu
$$4F(n)=\sum_{m\le \sqrt{n}} W(m)\left(Q\!\left(\left\lfloor\frac{n}{m^2}\right\rfloor\right)-1\right).$$
Damit ist die Struktur der Programme bereits klar: erst die arithmetischen Gewichte \(W(m)\) berechnen, dann die Kreiszählungen \(Q(\cdot)\) auswerten und beides verknüpfen.
Die Funktion \(W\) ist multiplikativ. Es genügt daher, ihre Werte für Primzahlpotenzen zu kennen. Diese folgen direkt aus den drei möglichen Verhalten rationaler Primzahlen in \(\mathbb{Z}[i]\).
Für \(p=2\) liegt Verzweigung vor:
$$2=-i(1+i)^2.$$
Über \(2\) gibt es also genau eine gaußsche Primzahl mit Norm \(2\). Der lokale Faktor ist daher
$$1-t.$$
Für \(p\equiv 1\pmod{4}\) zerfällt \(p\) in zwei nichtassoziierte konjugierte gaußsche Primzahlen mit Norm \(p\). Die quadratfreien lokalen Möglichkeiten sind: keine wählen, genau eine wählen oder beide wählen. Deshalb ist das lokale Polynom
$$\left(1-t\right)^2=1-2t+t^2.$$
Für \(p\equiv 3\pmod{4}\) bleibt \(p\) gaußsch prim, besitzt aber Norm \(p^2\). Somit ist der lokale Faktor
$$1-t^2.$$
Aus den Koeffizienten liest man ab
$$W(p^e)= \begin{cases} -1, & p=2,\ e=1,\\ -2, & p\equiv 1\pmod{4},\ e=1,\\ 1, & p\equiv 1\pmod{4},\ e=2,\\ -1, & p\equiv 3\pmod{4},\ e=2,\\ 0, & \text{otherwise}. \end{cases}$$
Durch Multiplikativität ist damit \(W(m)\) für jedes \(m\) bestimmt.
Der geometrische Teil ist die Kreisfunktion \(Q(x)\). Die Implementierungen benutzen die exakte Identität
$$Q(x)=1+4\sum_{a=0}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\sqrt{x-a^2}\right\rfloor.$$
Für jedes nichtnegative \(a\) liefert der innere Abrundungsausdruck den größten zulässigen Wert von \(b>0\). Wenn \(a\) wächst, kann dieser Maximalwert von \(b\) nie ansteigen. Deshalb genügt ein monotoner Zwei-Zeiger-Durchlauf, um die gesamte Summe in \(O(\sqrt{x})\) Zeit auszuwerten.
Die Größe
$$k=\left\lfloor\frac{n}{m^2}\right\rfloor$$
bleibt auf ganzen Intervallen aufeinanderfolgender \(m\)-Werte konstant. Beginnt ein Block bei \(\ell\), dann ist sein rechter Rand
$$r=\left\lfloor\sqrt{\frac{n}{k}}\right\rfloor.$$
Berechnet man zusätzlich die Präfixsummen
$$P(t)=\sum_{m\le t}W(m),$$
dann ist der Beitrag des gesamten Blocks gleich
$$\left(P(r)-P(\ell-1)\right)\left(Q(k)-1\right).$$
Genau hier entsteht der Geschwindigkeitsgewinn: Die teure Kreiszählung wird nur einmal pro verschiedenem Quotienten ausgewertet und nicht einmal pro \(m\).
Hier ist \(\lfloor\sqrt{10}\rfloor=3\), also kommen nur \(m=1,2,3\) infrage. Aus der Primzahlpotenz-Regel folgt
$$W(1)=1,\qquad W(2)=-1,\qquad W(3)=0.$$
Nun benötigt man nur zwei Kreiszählungen:
$$Q(10)=37,\qquad Q(2)=9.$$
Damit ergibt sich
$$4F(10)=1\cdot(37-1)+(-1)\cdot(9-1)+0\cdot(Q(1)-1)=36-8=28,$$
also
$$F(10)=7.$$
Dieser Wert stimmt mit dem kleinen Kontrollpunkt der Implementierungen überein.
Die C++-, Python- und Java-Implementierungen folgen demselben Plan. Zuerst setzen sie \(M=\lfloor\sqrt{n}\rfloor\) und bauen eine Tabelle der kleinsten Primfaktoren bis \(M\) auf. Damit lässt sich jedes \(m\le M\) schnell faktorisieren, und aus den lokalen Primzahlpotenz-Regeln erhält man das multiplikative Gewicht \(W(m)\).
Anschließend werden Präfixsummen von \(W\) aufgebaut, sodass jede Blocksumme \(\sum_{\ell\le m\le r}W(m)\) in konstanter Zeit verfügbar ist. Danach zerlegen die Programme das Intervall \(1\le m\le M\) in maximale Blöcke, auf denen \(\lfloor n/m^2\rfloor\) konstant ist. Für jeden Block wird \(Q(k)\) einmal berechnet, mit der entsprechenden Gewichtssumme multipliziert, alles addiert und ganz am Ende durch \(4\) geteilt.
Die C++-Version wertet die voneinander unabhängigen Kreiszählungen blockweise parallel aus. Die Python- und Java-Versionen führen dieselbe Mathematik sequentiell aus. Die zugrunde liegende Formel ist jedoch in allen drei Sprachen identisch.
Mit \(M=\lfloor\sqrt{n}\rfloor\) benötigen die kleinsten-Primfaktor-Daten und die Gewichte \(W(m)\) \(O(M)\) Speicher und annähernd lineare Zeit in \(M\). Auch die Präfixsummen kosten \(O(M)\).
Die Anzahl verschiedener Werte von \(\lfloor n/m^2\rfloor\) ist deutlich kleiner als \(M\); sie liegt bei \(O(n^{1/3})\). Für jeden solchen Wert wird genau eine Kreiszählung \(Q(k)\) benötigt, und diese kostet \(O(\sqrt{k})\) Zeit bei konstantem Zusatzspeicher. In der Praxis dominieren deshalb die geometrischen Auswertungen, sodass Quotientengruppierung und in C++ zusätzlich Parallelisierung die wesentlichen Optimierungen sind.
\(F(n)\), normu
$$N(z)=a^2+b^2\le n$$
olan, sıfırdan farklı kareözgür Gaussian tamsayıların \(z=a+bi\) sayısını göstersin. Birbirinden yalnızca \(\{\pm 1,\pm i\}\) birimleriyle çarpılma farkı olan değerler aynı nesne kabul edilir. Amaç \(F(10^{14})\) değerini hesaplamaktır.
\(a^2+b^2\le n\) diskindeki tüm kafes noktalarını doğrudan taramak çok pahalıdır. Çözüm bunun yerine kareözgürlük koşulunu \(\mathbb{Z}[i]\) içinde Möbius benzeri bir dahil et-çıkar formülüne dönüştürür, sonra Gaussian bölenleri normlarına göre bir araya toplar ve geriye kalan daire içi nokta sayımlarını eşit bölüm gruplarıyla hızlandırır.
Önce şu kafes nokta sayım fonksiyonunu tanımlayalım:
$$Q(x)=\#\{(a,b)\in\mathbb{Z}^2:a^2+b^2\le x\}.$$
Bu fonksiyon normu en çok \(x\) olan bütün Gaussian tamsayıları, yani \(0\)'ı da dahil ederek sayar. Son formülde orijini çıkarmak için \(Q(x)-1\) kullanılacaktır.
Tek çarpanlara ayırma bölgesi olan \(\mathbb{Z}[i]\) içinde kareözgürlük göstergesi, sıradan tamsayılardakiyle aynı biçimdedir:
$$\mathbf{1}_{\mathrm{sqfree}}(z)=\sum_{d^2\mid z}\mu_{\mathbb{Z}[i]}(d).$$
Burada \(\mu_{\mathbb{Z}[i]}(d)\) Gaussian Möbius fonksiyonudur. Eğer \(d\), bir Gaussian asalın karesiyle bölünüyorsa değer \(0\) olur; aksi halde \(d\), birimlere göre \(k\) farklı Gaussian asalın çarpımıysa değer \((-1)^k\) olur.
Bu ifadeyi normu en çok \(n\) olan sıfırdan farklı bütün Gaussian tamsayılar üzerinde toplarsak
$$4F(n)=\sum_{\substack{z\in\mathbb{Z}[i]\\0<N(z)\le n}}\mathbf{1}_{\mathrm{sqfree}}(z)=\sum_{d\ne 0}\mu_{\mathbb{Z}[i]}(d)\left(Q\!\left(\left\lfloor\frac{n}{N(d)^2}\right\rfloor\right)-1\right)$$
elde edilir. Baştaki \(4\) katsayısı, sıfırdan farklı her Gaussian tamsayının dört eşdeğer birim çarpanına sahip olmasından gelir.
Geometrik kısımda yalnızca \(N(d)\) normu önemlidir. Bu yüzden aynı norma sahip bütün Gaussian bölenleri tek bir ağırlık altında toplarız:
$$W(m)=\sum_{N(d)=m}\mu_{\mathbb{Z}[i]}(d).$$
Böylece ana özdeşlik
$$4F(n)=\sum_{m\le \sqrt{n}} W(m)\left(Q\!\left(\left\lfloor\frac{n}{m^2}\right\rfloor\right)-1\right)$$
biçimine iner. Programların iskeleti tam olarak budur: önce aritmetik ağırlık \(W(m)\), sonra daire içi nokta sayısı \(Q(\cdot)\), sonunda ikisinin birleşimi.
\(W\) çarpımsal bir fonksiyondur; dolayısıyla \(W(p^e)\) değerlerini bilmek yeterlidir. Bunlar, rasyonel asalların \(\mathbb{Z}[i]\) içindeki üç farklı davranışından doğrudan çıkar.
\(p=2\) için asal ramified olur:
$$2=-i(1+i)^2.$$
Yani \(2\)'nin üzerinde normu \(2\) olan yalnızca bir Gaussian asal vardır. Bu durumda yerel çarpan
$$1-t$$
olur.
\(p\equiv 1\pmod{4}\) için asal iki eşlenik ve birbirine associate olmayan Gaussian asala ayrılır; ikisinin de normu \(p\)'dir. Yerel kareözgür seçenekler: hiçbirini seçmemek, ikisinden birini seçmek veya ikisini birden seçmektir. Bu yüzden yerel polinom
$$\left(1-t\right)^2=1-2t+t^2$$
şeklindedir.
\(p\equiv 3\pmod{4}\) için asal Gaussian asal olarak kalır, fakat normu \(p^2\) olur. Bu durumda yerel çarpan
$$1-t^2$$
elde edilir. Katsayıları okuyunca
$$W(p^e)= \begin{cases} -1, & p=2,\ e=1,\\ -2, & p\equiv 1\pmod{4},\ e=1,\\ 1, & p\equiv 1\pmod{4},\ e=2,\\ -1, & p\equiv 3\pmod{4},\ e=2,\\ 0, & \text{otherwise}. \end{cases}$$
Çarpımsallık sayesinde artık her \(m\) için \(W(m)\) belirlenmiş olur.
Geometrik kısım \(Q(x)\) fonksiyonudur. Uygulamalar şu tam özdeşliği kullanır:
$$Q(x)=1+4\sum_{a=0}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\sqrt{x-a^2}\right\rfloor.$$
Her sabit \(a\ge 0\) için içteki taban fonksiyonu, \(b>0\) için izin verilen en büyük değeri verir. \(a\) arttıkça bu en büyük \(b\) değeri asla artmaz. Bu nedenle monoton iki işaretçili tarama toplamı \(O(\sqrt{x})\) sürede hesaplar.
Ana toplamda görülen
$$k=\left\lfloor\frac{n}{m^2}\right\rfloor$$
ifadesi, art arda gelen birçok \(m\) için sabit kalır. Bir blok \(\ell\) ile başlıyorsa sağ uç
$$r=\left\lfloor\sqrt{\frac{n}{k}}\right\rfloor$$
olur. Ayrıca ön toplamları
$$P(t)=\sum_{m\le t}W(m)$$
tanımlarsak, tüm blok katkısı
$$\left(P(r)-P(\ell-1)\right)\left(Q(k)-1\right)$$
şeklinde tek adımda alınır. En pahalı daire sayımı böylece her \(m\) için değil, yalnızca her farklı bölüm değeri için bir kez yapılır.
Burada \(\lfloor\sqrt{10}\rfloor=3\) olduğundan yalnızca \(m=1,2,3\) katkı verebilir. Asal kuvvet kuralından
$$W(1)=1,\qquad W(2)=-1,\qquad W(3)=0$$
elde edilir. Gerekli iki daire sayımı ise
$$Q(10)=37,\qquad Q(2)=9$$
olur. Dolayısıyla
$$4F(10)=1\cdot(37-1)+(-1)\cdot(9-1)+0\cdot(Q(1)-1)=36-8=28,$$
yani
$$F(10)=7.$$
Bu küçük değer, C++, Python ve Java uygulamalarının kontrol noktasıyla aynıdır.
C++, Python ve Java uygulamaları aynı planı izler. Önce \(M=\lfloor\sqrt{n}\rfloor\) alınır ve \(M\)'e kadar en küçük asal çarpan tablosu kurulur. Böylece her \(m\le M\) hızlıca çarpanlara ayrılır ve yukarıdaki asal kuvvet kurallarından çarpımsal ağırlık \(W(m)\) hesaplanır.
Ardından \(W\) için ön toplamlar hazırlanır; böylece \(\sum_{\ell\le m\le r}W(m)\) her blok için sabit zamanda bulunur. Sonra \(1\le m\le M\) aralığı, \(\lfloor n/m^2\rfloor\) sabit kalacak şekilde en büyük bloklara ayrılır. Her blok için \(Q(k)\) bir kez hesaplanır, uygun ağırlık toplamıyla çarpılır, bütün katkılar toplanır ve en sonda sonuç \(4\)'e bölünür.
C++ sürümü, birbirinden bağımsız bloklar için gereken daire sayımlarını paralel yürütür. Python ve Java sürümleri aynı matematiği sıralı olarak uygular. Üç dilde de hesaplanan formül aynıdır.
\(M=\lfloor\sqrt{n}\rfloor\) olmak üzere, en küçük asal çarpan bilgisi ve \(W(m)\) ağırlıkları \(O(M)\) bellek ve \(M\) açısından yaklaşık doğrusal zaman ister. Ön toplamlar da \(O(M)\) maliyetlidir.
\(\lfloor n/m^2\rfloor\) ifadesinin farklı değer sayısı \(M\)'den çok daha küçüktür; büyüklük mertebesi \(O(n^{1/3})\) olur. Her farklı değer için bir daire sayımı \(Q(k)\) gerekir ve bu işlem \(O(\sqrt{k})\) zamanda, sabit ek bellekle yapılır. Pratikte asıl maliyet geometrik sayımlardır; bu yüzden bölüm gruplama ve C++ tarafında paralelleştirme temel hızlandırmalardır.
Sea \(F(n)\) el número de enteros gaussianos no nulos y libres de cuadrados \(z=a+bi\) cuya norma satisface
$$N(z)=a^2+b^2\le n,$$
contados módulo las unidades \(\{\pm 1,\pm i\}\). Es decir, dos valores que solo difieren por multiplicación por una unidad representan la misma clase. La meta es calcular \(F(10^{14})\).
Enumerar directamente todos los puntos de rejilla del disco \(a^2+b^2\le n\) es inviable. La estrategia correcta es traducir la condición de ser libre de cuadrados a una suma de inclusión-exclusión con la función de Möbius gaussiana, reagrupar los divisores según su norma racional y acelerar las cuentas geométricas agrupando los valores repetidos de \(\lfloor n/m^2\rfloor\).
Definimos la función de conteo de puntos de rejilla
$$Q(x)=\#\{(a,b)\in\mathbb{Z}^2:a^2+b^2\le x\}.$$
Esta cantidad cuenta todos los enteros gaussianos de norma a lo sumo \(x\), incluyendo \(0\). Por eso la fórmula final usa \(Q(x)-1\) para eliminar el origen.
En el dominio de factorización única \(\mathbb{Z}[i]\), el indicador de ser libre de cuadrados tiene la forma
$$\mathbf{1}_{\mathrm{sqfree}}(z)=\sum_{d^2\mid z}\mu_{\mathbb{Z}[i]}(d).$$
Aquí \(\mu_{\mathbb{Z}[i]}(d)\) es la función de Möbius en los enteros gaussianos: vale \(0\) si \(d\) es divisible por el cuadrado de algún primo gaussiano, y si no vale \((-1)^k\) cuando \(d\) es producto de \(k\) primos gaussianos distintos, salvo unidades.
Al sumar sobre todos los enteros gaussianos no nulos con norma a lo sumo \(n\), obtenemos
$$4F(n)=\sum_{\substack{z\in\mathbb{Z}[i]\\0<N(z)\le n}}\mathbf{1}_{\mathrm{sqfree}}(z)=\sum_{d\ne 0}\mu_{\mathbb{Z}[i]}(d)\left(Q\!\left(\left\lfloor\frac{n}{N(d)^2}\right\rfloor\right)-1\right).$$
El factor \(4\) aparece porque cada entero gaussiano no nulo tiene exactamente cuatro asociados.
En la parte geométrica solo importa la norma \(N(d)\). Por tanto, reunimos todos los divisores gaussianos con la misma norma y definimos
$$W(m)=\sum_{N(d)=m}\mu_{\mathbb{Z}[i]}(d).$$
La identidad anterior se transforma entonces en
$$4F(n)=\sum_{m\le \sqrt{n}} W(m)\left(Q\!\left(\left\lfloor\frac{n}{m^2}\right\rfloor\right)-1\right).$$
Esta igualdad explica toda la arquitectura del algoritmo: construir el peso aritmético \(W(m)\), calcular \(Q(\cdot)\) y combinar ambos términos.
La función \(W\) es multiplicativa, así que basta determinar \(W(p^e)\). Eso depende del comportamiento de cada primo racional dentro de \(\mathbb{Z}[i]\).
Para \(p=2\), el primo está ramificado:
$$2=-i(1+i)^2.$$
Solo hay un primo gaussiano por encima de \(2\), con norma \(2\), así que el factor local es
$$1-t.$$
Para \(p\equiv 1\pmod{4}\), el primo se descompone como producto de dos primos gaussianos conjugados y no asociados, ambos de norma \(p\). Las posibilidades libres de cuadrados son: no elegir ninguno, elegir uno de los dos o elegir ambos. Por eso el polinomio local es
$$\left(1-t\right)^2=1-2t+t^2.$$
Para \(p\equiv 3\pmod{4}\), el primo sigue siendo gaussiano-primo, pero su norma es \(p^2\). Entonces el factor local es
$$1-t^2.$$
Al leer los coeficientes se obtiene
$$W(p^e)= \begin{cases} -1, & p=2,\ e=1,\\ -2, & p\equiv 1\pmod{4},\ e=1,\\ 1, & p\equiv 1\pmod{4},\ e=2,\\ -1, & p\equiv 3\pmod{4},\ e=2,\\ 0, & \text{otherwise}. \end{cases}$$
La multiplicatividad determina después \(W(m)\) para cualquier \(m\).
La parte geométrica es la función \(Q(x)\). Las implementaciones usan la identidad exacta
$$Q(x)=1+4\sum_{a=0}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\sqrt{x-a^2}\right\rfloor.$$
Para cada \(a\ge 0\), el término interior da el mayor \(b>0\) permitido. Cuando \(a\) aumenta, ese máximo de \(b\) nunca aumenta. Así, un recorrido monótono con dos punteros evalúa toda la suma en \(O(\sqrt{x})\).
La cantidad
$$k=\left\lfloor\frac{n}{m^2}\right\rfloor$$
permanece constante en bloques de valores consecutivos de \(m\). Si un bloque comienza en \(\ell\), su extremo derecho es
$$r=\left\lfloor\sqrt{\frac{n}{k}}\right\rfloor.$$
Si además construimos las sumas prefijas
$$P(t)=\sum_{m\le t}W(m),$$
la contribución de todo el bloque es
$$\left(P(r)-P(\ell-1)\right)\left(Q(k)-1\right).$$
Ese es el acelerador decisivo: la cuenta de puntos del círculo se hace una sola vez por cociente distinto y no una vez por cada \(m\).
Aquí \(\lfloor\sqrt{10}\rfloor=3\), de modo que solo pueden contribuir \(m=1,2,3\). A partir de la regla local se obtiene
$$W(1)=1,\qquad W(2)=-1,\qquad W(3)=0.$$
Solo hacen falta dos cuentas geométricas:
$$Q(10)=37,\qquad Q(2)=9.$$
Por tanto,
$$4F(10)=1\cdot(37-1)+(-1)\cdot(9-1)+0\cdot(Q(1)-1)=36-8=28,$$
y en consecuencia
$$F(10)=7.$$
Ese resultado coincide con el pequeño punto de control reproducido por las implementaciones.
Las implementaciones en C++, Python y Java siguen el mismo esquema. Primero fijan \(M=\lfloor\sqrt{n}\rfloor\) y construyen una tabla de menor factor primo hasta \(M\). Con ella factorizan cada entero \(m\le M\) y calculan el peso multiplicativo \(W(m)\) a partir de las reglas de potencias primas.
Después forman las sumas prefijas de \(W\), de modo que cualquier suma por bloques \(\sum_{\ell\le m\le r}W(m)\) se obtiene en tiempo constante. A continuación particionan \(1\le m\le M\) en intervalos máximos donde \(\lfloor n/m^2\rfloor\) es constante. Para cada intervalo calculan \(Q(k)\), multiplican por la suma de pesos correspondiente, acumulan todas las contribuciones y al final dividen entre \(4\).
La versión en C++ evalúa en paralelo los conteos geométricos independientes de los distintos bloques. Las versiones en Python y Java realizan las mismas operaciones de manera secuencial. La matemática subyacente es exactamente la misma en los tres casos.
Si \(M=\lfloor\sqrt{n}\rfloor\), la preparación de los menores factores primos y de los pesos \(W(m)\) usa \(O(M)\) memoria y tiempo casi lineal en \(M\). Las sumas prefijas también cuestan \(O(M)\).
El número de valores distintos de \(\lfloor n/m^2\rfloor\) es mucho menor que \(M\); de hecho es \(O(n^{1/3})\). Cada uno requiere una sola evaluación de \(Q(k)\), y esa evaluación tarda \(O(\sqrt{k})\) con memoria extra constante. En la práctica, la parte dominante es el conteo geométrico, por lo que la agrupación por cocientes y, en C++, la paralelización son las optimizaciones decisivas.
设 \(F(n)\) 表示满足
$$N(z)=a^2+b^2\le n$$
的非零平方自由高斯整数 \(z=a+bi\) 的个数，其中只差一个单位 \(\{\pm 1,\pm i\}\) 的两个数视为同一个对象。题目要求计算 \(F(10^{14})\)。
如果直接枚举圆盘 \(a^2+b^2\le n\) 内的所有格点，规模会完全失控。真正可行的做法是先把“平方自由”写成 \(\mathbb{Z}[i]\) 中的 Möbius 型反演公式，再把高斯除子的求和按普通整数范数归并，最后利用 \(\lfloor n/m^2\rfloor\) 在长区间上不变这一性质，减少圆内格点计数的次数。
先定义格点计数函数
$$Q(x)=\#\{(a,b)\in\mathbb{Z}^2:a^2+b^2\le x\}.$$
它表示范数不超过 \(x\) 的全部高斯整数个数，连原点 \(0\) 也算进去。因此最后的公式会使用 \(Q(x)-1\)，以便把原点剔除掉。
在唯一分解整环 \(\mathbb{Z}[i]\) 中，平方自由指示函数与普通整数情形完全类似：
$$\mathbf{1}_{\mathrm{sqfree}}(z)=\sum_{d^2\mid z}\mu_{\mathbb{Z}[i]}(d).$$
这里 \(\mu_{\mathbb{Z}[i]}(d)\) 是高斯整数上的 Möbius 函数。如果 \(d\) 被某个高斯素数的平方整除，它就等于 \(0\)；否则若 \(d\) 是 \(k\) 个互不相同的高斯素数的乘积（忽略单位因子），它就等于 \((-1)^k\)。
把这个指示函数对所有满足 \(0<N(z)\le n\) 的高斯整数求和，就得到
$$4F(n)=\sum_{\substack{z\in\mathbb{Z}[i]\\0<N(z)\le n}}\mathbf{1}_{\mathrm{sqfree}}(z)=\sum_{d\ne 0}\mu_{\mathbb{Z}[i]}(d)\left(Q\!\left(\left\lfloor\frac{n}{N(d)^2}\right\rfloor\right)-1\right).$$
前面的 \(4\) 来自单位群 \(\{\pm 1,\pm i\}\)：每个非零高斯整数恰好有四个 associate，因此若按 associate 类计数，就要在最后除以 \(4\)。
上式中几何部分只依赖于 \(N(d)\)，而不依赖于 \(d\) 的具体形状。所以可以把所有范数相同的高斯除子合并，定义
$$W(m)=\sum_{N(d)=m}\mu_{\mathbb{Z}[i]}(d).$$
这样就得到核心公式
$$4F(n)=\sum_{m\le \sqrt{n}} W(m)\left(Q\!\left(\left\lfloor\frac{n}{m^2}\right\rfloor\right)-1\right).$$
这一步非常关键，因为原本是在 \(\mathbb{Z}[i]\) 中对高斯除子求和，现在变成了对普通整数 \(m\) 求和。程序的整体结构也就一目了然了：先构造权函数 \(W(m)\)，再计算圆内格点数 \(Q(\cdot)\)，最后把两者组合起来。
函数 \(W\) 是乘法性的，所以只需求出每个素数幂 \(p^e\) 上的取值 \(W(p^e)\)。这些取值直接来自普通素数在 \(\mathbb{Z}[i]\) 中的三种行为。
第一种是 \(p=2\)。它在高斯整数中是分歧的，并满足
$$2=-i(1+i)^2.$$
因此在 \(2\) 上方只有一个高斯素数，范数为 \(2\)。对应的局部生成多项式是
$$1-t.$$
第二种是 \(p\equiv 1\pmod{4}\)。这时 \(p\) 分裂成两个互不 associate 的共轭高斯素数，而且它们的范数都等于 \(p\)。局部平方自由选择只有三种：一个都不选，选其中一个，或者两个都选。所以局部多项式是
$$\left(1-t\right)^2=1-2t+t^2.$$
第三种是 \(p\equiv 3\pmod{4}\)。这种素数在 \(\mathbb{Z}[i]\) 中保持高斯素数，但它的范数变成 \(p^2\)。于是局部多项式为
$$1-t^2.$$
从这些系数直接读出
$$W(p^e)= \begin{cases} -1, & p=2,\ e=1,\\ -2, & p\equiv 1\pmod{4},\ e=1,\\ 1, & p\equiv 1\pmod{4},\ e=2,\\ -1, & p\equiv 3\pmod{4},\ e=2,\\ 0, & \text{otherwise}. \end{cases}$$
然后由乘法性即可得到任意 \(m\) 的权重 \(W(m)\)。这正是实现中最核心的数论部分。
几何部分由 \(Q(x)\) 负责，程序采用的精确公式是
$$Q(x)=1+4\sum_{a=0}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\sqrt{x-a^2}\right\rfloor.$$
对每个固定的 \(a\ge 0\)，内部的取整平方根给出满足 \(a^2+b^2\le x\) 的最大正整数 \(b\)。随着 \(a\) 增大，这个最大 \(b\) 只会下降，不会回升，因此可以用单调双指针扫描整条边界，在 \(O(\sqrt{x})\) 时间内完成整个求和，而不必对每个 \(a\) 单独二分或重复查找。
在主和式里，
$$k=\left\lfloor\frac{n}{m^2}\right\rfloor$$
对一整段连续的 \(m\) 都保持不变。如果当前区间从 \(\ell\) 开始，那么它的右端点就是
$$r=\left\lfloor\sqrt{\frac{n}{k}}\right\rfloor.$$
再预处理前缀和
$$P(t)=\sum_{m\le t}W(m),$$
那么整个区间的贡献可一次写成
$$\left(P(r)-P(\ell-1)\right)\left(Q(k)-1\right).$$
这一步是提速的关键。原本若按 \(m=1,2,3,\dots\) 逐项计算，会重复计算很多相同的 \(Q(k)\)；现在每个不同的商值只需要算一次。
此时 \(\lfloor\sqrt{10}\rfloor=3\)，所以只可能出现 \(m=1,2,3\)。按照上面的局部规则，有
$$W(1)=1,\qquad W(2)=-1,\qquad W(3)=0.$$
接下来只需要两个圆盘计数：
$$Q(10)=37,\qquad Q(2)=9.$$
因此
$$4F(10)=1\cdot(37-1)+(-1)\cdot(9-1)+0\cdot(Q(1)-1)=36-8=28,$$
从而得到
$$F(10)=7.$$
这个结果和实现中的小规模校验值完全一致，也说明公式和程序的逻辑是吻合的。
C++、Python 和 Java 三个实现都遵循同一条路线。第一步是取 \(M=\lfloor\sqrt{n}\rfloor\)，并建立到 \(M\) 为止的最小素因子表。这样每个 \(m\le M\) 都能被快速分解，再根据上面的素数幂规则得到乘法性权重 \(W(m)\)。
第二步是构造 \(W\) 的前缀和，这样任意区间 \([\ell,r]\) 上的权重和都能在常数时间内取出。第三步是把 \(1\le m\le M\) 分解成若干个最大区间，使得 \(\lfloor n/m^2\rfloor\) 在每个区间内保持不变。对每个区间只计算一次 \(Q(k)\)，再乘上该区间的权重和并累加。最后把总和除以 \(4\)，就得到按 associate 类计数的答案。
C++ 实现还把不同区间所需的圆盘计数并行计算，因为这些计算彼此独立。Python 和 Java 实现则按顺序执行同样的数学步骤。三种语言只是实现方式不同，算法本身完全一致。
记 \(M=\lfloor\sqrt{n}\rfloor\)。最小素因子表、权重数组和前缀和数组都需要 \(O(M)\) 内存，预处理时间相对于 \(M\) 来说接近线性。
\(\lfloor n/m^2\rfloor\) 的不同取值个数远小于 \(M\)，数量级为 \(O(n^{1/3})\)。每个不同取值只需一次 \(Q(k)\) 计算，而单次圆盘计数的时间是 \(O(\sqrt{k})\)，额外空间为常数。因此实际运行时间主要由几何计数主导，而区间分组以及 C++ 版本中的并行化正是最重要的优化来源。
Пусть \(F(n)\) обозначает число ненулевых квадратсвободных гауссовых целых \(z=a+bi\), для которых
$$N(z)=a^2+b^2\le n,$$
причём числа, отличающиеся только умножением на единицу из \(\{\pm 1,\pm i\}\), считаются одной и той же сущностью. Требуется найти \(F(10^{14})\).
Прямой перебор всех узлов решётки внутри круга \(a^2+b^2\le n\) слишком дорог. Поэтому решение сначала выражает квадратсвободность через Möbius-подобную формулу включения-исключения в \(\mathbb{Z}[i]\), затем сводит сумму по гауссовым делителям к обычной сумме по их нормам и, наконец, ускоряет геометрическую часть, группируя одинаковые значения \(\lfloor n/m^2\rfloor\).
Введём функцию подсчёта решётчатых точек
$$Q(x)=\#\{(a,b)\in\mathbb{Z}^2:a^2+b^2\le x\}.$$
Она считает все гауссовы целые с нормой не больше \(x\), включая \(0\). Поэтому в окончательной формуле используется \(Q(x)-1\), чтобы исключить начало координат.
В области уникального разложения \(\mathbb{Z}[i]\) индикатор квадратсвободности имеет тот же вид, что и над обычными целыми:
$$\mathbf{1}_{\mathrm{sqfree}}(z)=\sum_{d^2\mid z}\mu_{\mathbb{Z}[i]}(d).$$
Здесь \(\mu_{\mathbb{Z}[i]}(d)\) равна \(0\), если \(d\) делится на квадрат гауссова простого, и равна \((-1)^k\), если \(d\) является произведением \(k\) различных гауссовых простых с точностью до единиц.
Суммируя по всем ненулевым гауссовым целым с нормой не больше \(n\), получаем
$$4F(n)=\sum_{\substack{z\in\mathbb{Z}[i]\\0<N(z)\le n}}\mathbf{1}_{\mathrm{sqfree}}(z)=\sum_{d\ne 0}\mu_{\mathbb{Z}[i]}(d)\left(Q\!\left(\left\lfloor\frac{n}{N(d)^2}\right\rfloor\right)-1\right).$$
Множитель \(4\) появляется потому, что у каждого ненулевого гауссова целого ровно четыре ассоциированных значения.
В геометрической части важна только норма \(N(d)\), а не сам делитель. Поэтому объединим все гауссовы делители с одной и той же нормой и определим
$$W(m)=\sum_{N(d)=m}\mu_{\mathbb{Z}[i]}(d).$$
Тогда основная формула переписывается так:
$$4F(n)=\sum_{m\le \sqrt{n}} W(m)\left(Q\!\left(\left\lfloor\frac{n}{m^2}\right\rfloor\right)-1\right).$$
Именно это равенство объясняет устройство программы: сначала строится арифметический вес \(W(m)\), затем вычисляется функция круга \(Q(\cdot)\), а потом эти части объединяются.
Функция \(W\) мультипликативна, поэтому достаточно знать значения \(W(p^e)\). Они напрямую следуют из того, как рациональные простые ведут себя в \(\mathbb{Z}[i]\).
Для \(p=2\) имеет место ветвление:
$$2=-i(1+i)^2.$$
Над числом \(2\) есть только один гауссов простой с нормой \(2\), так что локальный многочлен равен
$$1-t.$$
Для \(p\equiv 1\pmod{4}\) простой распадается на два сопряжённых неассоциированных гауссовых простых нормы \(p\). Локально квадратсвободные варианты таковы: не брать ни один, взять один из двух или взять оба. Поэтому локальный многочлен имеет вид
$$\left(1-t\right)^2=1-2t+t^2.$$
Для \(p\equiv 3\pmod{4}\) простой остаётся гауссово простым, но его норма становится равной \(p^2\). Отсюда локальный фактор
$$1-t^2.$$
По коэффициентам получаем
$$W(p^e)= \begin{cases} -1, & p=2,\ e=1,\\ -2, & p\equiv 1\pmod{4},\ e=1,\\ 1, & p\equiv 1\pmod{4},\ e=2,\\ -1, & p\equiv 3\pmod{4},\ e=2,\\ 0, & \text{otherwise}. \end{cases}$$
Мультипликативность затем определяет \(W(m)\) для любого \(m\).
Геометрическая часть задаётся функцией \(Q(x)\). Реализация использует точную формулу
$$Q(x)=1+4\sum_{a=0}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\sqrt{x-a^2}\right\rfloor.$$
Для каждого \(a\ge 0\) внутренний пол даёт максимальный положительный \(b\), удовлетворяющий неравенству \(a^2+b^2\le x\). При увеличении \(a\) это максимальное \(b\) никогда не возрастает, поэтому вся сумма вычисляется монотонным двухуказательным проходом за \(O(\sqrt{x})\).
В главной сумме величина
$$k=\left\lfloor\frac{n}{m^2}\right\rfloor$$
остаётся постоянной на целых блоках подряд идущих значений \(m\). Если текущий блок начинается с \(\ell\), его правая граница равна
$$r=\left\lfloor\sqrt{\frac{n}{k}}\right\rfloor.$$
Если заранее построить префиксные суммы
$$P(t)=\sum_{m\le t}W(m),$$
то вклад всего блока равен
$$\left(P(r)-P(\ell-1)\right)\left(Q(k)-1\right).$$
Это и даёт ускорение: дорогое вычисление числа решётчатых точек делается один раз для каждого различного значения частного, а не один раз на каждый \(m\).
Здесь \(\lfloor\sqrt{10}\rfloor=3\), поэтому возможны только \(m=1,2,3\). Из локального правила получаем
$$W(1)=1,\qquad W(2)=-1,\qquad W(3)=0.$$
Нужны лишь две геометрические величины:
$$Q(10)=37,\qquad Q(2)=9.$$
Следовательно,
$$4F(10)=1\cdot(37-1)+(-1)\cdot(9-1)+0\cdot(Q(1)-1)=36-8=28,$$
и потому
$$F(10)=7.$$
Это совпадает с малым контрольным значением, которое воспроизводят реализации.
Реализации на C++, Python и Java следуют одному и тому же плану. Сначала берётся \(M=\lfloor\sqrt{n}\rfloor\) и строится таблица наименьших простых делителей до \(M\). Это позволяет быстро факторизовать каждое \(m\le M\) и вычислить мультипликативный вес \(W(m)\) по локальным правилам для простых степеней.
Затем строятся префиксные суммы \(W\), благодаря чему любая сумма по блоку \(\sum_{\ell\le m\le r}W(m)\) получается за константное время. После этого интервал \(1\le m\le M\) разбивается на максимальные блоки, на которых \(\lfloor n/m^2\rfloor\) постоянно. Для каждого блока вычисляется \(Q(k)\), умножается на соответствующую сумму весов, все вклады складываются, и в самом конце результат делится на \(4\).
В версии на C++ независимые вычисления функции круга для разных блоков выполняются параллельно. В версиях на Python и Java та же арифметика выполняется последовательно. Математическая формула во всех трёх языках одна и та же.
Если \(M=\lfloor\sqrt{n}\rfloor\), то данные о наименьших простых делителях, веса \(W(m)\) и префиксные суммы требуют \(O(M)\) памяти и почти линейного времени по \(M\).
Число различных значений \(\lfloor n/m^2\rfloor\) намного меньше, чем \(M\); оно имеет порядок \(O(n^{1/3})\). Для каждого такого значения нужна ровно одна оценка \(Q(k)\), а она стоит \(O(\sqrt{k})\) времени при постоянной дополнительной памяти. На практике доминирует именно геометрическая часть, поэтому группировка одинаковых частных и параллелизм в C++ дают основной выигрыш.
ليكن \(F(n)\) عدد الأعداد الصحيحة الغاوسية غير الصفرية والخالية من المربعات \(z=a+bi\) التي تحقق
$$N(z)=a^2+b^2\le n,$$
مع اعتبار العددين متماثلين إذا اختلفا فقط بضرب أحدهما في وحدة من \(\{\pm 1,\pm i\}\). المطلوب هو حساب \(F(10^{14})\).
العد المباشر لكل نقاط الشبكة داخل القرص \(a^2+b^2\le n\) غير عملي إطلاقًا. لذلك يعتمد الحل على تحويل شرط الخلو من المربعات إلى صيغة إدراج واستبعاد باستخدام دالة Möbius على \(\mathbb{Z}[i]\)، ثم ضغط الجمع على القواسم الغاوسية بحسب معيارها الصحيح، وأخيرًا تقليل عدد مرات عدّ النقاط داخل الدوائر عبر تجميع القيم المتساوية لـ \(\lfloor n/m^2\rfloor\).
نعرّف أولًا دالة عدّ نقاط الشبكة
$$Q(x)=\#\{(a,b)\in\mathbb{Z}^2:a^2+b^2\le x\}.$$
هذه الدالة تعدّ جميع الأعداد الغاوسية ذات المعيار الذي لا يتجاوز \(x\)، بما في ذلك الصفر. ولهذا تظهر في الصيغة النهائية العبارة \(Q(x)-1\) حتى نستبعد الأصل.
في الحلقة ذات التحليل الفريد \(\mathbb{Z}[i]\)، يأخذ مؤشر الخلو من المربعات الصورة نفسها المعروفة على الأعداد الصحيحة العادية:
$$\mathbf{1}_{\mathrm{sqfree}}(z)=\sum_{d^2\mid z}\mu_{\mathbb{Z}[i]}(d).$$
حيث \(\mu_{\mathbb{Z}[i]}(d)\) هي دالة Möbius الغاوسية. فهي تساوي \(0\) إذا كان \(d\) قابلًا للقسمة على مربع أولي غاوسي، وتساوي \((-1)^k\) إذا كان \(d\) حاصل ضرب \(k\) أوليات غاوسية مختلفة مع تجاهل الوحدات.
بجمع هذا المؤشر على جميع الأعداد الغاوسية غير الصفرية ذات المعيار حتى \(n\)، نحصل على
$$4F(n)=\sum_{\substack{z\in\mathbb{Z}[i]\\0<N(z)\le n}}\mathbf{1}_{\mathrm{sqfree}}(z)=\sum_{d\ne 0}\mu_{\mathbb{Z}[i]}(d)\left(Q\!\left(\left\lfloor\frac{n}{N(d)^2}\right\rfloor\right)-1\right).$$
العامل \(4\) يظهر لأن كل عدد غاوسي غير صفري يملك أربع قيم associate ناتجة عن الضرب في الوحدات الأربع.
الجزء الهندسي يعتمد فقط على \(N(d)\)، لا على القاسم الغاوسي نفسه. لذلك نجمع كل القواسم التي لها المعيار نفسه في وزن واحد، ونعرف
$$W(m)=\sum_{N(d)=m}\mu_{\mathbb{Z}[i]}(d).$$
عندها تصبح الصيغة الرئيسية
$$4F(n)=\sum_{m\le \sqrt{n}} W(m)\left(Q\!\left(\left\lfloor\frac{n}{m^2}\right\rfloor\right)-1\right).$$
هذه هي البنية الحقيقية للخوارزمية: حساب الوزن الحسابي \(W(m)\)، ثم حساب \(Q(\cdot)\)، ثم دمجهما في مجموع واحد.
الدالة \(W\) ضربيّة، لذا يكفي معرفة \(W(p^e)\). وهذه القيم تأتي مباشرة من الأنماط الثلاثة المعروفة لتفكك الأوليات الصحيحة داخل \(\mathbb{Z}[i]\).
بالنسبة إلى \(p=2\)، يوجد تشعب:
$$2=-i(1+i)^2.$$
أي أن فوق \(2\) يوجد أولي غاوسي واحد فقط معياره \(2\)، ولذلك يكون العامل المحلي
$$1-t.$$
أما إذا كان \(p\equiv 1\pmod{4}\)، فإن \(p\) ينشطر إلى أوليين غاوسيين مترافقين وغير associate، وكلاهما معياره \(p\). حينئذ تكون الخيارات المحلية الخالية من المربعات هي: عدم اختيار أي منهما، أو اختيار واحد منهما، أو اختيارهما معًا. لذا يكون كثير الحدود المحلي
$$\left(1-t\right)^2=1-2t+t^2.$$
وإذا كان \(p\equiv 3\pmod{4}\)، فإنه يبقى أوليًا غاوسيًا، لكن معياره يصبح \(p^2\). ومن ثم يكون العامل المحلي
$$1-t^2.$$
وبقراءة المعاملات نحصل على
$$W(p^e)= \begin{cases} -1, & p=2,\ e=1,\\ -2, & p\equiv 1\pmod{4},\ e=1,\\ 1, & p\equiv 1\pmod{4},\ e=2,\\ -1, & p\equiv 3\pmod{4},\ e=2,\\ 0, & \text{otherwise}. \end{cases}$$
وبعد ذلك تحدد الضربيّة قيمة \(W(m)\) لكل \(m\) طبيعي.
الجزء الهندسي هو الدالة \(Q(x)\)، وتستخدم التطبيقات الهوية الدقيقة
$$Q(x)=1+4\sum_{a=0}^{\lfloor\sqrt{x}\rfloor}\left\lfloor\sqrt{x-a^2}\right\rfloor.$$
لكل \(a\ge 0\)، تعطي القيمة الداخلية أكبر \(b>0\) يمكن أن يحقق \(a^2+b^2\le x\). ومع ازدياد \(a\) فإن هذا الحد الأعلى لـ \(b\) لا يزداد أبدًا، ولذلك يكفي مسح أحادي الاتجاه بمؤشرين لحساب المجموع كله في \(O(\sqrt{x})\) زمنًا.
في المجموع الرئيسي تظهر الكمية
$$k=\left\lfloor\frac{n}{m^2}\right\rfloor$$
وهي تبقى ثابتة على فترات كاملة من قيم \(m\) المتتالية. إذا بدأنا من \(\ell\)، فإن نهاية الفترة هي
$$r=\left\lfloor\sqrt{\frac{n}{k}}\right\rfloor.$$
ومع مجاميع بادئة
$$P(t)=\sum_{m\le t}W(m),$$
تصبح مساهمة الفترة كلها
$$\left(P(r)-P(\ell-1)\right)\left(Q(k)-1\right).$$
وهنا يكمن التسريع الأساسي: العد الهندسي المكلف لا يُحسب لكل \(m\)، بل مرة واحدة فقط لكل قيمة متميزة من قيم خارج القسمة.
لدينا هنا \(\lfloor\sqrt{10}\rfloor=3\)، لذا لا يمكن أن تسهم إلا القيم \(m=1,2,3\). ومن قاعدة القوى الأولية نجد
$$W(1)=1,\qquad W(2)=-1,\qquad W(3)=0.$$
ثم نحتاج إلى قيمتين فقط من دالة العد الهندسي:
$$Q(10)=37,\qquad Q(2)=9.$$
وبالتالي
$$4F(10)=1\cdot(37-1)+(-1)\cdot(9-1)+0\cdot(Q(1)-1)=36-8=28,$$
ومن ثم
$$F(10)=7.$$
وهذه القيمة تطابق نقطة التحقق الصغيرة التي تنتجها التطبيقات.
تتبع تطبيقات C++ وPython وJava الخطة نفسها. أولًا تؤخذ \(M=\lfloor\sqrt{n}\rfloor\) ثم يُبنى جدول أصغر عامل أولي حتى \(M\). بهذا يمكن تحليل كل \(m\le M\) بسرعة إلى عوامله الأولية، ثم حساب الوزن الضربي \(W(m)\) من قواعد القوى الأولية السابقة.
بعد ذلك تُبنى المجاميع البادئة لـ \(W\)، فتُستخرج أي محصلة على فترة \([\ell,r]\) في زمن ثابت. ثم تُقسَّم القيم \(1\le m\le M\) إلى فترات قصوى يكون فيها \(\lfloor n/m^2\rfloor\) ثابتًا. لكل فترة تُحسب \(Q(k)\) مرة واحدة، ثم تُضرب في مجموع الأوزان المناسب، وتُجمَع كل المساهمات، وفي النهاية يُقسَم الناتج على \(4\).
تنفذ نسخة C++ حسابات الدوائر المستقلة الخاصة بالفترات المختلفة على التوازي، بينما تنفذ نسختا Python وJava العمليات نفسها بصورة متسلسلة. الصيغة الرياضية التي تُحسب واحدة في اللغات الثلاث.
إذا كان \(M=\lfloor\sqrt{n}\rfloor\)، فإن إعداد جدول أصغر عامل أولي مع أوزان \(W(m)\) ومجاميعها البادئة يحتاج إلى ذاكرة \(O(M)\) وزمن شبه خطي في \(M\).
عدد القيم المختلفة لـ \(\lfloor n/m^2\rfloor\) أصغر بكثير من \(M\)، ورتبته \(O(n^{1/3})\). وكل قيمة مختلفة تتطلب حسابًا واحدًا فقط لـ \(Q(k)\)، وهذا الحساب يكلف \(O(\sqrt{k})\) زمنًا مع ذاكرة إضافية ثابتة. عمليًا يهيمن العد الهندسي على زمن التنفيذ، ولهذا فإن تجميع الفترات، ومعه التنفيذ المتوازي في C++، هما أهم عنصرين في التسريع.