The solution works with a canonical parametrization by primitive integer pairs \((p,q)\) with \(q\ge 0\) and \(\gcd(p,q)=1\). In that parametrization, the size attached to a solution is the discriminant-\(13\) binary quadratic form
$$z=\left|3p^2-7pq+3q^2\right|,$$
and we must count all admissible primitive pairs for which \(z\le N\). A second quadratic coordinate must stay positive as well, so the ratio \(p/q\) is not free: only three regions survive after the sign and symmetry reductions.
The isolated case \(q=0\) gives \(z=3\). Every other contribution comes from \(q\ge 1\), which is why the implementations loop over \(q\) and count valid \(k\)-intervals instead of scanning the original search space directly.
Introduce the two quadratic forms
$$Q(p,q)=3p^2-7pq+3q^2,\qquad X(p,q)=p^2-6pq+6q^2.$$
The counted size is \(|Q(p,q)|\), while admissibility requires \(X(p,q)>0\), \(|p|>q\), and \(\gcd(p,q)=1\). The code turns those conditions into interval counts in a new variable \(k\).
Because \(|p|>q\), every admissible pair with \(q>0\) can be written in exactly one of the forms
$$p=q+k,\qquad p=-q-k.$$
with \(k\ge 1\). Coprimality is preserved under this substitution:
$$\gcd(q,q+k)=\gcd(q,k),\qquad \gcd(q,-q-k)=\gcd(q,k).$$
So after fixing \(q\), the arithmetic condition becomes simply \(\gcd(q,k)=1\).
For \(p=q+k\), the auxiliary form becomes
$$X(q+k,q)=q^2-4qk+k^2=(k-(2-\sqrt3)q)(k-(2+\sqrt3)q).$$
If we write
$$\alpha=2-\sqrt3,\qquad \beta=2+\sqrt3,$$
then \(X>0\) forces either \(k<\alpha q\) or \(k>\beta q\). These are the two branches called region A and region D:
$$z_A(q,k)=-(Q(q+k,q))=q^2+qk-3k^2,$$
$$z_D(q,k)=Q(q+k,q)=3k^2-qk-q^2.$$
For \(p=-q-k\), we get
$$X(-q-k,q)=13q^2+8qk+k^2>0,$$
so positivity is automatic, and the size becomes
$$z_C(q,k)=Q(-q-k,q)=13q^2+13qk+3k^2.$$
Thus the entire count is the sum of region A, region C, region D, and the single base solution \(z=3\).
For fixed \(q\), define
$$C_q(K)=\#\{1\le k\le K:\gcd(q,k)=1\}.$$
Using inclusion-exclusion over the prime divisors of \(q\),
$$C_q(K)=\sum_{d\mid q}\mu(d)\left\lfloor\frac{K}{d}\right\rfloor.$$
Only squarefree divisors matter, so after factoring \(q\) once, the implementations generate all \(2^{\omega(q)}\) squarefree divisors together with their Möbius signs. Any interval \([L,U]\) is then counted by
$$C_q(U)-C_q(L-1).$$
The discriminant-\(13\) form has a simple congruence:
$$Q(p,q)\equiv 3(p+q)^2 \pmod{13}.$$
Therefore
$$13\mid Q(p,q)\iff p\equiv -q \pmod{13}.$$
The implementations keep only the canonical branch with \(13\nmid Q(p,q)\), so one residue class must be removed:
$$p=q+k \Rightarrow k\equiv -2q \pmod{13},$$
$$p=-q-k \Rightarrow k\equiv 0 \pmod{13}.$$
When \(q\equiv 0\pmod{13}\), coprimality already forbids \(k\equiv 0\pmod{13}\), so there is no extra subtraction. Otherwise the forbidden class is counted with the same Möbius table, but now on an arithmetic progression modulo \(13\).
Each region contributes an interval in \(k\).
Region A starts with \(1\le k\le \lfloor \alpha q\rfloor\). The polynomial \(z_A(q,k)=q^2+qk-3k^2\) is concave, so the inequality \(z_A\le N\) can fail only on a middle interval. Solving \(z_A=N\) gives
$$k=\frac{q\pm\sqrt{13q^2-12N}}{6}.$$
If \(13q^2\le 12N\), the whole region survives; otherwise the integers between the two roots must be removed.
Region C is monotone increasing, so its upper bound comes from
$$13q^2+13qk+3k^2\le N\quad\Rightarrow\quad k\le \frac{\sqrt{13q^2+12N}-13q}{6}.$$
Region D is also monotone once \(k>\beta q\), so it uses
$$k\ge \lfloor \beta q\rfloor+1,\qquad k\le \frac{\sqrt{13q^2+12N}+q}{6}.$$
The implementations compute these bounds from integer square roots and then adjust the endpoints by a few exact checks, eliminating rounding errors.
Here \(\lfloor \alpha\rfloor=0\), so region A is empty.
Region C satisfies
$$13+13k+3k^2\le 100,$$
which gives \(k=1,2,3\), hence the contributions \(29,51,79\).
Region D starts at \(k=\lfloor \beta\rfloor+1=4\). There we need
$$3k^2-k-1\le 100,$$
so \(k=4,5\), giving \(43\) and \(69\).
For \(q=1\), the forbidden classes are \(k\equiv 11\pmod{13}\) in regions A and D, and \(k\equiv 0\pmod{13}\) in region C, so no subtraction occurs in this layer. Together with the isolated case \(q=0\), this shows exactly how the counting proceeds.
The C++, Python, and Java implementations first set \(Q_{\max}=\lfloor\sqrt N\rfloor\) and build a smallest-prime-factor table up to that limit. This lets them factor every \(q\) quickly and generate all squarefree divisors needed for the Möbius sums.
For each \(q\), the implementation evaluates the three regions separately. Region A is handled as a prefix count up to \(\lfloor\alpha q\rfloor\) minus the middle interval where \(z_A>N\). Regions C and D are monotone, so they are counted with direct prefix differences. In every case the coprime filter and, when necessary, the forbidden residue class modulo \(13\) are applied through the same divisor table.
The C++ implementation additionally splits the outer \(q\)-range across worker threads and accumulates partial sums in parallel. The Python and Java implementations perform the same arithmetic serially. After all \(q\ge 1\) contributions are finished, the isolated primitive case \(z=3\) is added once \(N\ge 3\).
Let \(Q=\lfloor\sqrt N\rfloor\). Building the smallest-prime-factor sieve costs \(O(Q\log\log Q)\) time and \(O(Q)\) memory. For a fixed \(q\), the number of squarefree divisors is \(2^{\omega(q)}\), and each region uses only a constant number of sums over that table.
So the total running time is
$$O\left(Q\log\log Q+\sum_{q\le Q}2^{\omega(q)}\right),$$
which is close to \(O(Q\log Q)\) on average and behaves near-linearly in practice. Memory usage stays \(O(Q)\), with only small per-\(q\) temporary arrays besides the sieve.
Die Lösung arbeitet mit einer kanonischen Parametrisierung durch primitive ganzzahlige Paare \((p,q)\) mit \(q\ge 0\) und \(\gcd(p,q)=1\). In dieser Darstellung ist die zu zählende Größe die binäre quadratische Form mit Diskriminante \(13\)
$$z=\left|3p^2-7pq+3q^2\right|,$$
und gezählt werden alle zulässigen primitiven Paare mit \(z\le N\). Zusätzlich muss eine zweite quadratische Koordinate positiv bleiben; genau diese Bedingung zerlegt den Quotienten \(p/q\) in einige wenige relevante Bereiche.
Der Sonderfall \(q=0\) liefert \(z=3\). Alle übrigen Beiträge kommen von \(q\ge 1\), weshalb die Implementierungen \(q\)-Schichten und passende \(k\)-Intervalle zählen, statt den ursprünglichen Suchraum direkt zu durchsuchen.
Wir führen die beiden quadratischen Formen ein:
$$Q(p,q)=3p^2-7pq+3q^2,\qquad X(p,q)=p^2-6pq+6q^2.$$
Die gezählte Größe ist \(|Q(p,q)|\), während die Zulässigkeit \(X(p,q)>0\), \(|p|>q\) und \(\gcd(p,q)=1\) verlangt. Der Code übersetzt diese Bedingungen in Intervallzählungen für eine neue Variable \(k\).
Wegen \(|p|>q\) kann jedes zulässige Paar mit \(q>0\) eindeutig in einer der beiden Formen geschrieben werden:
$$p=q+k,\qquad p=-q-k.$$
mit \(k\ge 1\). Die Teilerfremdheit bleibt dabei erhalten:
$$\gcd(q,q+k)=\gcd(q,k),\qquad \gcd(q,-q-k)=\gcd(q,k).$$
Für festes \(q\) reduziert sich die arithmetische Bedingung also auf \(\gcd(q,k)=1\).
Für \(p=q+k\) gilt
$$X(q+k,q)=q^2-4qk+k^2=(k-(2-\sqrt3)q)(k-(2+\sqrt3)q).$$
Mit
$$\alpha=2-\sqrt3,\qquad \beta=2+\sqrt3$$
erzwingt \(X>0\) entweder \(k<\alpha q\) oder \(k>\beta q\). Das sind die beiden Zweige A und D:
$$z_A(q,k)=-(Q(q+k,q))=q^2+qk-3k^2,$$
$$z_D(q,k)=Q(q+k,q)=3k^2-qk-q^2.$$
Für \(p=-q-k\) erhält man
$$X(-q-k,q)=13q^2+8qk+k^2>0,$$
also automatisch Positivität, und die Größe wird zu
$$z_C(q,k)=Q(-q-k,q)=13q^2+13qk+3k^2.$$
Damit besteht die gesamte Zählung aus Region A, Region C, Region D und dem einzelnen Basiswert \(z=3\).
Für festes \(q\) definieren wir
$$C_q(K)=\#\{1\le k\le K:\gcd(q,k)=1\}.$$
Mit Inklusion-Exklusion über die Primteiler von \(q\) folgt
$$C_q(K)=\sum_{d\mid q}\mu(d)\left\lfloor\frac{K}{d}\right\rfloor.$$
Nur quadratfreie Teiler sind relevant. Nach einer Faktorisierung von \(q\) erzeugen die Implementierungen daher alle \(2^{\omega(q)}\) quadratfreien Teiler samt Möbius-Vorzeichen. Ein Intervall \([L,U]\) wird dann durch
$$C_q(U)-C_q(L-1)$$
gezählt.
Für die Form mit Diskriminante \(13\) gilt die Kongruenz
$$Q(p,q)\equiv 3(p+q)^2 \pmod{13}.$$
Daher gilt
$$13\mid Q(p,q)\iff p\equiv -q \pmod{13}.$$
Die Implementierungen behalten nur den kanonischen Zweig mit \(13\nmid Q(p,q)\). Also muss genau eine Restklasse ausgeschlossen werden:
$$p=q+k \Rightarrow k\equiv -2q \pmod{13},$$
$$p=-q-k \Rightarrow k\equiv 0 \pmod{13}.$$
Falls \(q\equiv 0\pmod{13}\), verbietet die Teilerfremdheit bereits \(k\equiv 0\pmod{13}\); dann ist keine zusätzliche Subtraktion nötig. Andernfalls zählt man die verbotene Klasse mit derselben Möbius-Tabelle, jetzt aber auf einer arithmetischen Progression modulo \(13\).
Jede Region liefert ein Intervall in \(k\).
Region A beginnt mit \(1\le k\le \lfloor \alpha q\rfloor\). Das Polynom \(z_A(q,k)=q^2+qk-3k^2\) ist konkav, daher kann die Ungleichung \(z_A\le N\) nur in einem mittleren Intervall verletzt werden. Aus \(z_A=N\) folgen die Wurzeln
$$k=\frac{q\pm\sqrt{13q^2-12N}}{6}.$$
Ist \(13q^2\le 12N\), überlebt die ganze Region; sonst müssen die ganzen Zahlen zwischen diesen beiden Wurzeln entfernt werden.
Region C ist monoton wachsend, also kommt ihre obere Schranke aus
$$13q^2+13qk+3k^2\le N\quad\Rightarrow\quad k\le \frac{\sqrt{13q^2+12N}-13q}{6}.$$
Region D ist für \(k>\beta q\) ebenfalls monoton, also benutzt man
$$k\ge \lfloor \beta q\rfloor+1,\qquad k\le \frac{\sqrt{13q^2+12N}+q}{6}.$$
Die Implementierungen bestimmen diese Grenzen mit ganzzahligen Quadratwurzeln und korrigieren sie anschließend durch wenige exakte Tests, damit keine Rundungsfehler übrig bleiben.
Hier ist \(\lfloor \alpha\rfloor=0\), also ist Region A leer.
Für Region C gilt
$$13+13k+3k^2\le 100,$$
also \(k=1,2,3\) und damit die Beiträge \(29,51,79\).
Region D startet bei \(k=\lfloor \beta\rfloor+1=4\). Dort braucht man
$$3k^2-k-1\le 100,$$
also \(k=4,5\) und damit \(43\) sowie \(69\).
Für \(q=1\) lauten die verbotenen Restklassen \(k\equiv 11\pmod{13}\) in A und D sowie \(k\equiv 0\pmod{13}\) in C, daher wird in dieser Schicht nichts abgezogen. Zusammen mit dem isolierten Fall \(q=0\) sieht man hier exakt, wie die Zählung arbeitet.
Die C++-, Python- und Java-Implementierungen setzen zunächst \(Q_{\max}=\lfloor\sqrt N\rfloor\) und bauen bis zu dieser Grenze eine Tabelle des kleinsten Primteilers auf. Damit lässt sich jedes \(q\) schnell faktorisieren, und alle quadratfreien Teiler für die Möbius-Summen können direkt erzeugt werden.
Für jedes \(q\) werden die Regionen A, C und D getrennt ausgewertet. Region A wird als Präfixzählung bis \(\lfloor\alpha q\rfloor\) behandelt, von der anschließend das mittlere Intervall mit \(z_A>N\) abgezogen wird. Die Regionen C und D sind monoton und werden daher über direkte Präfixdifferenzen gezählt. In allen Fällen laufen der Teilerfremdheitsfilter und, falls nötig, die verbotene Restklasse modulo \(13\) über dieselbe Divisorentabelle.
Die C++-Implementierung teilt den äußeren \(q\)-Bereich zusätzlich auf mehrere Threads auf und summiert partielle Ergebnisse parallel. Python und Java führen dieselbe Arithmetik seriell aus. Nach allen Beiträgen mit \(q\ge 1\) wird der isolierte primitive Fall \(z=3\) genau einmal addiert, sobald \(N\ge 3\) ist.
Sei \(Q=\lfloor\sqrt N\rfloor\). Das SPF-Sieb kostet \(O(Q\log\log Q)\) Zeit und \(O(Q)\) Speicher. Für festes \(q\) ist die Anzahl quadratfreier Teiler gleich \(2^{\omega(q)}\), und jede Region benötigt nur konstant viele Summen über diese Tabelle.
Damit ergibt sich insgesamt die Laufzeit
$$O\left(Q\log\log Q+\sum_{q\le Q}2^{\omega(q)}\right),$$
also im Mittel ungefähr \(O(Q\log Q)\) und praktisch sehr nahe an linear. Der Speicherverbrauch bleibt \(O(Q)\); zusätzlich gibt es nur kleine temporäre Felder pro \(q\).
Çözüm, \((p,q)\) biçimindeki asal aralarında tamsayı çiftlerinin kanonik bir parametrelemesini kullanır; burada \(q\ge 0\) ve \(\gcd(p,q)=1\) şartı vardır. Bu parametrelemede sayılan büyüklük, diskriminantı \(13\) olan ikili kuadratik formdur:
$$z=\left|3p^2-7pq+3q^2\right|.$$
Amaç, \(z\le N\) koşulunu sağlayan tüm geçerli asal çiftleri saymaktır. Ayrıca ikinci bir kuadratik koordinatın da pozitif kalması gerekir; bu da \(p/q\) oranını birkaç dar bölgeye ayırır ve tüm aramayı ciddi biçimde küçültür.
\(q=0\) olan tekil durum \(z=3\) verir. Geri kalan tüm katkılar \(q\ge 1\) üzerinden geldiği için, uygulamalar doğrudan özgün arama uzayını taramak yerine her \(q\) için uygun \(k\)-aralıklarını sayar.
İki kuadratik form tanımlayalım:
$$Q(p,q)=3p^2-7pq+3q^2,\qquad X(p,q)=p^2-6pq+6q^2.$$
Sayılan büyüklük \(|Q(p,q)|\)'dir. Geçerlilik ise \(X(p,q)>0\), \(|p|>q\) ve \(\gcd(p,q)=1\) ister. Kod bu şartları yeni bir \(k\) değişkeni üzerinde aralık sayımına dönüştürür.
\(|p|>q\) olduğundan, \(q>0\) için her geçerli çift tam olarak şu iki biçimden birine sahiptir:
$$p=q+k,\qquad p=-q-k.$$
burada \(k\ge 1\). Aralarında asallık bu dönüşüm altında korunur:
$$\gcd(q,q+k)=\gcd(q,k),\qquad \gcd(q,-q-k)=\gcd(q,k).$$
Dolayısıyla sabit bir \(q\) için aritmetik koşul yalnızca \(\gcd(q,k)=1\) olur.
\(p=q+k\) için yardımcı form
$$X(q+k,q)=q^2-4qk+k^2=(k-(2-\sqrt3)q)(k-(2+\sqrt3)q)$$
şeklinde çarpanlara ayrılır. Şimdi
$$\alpha=2-\sqrt3,\qquad \beta=2+\sqrt3$$
dersek, \(X>0\) ancak \(k<\alpha q\) ya da \(k>\beta q\) iken mümkündür. Böylece A ve D bölgeleri elde edilir:
$$z_A(q,k)=-(Q(q+k,q))=q^2+qk-3k^2,$$
$$z_D(q,k)=Q(q+k,q)=3k^2-qk-q^2.$$
\(p=-q-k\) için ise
$$X(-q-k,q)=13q^2+8qk+k^2>0$$
otomatik olarak pozitiftir ve sayılan büyüklük
$$z_C(q,k)=Q(-q-k,q)=13q^2+13qk+3k^2$$
olur. Sonuç olarak toplam sayım, A bölgesi, C bölgesi, D bölgesi ve tekil temel değer \(z=3\)'ün toplamıdır.
Sabit bir \(q\) için
$$C_q(K)=\#\{1\le k\le K:\gcd(q,k)=1\}$$
tanımını yapalım. \(q\)'nun asal bölenleri üzerinde dahil et-çıkar uygulanınca
$$C_q(K)=\sum_{d\mid q}\mu(d)\left\lfloor\frac{K}{d}\right\rfloor$$
elde edilir. Burada sadece kareden arınmış bölenler önemlidir. Bu yüzden uygulamalar önce \(q\)'yu çarpanlarına ayırır, sonra da tüm \(2^{\omega(q)}\) kareden arınmış böleni ve Möbius işaretini üretir. Her \([L,U]\) aralığı
$$C_q(U)-C_q(L-1)$$
ile sayılır.
Diskriminantı \(13\) olan bu form için basit bir kongruans vardır:
$$Q(p,q)\equiv 3(p+q)^2 \pmod{13}.$$
Dolayısıyla
$$13\mid Q(p,q)\iff p\equiv -q \pmod{13}.$$
Uygulamalar yalnızca \(13\nmid Q(p,q)\) olan kanonik dalı tutar. Bu nedenle tam bir kalıntı sınıfı çıkarılır:
$$p=q+k \Rightarrow k\equiv -2q \pmod{13},$$
$$p=-q-k \Rightarrow k\equiv 0 \pmod{13}.$$
Eğer \(q\equiv 0\pmod{13}\) ise, \(\gcd(q,k)=1\) şartı zaten \(k\equiv 0\pmod{13}\) olasılığını yasaklar; dolayısıyla ek bir düzeltmeye ihtiyaç kalmaz. Aksi halde yasak kalıntı sınıfı da yine aynı Möbius tablosu ile, bu kez mod \(13\) aritmetik dizisi üzerinde sayılır.
Her bölge, \(k\) için bir aralık üretir.
A bölgesi başlangıçta \(1\le k\le \lfloor\alpha q\rfloor\) aralığıdır. \(z_A(q,k)=q^2+qk-3k^2\) aşağı bakan bir paraboldür; bu yüzden \(z_A\le N\) eşitsizliği yalnızca ortadaki bir aralıkta bozulabilir. \(z_A=N\) denkleminin kökleri
$$k=\frac{q\pm\sqrt{13q^2-12N}}{6}$$
olur. Eğer \(13q^2\le 12N\) ise bölgenin tamamı korunur; aksi halde bu iki kök arasındaki tamsayılar çıkarılır.
C bölgesi monoton artandır; üst sınır
$$13q^2+13qk+3k^2\le N\quad\Rightarrow\quad k\le \frac{\sqrt{13q^2+12N}-13q}{6}$$
ile gelir.
D bölgesi de \(k>\beta q\) sonrasında monotondur; burada
$$k\ge \lfloor\beta q\rfloor+1,\qquad k\le \frac{\sqrt{13q^2+12N}+q}{6}$$
kullanılır. Uygulamalar bu sınırları önce tam sayı karekökü ile hesaplar, sonra birkaç kesin kontrolle düzeltir; böylece yuvarlama hatası kalmaz.
Bu durumda \(\lfloor\alpha\rfloor=0\) olduğundan A bölgesi boştur.
C bölgesinde
$$13+13k+3k^2\le 100$$
eşitsizliği \(k=1,2,3\) verir; dolayısıyla katkılar \(29,51,79\) olur.
D bölgesi \(k=\lfloor\beta\rfloor+1=4\) noktasından başlar. Burada
$$3k^2-k-1\le 100$$
eşitsizliği \(k=4,5\) verir; katkılar \(43\) ve \(69\)'dur.
\(q=1\) için yasak sınıflar A ve D'de \(k\equiv 11\pmod{13}\), C'de ise \(k\equiv 0\pmod{13}\) olduğundan bu katmanda herhangi bir çıkarma yapılmaz. Tekil \(q=0\) durumu ile birlikte, bir \(q\)-katmanının nasıl sayıldığı açıkça görülür.
C++, Python ve Java uygulamaları önce \(Q_{\max}=\lfloor\sqrt N\rfloor\) değerini alır ve bu sınıra kadar en küçük asal bölen tablosu kurar. Böylece her \(q\) hızlıca çarpanlara ayrılır ve Möbius toplamlarında kullanılacak kareden arınmış bölenler doğrudan üretilir.
Ardından her \(q\) için A, C ve D bölgeleri ayrı ayrı değerlendirilir. A bölgesi, \(\lfloor\alpha q\rfloor\)'ya kadar olan ön-ek sayımından \(z_A>N\) olan orta aralık çıkarılarak elde edilir. C ve D bölgeleri monoton olduğundan doğrudan ön-ek farklarıyla sayılır. Her bölgede hem aralarında asallık filtresi hem de gerekirse mod \(13\) yasak kalıntı sınıfı aynı bölen tablosu üstünden uygulanır.
C++ sürümü dıştaki \(q\) aralığını iş parçacıklarına bölerek kısmı toplamları paralel toplar. Python ve Java aynı aritmetiği seri biçimde yürütür. Tüm \(q\ge 1\) katkıları bittikten sonra, \(N\ge 3\) ise izole temel durum \(z=3\) bir kez eklenir.
\(Q=\lfloor\sqrt N\rfloor\) olsun. En küçük asal bölen süzgeci \(O(Q\log\log Q)\) zaman ve \(O(Q)\) bellek ister. Sabit bir \(q\) için kareden arınmış bölen sayısı \(2^{\omega(q)}\)'dir ve her bölge bu tablo üzerinde sabit sayıda toplama yapar.
Böylece toplam çalışma süresi
$$O\left(Q\log\log Q+\sum_{q\le Q}2^{\omega(q)}\right)$$
olur; ortalama davranış yaklaşık \(O(Q\log Q)\), pratikte ise \(\sqrt N\)'de neredeyse lineerdir. Bellek kullanımı \(O(Q)\) olarak kalır; süzgeç dışında yalnızca küçük geçici diziler gerekir.
La solución usa una parametrización canónica mediante pares enteros primitivos \((p,q)\) con \(q\ge 0\) y \(\gcd(p,q)=1\). En esa parametrización, el tamaño que se cuenta está dado por la forma cuadrática binaria de discriminante \(13\)
$$z=\left|3p^2-7pq+3q^2\right|,$$
y debemos contar todos los pares primitivos admisibles con \(z\le N\). Además, otra coordenada cuadrática debe mantenerse positiva, y justamente esa condición divide la razón \(p/q\) en unas pocas regiones viables.
El caso aislado \(q=0\) produce \(z=3\). Todas las demás contribuciones vienen de \(q\ge 1\), por lo que las implementaciones recorren capas en \(q\) y cuentan intervalos válidos de \(k\) en lugar de explorar directamente el espacio original.
Introducimos las dos formas cuadráticas
$$Q(p,q)=3p^2-7pq+3q^2,\qquad X(p,q)=p^2-6pq+6q^2.$$
La magnitud contada es \(|Q(p,q)|\), mientras que la admisibilidad exige \(X(p,q)>0\), \(|p|>q\) y \(\gcd(p,q)=1\). El programa transforma esas condiciones en conteos por intervalos de una nueva variable \(k\).
Como \(|p|>q\), todo par admisible con \(q>0\) puede escribirse de manera única como
$$p=q+k,\qquad p=-q-k.$$
con \(k\ge 1\). La coprimalidad se conserva:
$$\gcd(q,q+k)=\gcd(q,k),\qquad \gcd(q,-q-k)=\gcd(q,k).$$
Así, para \(q\) fijo, la condición aritmética se reduce a \(\gcd(q,k)=1\).
Si \(p=q+k\), la forma auxiliar vale
$$X(q+k,q)=q^2-4qk+k^2=(k-(2-\sqrt3)q)(k-(2+\sqrt3)q).$$
Definiendo
$$\alpha=2-\sqrt3,\qquad \beta=2+\sqrt3,$$
la condición \(X>0\) obliga a \(k<\alpha q\) o bien a \(k>\beta q\). Eso produce las regiones A y D:
$$z_A(q,k)=-(Q(q+k,q))=q^2+qk-3k^2,$$
$$z_D(q,k)=Q(q+k,q)=3k^2-qk-q^2.$$
Si \(p=-q-k\), entonces
$$X(-q-k,q)=13q^2+8qk+k^2>0,$$
de modo que la positividad es automática y el tamaño queda
$$z_C(q,k)=Q(-q-k,q)=13q^2+13qk+3k^2.$$
Por tanto, el conteo total es la suma de las regiones A, C, D y el único caso base \(z=3\).
Para un \(q\) fijo, definimos
$$C_q(K)=\#\{1\le k\le K:\gcd(q,k)=1\}.$$
Aplicando inclusión-exclusión sobre los primos que dividen a \(q\), se obtiene
$$C_q(K)=\sum_{d\mid q}\mu(d)\left\lfloor\frac{K}{d}\right\rfloor.$$
Solo importan los divisores libres de cuadrados. Por eso, una vez factorizado \(q\), las implementaciones generan los \(2^{\omega(q)}\) divisores libres de cuadrados junto con sus signos de Möbius. Entonces cualquier intervalo \([L,U]\) se cuenta mediante
$$C_q(U)-C_q(L-1).$$
La forma de discriminante \(13\) satisface
$$Q(p,q)\equiv 3(p+q)^2 \pmod{13}.$$
Así,
$$13\mid Q(p,q)\iff p\equiv -q \pmod{13}.$$
Las implementaciones conservan solo la rama canónica con \(13\nmid Q(p,q)\), así que hay que restar exactamente una clase residual:
$$p=q+k \Rightarrow k\equiv -2q \pmod{13},$$
$$p=-q-k \Rightarrow k\equiv 0 \pmod{13}.$$
Si \(q\equiv 0\pmod{13}\), la coprimalidad ya impide \(k\equiv 0\pmod{13}\), de modo que no hace falta corrección adicional. En caso contrario, la clase prohibida se cuenta con la misma tabla de Möbius, ahora sobre una progresión aritmética módulo \(13\).
Cada región produce un intervalo en \(k\).
La región A comienza con \(1\le k\le \lfloor\alpha q\rfloor\). El polinomio \(z_A(q,k)=q^2+qk-3k^2\) es cóncavo, así que la desigualdad \(z_A\le N\) solo puede fallar en un intervalo central. Resolviendo \(z_A=N\), se obtiene
$$k=\frac{q\pm\sqrt{13q^2-12N}}{6}.$$
Si \(13q^2\le 12N\), toda la región es válida; en caso contrario, se eliminan los enteros comprendidos entre esas dos raíces.
La región C es creciente, así que su cota superior viene de
$$13q^2+13qk+3k^2\le N\quad\Rightarrow\quad k\le \frac{\sqrt{13q^2+12N}-13q}{6}.$$
La región D también es creciente una vez que \(k>\beta q\), por lo que usa
$$k\ge \lfloor\beta q\rfloor+1,\qquad k\le \frac{\sqrt{13q^2+12N}+q}{6}.$$
Las implementaciones calculan primero esas fronteras con raíces cuadradas enteras y luego las corrigen con unas pocas comprobaciones exactas, eliminando cualquier error de redondeo.
Aquí \(\lfloor\alpha\rfloor=0\), por lo que la región A está vacía.
En la región C tenemos
$$13+13k+3k^2\le 100,$$
que da \(k=1,2,3\), y por tanto las contribuciones \(29,51,79\).
La región D empieza en \(k=\lfloor\beta\rfloor+1=4\). Allí se necesita
$$3k^2-k-1\le 100,$$
de modo que \(k=4,5\), produciendo \(43\) y \(69\).
Para \(q=1\), las clases prohibidas son \(k\equiv 11\pmod{13}\) en A y D, y \(k\equiv 0\pmod{13}\) en C, así que en esta capa no se resta nada. Junto con el caso aislado \(q=0\), esto muestra exactamente cómo se organiza el conteo.
Las implementaciones en C++, Python y Java fijan primero \(Q_{\max}=\lfloor\sqrt N\rfloor\) y construyen una tabla de factor primo mínimo hasta ese límite. Eso permite factorizar cada \(q\) con rapidez y generar todos los divisores libres de cuadrados necesarios para las sumas de Möbius.
Luego cada \(q\) se procesa por separado en las regiones A, C y D. La región A se trata como un conteo prefijo hasta \(\lfloor\alpha q\rfloor\), del que se resta después el intervalo central donde \(z_A>N\). Las regiones C y D son monótonas, así que se cuentan con diferencias de prefijos. En todos los casos, el filtro de coprimalidad y, cuando hace falta, la clase residual prohibida módulo \(13\) se aplican usando la misma tabla de divisores.
La implementación en C++ además divide el rango exterior de \(q\) entre varios hilos y suma los totales parciales en paralelo. Las versiones en Python y Java realizan la misma aritmética de forma secuencial. Al final se añade una sola vez el caso primitivo aislado \(z=3\) cuando \(N\ge 3\).
Sea \(Q=\lfloor\sqrt N\rfloor\). Construir la criba de factores primos mínimos cuesta \(O(Q\log\log Q)\) tiempo y \(O(Q)\) memoria. Para un \(q\) fijo, el número de divisores libres de cuadrados es \(2^{\omega(q)}\), y cada región usa solo un número constante de sumas sobre esa tabla.
Por tanto, el tiempo total es
$$O\left(Q\log\log Q+\sum_{q\le Q}2^{\omega(q)}\right),$$
que en promedio está cerca de \(O(Q\log Q)\) y en la práctica se comporta casi linealmente en \(\sqrt N\). La memoria permanece en \(O(Q)\), con arreglos temporales pequeños además de la criba.
这道题的实现并不是直接在原始对象空间里暴力搜索，而是先把可行解改写成一组规范化的整数参数 \((p,q)\)，其中 \(q\ge 0\)、\(\gcd(p,q)=1\)。在这个参数化里，决定“大小”的核心量是判别式为 \(13\) 的二元二次型
$$z=\left|3p^2-7pq+3q^2\right|,$$
目标就是统计所有满足 \(z\le N\) 的规范原始参数对。除此之外，还有另一个二次型必须保持正值，因此 \(p/q\) 不能任意取值，最终只剩下三个有效区域需要计数。
当 \(q=0\) 时，会得到唯一的基例 \(z=3\)。除这个特例外，所有贡献都来自 \(q\ge 1\)，所以程序的主循环是按 \(q\) 分层，再在每一层里统计满足条件的 \(k\) 区间。
先引入两个二次型：
$$Q(p,q)=3p^2-7pq+3q^2,\qquad X(p,q)=p^2-6pq+6q^2.$$
真正计数的大小是 \(|Q(p,q)|\)，而可行性还要求 \(X(p,q)>0\)、\(|p|>q\)、以及 \(\gcd(p,q)=1\)。实现的关键，就是把这些条件系统地转写成关于一个新变量 \(k\) 的区间计数问题。
由于必须满足 \(|p|>q\)，所以对任意 \(q>0\) 的可行点，\(p\) 必然唯一地落在下面两种形式之一：
$$p=q+k,\qquad p=-q-k.$$
其中 \(k\ge 1\)。这个替换最大的好处是，互素条件会立刻简化：
$$\gcd(q,q+k)=\gcd(q,k),\qquad \gcd(q,-q-k)=\gcd(q,k).$$
也就是说，对固定的 \(q\) 来说，我们只需要统计满足 \(\gcd(q,k)=1\) 的正整数 \(k\)。
先看 \(p=q+k\) 这支。此时辅助二次型变成
$$X(q+k,q)=q^2-4qk+k^2=(k-(2-\sqrt3)q)(k-(2+\sqrt3)q).$$
设
$$\alpha=2-\sqrt3,\qquad \beta=2+\sqrt3,$$
则 \(X>0\) 当且仅当 \(k<\alpha q\) 或 \(k>\beta q\)。于是得到两个分支：
$$z_A(q,k)=-(Q(q+k,q))=q^2+qk-3k^2,$$
$$z_D(q,k)=Q(q+k,q)=3k^2-qk-q^2.$$
再看 \(p=-q-k\) 这一支。此时
$$X(-q-k,q)=13q^2+8qk+k^2>0,$$
永远自动成立，而被计数的大小变成
$$z_C(q,k)=Q(-q-k,q)=13q^2+13qk+3k^2.$$
因此，全部答案由 A 区、C 区、D 区，再加上孤立的基例 \(z=3\) 组成。
对固定的 \(q\)，定义
$$C_q(K)=\#\{1\le k\le K:\gcd(q,k)=1\}.$$
对 \(q\) 的质因子使用容斥，得到
$$C_q(K)=\sum_{d\mid q}\mu(d)\left\lfloor\frac{K}{d}\right\rfloor.$$
这里真正有用的只有无平方因子除数，所以实现先分解 \(q\)，再生成全部 \(2^{\omega(q)}\) 个无平方因子除数及其 Möbius 符号。于是任意区间 \([L,U]\) 内的互素个数都可以写成
$$C_q(U)-C_q(L-1).$$
这就是三种区域里所有前缀计数和区间差分的统一基础。
判别式为 \(13\) 的这个二次型有一个非常关键的同余关系：
$$Q(p,q)\equiv 3(p+q)^2 \pmod{13}.$$
因此
$$13\mid Q(p,q)\iff p\equiv -q \pmod{13}.$$
实现保留的是 \(13\nmid Q(p,q)\) 的规范分支，所以每层都要减掉一整类坏余数：
$$p=q+k \Rightarrow k\equiv -2q \pmod{13},$$
$$p=-q-k \Rightarrow k\equiv 0 \pmod{13}.$$
当 \(q\equiv 0\pmod{13}\) 时，互素条件已经自动排除了 \(k\equiv 0\pmod{13}\)，所以不需要额外修正。否则，就用同一张 Möbius 除数表在模 \(13\) 的等差数列上再做一次计数，把坏类减掉。
每个区域最终都会化成 \(k\) 的一个区间。
A 区的初始范围是 \(1\le k\le \lfloor\alpha q\rfloor\)。由于 \(z_A(q,k)=q^2+qk-3k^2\) 是开口向下的抛物线，不满足 \(z_A\le N\) 的部分只可能出现在中间一段。解方程 \(z_A=N\) 可得根
$$k=\frac{q\pm\sqrt{13q^2-12N}}{6}.$$
如果 \(13q^2\le 12N\)，说明整段都合法；否则就把两根之间的整数全部扣掉。
C 区单调递增，因此其上界由
$$13q^2+13qk+3k^2\le N\quad\Rightarrow\quad k\le \frac{\sqrt{13q^2+12N}-13q}{6}$$
给出。
D 区在 \(k>\beta q\) 之后也单调递增，所以使用
$$k\ge \lfloor\beta q\rfloor+1,\qquad k\le \frac{\sqrt{13q^2+12N}+q}{6}.$$
实现里先用整数平方根给出这些边界，再通过极少量的精确检查把端点调到完全正确的位置，从而避免浮点误差带来的 off-by-one 问题。
此时 \(\lfloor\alpha\rfloor=0\)，所以 A 区为空。
C 区满足
$$13+13k+3k^2\le 100,$$
得到 \(k=1,2,3\)，对应的贡献是 \(29,51,79\)。
D 区从 \(k=\lfloor\beta\rfloor+1=4\) 开始，此时需要
$$3k^2-k-1\le 100,$$
所以 \(k=4,5\)，对应 \(43\) 和 \(69\)。
对 \(q=1\) 来说，A、D 区的禁用余数类是 \(k\equiv 11\pmod{13}\)，C 区则是 \(k\equiv 0\pmod{13}\)，这一层里都不会命中，所以没有额外扣除。再加上 \(q=0\) 的基例，就可以清楚看到程序是如何逐层累加答案的。
C++、Python 和 Java 实现都会先取 \(Q_{\max}=\lfloor\sqrt N\rfloor\)，并建立到这个上界为止的最小质因子表。这样每个 \(q\) 都能被快速分解，同时可以直接生成 Möbius 求和所需的全部无平方因子除数。
接下来，对每个 \(q\) 分别处理 A、C、D 三个区域。A 区的做法是先数到 \(\lfloor\alpha q\rfloor\) 为止的前缀，再减去其中 \(z_A>N\) 的中间坏区间。C 和 D 区是单调的，因此可以直接做前缀差。无论哪个区域，互素过滤和必要时的模 \(13\) 禁用余数类，都是通过同一组除数与 Möbius 符号来完成的。
C++ 版本还会把外层的 \(q\) 范围切分给多个线程并行累加部分和；Python 和 Java 版本则按相同数学步骤串行执行。所有 \(q\ge 1\) 的贡献处理完之后，如果 \(N\ge 3\)，还要再加上孤立的原始基例 \(z=3\)。
记 \(Q=\lfloor\sqrt N\rfloor\)。建立最小质因子筛需要 \(O(Q\log\log Q)\) 时间和 \(O(Q)\) 空间。对固定的 \(q\)，无平方因子除数的个数是 \(2^{\omega(q)}\)，而每个区域只会对这张表做常数次求和。
因此总时间复杂度为
$$O\left(Q\log\log Q+\sum_{q\le Q}2^{\omega(q)}\right),$$
平均意义下接近 \(O(Q\log Q)\)，实践中对 \(\sqrt N\) 来说几乎是线性的。空间复杂度保持为 \(O(Q)\)，除了筛表之外，每次 \(q\) 只需要很小的临时数组。
Решение использует каноническую параметризацию примитивными целочисленными парами \((p,q)\), где \(q\ge 0\) и \(\gcd(p,q)=1\). В этой записи размер, который нужно считать, задается бинарной квадратичной формой с дискриминантом \(13\)
$$z=\left|3p^2-7pq+3q^2\right|,$$
и требуется посчитать все допустимые примитивные пары с \(z\le N\). При этом еще одна квадратичная координата должна быть положительной, поэтому отношение \(p/q\) распадается лишь на несколько разрешенных областей.
Изолированный случай \(q=0\) дает \(z=3\). Все остальные вклады приходят из \(q\ge 1\), поэтому реализации перебирают слои по \(q\) и считают допустимые интервалы для \(k\), а не исходное пространство напрямую.
Введем две квадратичные формы:
$$Q(p,q)=3p^2-7pq+3q^2,\qquad X(p,q)=p^2-6pq+6q^2.$$
Считаемое значение равно \(|Q(p,q)|\), а допустимость требует \(X(p,q)>0\), \(|p|>q\) и \(\gcd(p,q)=1\). Код переводит все эти условия в интервальные подсчеты по новой переменной \(k\).
Так как \(|p|>q\), каждую допустимую пару при \(q>0\) можно единственным образом представить в одном из двух видов:
$$p=q+k,\qquad p=-q-k.$$
где \(k\ge 1\). При этом взаимная простота сохраняется:
$$\gcd(q,q+k)=\gcd(q,k),\qquad \gcd(q,-q-k)=\gcd(q,k).$$
Следовательно, для фиксированного \(q\) арифметическое условие сводится к \(\gcd(q,k)=1\).
Если \(p=q+k\), то вспомогательная форма принимает вид
$$X(q+k,q)=q^2-4qk+k^2=(k-(2-\sqrt3)q)(k-(2+\sqrt3)q).$$
Обозначив
$$\alpha=2-\sqrt3,\qquad \beta=2+\sqrt3,$$
получаем, что \(X>0\) только при \(k<\alpha q\) либо \(k>\beta q\). Это и есть области A и D:
$$z_A(q,k)=-(Q(q+k,q))=q^2+qk-3k^2,$$
$$z_D(q,k)=Q(q+k,q)=3k^2-qk-q^2.$$
Если же \(p=-q-k\), то
$$X(-q-k,q)=13q^2+8qk+k^2>0,$$
то есть положительность выполняется автоматически, а размер равен
$$z_C(q,k)=Q(-q-k,q)=13q^2+13qk+3k^2.$$
Значит, полный ответ состоит из вкладов областей A, C, D и единственного базового случая \(z=3\).
Для фиксированного \(q\) определим
$$C_q(K)=\#\{1\le k\le K:\gcd(q,k)=1\}.$$
По принципу включения-исключения по простым делителям \(q\) имеем
$$C_q(K)=\sum_{d\mid q}\mu(d)\left\lfloor\frac{K}{d}\right\rfloor.$$
Нужны только квадратсвободные делители, поэтому после факторизации \(q\) реализации строят все \(2^{\omega(q)}\) такие делители вместе со знаками Мёбиуса. Тогда число допустимых \(k\) на любом интервале \([L,U]\) равно
$$C_q(U)-C_q(L-1).$$
Для формы с дискриминантом \(13\) выполняется сравнение
$$Q(p,q)\equiv 3(p+q)^2 \pmod{13}.$$
Отсюда следует
$$13\mid Q(p,q)\iff p\equiv -q \pmod{13}.$$
Реализации оставляют только каноническую ветвь с \(13\nmid Q(p,q)\), поэтому приходится исключать ровно один класс вычетов:
$$p=q+k \Rightarrow k\equiv -2q \pmod{13},$$
$$p=-q-k \Rightarrow k\equiv 0 \pmod{13}.$$
Если \(q\equiv 0\pmod{13}\), то взаимная простота уже запрещает \(k\equiv 0\pmod{13}\), и дополнительное вычитание не требуется. В остальных случаях запрещенный класс считается по той же таблице Мёбиуса, но уже на арифметической прогрессии по модулю \(13\).
Каждая область дает интервал по \(k\).
Область A сначала ограничена как \(1\le k\le \lfloor\alpha q\rfloor\). Полином \(z_A(q,k)=q^2+qk-3k^2\) вогнутый, поэтому неравенство \(z_A\le N\) может нарушаться только на среднем отрезке. Из уравнения \(z_A=N\) получаем корни
$$k=\frac{q\pm\sqrt{13q^2-12N}}{6}.$$
Если \(13q^2\le 12N\), вся область допустима; иначе нужно вычесть целые числа между этими корнями.
Область C монотонно возрастает, значит ее верхняя граница следует из
$$13q^2+13qk+3k^2\le N\quad\Rightarrow\quad k\le \frac{\sqrt{13q^2+12N}-13q}{6}.$$
Область D также монотонна после \(k>\beta q\), поэтому используется
$$k\ge \lfloor\beta q\rfloor+1,\qquad k\le \frac{\sqrt{13q^2+12N}+q}{6}.$$
В реализациях эти границы сначала вычисляются через целочисленный квадратный корень, а затем поправляются несколькими точными проверками, чтобы избежать ошибок округления.
Здесь \(\lfloor\alpha\rfloor=0\), поэтому область A пуста.
Для области C нужно решить
$$13+13k+3k^2\le 100,$$
откуда \(k=1,2,3\), а значит, вклад равен \(29,51,79\).
Область D начинается с \(k=\lfloor\beta\rfloor+1=4\). Там требуется
$$3k^2-k-1\le 100,$$
что дает \(k=4,5\), то есть значения \(43\) и \(69\).
При \(q=1\) запрещенные классы равны \(k\equiv 11\pmod{13}\) в A и D и \(k\equiv 0\pmod{13}\) в C, поэтому на этом слое ничего не вычитается. Вместе с базовым случаем \(q=0\) это хорошо показывает механику подсчета.
Реализации на C++, Python и Java сначала задают \(Q_{\max}=\lfloor\sqrt N\rfloor\) и строят таблицу наименьшего простого делителя до этого предела. Благодаря этому каждое \(q\) быстро факторизуется, а все квадратсвободные делители для сумм Мёбиуса порождаются сразу.
Далее для каждого \(q\) отдельно считаются области A, C и D. Область A обрабатывается как префиксный подсчет до \(\lfloor\alpha q\rfloor\) с последующим вычитанием среднего интервала, где \(z_A>N\). Области C и D монотонны, поэтому считаются через разность префиксов. И фильтр взаимной простоты, и при необходимости запрещенный класс по модулю \(13\) используют одну и ту же таблицу делителей.
В версии на C++ внешний диапазон по \(q\) дополнительно делится между потоками, и частичные суммы складываются параллельно. Версии на Python и Java выполняют ту же арифметику последовательно. После завершения всех вкладов с \(q\ge 1\) единственный базовый примитивный случай \(z=3\) добавляется один раз при \(N\ge 3\).
Пусть \(Q=\lfloor\sqrt N\rfloor\). Построение SPF-решета занимает \(O(Q\log\log Q)\) времени и \(O(Q)\) памяти. Для фиксированного \(q\) число квадратсвободных делителей равно \(2^{\omega(q)}\), и каждая область использует лишь константное число сумм по этой таблице.
Итоговая временная сложность равна
$$O\left(Q\log\log Q+\sum_{q\le Q}2^{\omega(q)}\right),$$
что в среднем близко к \(O(Q\log Q)\) и на практике ведет себя почти линейно по \(\sqrt N\). Память остается \(O(Q)\), помимо решета нужны только небольшие временные массивы.
تعتمد الفكرة على تمثيلٍ قانوني للحلول بواسطة أزواج صحيحة أولية نسبيًا \((p,q)\) بحيث \(q\ge 0\) و \(\gcd(p,q)=1\). في هذا التمثيل تكون الكمية التي نعدها معطاةً بالشكل التربيعي الثنائي ذي المميز \(13\):
$$z=\left|3p^2-7pq+3q^2\right|.$$
والمطلوب هو عدّ جميع الأزواج البدائية المقبولة التي تحقق \(z\le N\). لكن توجد أيضًا صيغة تربيعية ثانية يجب أن تبقى موجبة، ولهذا لا تكون النسبة \(p/q\) حرة بالكامل، بل تنقسم إلى ثلاث مناطق فعالة فقط.
الحالة المنفردة \(q=0\) تعطي \(z=3\). وكل مساهمة أخرى تأتي من \(q\ge 1\)، ولذلك تدور التطبيقات على قيم \(q\) ثم تعدّ فترات \(k\) الصحيحة الصالحة بدلًا من فحص فضاء البحث الأصلي مباشرة.
لنعرّف الشكلين التربيعيين
$$Q(p,q)=3p^2-7pq+3q^2,\qquad X(p,q)=p^2-6pq+6q^2.$$
القيمة التي تدخل في العد هي \(|Q(p,q)|\)، بينما القبول يتطلب \(X(p,q)>0\) و \(|p|>q\) و \(\gcd(p,q)=1\). يقوم الكود بتحويل هذه الشروط إلى عدٍّ على فترات في متغير جديد هو \(k\).
بما أن \(|p|>q\)، فإن كل زوج مقبول مع \(q>0\) يمكن كتابته بصورة وحيدة على إحدى الصورتين
$$p=q+k,\qquad p=-q-k.$$
حيث \(k\ge 1\). كما أن شرط الأولية النسبية يبقى كما هو:
$$\gcd(q,q+k)=\gcd(q,k),\qquad \gcd(q,-q-k)=\gcd(q,k).$$
إذًا، بعد تثبيت \(q\)، يصبح الشرط الحسابي الوحيد هو \(\gcd(q,k)=1\).
إذا كان \(p=q+k\)، فإن الصيغة المساعدة تصبح
$$X(q+k,q)=q^2-4qk+k^2=(k-(2-\sqrt3)q)(k-(2+\sqrt3)q).$$
وبوضع
$$\alpha=2-\sqrt3,\qquad \beta=2+\sqrt3,$$
فإن الشرط \(X>0\) يفرض إما \(k<\alpha q\) أو \(k>\beta q\). ومن هنا تظهر المنطقتان A و D:
$$z_A(q,k)=-(Q(q+k,q))=q^2+qk-3k^2,$$
$$z_D(q,k)=Q(q+k,q)=3k^2-qk-q^2.$$
أما إذا كان \(p=-q-k\)، فنحصل على
$$X(-q-k,q)=13q^2+8qk+k^2>0,$$
أي إن شرط الإيجابية يتحقق تلقائيًا، والقيمة المعدودة تصبح
$$z_C(q,k)=Q(-q-k,q)=13q^2+13qk+3k^2.$$
إذن الجواب الكلي هو مجموع مساهمات المناطق A و C و D، بالإضافة إلى الحالة الأساسية الوحيدة \(z=3\).
لـ \(q\) ثابت، نعرّف
$$C_q(K)=\#\{1\le k\le K:\gcd(q,k)=1\}.$$
وباستخدام مبدأ الشمول والاستبعاد على القواسم الأولية لـ \(q\)، نحصل على
$$C_q(K)=\sum_{d\mid q}\mu(d)\left\lfloor\frac{K}{d}\right\rfloor.$$
القواسم المربعة الحرة فقط هي التي تؤثر، لذلك بعد تحليل \(q\) إلى عوامله الأولية تولد التطبيقات كل القواسم المربعة الحرة وعددها \(2^{\omega(q)}\) مع إشارات Möbius الخاصة بها. وبعد ذلك يمكن عدّ أي فترة \([L,U]\) بالصيغة
$$C_q(U)-C_q(L-1).$$
الشكل ذو المميز \(13\) يحقق العلاقة
$$Q(p,q)\equiv 3(p+q)^2 \pmod{13}.$$
ومن ثم
$$13\mid Q(p,q)\iff p\equiv -q \pmod{13}.$$
التطبيقات تبقي فقط الفرع القانوني الذي لا يقبل القسمة على \(13\)، ولذلك يجب حذف فئة باقية واحدة تمامًا:
$$p=q+k \Rightarrow k\equiv -2q \pmod{13},$$
$$p=-q-k \Rightarrow k\equiv 0 \pmod{13}.$$
إذا كان \(q\equiv 0\pmod{13}\)، فإن شرط الأولية النسبية يمنع أصلًا \(k\equiv 0\pmod{13}\)، فلا نحتاج إلى طرح إضافي. أما في غير ذلك، فتُعدّ الفئة الممنوعة باستخدام جدول Möbius نفسه ولكن على متتالية حسابية modulo \(13\).
كل منطقة تعطي فترة في \(k\).
المنطقة A تبدأ من \(1\le k\le \lfloor\alpha q\rfloor\). وبما أن \(z_A(q,k)=q^2+qk-3k^2\) قطع مكافئ مفتوح إلى الأسفل، فإن مخالفة الشرط \(z_A\le N\) لا يمكن أن تقع إلا في فترة وسطى. ومن حل المعادلة \(z_A=N\) نحصل على
$$k=\frac{q\pm\sqrt{13q^2-12N}}{6}.$$
إذا كان \(13q^2\le 12N\)، فكل المنطقة تبقى صالحة، وإلا يجب حذف الأعداد الصحيحة الواقعة بين الجذرين.
المنطقة C متزايدة رتيبًا، لذا يأتي حدها الأعلى من
$$13q^2+13qk+3k^2\le N\quad\Rightarrow\quad k\le \frac{\sqrt{13q^2+12N}-13q}{6}.$$
والمنطقة D تصبح متزايدة أيضًا بعد \(k>\beta q\)، ولذلك تستخدم
$$k\ge \lfloor\beta q\rfloor+1,\qquad k\le \frac{\sqrt{13q^2+12N}+q}{6}.$$
في التنفيذ تُحسب هذه الحدود أولًا بواسطة الجذر التربيعي الصحيح، ثم تُصحح بعدد قليل من الاختبارات الدقيقة حتى لا تبقى أي أخطاء تقريبية.
هنا لدينا \(\lfloor\alpha\rfloor=0\)، ولذلك تكون المنطقة A فارغة.
في المنطقة C نحل
$$13+13k+3k^2\le 100,$$
فنحصل على \(k=1,2,3\)، أي القيم \(29,51,79\).
المنطقة D تبدأ من \(k=\lfloor\beta\rfloor+1=4\). وهناك نحتاج إلى
$$3k^2-k-1\le 100,$$
ومن ثم \(k=4,5\)، فتكون القيم \(43\) و \(69\).
عندما \(q=1\)، تكون الفئات الممنوعة هي \(k\equiv 11\pmod{13}\) في A و D، و \(k\equiv 0\pmod{13}\) في C، ولذلك لا يحدث أي طرح في هذه الطبقة. ومع إضافة الحالة الأساسية \(q=0\)، يصبح مسار العد واضحًا بالكامل.
تبدأ تطبيقات C++ و Python و Java بحساب \(Q_{\max}=\lfloor\sqrt N\rfloor\) ثم بناء جدول لأصغر عامل أولي حتى هذا الحد. وبهذا يصبح تحليل كل \(q\) إلى عوامله سريعًا، ويمكن توليد جميع القواسم المربعة الحرة اللازمة في مجاميع Möbius مباشرة.
بعد ذلك تعالج كل قيمة \(q\) المناطق A و C و D بصورة مستقلة. في المنطقة A يتم عدّ كل ما يصل إلى \(\lfloor\alpha q\rfloor\) ثم طرح الفترة الوسطى التي تحقق \(z_A>N\). أما المنطقتان C و D فهما رتيبتان، ولذلك تُحسبان بفروق المجاميع التراكمية. وفي جميع الحالات يُطبّق شرط الأولية النسبية، ومعه عند الحاجة حذف فئة البواقي المحظورة modulo \(13\)، باستخدام جدول القواسم نفسه.
تنفذ نسخة C++ هذه الطبقة الخارجية على عدة خيوط وتجمع المجاميع الجزئية بالتوازي. أما نسختا Python و Java فتطبقان الحساب نفسه بصورة متسلسلة. وبعد إنهاء جميع مساهمات \(q\ge 1\)، تُضاف الحالة البدائية المعزولة \(z=3\) مرة واحدة إذا كان \(N\ge 3\).
لنكتب \(Q=\lfloor\sqrt N\rfloor\). بناء غربال أصغر عامل أولي يكلف \(O(Q\log\log Q)\) زمنًا و \(O(Q)\) ذاكرة. ولـ \(q\) ثابت، فإن عدد القواسم المربعة الحرة هو \(2^{\omega(q)}\)، وكل منطقة تستخدم عددًا ثابتًا فقط من المجاميع فوق هذه القائمة.
إذًا يكون الزمن الكلي
$$O\left(Q\log\log Q+\sum_{q\le Q}2^{\omega(q)}\right),$$
وهو قريب من \(O(Q\log Q)\) في المتوسط، ويتصرف عمليًا بشكل شبه خطي بالنسبة إلى \(\sqrt N\). أما الذاكرة فتبقى \(O(Q)\)، ولا توجد سوى مصفوفات مؤقتة صغيرة إضافة إلى الغربال.