We must count integer pairs \((a,b)\neq(0,0)\) such that
$$a^2+ab+41b^2\le N.$$
This binary quadratic form is positive definite, so all valid pairs lie in a bounded region. A naive two-dimensional scan would still be far too slow for \(N=10^{16}\). The key observation is that a simple linear change of variables turns the form into an axis-aligned ellipse with a parity condition, which makes the counting almost one-dimensional.
Let \(F(N)\) denote the number of nonzero integer pairs \((a,b)\) satisfying \(a^2+ab+41b^2\le N\). The implementations count \(F(N)\) by rewriting the quadratic form and then sweeping over one coordinate.
Introduce the substitution
$$x=2a+b,\qquad y=b.$$
Then
$$4(a^2+ab+41b^2)=(2a+b)^2+163b^2=x^2+163y^2.$$
So the original inequality is equivalent to
$$x^2+163y^2\le 4N.$$
The coefficient \(163\) is the absolute value of the discriminant \(1-4\cdot 41=-163\), so its appearance is natural rather than accidental.
The substitution is reversible:
$$a=\frac{x-y}{2},\qquad b=y.$$
Therefore \((a,b)\) is an integer pair exactly when \(x-y\) is even, which is the same as
$$x\equiv y\pmod{2}.$$
This gives the equivalent counting formula
$$F(N)=\#\left\{(x,y)\in\mathbb{Z}^2:x^2+163y^2\le 4N,\ x\equiv y\pmod{2}\right\}-1.$$
The final \(-1\) removes the origin, which is the unique lattice point producing the value \(0\).
For any admissible point,
$$163y^2\le 4N.$$
Hence
$$|y|\le Y=\left\lfloor\sqrt{\frac{4N}{163}}\right\rfloor.$$
So there are only \(2Y+1\) possible values of \(y\). This collapses the search from a two-variable problem to a one-variable sweep over horizontal slices of the ellipse.
Fix one integer \(y\) with \(|y|\le Y\), and define
$$R(y)=4N-163y^2,\qquad M(y)=\left\lfloor\sqrt{R(y)}\right\rfloor.$$
Then \(x^2\le R(y)\) is equivalent to
$$-M(y)\le x\le M(y).$$
Without the parity restriction, this interval contains \(2M(y)+1\) integers. We only want those whose parity matches \(y\).
The interval \([-M,M]\) is symmetric around \(0\). If \(M\) is even, it contains \(M+1\) even integers and \(M\) odd integers. If \(M\) is odd, the roles reverse. Therefore the number of valid \(x\)-values for a fixed \(y\) is
$$C(y)=\begin{cases} M(y)+1,& y\equiv M(y)\pmod{2},\\ M(y),& y\not\equiv M(y)\pmod{2}. \end{cases}$$
This is the crucial simplification: one integer square root and one parity test replace a full inner enumeration of \(x\).
Adding the contribution of every permitted \(y\) gives
$$F(N)=\sum_{y=-Y}^{Y} C(y)-1.$$
Because \(M(y)=M(-y)\), the picture is symmetric. The implementations still loop over the full range directly, which keeps the logic simple and avoids any special handling at \(y=0\).
For \(N=1000\),
$$Y=\left\lfloor\sqrt{\frac{4000}{163}}\right\rfloor=4.$$
Now evaluate the slices:
$$\begin{aligned} y=0 &:\quad R=4000,\quad M=63,\quad C=63,\\ y=\pm 1 &:\quad R=3837,\quad M=61,\quad C=62,\\ y=\pm 2 &:\quad R=3348,\quad M=57,\quad C=57,\\ y=\pm 3 &:\quad R=2533,\quad M=50,\quad C=50,\\ y=\pm 4 &:\quad R=1392,\quad M=37,\quad C=37. \end{aligned}$$
Therefore
$$F(1000)=63+2(62+57+50+37)-1=474.$$
This matches the count returned by the implementation and confirms the transformation and parity formula.
The C++, Python, and Java implementations follow the same mathematical plan. First they compute an exact floor square root. Python uses the standard exact integer routine, while the C++ and Java versions start from a floating approximation and then correct it upward or downward until the exact integer floor is reached. That correction step matters because the input is large enough that a one-unit error near a square boundary would change the final count.
After that, the implementation finds the largest permitted absolute value of the outer coordinate, iterates through every integer in that range, forms the remaining budget \(4N-163y^2\), takes its integer square root, and applies the parity formula above to count the matching values of the inner coordinate. Each loop iteration contributes one horizontal slice of the ellipse, and the origin is removed once at the end.
The loop length is
$$2\left\lfloor\sqrt{\frac{4N}{163}}\right\rfloor+1=\Theta(\sqrt N).$$
Each iteration performs constant-time arithmetic plus one integer square root, so the total running time is \(O(\sqrt N)\). Only a handful of integer variables are stored, so the memory usage is \(O(1)\).
Gesucht ist die Anzahl der ganzzahligen Paare \((a,b)\neq(0,0)\) mit
$$a^2+ab+41b^2\le N.$$
Diese binäre quadratische Form ist positiv definit, daher liegen alle zulässigen Punkte in einem beschränkten Bereich. Eine naive Suche über beide Koordinaten wäre für \(N=10^{16}\) dennoch viel zu langsam. Der entscheidende Schritt ist eine lineare Variablentransformation, die die Form in eine achsenparallele Ellipse mit zusätzlicher Paritätsbedingung umschreibt.
Sei \(F(N)\) die Anzahl der von null verschiedenen ganzzahligen Paare \((a,b)\) mit \(a^2+ab+41b^2\le N\). Die Implementierungen berechnen \(F(N)\), indem sie die Form umschreiben und dann nur noch über eine Koordinate iterieren.
Wir setzen
$$x=2a+b,\qquad y=b.$$
Dann gilt
$$4(a^2+ab+41b^2)=(2a+b)^2+163b^2=x^2+163y^2.$$
Die ursprüngliche Ungleichung ist also äquivalent zu
$$x^2+163y^2\le 4N.$$
Die Zahl \(163\) ist dabei nicht zufällig: Sie ist der Betrag der Diskriminante \(1-4\cdot 41=-163\).
Die Transformation ist umkehrbar:
$$a=\frac{x-y}{2},\qquad b=y.$$
Damit ist \((a,b)\) genau dann ganzzahlig, wenn \(x-y\) gerade ist, also
$$x\equiv y\pmod{2}.$$
Somit erhalten wir die äquivalente Zählformel
$$F(N)=\#\left\{(x,y)\in\mathbb{Z}^2:x^2+163y^2\le 4N,\ x\equiv y\pmod{2}\right\}-1.$$
Das abschließende \(-1\) entfernt den Ursprung, also den einzigen Gitterpunkt mit Wert \(0\).
Für jeden zulässigen Punkt muss gelten
$$163y^2\le 4N.$$
Daraus folgt
$$|y|\le Y=\left\lfloor\sqrt{\frac{4N}{163}}\right\rfloor.$$
Es sind also nur \(2Y+1\) Werte für \(y\) möglich. Genau dadurch wird aus einem zweidimensionalen Suchproblem ein eindimensionaler Sweep über horizontale Schnitte der Ellipse.
Fixiere ein \(y\) mit \(|y|\le Y\) und definiere
$$R(y)=4N-163y^2,\qquad M(y)=\left\lfloor\sqrt{R(y)}\right\rfloor.$$
Dann ist \(x^2\le R(y)\) äquivalent zu
$$-M(y)\le x\le M(y).$$
Ohne Paritätsbedingung enthält dieses Intervall \(2M(y)+1\) ganze Zahlen. Gesucht sind aber nur diejenigen mit derselben Parität wie \(y\).
Das Intervall \([-M,M]\) ist um \(0\) symmetrisch. Ist \(M\) gerade, enthält es \(M+1\) gerade und \(M\) ungerade Zahlen. Ist \(M\) ungerade, vertauschen sich diese Rollen. Deshalb ist die Anzahl zulässiger \(x\)-Werte für festes \(y\)
$$C(y)=\begin{cases} M(y)+1,& y\equiv M(y)\pmod{2},\\ M(y),& y\not\equiv M(y)\pmod{2}. \end{cases}$$
Das ist die zentrale Vereinfachung: Eine einzige ganzzahlige Quadratwurzel und ein Paritätstest ersetzen die vollständige innere Enumeration über \(x\).
Summiert man die Beiträge aller zulässigen \(y\), erhält man
$$F(N)=\sum_{y=-Y}^{Y} C(y)-1.$$
Da \(M(y)=M(-y)\) gilt, ist die Figur symmetrisch. Die Implementierungen laufen trotzdem direkt über den gesamten Bereich, weil das den Code einfacher hält und keine Sonderbehandlung für \(y=0\) braucht.
Für \(N=1000\) gilt
$$Y=\left\lfloor\sqrt{\frac{4000}{163}}\right\rfloor=4.$$
Die einzelnen Schnitte sind dann
$$\begin{aligned} y=0 &:\quad R=4000,\quad M=63,\quad C=63,\\ y=\pm 1 &:\quad R=3837,\quad M=61,\quad C=62,\\ y=\pm 2 &:\quad R=3348,\quad M=57,\quad C=57,\\ y=\pm 3 &:\quad R=2533,\quad M=50,\quad C=50,\\ y=\pm 4 &:\quad R=1392,\quad M=37,\quad C=37. \end{aligned}$$
Also
$$F(1000)=63+2(62+57+50+37)-1=474.$$
Genau diesen Wert liefert auch die Implementierung; damit sind Transformation und Paritätsformel überprüft.
Die C++-, Python- und Java-Implementierungen folgen demselben mathematischen Ablauf. Zuerst wird jeweils eine exakte ganzzahlige Quadratwurzel berechnet. Python verwendet dafür die exakte Standardbibliothek, während C++ und Java von einer Gleitkomma-Näherung ausgehen und diese anschließend nach oben oder unten korrigieren, bis der exakte ganzzahlige Bodenwert erreicht ist. Dieser Korrekturschritt ist wichtig, weil bei so großen Eingaben schon ein Fehler von \(1\) an einer Quadratgrenze die Gesamtanzahl verfälschen würde.
Anschließend bestimmt die Implementierung den größten möglichen Betrag der äußeren Koordinate, iteriert über alle ganzzahligen Werte in diesem Bereich, bildet das Restbudget \(4N-163y^2\), zieht daraus die ganzzahlige Quadratwurzel und wendet die obige Paritätsformel an. Jeder Schleifendurchlauf zählt genau einen horizontalen Schnitt der Ellipse; am Ende wird der Ursprung einmal entfernt.
Die Schleifenlänge ist
$$2\left\lfloor\sqrt{\frac{4N}{163}}\right\rfloor+1=\Theta(\sqrt N).$$
Pro Iteration fallen konstante Arithmetik und eine ganzzahlige Quadratwurzel an. Damit beträgt die Gesamtlaufzeit \(O(\sqrt N)\). Gespeichert werden nur wenige Ganzzahlen, also ist der Speicherbedarf \(O(1)\).
Sayılması gereken şey,
$$a^2+ab+41b^2\le N$$
eşitsizliğini sağlayan ve \((a,b)\neq(0,0)\) olan tamsayı çiftleridir. Bu ikili kuadratik form pozitif tanımlıdır; dolayısıyla geçerli noktalar sınırlı bir bölgede kalır. Yine de \(N=10^{16}\) için iki değişken üzerinde doğrudan tarama çok yavaştır. Çözüm, formu eksenlere hizalı bir elipse ve ek bir parite koşuluna dönüştüren lineer değişken dönüşümünden yararlanır.
\(F(N)\), \(a^2+ab+41b^2\le N\) koşulunu sağlayan sıfırdan farklı tamsayı çiftlerinin sayısı olsun. Uygulamalar \(F(N)\)'i, kuadratik formu yeniden yazarak ve sonra tek bir koordinat üzerinde gezerek hesaplar.
Şu dönüşümü yapalım:
$$x=2a+b,\qquad y=b.$$
O zaman
$$4(a^2+ab+41b^2)=(2a+b)^2+163b^2=x^2+163y^2.$$
Böylece başlangıçtaki eşitsizlik
$$x^2+163y^2\le 4N$$
şeklini alır. Buradaki \(163\) sayısı tesadüf değildir; diskriminant \(1-4\cdot 41=-163\) olduğundan doğal olarak ortaya çıkar.
Dönüşüm tersinir:
$$a=\frac{x-y}{2},\qquad b=y.$$
Demek ki \((a,b)\) ancak ve ancak \(x-y\) çift olduğunda tam sayıdır; bu da
$$x\equiv y\pmod{2}$$
ile aynıdır. Böylece eşdeğer sayım formülü
$$F(N)=\#\left\{(x,y)\in\mathbb{Z}^2:x^2+163y^2\le 4N,\ x\equiv y\pmod{2}\right\}-1$$
olur. Sondaki \(-1\), değeri \(0\) yapan tek nokta olan orijini çıkarır.
Geçerli her nokta için
$$163y^2\le 4N$$
olmalıdır. Buradan
$$|y|\le Y=\left\lfloor\sqrt{\frac{4N}{163}}\right\rfloor$$
elde edilir. Yani yalnızca \(2Y+1\) farklı \(y\) değeri incelenir. Böylece iki değişkenli arama, elipsin yatay kesitleri üzerinde tek boyutlu bir dolaşıma indirgenir.
\(|y|\le Y\) olacak biçimde bir \(y\) sabitleyelim ve
$$R(y)=4N-163y^2,\qquad M(y)=\left\lfloor\sqrt{R(y)}\right\rfloor$$
tanımlarını yapalım. O zaman \(x^2\le R(y)\) koşulu
$$-M(y)\le x\le M(y)$$
aralığına eşdeğerdir. Parite koşulu olmasa bu aralıkta \(2M(y)+1\) tamsayı olurdu. Biz sadece \(y\) ile aynı pariteye sahip olanları istiyoruz.
\([-M,M]\) aralığı \(0\) etrafında simetriktir. \(M\) çiftse \(M+1\) çift ve \(M\) tek sayı vardır. \(M\) tekse bu roller yer değiştirir. Dolayısıyla sabit bir \(y\) için geçerli \(x\) sayısı
$$C(y)=\begin{cases} M(y)+1,& y\equiv M(y)\pmod{2},\\ M(y),& y\not\equiv M(y)\pmod{2}. \end{cases}$$
Kritik sadeleşme budur: İç döngüyle tüm \(x\) değerlerini denemek yerine tek bir tamsayı karekök ve bir parite testi yeterlidir.
Tüm uygun \(y\) değerlerinin katkıları toplandığında
$$F(N)=\sum_{y=-Y}^{Y} C(y)-1$$
elde edilir. \(M(y)=M(-y)\) olduğu için şekil simetriktir; yine de uygulamalar tüm aralığı doğrudan dolaşır. Bu yaklaşım kodu kısa tutar ve \(y=0\) için özel durum yazmayı gerektirmez.
\(N=1000\) için
$$Y=\left\lfloor\sqrt{\frac{4000}{163}}\right\rfloor=4$$
olur. Kesitler tek tek şöyledir:
$$\begin{aligned} y=0 &:\quad R=4000,\quad M=63,\quad C=63,\\ y=\pm 1 &:\quad R=3837,\quad M=61,\quad C=62,\\ y=\pm 2 &:\quad R=3348,\quad M=57,\quad C=57,\\ y=\pm 3 &:\quad R=2533,\quad M=50,\quad C=50,\\ y=\pm 4 &:\quad R=1392,\quad M=37,\quad C=37. \end{aligned}$$
Dolayısıyla
$$F(1000)=63+2(62+57+50+37)-1=474.$$
Bu değer, uygulamanın ürettiği küçük kontrol sonucuyla aynıdır ve dönüşümün doğru kurulduğunu gösterir.
C++, Python ve Java uygulamaları aynı matematiksel planı izler. Önce tam olarak \(\lfloor\sqrt{n}\rfloor\) değerini hesaplarlar. Python bunu standart kütüphanedeki kesin yordamla yapar; C++ ve Java ise kayan noktalı bir başlangıç tahmini alıp bunu yukarı veya aşağı düzelterek tam tamsayı sonuca ulaşır. Bu düzeltme önemlidir, çünkü giriş çok büyük olduğunda kare sınırlarında oluşacak tek bir hata bile toplam sayımı bozabilir.
Daha sonra uygulama dış koordinatın alabileceği en büyük mutlak değeri bulur, bu aralıktaki tüm tam sayılar üzerinde dolaşır, kalan bütçe \(4N-163y^2\)'yi hesaplar, bunun tamsayı karekökünü alır ve yukarıdaki parite formülünü uygular. Her yineleme elipsin tek bir yatay kesitine karşılık gelir; döngü bitince orijin bir kez çıkarılır.
Döngü uzunluğu
$$2\left\lfloor\sqrt{\frac{4N}{163}}\right\rfloor+1=\Theta(\sqrt N)$$
kadardır. Her adımda sabit sayıda aritmetik işlem ve bir tamsayı karekök vardır. Bu yüzden toplam süre \(O(\sqrt N)\), kullanılan bellek ise \(O(1)\) olur.
Debemos contar los pares enteros \((a,b)\neq(0,0)\) tales que
$$a^2+ab+41b^2\le N.$$
Esta forma cuadrática binaria es definida positiva, de modo que todos los pares válidos quedan dentro de una región acotada. Aun así, un barrido ingenuo en dos dimensiones sería demasiado costoso para \(N=10^{16}\). La idea decisiva es aplicar un cambio lineal de variables que transforma la forma en una elipse alineada con los ejes, junto con una condición de paridad.
Sea \(F(N)\) el número de pares enteros no nulos \((a,b)\) que satisfacen \(a^2+ab+41b^2\le N\). Las implementaciones calculan \(F(N)\) reescribiendo la forma y recorriendo luego una sola coordenada.
Introducimos el cambio
$$x=2a+b,\qquad y=b.$$
Entonces
$$4(a^2+ab+41b^2)=(2a+b)^2+163b^2=x^2+163y^2.$$
Por tanto, la desigualdad original equivale a
$$x^2+163y^2\le 4N.$$
El número \(163\) aparece de forma natural porque es el valor absoluto del discriminante \(1-4\cdot 41=-163\).
La transformación es invertible:
$$a=\frac{x-y}{2},\qquad b=y.$$
Así, \((a,b)\) es entero exactamente cuando \(x-y\) es par, es decir, cuando
$$x\equiv y\pmod{2}.$$
De este modo obtenemos la fórmula equivalente
$$F(N)=\#\left\{(x,y)\in\mathbb{Z}^2:x^2+163y^2\le 4N,\ x\equiv y\pmod{2}\right\}-1.$$
El \(-1\) final elimina el origen, que es el único punto con valor cuadrático igual a \(0\).
Todo punto admisible debe cumplir
$$163y^2\le 4N.$$
Luego
$$|y|\le Y=\left\lfloor\sqrt{\frac{4N}{163}}\right\rfloor.$$
Por lo tanto solo hay \(2Y+1\) valores posibles de \(y\). Esto reduce el problema de una búsqueda bidimensional a un barrido unidimensional sobre secciones horizontales de la elipse.
Fijemos un entero \(y\) con \(|y|\le Y\) y definamos
$$R(y)=4N-163y^2,\qquad M(y)=\left\lfloor\sqrt{R(y)}\right\rfloor.$$
Entonces \(x^2\le R(y)\) equivale a
$$-M(y)\le x\le M(y).$$
Sin la restricción de paridad, ese intervalo contiene \(2M(y)+1\) enteros. Solo interesan los que tienen la misma paridad que \(y\).
El intervalo \([-M,M]\) es simétrico respecto de \(0\). Si \(M\) es par, contiene \(M+1\) enteros pares y \(M\) impares. Si \(M\) es impar, sucede al revés. Por tanto, el número de valores válidos de \(x\) para un \(y\) fijo es
$$C(y)=\begin{cases} M(y)+1,& y\equiv M(y)\pmod{2},\\ M(y),& y\not\equiv M(y)\pmod{2}. \end{cases}$$
Aquí está la simplificación clave: una sola raíz cuadrada entera y una comprobación de paridad sustituyen a un bucle interno completo sobre \(x\).
Al sumar la contribución de todos los \(y\) permitidos obtenemos
$$F(N)=\sum_{y=-Y}^{Y} C(y)-1.$$
Como \(M(y)=M(-y)\), la figura es simétrica. Aun así, las implementaciones recorren directamente todo el intervalo, lo cual mantiene el código sencillo y evita casos especiales en \(y=0\).
Para \(N=1000\),
$$Y=\left\lfloor\sqrt{\frac{4000}{163}}\right\rfloor=4.$$
Las secciones quedan así:
$$\begin{aligned} y=0 &:\quad R=4000,\quad M=63,\quad C=63,\\ y=\pm 1 &:\quad R=3837,\quad M=61,\quad C=62,\\ y=\pm 2 &:\quad R=3348,\quad M=57,\quad C=57,\\ y=\pm 3 &:\quad R=2533,\quad M=50,\quad C=50,\\ y=\pm 4 &:\quad R=1392,\quad M=37,\quad C=37. \end{aligned}$$
Por consiguiente,
$$F(1000)=63+2(62+57+50+37)-1=474.$$
Ese resultado coincide con el valor que produce la implementación y valida tanto la transformación como la fórmula de paridad.
Las implementaciones en C++, Python y Java siguen exactamente el mismo esquema matemático. Primero calculan una raíz cuadrada entera exacta. Python utiliza la rutina exacta de la biblioteca estándar; C++ y Java parten de una aproximación en coma flotante y luego la corrigen hacia arriba o hacia abajo hasta alcanzar el piso entero exacto. Esa corrección es importante porque, con entradas tan grandes, un error de una sola unidad cerca de un cuadrado perfecto alteraría la cuenta final.
Después, la implementación obtiene el máximo valor absoluto posible de la coordenada exterior, recorre todos los enteros de ese rango, forma el presupuesto restante \(4N-163y^2\), toma su raíz cuadrada entera y aplica la fórmula de paridad anterior. Cada iteración aporta una sección horizontal de la elipse, y al final se resta el origen una sola vez.
La longitud del bucle es
$$2\left\lfloor\sqrt{\frac{4N}{163}}\right\rfloor+1=\Theta(\sqrt N).$$
Cada iteración realiza aritmética de tiempo constante y una raíz cuadrada entera, así que el tiempo total es \(O(\sqrt N)\). Solo se almacenan unas pocas variables enteras, por lo que la memoria es \(O(1)\).
我们要统计满足
$$a^2+ab+41b^2\le N$$
且 \((a,b)\neq(0,0)\) 的所有整数对。这个二元二次型是正定的，所以可行点只会落在一个有界区域内。但当 \(N=10^{16}\) 时，直接对两个变量做双重枚举仍然太慢。实现采用的核心技巧是做一次线性变量替换，把原式改写成一个与坐标轴对齐的椭圆不等式，同时附带一个同奇偶条件。
记 \(F(N)\) 为所有满足 \(a^2+ab+41b^2\le N\) 的非零整数对 \((a,b)\) 的数量。实现的思路是先重写这个二次型，再只对一个坐标做扫描。
作变量替换
$$x=2a+b,\qquad y=b.$$
于是有
$$4(a^2+ab+41b^2)=(2a+b)^2+163b^2=x^2+163y^2.$$
因此原来的不等式等价于
$$x^2+163y^2\le 4N.$$
这里出现的 \(163\) 并不是偶然，它正好是判别式 \(1-4\cdot 41=-163\) 的绝对值。
这个变换可以反解：
$$a=\frac{x-y}{2},\qquad b=y.$$
所以 \((a,b)\) 为整数对，当且仅当 \(x-y\) 为偶数，也就是
$$x\equiv y\pmod{2}.$$
于是计数问题完全等价于
$$F(N)=\#\left\{(x,y)\in\mathbb{Z}^2:x^2+163y^2\le 4N,\ x\equiv y\pmod{2}\right\}-1.$$
最后的 \(-1\) 用来去掉原点，因为只有原点会让二次型取值为 \(0\)。
任意可行点都必须满足
$$163y^2\le 4N.$$
因此
$$|y|\le Y=\left\lfloor\sqrt{\frac{4N}{163}}\right\rfloor.$$
也就是说，外层只需要考察 \(2Y+1\) 个整数 \(y\)。这一步把原本的二维搜索降成了沿椭圆水平切片的一维扫描。
固定一个满足 \(|y|\le Y\) 的整数 \(y\)，定义
$$R(y)=4N-163y^2,\qquad M(y)=\left\lfloor\sqrt{R(y)}\right\rfloor.$$
那么条件 \(x^2\le R(y)\) 等价于
$$-M(y)\le x\le M(y).$$
如果暂时忽略奇偶限制，这个区间一共有 \(2M(y)+1\) 个整数。真正需要保留的，是与 \(y\) 同奇偶的那些 \(x\)。
区间 \([-M,M]\) 关于 \(0\) 对称。如果 \(M\) 是偶数，那么其中有 \(M+1\) 个偶数和 \(M\) 个奇数；如果 \(M\) 是奇数，情况正好反过来。所以固定 \(y\) 时，可行的 \(x\) 个数为
$$C(y)=\begin{cases} M(y)+1,& y\equiv M(y)\pmod{2},\\ M(y),& y\not\equiv M(y)\pmod{2}. \end{cases}$$
这就是算法的关键简化：不再需要对区间中的所有 \(x\) 逐个测试，只需一次整数平方根和一次奇偶判断即可。
把每个合法 \(y\) 对应的切片贡献加起来，就得到
$$F(N)=\sum_{y=-Y}^{Y} C(y)-1.$$
由于 \(M(y)=M(-y)\)，几何图形本身是对称的。不过实现仍然直接遍历整个区间，这样逻辑更直接，也不需要对 \(y=0\) 单独写分支。
当 \(N=1000\) 时，
$$Y=\left\lfloor\sqrt{\frac{4000}{163}}\right\rfloor=4.$$
逐层计算可得
$$\begin{aligned} y=0 &:\quad R=4000,\quad M=63,\quad C=63,\\ y=\pm 1 &:\quad R=3837,\quad M=61,\quad C=62,\\ y=\pm 2 &:\quad R=3348,\quad M=57,\quad C=57,\\ y=\pm 3 &:\quad R=2533,\quad M=50,\quad C=50,\\ y=\pm 4 &:\quad R=1392,\quad M=37,\quad C=37. \end{aligned}$$
因此
$$F(1000)=63+2(62+57+50+37)-1=474.$$
这个结果与实现给出的校验值一致，说明变换和计数公式都是正确的。
C++、Python 和 Java 实现遵循完全相同的数学流程。第一步都是计算精确的整数平方根。Python 直接使用标准库中的精确整数例程；C++ 和 Java 则先取一个浮点近似值，再向上或向下微调，直到得到真正的整数下取整结果。由于输入规模非常大，如果在接近平方数边界的地方错一位，最终答案就会被破坏，所以这一步必须做得非常严谨。
随后，实现先求出外层坐标可能达到的最大绝对值，再遍历这个范围内的每一个整数，计算剩余预算 \(4N-163y^2\)，对它取整数平方根，然后套用上面的奇偶闭式公式。每次循环都精确对应椭圆上的一条水平切片，全部累加之后再减去原点一次即可。
循环次数为
$$2\left\lfloor\sqrt{\frac{4N}{163}}\right\rfloor+1=\Theta(\sqrt N).$$
每次迭代只做常数次整数运算和一次整数平方根，所以总时间复杂度是 \(O(\sqrt N)\)。算法只保存少量整数变量，因此空间复杂度是 \(O(1)\)。
Нужно подсчитать количество целочисленных пар \((a,b)\neq(0,0)\), для которых
$$a^2+ab+41b^2\le N.$$
Эта бинарная квадратичная форма положительно определена, поэтому все допустимые пары лежат в ограниченной области. Однако для \(N=10^{16}\) прямой перебор по двум координатам всё равно слишком дорог. Ключевой шаг решения состоит в линейной замене переменных, которая превращает форму в эллипс, выровненный по осям, с дополнительным условием чётности.
Обозначим через \(F(N)\) число ненулевых целочисленных пар \((a,b)\), удовлетворяющих неравенству \(a^2+ab+41b^2\le N\). Реализации вычисляют \(F(N)\), переписывая форму и затем проходя только по одной координате.
Сделаем замену
$$x=2a+b,\qquad y=b.$$
Тогда
$$4(a^2+ab+41b^2)=(2a+b)^2+163b^2=x^2+163y^2.$$
Следовательно, исходное неравенство эквивалентно условию
$$x^2+163y^2\le 4N.$$
Число \(163\) здесь возникает естественно: это модуль дискриминанта \(1-4\cdot 41=-163\).
Замена обратима:
$$a=\frac{x-y}{2},\qquad b=y.$$
Значит, пара \((a,b)\) целочисленна тогда и только тогда, когда \(x-y\) чётно, то есть
$$x\equiv y\pmod{2}.$$
Отсюда получаем эквивалентную формулу подсчёта
$$F(N)=\#\left\{(x,y)\in\mathbb{Z}^2:x^2+163y^2\le 4N,\ x\equiv y\pmod{2}\right\}-1.$$
Последнее \(-1\) убирает начало координат, единственную точку, в которой квадратичная форма равна \(0\).
Для любой допустимой точки необходимо
$$163y^2\le 4N.$$
Поэтому
$$|y|\le Y=\left\lfloor\sqrt{\frac{4N}{163}}\right\rfloor.$$
То есть нужно рассмотреть только \(2Y+1\) значений \(y\). Именно здесь двумерный перебор сокращается до одномерного прохода по горизонтальным сечениям эллипса.
Зафиксируем целое \(y\) с \(|y|\le Y\) и введём
$$R(y)=4N-163y^2,\qquad M(y)=\left\lfloor\sqrt{R(y)}\right\rfloor.$$
Тогда условие \(x^2\le R(y)\) равносильно
$$-M(y)\le x\le M(y).$$
Без ограничения по чётности этот отрезок содержит \(2M(y)+1\) целых чисел. Нам нужны только те \(x\), у которых чётность совпадает с чётностью \(y\).
Интервал \([-M,M]\) симметричен относительно нуля. Если \(M\) чётно, в нём \(M+1\) чётных чисел и \(M\) нечётных. Если \(M\) нечётно, всё наоборот. Поэтому число допустимых значений \(x\) при фиксированном \(y\) равно
$$C(y)=\begin{cases} M(y)+1,& y\equiv M(y)\pmod{2},\\ M(y),& y\not\equiv M(y)\pmod{2}. \end{cases}$$
В этом и состоит ключевое упрощение: вместо полного внутреннего перебора нужен только один целочисленный квадратный корень и одна проверка чётности.
Складывая вклад всех допустимых значений \(y\), получаем
$$F(N)=\sum_{y=-Y}^{Y} C(y)-1.$$
Так как \(M(y)=M(-y)\), область симметрична. Тем не менее реализации просто проходят по всему диапазону целиком: так код остаётся короче и не требует отдельной обработки случая \(y=0\).
При \(N=1000\)
$$Y=\left\lfloor\sqrt{\frac{4000}{163}}\right\rfloor=4.$$
Сечения имеют вид
$$\begin{aligned} y=0 &:\quad R=4000,\quad M=63,\quad C=63,\\ y=\pm 1 &:\quad R=3837,\quad M=61,\quad C=62,\\ y=\pm 2 &:\quad R=3348,\quad M=57,\quad C=57,\\ y=\pm 3 &:\quad R=2533,\quad M=50,\quad C=50,\\ y=\pm 4 &:\quad R=1392,\quad M=37,\quad C=37. \end{aligned}$$
Следовательно,
$$F(1000)=63+2(62+57+50+37)-1=474.$$
Это совпадает со значением, которое выдаёт реализация, и тем самым подтверждает корректность преобразования и формулы чётности.
Реализации на C++, Python и Java устроены одинаково с математической точки зрения. Сначала они вычисляют точный целый квадратный корень. В Python для этого используется точная стандартная процедура, а в C++ и Java берётся приближение через вещественный корень, после чего результат корректируется вверх или вниз, пока не получится точный целый пол. Это важно, потому что при таком размере входа даже ошибка на единицу рядом с квадратной границей изменит итоговый ответ.
Далее реализация находит максимально возможное по модулю значение внешней координаты, проходит по всем целым в этом диапазоне, вычисляет остаток \(4N-163y^2\), берёт от него целый квадратный корень и применяет приведённую выше формулу чётности. Каждая итерация соответствует одному горизонтальному сечению эллипса, а после завершения цикла начало координат вычитается один раз.
Длина цикла равна
$$2\left\lfloor\sqrt{\frac{4N}{163}}\right\rfloor+1=\Theta(\sqrt N).$$
На каждой итерации выполняются константное число арифметических операций и один целочисленный квадратный корень. Поэтому общее время работы равно \(O(\sqrt N)\), а память — \(O(1)\).
نريد عدَّ الأزواج الصحيحة \((a,b)\neq(0,0)\) التي تحقق
$$a^2+ab+41b^2\le N.$$
هذه صيغة تربيعية ثنائية موجبة التعريف، ولذلك فإن جميع الأزواج الممكنة تقع داخل منطقة محدودة. لكن الفحص المباشر على متغيرين يظل بطيئًا جدًا عندما يكون \(N=10^{16}\). الفكرة الأساسية في الحل هي إجراء تحويل خطي للمتغيرات يحول الصيغة إلى متباينة إهليلجية مصطفة مع المحاور، مع شرط إضافي على التماثل الزوجي والفردي.
لنرمز بـ \(F(N)\) إلى عدد الأزواج الصحيحة غير الصفرية \((a,b)\) التي تحقق \(a^2+ab+41b^2\le N\). تحسب التطبيقات هذا المقدار بإعادة كتابة الصيغة ثم الدوران على إحداثي واحد فقط.
نأخذ التحويل
$$x=2a+b,\qquad y=b.$$
وعندئذ نحصل على
$$4(a^2+ab+41b^2)=(2a+b)^2+163b^2=x^2+163y^2.$$
إذًا تصبح المتباينة الأصلية مكافئة لـ
$$x^2+163y^2\le 4N.$$
وظهور العدد \(163\) طبيعي هنا، لأنه القيمة المطلقة للمميِّز \(1-4\cdot 41=-163\).
التحويل قابل للعكس:
$$a=\frac{x-y}{2},\qquad b=y.$$
ولذلك يكون \((a,b)\) زوجًا صحيحًا بالضبط عندما يكون \(x-y\) عددًا زوجيًا، أي عندما
$$x\equiv y\pmod{2}.$$
ومن ثم يصبح العد مكافئًا للصيغة
$$F(N)=\#\left\{(x,y)\in\mathbb{Z}^2:x^2+163y^2\le 4N,\ x\equiv y\pmod{2}\right\}-1.$$
والطرح الأخير بمقدار \(1\) يحذف نقطة الأصل، لأنها النقطة الوحيدة التي تعطي القيمة \(0\).
لكل نقطة مسموحة يجب أن يتحقق
$$163y^2\le 4N.$$
ومن هنا نحصل على
$$|y|\le Y=\left\lfloor\sqrt{\frac{4N}{163}}\right\rfloor.$$
أي إننا لا نحتاج إلا إلى فحص \(2Y+1\) قيمة صحيحة للمتغير \(y\). بهذه الخطوة يتحول البحث من مسألة ثنائية الأبعاد إلى مسح أحادي البعد عبر المقاطع الأفقية للإهليلج.
نثبت عددًا صحيحًا \(y\) بحيث \(|y|\le Y\)، ونعرف
$$R(y)=4N-163y^2,\qquad M(y)=\left\lfloor\sqrt{R(y)}\right\rfloor.$$
عندئذ يكون الشرط \(x^2\le R(y)\) مكافئًا لـ
$$-M(y)\le x\le M(y).$$
من دون شرط التماثل الزوجي والفردي يحتوي هذا المجال على \(2M(y)+1\) عددًا صحيحًا. لكننا نحتفظ فقط بالقيم التي تملك نفس زوجية \(y\).
الفترة \([-M,M]\) متناظرة حول الصفر. إذا كان \(M\) زوجيًا فهي تحتوي على \(M+1\) أعداد زوجية و\(M\) أعداد فردية، وإذا كان \(M\) فرديًا تنعكس الأدوار. لذلك فإن عدد قيم \(x\) المسموح بها عند تثبيت \(y\) يساوي
$$C(y)=\begin{cases} M(y)+1,& y\equiv M(y)\pmod{2},\\ M(y),& y\not\equiv M(y)\pmod{2}. \end{cases}$$
وهذا هو التبسيط الحاسم في الحل: بدل تعداد كل قيم \(x\) داخليًا، تكفي جذور تربيعية صحيحة واختبار زوجية واحد.
بعد جمع مساهمة كل قيمة ممكنة لـ \(y\) نحصل على
$$F(N)=\sum_{y=-Y}^{Y} C(y)-1.$$
ولأن \(M(y)=M(-y)\)، فالشكل الهندسي متناظر. ومع ذلك تنفذ التطبيقات حلقة مباشرة على المجال كله، لأن ذلك يجعل المنطق أبسط ولا يحتاج إلى معالجة خاصة عند \(y=0\).
عندما يكون \(N=1000\)، فإن
$$Y=\left\lfloor\sqrt{\frac{4000}{163}}\right\rfloor=4.$$
وتصبح المقاطع كما يلي:
$$\begin{aligned} y=0 &:\quad R=4000,\quad M=63,\quad C=63,\\ y=\pm 1 &:\quad R=3837,\quad M=61,\quad C=62,\\ y=\pm 2 &:\quad R=3348,\quad M=57,\quad C=57,\\ y=\pm 3 &:\quad R=2533,\quad M=50,\quad C=50,\\ y=\pm 4 &:\quad R=1392,\quad M=37,\quad C=37. \end{aligned}$$
إذن
$$F(1000)=63+2(62+57+50+37)-1=474.$$
وهذه هي القيمة نفسها التي تنتجها التطبيقات، ولذلك فهي تحقق جيد من صحة التحويل وصيغة العد.
تتبع تطبيقات ++C وPython وJava الخطة الرياضية نفسها. أولًا تُحسب قيمة الجذر التربيعي الصحيح بدقة. في Python يتم ذلك عبر دالة معيارية دقيقة، أما في ++C وJava فتبدأ العملية بتقريب عشري ثم يُصحح الناتج صعودًا أو هبوطًا حتى يصل إلى القيمة الصحيحة المكافئة لـ \(\lfloor\sqrt{n}\rfloor\). هذه الدقة مهمة جدًا، لأن خطأً بمقدار \(1\) قرب حد مربع كامل يكفي لتغيير الحصيلة النهائية عندما تكون القيم كبيرة.
بعد ذلك تحدد الشيفرة أكبر قيمة مطلقة ممكنة للإحداثي الخارجي، ثم تمر على كل عدد صحيح في هذا المجال، وتحسب المتبقي \(4N-163y^2\)، وتأخذ جذره التربيعي الصحيح، ثم تطبق صيغة الزوجية السابقة لتحصل على عدد قيم الإحداثي الداخلي الموافقة. كل دورة تمثل مقطعًا أفقيًا واحدًا من الإهليلج، وفي النهاية تُطرح نقطة الأصل مرة واحدة.
طول الحلقة يساوي
$$2\left\lfloor\sqrt{\frac{4N}{163}}\right\rfloor+1=\Theta(\sqrt N).$$
وفي كل دورة توجد عمليات حسابية ثابتة بالإضافة إلى جذر تربيعي صحيح واحد، لذا يكون الزمن الكلي \(O(\sqrt N)\)، بينما تبقى الذاكرة \(O(1)\) لأن الخوارزمية تخزن عددًا قليلًا فقط من المتغيرات الصحيحة.