Project Euler 372 defines \(R(M,N)\) as the number of lattice points \((x,y)\) such that \(M \lt x \le N\), \(M \lt y \le N\), and
$$\left\lfloor \frac{y^2}{x^2} \right\rfloor$$
is odd. The implementation uses the parameter names m and n, so below we write
\(R(m,n)\) and set \(\ell=m+1\). We must therefore count all integer pairs
$$\ell \le x \le n,\qquad \ell \le y \le n,$$
for which the squared slope \((y/x)^2\) falls into an odd strip \([k,k+1)\).
For every odd integer \(k \ge 1\), define
$$A_k=\left\{(x,y)\in \mathbb{Z}^2:\ell \le x,y \le n,\ k \le \frac{y^2}{x^2} \lt k+1\right\}.$$
The sets \(A_k\) are disjoint, and every admissible lattice point belongs to exactly one such set because the floor value is unique. Therefore
$$R(m,n)=\sum_{\substack{k \ge 1\\k\text{ odd}}} |A_k|.$$
Since \(x \ge \ell\) and \(y \le n\), we have
$$\frac{y^2}{x^2}\le \frac{n^2}{\ell^2},$$
so only odd \(k\) up to
$$K_{\max}=\left\lfloor\frac{n^2}{\ell^2}\right\rfloor$$
can contribute. This is the outer loop in all three solution files.
Fix an odd \(k\) and an \(x\). The condition
$$k \le \frac{y^2}{x^2} \lt k+1$$
is equivalent to
$$x\sqrt{k}\le y \lt x\sqrt{k+1}.$$
Because \(y\) is an integer and we also need \(y \le n\), the number of admissible \(y\)-values is
$$C_k(x)=\max\!\left(0,\ \min\!\bigl(n+1,\lceil x\sqrt{k+1}\rceil\bigr)-\lceil x\sqrt{k}\rceil\right).$$
The form with \(n+1\) is convenient because the count of integers in a half-open interval \([a,b)\cap \mathbb{Z}\) is \(\lceil b\rceil-\lceil a\rceil\), and truncating at \(y \le n\) means replacing the upper endpoint by \(n+1\).
The lower bound \(y \ge \lceil x\sqrt{k}\rceil\) is impossible once \(x\sqrt{k} \gt n\). Therefore define
$$u_k=\left\lfloor\frac{n}{\sqrt{k}}\right\rfloor.$$
Only \(x \le u_k\) can contribute.
The upper endpoint changes behavior when \(x\sqrt{k+1}\) crosses \(n+1\). Define
$$v_k=\left\lfloor\frac{n+1}{\sqrt{k+1}}\right\rfloor.$$
Then we have two cases:
$$C_k(x)=\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil \qquad (\ell \le x \le \min(u_k,v_k)),$$
$$C_k(x)=(n+1)-\lceil x\sqrt{k}\rceil \qquad (\max(\ell,v_k+1)\le x \le u_k).$$
This is exactly why the code splits each odd \(k\) into the intervals \([\ell,\min(u_k,v_k)]\) and \([\max(\ell,v_k+1),u_k]\).
Let
$$T_d(a,b)=\sum_{x=a}^{b}\lceil x\sqrt{d}\rceil.$$
For one odd \(k\), set
$$x_1=\min(u_k,v_k),\qquad x_2=\max(\ell,v_k+1).$$
Then the full contribution of \(k\) is
$$\sum_{x=\ell}^{x_1}\bigl(\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil\bigr) +\sum_{x=x_2}^{u_k}\bigl((n+1)-\lceil x\sqrt{k}\rceil\bigr),$$
that is,
$$T_{k+1}(\ell,x_1)-T_k(\ell,x_1)+(u_k-x_2+1)(n+1)-T_k(x_2,u_k).$$
So the geometric lattice-point problem has been converted into fast evaluation of interval sums of the form \(T_d(a,b)\).
If \(d=s^2\) is a perfect square, then \(\lceil x\sqrt{d}\rceil=sx\), so
$$T_d(a,b)=s\sum_{x=a}^{b}x=s\frac{(a+b)(b-a+1)}{2}.$$
The non-square case is the interesting one. Define the prefix sum
$$F_d(N)=\sum_{x=1}^{N}\lfloor x\sqrt{d}\rfloor.$$
For non-square \(d\), every \(x\sqrt{d}\) is irrational, hence \(\lceil x\sqrt{d}\rceil=\lfloor x\sqrt{d}\rfloor+1\), and therefore
$$T_d(a,b)=F_d(b)-F_d(a-1)+(b-a+1).$$
The helper class computes a more general quantity
$$F_{d,p,q}(N)=\sum_{i=1}^{N}\left\lfloor i\frac{\sqrt{d}+p}{q}\right\rfloor,$$
with the target value \(F_d(N)=F_{d,0,1}(N)\). Write
$$\alpha=\frac{\sqrt{d}+p}{q}=a+\beta,\qquad a=\lfloor \alpha \rfloor,\qquad 0 \lt \beta \lt 1.$$
The transformed parameters used by the code are
$$p_1=aq-p,\qquad q_1=\frac{d-p_1^2}{q}.$$
Then
$$\beta=\frac{\sqrt{d}-p_1}{q},\qquad \frac{1}{\beta}=\frac{\sqrt{d}+p_1}{q_1}.$$
If we set
$$r=\left\lfloor N\beta \right\rfloor=\left\lfloor \frac{N\sqrt{d}-Np_1}{q}\right\rfloor,$$
the complementary Beatty identity yields
$$\sum_{i=1}^{N}\lfloor i\beta\rfloor = Nr - F_{d,p_1,q_1}(r).$$
Since \(\lfloor i\alpha\rfloor = ai + \lfloor i\beta\rfloor\), we obtain the recurrence used verbatim by the implementation:
$$\boxed{F_{d,p,q}(N)=a\frac{N(N+1)}{2}+Nr-F_{d,p_1,q_1}(r).}$$
Memoization on the key \((d,p,q,N)\) makes repeated prefix queries cheap, which is essential because neighboring odd \(k\)-layers call the same surd sums again and again.
Take \(m=0\), \(n=5\), so \(\ell=1\) and \(K_{\max}=\lfloor 25/1\rfloor=25\). In practice only \(k=1\) contributes. For \(k=1\),
$$u_1=\left\lfloor\frac{5}{1}\right\rfloor=5,\qquad v_1=\left\lfloor\frac{6}{\sqrt{2}}\right\rfloor=4.$$
So
$$\sum_{x=1}^{4}\bigl(\lceil x\sqrt{2}\rceil-\lceil x\rceil\bigr) +\sum_{x=5}^{5}\bigl(6-\lceil x\rceil\bigr)$$
becomes
$$[(2-1)+(3-2)+(5-3)+(6-4)] + (6-5)=6+1=7.$$
The valid points are \((1,1),(2,2),(3,3),(3,4),(4,4),(4,5),(5,5)\), confirming the decomposition.
The C++, Python, and Java programs implement the same structure.
First, k_max = (n*n)/(l*l) bounds the odd layers. For each odd \(k\), the helpers
floor_div_by_sqrt / floor_div_sqrt / isqrt((t*t)/d) compute
\(\lfloor t/\sqrt{d}\rfloor\) exactly using integer square roots, avoiding floating-point boundary errors.
Next, sum_ceil_sqrt returns \(T_d(a,b)\). Perfect squares are handled by the closed arithmetic-series
formula. Otherwise the program asks BeattySumSqrt for prefix floor sums and converts them into ceiling
sums. The recursive helper starts from a floating estimate for
\(\left\lfloor\frac{a\sqrt{d}+b}{c}\right\rfloor\), then corrects the estimate by exact integer comparisons, so the
final result is mathematically exact.
Finally, the total contribution of each odd \(k\) is added to a big integer accumulator. The C++ version also checks the published checkpoints \(R(0,100)=3019\) and \(R(100,10000)=29750422\) before computing the final answer.
Let \(K_{\max}=\lfloor n^2/\ell^2\rfloor\). The outer enumeration touches only the odd integers up to \(K_{\max}\), so there are \(O(K_{\max})\) main iterations. Each iteration performs a constant number of exact boundary computations and up to three interval sums \(T_d(a,b)\).
When \(d\) is a square, the interval sum is \(O(1)\). Otherwise it is evaluated by the memoized quadratic-irrational recurrence above; the recursion depth is small in practice and repeated states are reused aggressively. A safe way to summarize the implementation is
$$O\!\bigl(K_{\max}\cdot T_{\mathrm{surd}}\bigr)\ \text{time},$$
where \(T_{\mathrm{surd}}\) is the amortized cost of one memoized Beatty query, with memory proportional to the number of cached states. For the actual parameters \(m=2\cdot 10^6\) and \(n=10^9\), this is fast enough because \(K_{\max}\approx 2.5\cdot 10^5\), not \(10^{18}\).
Project Euler 372 definiert \(R(M,N)\) als die Anzahl der Gitterpunkte \((x,y)\) mit \(M \lt x \le N\), \(M \lt y \le N\), für die
$$\left\lfloor \frac{y^2}{x^2} \right\rfloor$$
ungerade ist. Die Implementierungen verwenden die Parameternamen m und n, daher schreiben
wir im Folgenden \(R(m,n)\) und setzen \(\ell=m+1\). Gesucht sind also alle ganzzahligen Paare
$$\ell \le x \le n,\qquad \ell \le y \le n,$$
deren quadratische Steigung \((y/x)^2\) in einem ungeraden Streifen \([k,k+1)\) liegt.
Für jede ungerade ganze Zahl \(k \ge 1\) definieren wir
$$A_k=\left\{(x,y)\in \mathbb{Z}^2:\ell \le x,y \le n,\ k \le \frac{y^2}{x^2} \lt k+1\right\}.$$
Die Mengen \(A_k\) sind disjunkt, und jeder zulässige Gitterpunkt liegt genau in einer dieser Mengen, weil der Floor-Wert eindeutig ist. Also gilt
$$R(m,n)=\sum_{\substack{k \ge 1\\k\text{ ungerade}}} |A_k|.$$
Aus \(x \ge \ell\) und \(y \le n\) folgt
$$\frac{y^2}{x^2}\le \frac{n^2}{\ell^2},$$
also können nur ungerade \(k\) bis
$$K_{\max}=\left\lfloor\frac{n^2}{\ell^2}\right\rfloor$$
beitragen. Genau darüber läuft die äußere Schleife der Programme.
Fixiere ein ungerades \(k\) und ein \(x\). Dann ist
$$k \le \frac{y^2}{x^2} \lt k+1$$
äquivalent zu
$$x\sqrt{k}\le y \lt x\sqrt{k+1}.$$
Da \(y\) ganzzahlig sein und zusätzlich \(y \le n\) gelten muss, ist die Anzahl möglicher \(y\)-Werte
$$C_k(x)=\max\!\left(0,\ \min\!\bigl(n+1,\lceil x\sqrt{k+1}\rceil\bigr)-\lceil x\sqrt{k}\rceil\right).$$
Die Schreibweise mit \(n+1\) ist praktisch, weil die Anzahl ganzer Zahlen im halboffenen Intervall \([a,b)\cap \mathbb{Z}\) gleich \(\lceil b\rceil-\lceil a\rceil\) ist; die Zusatzbedingung \(y \le n\) ersetzt den oberen Endpunkt einfach durch \(n+1\).
Die untere Schranke \(y \ge \lceil x\sqrt{k}\rceil\) ist unmöglich, sobald \(x\sqrt{k} \gt n\). Daher definieren wir
$$u_k=\left\lfloor\frac{n}{\sqrt{k}}\right\rfloor.$$
Nur \(x \le u_k\) kann also überhaupt beitragen.
Der obere Endpunkt ändert sein Verhalten, sobald \(x\sqrt{k+1}\) die Grenze \(n+1\) überschreitet. Dazu setzen wir
$$v_k=\left\lfloor\frac{n+1}{\sqrt{k+1}}\right\rfloor.$$
Damit entstehen genau zwei Fälle:
$$C_k(x)=\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil \qquad (\ell \le x \le \min(u_k,v_k)),$$
$$C_k(x)=(n+1)-\lceil x\sqrt{k}\rceil \qquad (\max(\ell,v_k+1)\le x \le u_k).$$
Genau deshalb zerlegt der Code jeden ungeraden \(k\)-Beitrag in die Intervalle \([\ell,\min(u_k,v_k)]\) und \([\max(\ell,v_k+1),u_k]\).
Wir definieren
$$T_d(a,b)=\sum_{x=a}^{b}\lceil x\sqrt{d}\rceil.$$
Für ein fixes ungerades \(k\) setzen wir
$$x_1=\min(u_k,v_k),\qquad x_2=\max(\ell,v_k+1).$$
Dann ist der Gesamtbeitrag dieses \(k\)
$$\sum_{x=\ell}^{x_1}\bigl(\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil\bigr) +\sum_{x=x_2}^{u_k}\bigl((n+1)-\lceil x\sqrt{k}\rceil\bigr),$$
also
$$T_{k+1}(\ell,x_1)-T_k(\ell,x_1)+(u_k-x_2+1)(n+1)-T_k(x_2,u_k).$$
Damit wird das geometrische Gitterpunktproblem auf die schnelle Auswertung von Intervallsummen der Form \(T_d(a,b)\) zurückgeführt.
Ist \(d=s^2\) ein Quadrat, dann gilt \(\lceil x\sqrt{d}\rceil=sx\), also
$$T_d(a,b)=s\sum_{x=a}^{b}x=s\frac{(a+b)(b-a+1)}{2}.$$
Interessant ist der nichtquadratische Fall. Definiere die Präfixsumme
$$F_d(N)=\sum_{x=1}^{N}\lfloor x\sqrt{d}\rfloor.$$
Für nichtquadratisches \(d\) ist jedes \(x\sqrt{d}\) irrational, also \(\lceil x\sqrt{d}\rceil=\lfloor x\sqrt{d}\rfloor+1\). Damit folgt
$$T_d(a,b)=F_d(b)-F_d(a-1)+(b-a+1).$$
Die Hilfsklasse berechnet die allgemeinere Summe
$$F_{d,p,q}(N)=\sum_{i=1}^{N}\left\lfloor i\frac{\sqrt{d}+p}{q}\right\rfloor,$$
wobei \(F_d(N)=F_{d,0,1}(N)\) unser Ziel ist. Schreibe
$$\alpha=\frac{\sqrt{d}+p}{q}=a+\beta,\qquad a=\lfloor \alpha \rfloor,\qquad 0 \lt \beta \lt 1.$$
Die im Code verwendeten transformierten Parameter sind
$$p_1=aq-p,\qquad q_1=\frac{d-p_1^2}{q}.$$
Dann gilt
$$\beta=\frac{\sqrt{d}-p_1}{q},\qquad \frac{1}{\beta}=\frac{\sqrt{d}+p_1}{q_1}.$$
Setzen wir
$$r=\left\lfloor N\beta \right\rfloor=\left\lfloor \frac{N\sqrt{d}-Np_1}{q}\right\rfloor,$$
liefert die komplementäre Beatty-Identität
$$\sum_{i=1}^{N}\lfloor i\beta\rfloor = Nr - F_{d,p_1,q_1}(r).$$
Wegen \(\lfloor i\alpha\rfloor = ai + \lfloor i\beta\rfloor\) folgt exakt die Rekursion der Implementierung:
$$\boxed{F_{d,p,q}(N)=a\frac{N(N+1)}{2}+Nr-F_{d,p_1,q_1}(r).}$$
Durch Memoisierung über dem Schlüssel \((d,p,q,N)\) werden wiederholte Präfixanfragen billig, was entscheidend ist, weil benachbarte ungerade \(k\)-Schichten dieselben Surdensummen mehrfach benötigen.
Nehmen wir \(m=0\), \(n=5\). Dann ist \(\ell=1\) und \(K_{\max}=\lfloor 25/1\rfloor=25\). Tatsächlich trägt nur \(k=1\) etwas bei. Für \(k=1\) ist
$$u_1=\left\lfloor\frac{5}{1}\right\rfloor=5,\qquad v_1=\left\lfloor\frac{6}{\sqrt{2}}\right\rfloor=4.$$
Damit erhält man
$$\sum_{x=1}^{4}\bigl(\lceil x\sqrt{2}\rceil-\lceil x\rceil\bigr) +\sum_{x=5}^{5}\bigl(6-\lceil x\rceil\bigr),$$
also
$$[(2-1)+(3-2)+(5-3)+(6-4)] + (6-5)=6+1=7.$$
Die gültigen Punkte sind \((1,1),(2,2),(3,3),(3,4),(4,4),(4,5),(5,5)\), was die Zerlegung bestätigt.
Die Fassungen in C++, Python und Java setzen genau dieselbe Struktur um.
Zuerst beschränkt k_max = (n*n)/(l*l) die ungeraden Schichten. Für jedes ungerade \(k\) berechnen die
Hilfsfunktionen floor_div_by_sqrt / floor_div_sqrt /
isqrt((t*t)/d) den Wert \(\lfloor t/\sqrt{d}\rfloor\) exakt per ganzzahliger Quadratwurzel; dadurch
werden Floating-Point-Grenzfehler vermieden.
Danach liefert sum_ceil_sqrt die Größe \(T_d(a,b)\). Quadrate werden direkt per arithmetischer
Summenformel behandelt. Im nichtquadratischen Fall fragt das Programm bei BeattySumSqrt Präfixsummen
von Floors ab und verwandelt sie in Ceiling-Summen. Die rekursive Hilfsfunktion startet mit einer Gleitkomma-Schätzung
für \(\left\lfloor\frac{a\sqrt{d}+b}{c}\right\rfloor\) und korrigiert diese anschließend nur noch mit exakten
Ganzzahlvergleichen; das Endergebnis ist also mathematisch exakt.
Zum Schluss wird der Beitrag jedes ungeraden \(k\) zu einem großen Akkumulator addiert. Die C++-Version prüft vor der Endberechnung zusätzlich die veröffentlichten Kontrollwerte \(R(0,100)=3019\) und \(R(100,10000)=29750422\).
Sei \(K_{\max}=\lfloor n^2/\ell^2\rfloor\). Die äußere Enumeration betrachtet nur die ungeraden Zahlen bis \(K_{\max}\), also \(O(K_{\max})\) Hauptiterationen. Jede Iteration führt eine konstante Anzahl exakter Grenzberechnungen und bis zu drei Intervallsummen \(T_d(a,b)\) aus.
Ist \(d\) ein Quadrat, dann kostet die Intervallsumme \(O(1)\). Andernfalls wird sie durch die oben beschriebene, memoisierten Rekursion über quadratische Irrationale ausgewertet; die Rekursionstiefe ist in der Praxis klein, und wiederkehrende Zustände werden stark wiederverwendet. Eine sichere Zusammenfassung ist daher
$$O\!\bigl(K_{\max}\cdot T_{\mathrm{surd}}\bigr)\ \text{Zeit},$$
wobei \(T_{\mathrm{surd}}\) die amortisierten Kosten einer memoisierten Beatty-Anfrage bezeichnet; der Speicher ist proportional zur Anzahl der zwischengespeicherten Zustände. Für die echten Parameter \(m=2\cdot 10^6\) und \(n=10^9\) ist das praktikabel, weil \(K_{\max}\approx 2{,}5\cdot 10^5\) ist und nicht \(10^{18}\).
Project Euler 372, \(R(M,N)\) fonksiyonunu \(M \lt x \le N\), \(M \lt y \le N\) koşullarını sağlayan ve
$$\left\lfloor \frac{y^2}{x^2} \right\rfloor$$
değeri tek olan kafes noktalarının \((x,y)\) sayısı olarak tanımlar. Kod dosyaları parametre adları olarak
m ve n kullandığı için aşağıda \(R(m,n)\) yazacağız ve \(\ell=m+1\) tanımını
kullanacağız. Böylece saymamız gereken tam sayı çiftleri
$$\ell \le x \le n,\qquad \ell \le y \le n,$$
ve \((y/x)^2\) değeri tek numaralı bir şerit olan \([k,k+1)\) aralığına düşmelidir.
Her tek tamsayı \(k \ge 1\) için
$$A_k=\left\{(x,y)\in \mathbb{Z}^2:\ell \le x,y \le n,\ k \le \frac{y^2}{x^2} \lt k+1\right\}$$
kümesini tanımlayalım. Bu kümeler ayrışıktır ve her uygun nokta, floor değeri tek ve sabit olduğu için tam olarak bir kümeye düşer. Dolayısıyla
$$R(m,n)=\sum_{\substack{k \ge 1\\k\text{ tek}}} |A_k|$$
elde edilir. \(x \ge \ell\) ve \(y \le n\) olduğundan
$$\frac{y^2}{x^2}\le \frac{n^2}{\ell^2}$$
olur; bu yüzden yalnızca
$$K_{\max}=\left\lfloor\frac{n^2}{\ell^2}\right\rfloor$$
değerine kadar olan tek \(k\)'lar katkı verir. Çözüm dosyalarındaki dış döngü tam olarak budur.
Tek bir \(k\) ve sabit bir \(x\) seçelim. O zaman
$$k \le \frac{y^2}{x^2} \lt k+1$$
koşulu
$$x\sqrt{k}\le y \lt x\sqrt{k+1}$$
ile aynıdır. \(y\) tam sayı olmalı ve ayrıca \(y \le n\) sınırı vardır. Bu nedenle geçerli \(y\) sayısı
$$C_k(x)=\max\!\left(0,\ \min\!\bigl(n+1,\lceil x\sqrt{k+1}\rceil\bigr)-\lceil x\sqrt{k}\rceil\right)$$
olur. Buradaki \(n+1\) yazımı kullanışlıdır; çünkü \([a,b)\cap\mathbb{Z}\) içindeki tam sayı sayısı \(\lceil b\rceil-\lceil a\rceil\) olur ve \(y \le n\) şartı üst sınırı doğrudan \(n+1\) ile kesmek anlamına gelir.
Alt sınır \(y \ge \lceil x\sqrt{k}\rceil\), \(x\sqrt{k} \gt n\) olduğunda artık mümkün değildir. Bu yüzden
$$u_k=\left\lfloor\frac{n}{\sqrt{k}}\right\rfloor$$
tanımını yaparız. Yalnızca \(x \le u_k\) değerleri katkı verebilir.
Üst sınır ise \(x\sqrt{k+1}\) ifadesi \(n+1\) eşiğini geçtiğinde biçim değiştirir. Bunun için
$$v_k=\left\lfloor\frac{n+1}{\sqrt{k+1}}\right\rfloor$$
tanımını kullanırız. Böylece tam olarak iki durum vardır:
$$C_k(x)=\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil \qquad (\ell \le x \le \min(u_k,v_k)),$$
$$C_k(x)=(n+1)-\lceil x\sqrt{k}\rceil \qquad (\max(\ell,v_k+1)\le x \le u_k).$$
Koddaki her tek \(k\) için yapılan iki parçalı toplama işleminin nedeni tam olarak budur.
Şunu tanımlayalım:
$$T_d(a,b)=\sum_{x=a}^{b}\lceil x\sqrt{d}\rceil.$$
Sabit bir tek \(k\) için
$$x_1=\min(u_k,v_k),\qquad x_2=\max(\ell,v_k+1)$$
olsun. O zaman bu \(k\)'nın toplam katkısı
$$\sum_{x=\ell}^{x_1}\bigl(\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil\bigr) +\sum_{x=x_2}^{u_k}\bigl((n+1)-\lceil x\sqrt{k}\rceil\bigr)$$
ve dolayısıyla
$$T_{k+1}(\ell,x_1)-T_k(\ell,x_1)+(u_k-x_2+1)(n+1)-T_k(x_2,u_k)$$
şeklindedir. Böylece geometrik nokta sayımı, \(T_d(a,b)\) türünden aralık toplamlarını hızlı hesaba dönüştürülür.
Eğer \(d=s^2\) tam kare ise, \(\lceil x\sqrt{d}\rceil=sx\) olur ve
$$T_d(a,b)=s\sum_{x=a}^{b}x=s\frac{(a+b)(b-a+1)}{2}$$
elde edilir. Asıl ilginç durum tam kare olmayan \(d\) değerleridir. Şu önek toplamını tanımlayalım:
$$F_d(N)=\sum_{x=1}^{N}\lfloor x\sqrt{d}\rfloor.$$
\(d\) tam kare değilse her \(x\sqrt{d}\) irrasyoneldir; dolayısıyla \(\lceil x\sqrt{d}\rceil=\lfloor x\sqrt{d}\rfloor+1\). O halde
$$T_d(a,b)=F_d(b)-F_d(a-1)+(b-a+1)$$
olur. Yardımcı sınıf bunun daha genel halini hesaplar:
$$F_{d,p,q}(N)=\sum_{i=1}^{N}\left\lfloor i\frac{\sqrt{d}+p}{q}\right\rfloor,$$
ve bizim hedefimiz \(F_d(N)=F_{d,0,1}(N)\)'dir. Şimdi
$$\alpha=\frac{\sqrt{d}+p}{q}=a+\beta,\qquad a=\lfloor \alpha \rfloor,\qquad 0 \lt \beta \lt 1$$
yazalım. Kodun kullandığı dönüştürülmüş parametreler
$$p_1=aq-p,\qquad q_1=\frac{d-p_1^2}{q}$$
şeklindedir. Böylece
$$\beta=\frac{\sqrt{d}-p_1}{q},\qquad \frac{1}{\beta}=\frac{\sqrt{d}+p_1}{q_1}$$
olur. Ayrıca
$$r=\left\lfloor N\beta \right\rfloor=\left\lfloor \frac{N\sqrt{d}-Np_1}{q}\right\rfloor$$
tanımıyla, tamamlayıcı Beatty özdeşliği
$$\sum_{i=1}^{N}\lfloor i\beta\rfloor = Nr - F_{d,p_1,q_1}(r)$$
sonucunu verir. \(\lfloor i\alpha\rfloor = ai + \lfloor i\beta\rfloor\) olduğundan kodun birebir kullandığı bağıntı
$$\boxed{F_{d,p,q}(N)=a\frac{N(N+1)}{2}+Nr-F_{d,p_1,q_1}(r)}$$
olur. \((d,p,q,N)\) anahtarıyla yapılan memoization, aynı köklü toplamların tekrar tekrar hesaplanmasını önler.
\(m=0\), \(n=5\) alalım. O zaman \(\ell=1\) ve \(K_{\max}=\lfloor 25/1\rfloor=25\) olur; pratikte yalnızca \(k=1\) katkı verir. \(k=1\) için
$$u_1=\left\lfloor\frac{5}{1}\right\rfloor=5,\qquad v_1=\left\lfloor\frac{6}{\sqrt{2}}\right\rfloor=4$$
elde edilir. Katkı
$$\sum_{x=1}^{4}\bigl(\lceil x\sqrt{2}\rceil-\lceil x\rceil\bigr) +\sum_{x=5}^{5}\bigl(6-\lceil x\rceil\bigr)$$
yani
$$[(2-1)+(3-2)+(5-3)+(6-4)] + (6-5)=6+1=7$$
olur. Geçerli noktalar \((1,1),(2,2),(3,3),(3,4),(4,4),(4,5),(5,5)\) şeklindedir; ayrıştırma tam olarak bunu üretir.
C++, Python ve Java sürümleri aynı mantığı uygular.
İlk olarak k_max = (n*n)/(l*l) ile tek katmanların üst sınırı alınır. Her tek \(k\) için
floor_div_by_sqrt / floor_div_sqrt /
isqrt((t*t)/d) yardımcıları \(\lfloor t/\sqrt{d}\rfloor\) değerini tam olarak hesaplar; böylece kayan
nokta sınır hataları ortadan kalkar.
Daha sonra sum_ceil_sqrt, \(T_d(a,b)\) toplamını döndürür. Tam kareler doğrudan aritmetik dizi formülüyle
hesaplanır. Tam kare olmayan durumda program, BeattySumSqrt içinden önek floor toplamlarını alır ve
bunları ceiling toplamına çevirir. Özyinelemeli yardımcı, önce
\(\left\lfloor\frac{a\sqrt{d}+b}{c}\right\rfloor\) için bir kayan nokta tahmini kullanır, sonra sonucu yalnızca tam
sayı karşılaştırmaları ile düzeltir; sonuç bu yüzden tam doğrudur.
Son aşamada her tek \(k\)'nın katkısı büyük bir toplama eklenir. C++ sürümü ayrıca son cevaptan önce \(R(0,100)=3019\) ve \(R(100,10000)=29750422\) kontrol değerlerini de doğrular.
\(K_{\max}=\lfloor n^2/\ell^2\rfloor\) olsun. Dış döngü yalnızca \(K_{\max}\)'e kadar olan tek sayıları dolaşır; bu yüzden ana iterasyon sayısı \(O(K_{\max})\)'tir. Her iterasyonda sabit sayıda tam sınır hesabı ve en çok üç adet \(T_d(a,b)\) aralık toplamı vardır.
\(d\) tam kare ise bu toplam \(O(1)\) zamandır. Değilse, yukarıdaki memoization'lı kuadratik-irrasyonel bağıntı ile hesaplanır; özyineleme derinliği pratikte küçüktür ve tekrar eden durumlar yeniden kullanılır. Bu yüzden yöntem güvenli biçimde
$$O\!\bigl(K_{\max}\cdot T_{\mathrm{surd}}\bigr)\ \text{zaman}$$
olarak özetlenebilir; burada \(T_{\mathrm{surd}}\), memoize edilmiş bir Beatty sorgusunun amortize maliyetidir. Bellek ihtiyacı da önbelleğe alınan durum sayısıyla orantılıdır. Gerçek parametrelerde \(K_{\max}\approx 2{,}5\cdot 10^5\) olduğu için yöntem uygulanabilirdir.
Project Euler 372 define \(R(M,N)\) como el número de puntos de la red \((x,y)\) tales que \(M \lt x \le N\), \(M \lt y \le N\), y
$$\left\lfloor \frac{y^2}{x^2} \right\rfloor$$
es impar. Las implementaciones usan los nombres m y n, así que a continuación escribiremos
\(R(m,n)\) y \(\ell=m+1\). Debemos contar, por tanto, todos los pares enteros
$$\ell \le x \le n,\qquad \ell \le y \le n,$$
para los que la pendiente cuadrática \((y/x)^2\) cae dentro de una franja impar \([k,k+1)\).
Para cada entero impar \(k \ge 1\), definimos
$$A_k=\left\{(x,y)\in \mathbb{Z}^2:\ell \le x,y \le n,\ k \le \frac{y^2}{x^2} \lt k+1\right\}.$$
Los conjuntos \(A_k\) son disjuntos, y cada punto válido pertenece a exactamente uno de ellos porque el valor del piso es único. En consecuencia,
$$R(m,n)=\sum_{\substack{k \ge 1\\k\text{ impar}}} |A_k|.$$
Como \(x \ge \ell\) y \(y \le n\), tenemos
$$\frac{y^2}{x^2}\le \frac{n^2}{\ell^2},$$
así que solo pueden contribuir los \(k\) impares hasta
$$K_{\max}=\left\lfloor\frac{n^2}{\ell^2}\right\rfloor.$$
Ese es exactamente el límite del bucle exterior en los tres lenguajes.
Fijemos un \(k\) impar y un valor de \(x\). Entonces
$$k \le \frac{y^2}{x^2} \lt k+1$$
equivale a
$$x\sqrt{k}\le y \lt x\sqrt{k+1}.$$
Como \(y\) debe ser entero y además satisfacer \(y \le n\), el número de opciones es
$$C_k(x)=\max\!\left(0,\ \min\!\bigl(n+1,\lceil x\sqrt{k+1}\rceil\bigr)-\lceil x\sqrt{k}\rceil\right).$$
La forma con \(n+1\) es útil porque el número de enteros en un intervalo semiabierto \([a,b)\cap\mathbb{Z}\) es \(\lceil b\rceil-\lceil a\rceil\), y la restricción \(y \le n\) se incorpora simplemente sustituyendo el extremo superior por \(n+1\).
La cota inferior \(y \ge \lceil x\sqrt{k}\rceil\) deja de ser posible cuando \(x\sqrt{k} \gt n\). Por eso definimos
$$u_k=\left\lfloor\frac{n}{\sqrt{k}}\right\rfloor.$$
Solo \(x \le u_k\) puede aportar algo.
El comportamiento del extremo superior cambia cuando \(x\sqrt{k+1}\) cruza \(n+1\). Definimos entonces
$$v_k=\left\lfloor\frac{n+1}{\sqrt{k+1}}\right\rfloor.$$
Así aparecen dos casos exactos:
$$C_k(x)=\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil \qquad (\ell \le x \le \min(u_k,v_k)),$$
$$C_k(x)=(n+1)-\lceil x\sqrt{k}\rceil \qquad (\max(\ell,v_k+1)\le x \le u_k).$$
Por eso el código separa cada contribución impar de \(k\) en los intervalos \([\ell,\min(u_k,v_k)]\) y \([\max(\ell,v_k+1),u_k]\).
Definimos
$$T_d(a,b)=\sum_{x=a}^{b}\lceil x\sqrt{d}\rceil.$$
Para un \(k\) impar dado, sea
$$x_1=\min(u_k,v_k),\qquad x_2=\max(\ell,v_k+1).$$
La contribución total de ese \(k\) es
$$\sum_{x=\ell}^{x_1}\bigl(\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil\bigr) +\sum_{x=x_2}^{u_k}\bigl((n+1)-\lceil x\sqrt{k}\rceil\bigr),$$
es decir,
$$T_{k+1}(\ell,x_1)-T_k(\ell,x_1)+(u_k-x_2+1)(n+1)-T_k(x_2,u_k).$$
El problema geométrico de contar puntos de la red queda así convertido en evaluar con rapidez sumas del tipo \(T_d(a,b)\).
Si \(d=s^2\) es un cuadrado perfecto, entonces \(\lceil x\sqrt{d}\rceil=sx\), y por tanto
$$T_d(a,b)=s\sum_{x=a}^{b}x=s\frac{(a+b)(b-a+1)}{2}.$$
El caso interesante es cuando \(d\) no es cuadrado. Definimos la suma prefija
$$F_d(N)=\sum_{x=1}^{N}\lfloor x\sqrt{d}\rfloor.$$
Si \(d\) no es cuadrado, todo \(x\sqrt{d}\) es irracional y entonces \(\lceil x\sqrt{d}\rceil=\lfloor x\sqrt{d}\rfloor+1\). De ahí sale
$$T_d(a,b)=F_d(b)-F_d(a-1)+(b-a+1).$$
La clase auxiliar calcula una cantidad más general:
$$F_{d,p,q}(N)=\sum_{i=1}^{N}\left\lfloor i\frac{\sqrt{d}+p}{q}\right\rfloor,$$
de modo que \(F_d(N)=F_{d,0,1}(N)\) es el caso que nos interesa. Escribimos
$$\alpha=\frac{\sqrt{d}+p}{q}=a+\beta,\qquad a=\lfloor \alpha \rfloor,\qquad 0 \lt \beta \lt 1.$$
Los parámetros transformados usados por el programa son
$$p_1=aq-p,\qquad q_1=\frac{d-p_1^2}{q}.$$
Entonces
$$\beta=\frac{\sqrt{d}-p_1}{q},\qquad \frac{1}{\beta}=\frac{\sqrt{d}+p_1}{q_1}.$$
Si definimos
$$r=\left\lfloor N\beta \right\rfloor=\left\lfloor \frac{N\sqrt{d}-Np_1}{q}\right\rfloor,$$
la identidad complementaria de Beatty da
$$\sum_{i=1}^{N}\lfloor i\beta\rfloor = Nr - F_{d,p_1,q_1}(r).$$
Como \(\lfloor i\alpha\rfloor = ai + \lfloor i\beta\rfloor\), obtenemos exactamente la recurrencia implementada:
$$\boxed{F_{d,p,q}(N)=a\frac{N(N+1)}{2}+Nr-F_{d,p_1,q_1}(r).}$$
La memoización sobre la clave \((d,p,q,N)\) vuelve baratas las consultas repetidas, algo fundamental porque capas vecinas de \(k\) reutilizan una y otra vez las mismas sumas con radicales.
Tomemos \(m=0\) y \(n=5\). Entonces \(\ell=1\) y \(K_{\max}=\lfloor 25/1\rfloor=25\). En realidad solo contribuye \(k=1\). Para \(k=1\),
$$u_1=\left\lfloor\frac{5}{1}\right\rfloor=5,\qquad v_1=\left\lfloor\frac{6}{\sqrt{2}}\right\rfloor=4.$$
Por tanto, la cuenta es
$$\sum_{x=1}^{4}\bigl(\lceil x\sqrt{2}\rceil-\lceil x\rceil\bigr) +\sum_{x=5}^{5}\bigl(6-\lceil x\rceil\bigr),$$
esto es,
$$[(2-1)+(3-2)+(5-3)+(6-4)] + (6-5)=6+1=7.$$
Los puntos válidos son \((1,1),(2,2),(3,3),(3,4),(4,4),(4,5),(5,5)\), lo que confirma la descomposición.
Las versiones en C++, Python y Java siguen la misma estructura.
Primero, k_max = (n*n)/(l*l) acota las capas impares. Para cada \(k\) impar, las funciones auxiliares
floor_div_by_sqrt / floor_div_sqrt / isqrt((t*t)/d) calculan
\(\lfloor t/\sqrt{d}\rfloor\) exactamente mediante raíces cuadradas enteras, evitando errores de frontera por coma
flotante.
Después, sum_ceil_sqrt devuelve \(T_d(a,b)\). Los cuadrados perfectos se resuelven con la fórmula
cerrada de progresión aritmética. En el caso no cuadrado, el programa pide a BeattySumSqrt las sumas
prefijas de floors y las convierte en sumas de ceilings. La rutina recursiva parte de una estimación en coma
flotante para \(\left\lfloor\frac{a\sqrt{d}+b}{c}\right\rfloor\), y luego corrige esa estimación con comparaciones
enteras exactas, de modo que el resultado final es exacto.
Por último, la contribución de cada \(k\) impar se añade a un acumulador grande. La versión en C++ también comprueba los puntos de control publicados \(R(0,100)=3019\) y \(R(100,10000)=29750422\) antes de calcular la respuesta final.
Sea \(K_{\max}=\lfloor n^2/\ell^2\rfloor\). La enumeración exterior recorre solo los enteros impares hasta \(K_{\max}\), así que hay \(O(K_{\max})\) iteraciones principales. En cada una se realizan un número constante de cálculos exactos de frontera y hasta tres sumas de intervalo \(T_d(a,b)\).
Cuando \(d\) es cuadrado, la suma cuesta \(O(1)\). En otro caso, se evalúa mediante la recurrencia memoizada sobre irracionales cuadráticos descrita arriba; la profundidad recursiva es pequeña en la práctica y los estados repetidos se reutilizan mucho. Una forma segura de resumir el coste es
$$O\!\bigl(K_{\max}\cdot T_{\mathrm{surd}}\bigr)\ \text{tiempo},$$
donde \(T_{\mathrm{surd}}\) es el coste amortizado de una consulta Beatty memoizada. La memoria es proporcional al número de estados almacenados. Para los parámetros reales, \(K_{\max}\approx 2{,}5\cdot 10^5\), así que el método es perfectamente viable.
Project Euler 372 把 \(R(M,N)\) 定义为满足 \(M \lt x \le N\)、\(M \lt y \le N\) 且
$$\left\lfloor \frac{y^2}{x^2} \right\rfloor$$
为奇数的格点 \((x,y)\) 的个数。代码里参数名使用的是 m 和 n，因此下面记作
\(R(m,n)\)，并设 \(\ell=m+1\)。也就是说，我们要统计所有满足
$$\ell \le x \le n,\qquad \ell \le y \le n$$
且斜率平方 \((y/x)^2\) 落在某个奇数带 \([k,k+1)\) 中的整点。
对每个奇数 \(k \ge 1\)，定义
$$A_k=\left\{(x,y)\in \mathbb{Z}^2:\ell \le x,y \le n,\ k \le \frac{y^2}{x^2} \lt k+1\right\}.$$
这些 \(A_k\) 两两不交，而且每个合法点都恰好落入一个这样的集合，因为 floor 值唯一。因此
$$R(m,n)=\sum_{\substack{k \ge 1\\k\text{ 为奇数}}} |A_k|.$$
由于 \(x \ge \ell\) 且 \(y \le n\)，有
$$\frac{y^2}{x^2}\le \frac{n^2}{\ell^2},$$
所以只有不超过
$$K_{\max}=\left\lfloor\frac{n^2}{\ell^2}\right\rfloor$$
的奇数 \(k\) 才可能产生贡献。这正是实现中的外层循环上界。
固定一个奇数 \(k\) 和一个 \(x\)。条件
$$k \le \frac{y^2}{x^2} \lt k+1$$
等价于
$$x\sqrt{k}\le y \lt x\sqrt{k+1}.$$
\(y\) 必须是整数，同时还要满足 \(y \le n\)。因此合法的 \(y\) 个数为
$$C_k(x)=\max\!\left(0,\ \min\!\bigl(n+1,\lceil x\sqrt{k+1}\rceil\bigr)-\lceil x\sqrt{k}\rceil\right).$$
写成 \(n+1\) 的好处是：半开区间 \([a,b)\cap\mathbb{Z}\) 中整数的个数就是 \(\lceil b\rceil-\lceil a\rceil\)，而限制 \(y \le n\) 只需要把上端点截断为 \(n+1\)。
当 \(x\sqrt{k} \gt n\) 时，下界 \(y \ge \lceil x\sqrt{k}\rceil\) 已经不可能满足，所以定义
$$u_k=\left\lfloor\frac{n}{\sqrt{k}}\right\rfloor.$$
只有 \(x \le u_k\) 才可能有贡献。
上界的行为则在 \(x\sqrt{k+1}\) 穿过 \(n+1\) 时发生变化，因此再定义
$$v_k=\left\lfloor\frac{n+1}{\sqrt{k+1}}\right\rfloor.$$
于是恰好分成两种情况：
$$C_k(x)=\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil \qquad (\ell \le x \le \min(u_k,v_k)),$$
$$C_k(x)=(n+1)-\lceil x\sqrt{k}\rceil \qquad (\max(\ell,v_k+1)\le x \le u_k).$$
这就是代码为何把每个奇数 \(k\) 的贡献拆成 \([\ell,\min(u_k,v_k)]\) 与 \([\max(\ell,v_k+1),u_k]\) 两段来处理。
定义
$$T_d(a,b)=\sum_{x=a}^{b}\lceil x\sqrt{d}\rceil.$$
对固定奇数 \(k\)，令
$$x_1=\min(u_k,v_k),\qquad x_2=\max(\ell,v_k+1).$$
则该 \(k\) 的总贡献为
$$\sum_{x=\ell}^{x_1}\bigl(\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil\bigr) +\sum_{x=x_2}^{u_k}\bigl((n+1)-\lceil x\sqrt{k}\rceil\bigr),$$
也就是
$$T_{k+1}(\ell,x_1)-T_k(\ell,x_1)+(u_k-x_2+1)(n+1)-T_k(x_2,u_k).$$
这样一来，几何计数问题就被转化为高效计算 \(T_d(a,b)\) 这样的区间和。
若 \(d=s^2\) 是完全平方数，那么 \(\lceil x\sqrt{d}\rceil=sx\)，于是
$$T_d(a,b)=s\sum_{x=a}^{b}x=s\frac{(a+b)(b-a+1)}{2}.$$
真正关键的是 \(d\) 不是平方数的情况。定义前缀和
$$F_d(N)=\sum_{x=1}^{N}\lfloor x\sqrt{d}\rfloor.$$
对于非平方数 \(d\)，每个 \(x\sqrt{d}\) 都是无理数，因此 \(\lceil x\sqrt{d}\rceil=\lfloor x\sqrt{d}\rfloor+1\)。所以
$$T_d(a,b)=F_d(b)-F_d(a-1)+(b-a+1).$$
辅助类实际计算的是更一般的量
$$F_{d,p,q}(N)=\sum_{i=1}^{N}\left\lfloor i\frac{\sqrt{d}+p}{q}\right\rfloor,$$
其中我们的目标就是 \(F_d(N)=F_{d,0,1}(N)\)。写成
$$\alpha=\frac{\sqrt{d}+p}{q}=a+\beta,\qquad a=\lfloor \alpha \rfloor,\qquad 0 \lt \beta \lt 1.$$
代码中的变换参数为
$$p_1=aq-p,\qquad q_1=\frac{d-p_1^2}{q}.$$
于是
$$\beta=\frac{\sqrt{d}-p_1}{q},\qquad \frac{1}{\beta}=\frac{\sqrt{d}+p_1}{q_1}.$$
再令
$$r=\left\lfloor N\beta \right\rfloor=\left\lfloor \frac{N\sqrt{d}-Np_1}{q}\right\rfloor,$$
由互补 Beatty 恒等式得到
$$\sum_{i=1}^{N}\lfloor i\beta\rfloor = Nr - F_{d,p_1,q_1}(r).$$
由于 \(\lfloor i\alpha\rfloor = ai + \lfloor i\beta\rfloor\)，便得到实现中原封不动使用的递推式：
$$\boxed{F_{d,p,q}(N)=a\frac{N(N+1)}{2}+Nr-F_{d,p_1,q_1}(r).}$$
对键 \((d,p,q,N)\) 做记忆化后，重复的前缀查询就会非常便宜，而相邻的奇数层 \(k\) 正好会大量重用同一类根式求和。
取 \(m=0\)、\(n=5\)。此时 \(\ell=1\)，且 \(K_{\max}=\lfloor 25/1\rfloor=25\)。实际只有 \(k=1\) 有贡献。 对 \(k=1\)，有
$$u_1=\left\lfloor\frac{5}{1}\right\rfloor=5,\qquad v_1=\left\lfloor\frac{6}{\sqrt{2}}\right\rfloor=4.$$
因此总数为
$$\sum_{x=1}^{4}\bigl(\lceil x\sqrt{2}\rceil-\lceil x\rceil\bigr) +\sum_{x=5}^{5}\bigl(6-\lceil x\rceil\bigr),$$
即
$$[(2-1)+(3-2)+(5-3)+(6-4)] + (6-5)=6+1=7.$$
对应的合法点是 \((1,1),(2,2),(3,3),(3,4),(4,4),(4,5),(5,5)\)， 与分层公式完全一致。
C++、Python 和 Java 三个版本实现的是同一套结构。
首先用 k_max = (n*n)/(l*l) 确定奇数层上界。对每个奇数 \(k\)，辅助函数
floor_div_by_sqrt / floor_div_sqrt / isqrt((t*t)/d) 通过整数平方根精确计算
\(\lfloor t/\sqrt{d}\rfloor\)，从而避免浮点临界误差。
随后 sum_ceil_sqrt 负责返回 \(T_d(a,b)\)。若 \(d\) 是完全平方数，就直接套用等差求和公式；
否则程序向 BeattySumSqrt 查询前缀 floor 和，再转换成 ceiling 和。递归辅助函数先对
\(\left\lfloor\frac{a\sqrt{d}+b}{c}\right\rfloor\) 做一次浮点估计，再用精确整数比较把估计修正到真正的值，
因而最终结果仍然是严格正确的。
最后，把每个奇数 \(k\) 的贡献加到大整数累加器里。C++ 版本还会在正式求解前验证公开给出的检查点 \(R(0,100)=3019\) 和 \(R(100,10000)=29750422\)。
令 \(K_{\max}=\lfloor n^2/\ell^2\rfloor\)。外层枚举只遍历到 \(K_{\max}\) 为止的奇数，因此主循环是 \(O(K_{\max})\) 次。每次迭代只做常数个边界计算，以及最多三个区间和 \(T_d(a,b)\)。
若 \(d\) 为平方数，则区间和是 \(O(1)\)。否则它通过上面的二次无理数递推加记忆化来求值；递归深度在实践中较小， 而重复状态会被大量复用。安全的复杂度写法是
$$O\!\bigl(K_{\max}\cdot T_{\mathrm{surd}}\bigr)\ \text{时间},$$
其中 \(T_{\mathrm{surd}}\) 表示一次记忆化 Beatty 查询的摊还代价，内存则与缓存状态数量成正比。对实际参数而言， \(K_{\max}\approx 2.5\times 10^5\)，因此算法是可行的。
В задаче Project Euler 372 функция \(R(M,N)\) определяется как число решётчатых точек \((x,y)\), для которых \(M \lt x \le N\), \(M \lt y \le N\), и
$$\left\lfloor \frac{y^2}{x^2} \right\rfloor$$
нечётно. В коде используются имена параметров m и n, поэтому ниже мы пишем \(R(m,n)\) и
вводим \(\ell=m+1\). Значит, нужно посчитать все целочисленные пары
$$\ell \le x \le n,\qquad \ell \le y \le n,$$
для которых квадрат наклона \((y/x)^2\) попадает в нечётную полосу \([k,k+1)\).
Для каждого нечётного целого \(k \ge 1\) определим
$$A_k=\left\{(x,y)\in \mathbb{Z}^2:\ell \le x,y \le n,\ k \le \frac{y^2}{x^2} \lt k+1\right\}.$$
Множества \(A_k\) не пересекаются, и каждая допустимая точка попадает ровно в одно из них, поскольку значение floor единственно. Поэтому
$$R(m,n)=\sum_{\substack{k \ge 1\\k\text{ нечётно}}} |A_k|.$$
Так как \(x \ge \ell\) и \(y \le n\), получаем
$$\frac{y^2}{x^2}\le \frac{n^2}{\ell^2},$$
а значит, вклад возможен только для нечётных \(k\) не больше
$$K_{\max}=\left\lfloor\frac{n^2}{\ell^2}\right\rfloor.$$
Именно до этого предела идёт внешний цикл в программах.
Зафиксируем нечётное \(k\) и значение \(x\). Тогда условие
$$k \le \frac{y^2}{x^2} \lt k+1$$
эквивалентно
$$x\sqrt{k}\le y \lt x\sqrt{k+1}.$$
Так как \(y\) должно быть целым и дополнительно удовлетворять \(y \le n\), число допустимых значений равно
$$C_k(x)=\max\!\left(0,\ \min\!\bigl(n+1,\lceil x\sqrt{k+1}\rceil\bigr)-\lceil x\sqrt{k}\rceil\right).$$
Форма с \(n+1\) удобна потому, что количество целых чисел в полуинтервале \([a,b)\cap\mathbb{Z}\) равно \(\lceil b\rceil-\lceil a\rceil\), а ограничение \(y \le n\) просто заменяет верхнюю границу на \(n+1\).
Нижняя граница \(y \ge \lceil x\sqrt{k}\rceil\) становится невозможной, как только \(x\sqrt{k} \gt n\). Поэтому вводим
$$u_k=\left\lfloor\frac{n}{\sqrt{k}}\right\rfloor.$$
Только \(x \le u_k\) может дать вклад.
Поведение верхней границы меняется, когда \(x\sqrt{k+1}\) пересекает порог \(n+1\). Для этого определяем
$$v_k=\left\lfloor\frac{n+1}{\sqrt{k+1}}\right\rfloor.$$
Тогда существуют ровно два случая:
$$C_k(x)=\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil \qquad (\ell \le x \le \min(u_k,v_k)),$$
$$C_k(x)=(n+1)-\lceil x\sqrt{k}\rceil \qquad (\max(\ell,v_k+1)\le x \le u_k).$$
Именно поэтому код разбивает вклад каждого нечётного \(k\) на интервалы \([\ell,\min(u_k,v_k)]\) и \([\max(\ell,v_k+1),u_k]\).
Определим
$$T_d(a,b)=\sum_{x=a}^{b}\lceil x\sqrt{d}\rceil.$$
Для фиксированного нечётного \(k\) положим
$$x_1=\min(u_k,v_k),\qquad x_2=\max(\ell,v_k+1).$$
Тогда полный вклад этого \(k\) равен
$$\sum_{x=\ell}^{x_1}\bigl(\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil\bigr) +\sum_{x=x_2}^{u_k}\bigl((n+1)-\lceil x\sqrt{k}\rceil\bigr),$$
то есть
$$T_{k+1}(\ell,x_1)-T_k(\ell,x_1)+(u_k-x_2+1)(n+1)-T_k(x_2,u_k).$$
Таким образом, геометрическая задача о точках решётки сводится к быстрому вычислению сумм вида \(T_d(a,b)\).
Если \(d=s^2\) является точным квадратом, то \(\lceil x\sqrt{d}\rceil=sx\), и потому
$$T_d(a,b)=s\sum_{x=a}^{b}x=s\frac{(a+b)(b-a+1)}{2}.$$
Интересен случай, когда \(d\) не квадрат. Введём префиксную сумму
$$F_d(N)=\sum_{x=1}^{N}\lfloor x\sqrt{d}\rfloor.$$
Для неквадратного \(d\) каждое число \(x\sqrt{d}\) иррационально, значит \(\lceil x\sqrt{d}\rceil=\lfloor x\sqrt{d}\rfloor+1\), и поэтому
$$T_d(a,b)=F_d(b)-F_d(a-1)+(b-a+1).$$
Вспомогательный класс считает более общую величину
$$F_{d,p,q}(N)=\sum_{i=1}^{N}\left\lfloor i\frac{\sqrt{d}+p}{q}\right\rfloor,$$
где нужный нам случай есть \(F_d(N)=F_{d,0,1}(N)\). Запишем
$$\alpha=\frac{\sqrt{d}+p}{q}=a+\beta,\qquad a=\lfloor \alpha \rfloor,\qquad 0 \lt \beta \lt 1.$$
Преобразованные параметры в коде таковы:
$$p_1=aq-p,\qquad q_1=\frac{d-p_1^2}{q}.$$
Тогда
$$\beta=\frac{\sqrt{d}-p_1}{q},\qquad \frac{1}{\beta}=\frac{\sqrt{d}+p_1}{q_1}.$$
Если обозначить
$$r=\left\lfloor N\beta \right\rfloor=\left\lfloor \frac{N\sqrt{d}-Np_1}{q}\right\rfloor,$$
то дополнительная тождественность Битти даёт
$$\sum_{i=1}^{N}\lfloor i\beta\rfloor = Nr - F_{d,p_1,q_1}(r).$$
Так как \(\lfloor i\alpha\rfloor = ai + \lfloor i\beta\rfloor\), получаем ровно ту рекурсию, которая реализована в программе:
$$\boxed{F_{d,p,q}(N)=a\frac{N(N+1)}{2}+Nr-F_{d,p_1,q_1}(r).}$$
Мемоизация по ключу \((d,p,q,N)\) делает повторные запросы дешёвыми, а соседние нечётные слои \(k\) как раз многократно используют одни и те же суммы с квадратными иррациональностями.
Возьмём \(m=0\) и \(n=5\). Тогда \(\ell=1\), а \(K_{\max}=\lfloor 25/1\rfloor=25\). На практике вклад даёт только \(k=1\). Для \(k=1\) имеем
$$u_1=\left\lfloor\frac{5}{1}\right\rfloor=5,\qquad v_1=\left\lfloor\frac{6}{\sqrt{2}}\right\rfloor=4.$$
Следовательно, суммарный вклад равен
$$\sum_{x=1}^{4}\bigl(\lceil x\sqrt{2}\rceil-\lceil x\rceil\bigr) +\sum_{x=5}^{5}\bigl(6-\lceil x\rceil\bigr),$$
то есть
$$[(2-1)+(3-2)+(5-3)+(6-4)] + (6-5)=6+1=7.$$
Подходящие точки: \((1,1),(2,2),(3,3),(3,4),(4,4),(4,5),(5,5)\), что полностью подтверждает разбиение.
Версии на C++, Python и Java реализуют одну и ту же схему.
Сначала k_max = (n*n)/(l*l) ограничивает число нечётных слоёв. Для каждого нечётного \(k\) вспомогательные
функции floor_div_by_sqrt / floor_div_sqrt /
isqrt((t*t)/d) точно вычисляют \(\lfloor t/\sqrt{d}\rfloor\) с помощью целочисленного квадратного корня,
избегая пограничных ошибок плавающей точки.
Затем sum_ceil_sqrt возвращает \(T_d(a,b)\). Для точных квадратов используется явная формула суммы
арифметической прогрессии. В неквадратном случае программа запрашивает префиксные суммы floor у класса
BeattySumSqrt и превращает их в суммы ceiling. Рекурсивный помощник стартует с вещественной оценки для
\(\left\lfloor\frac{a\sqrt{d}+b}{c}\right\rfloor\), а затем исправляет её только точными целочисленными сравнениями,
так что итоговый результат остаётся строго точным.
После этого вклад каждого нечётного \(k\) добавляется в большой аккумулятор. Вариант на C++ дополнительно проверяет опубликованные контрольные значения \(R(0,100)=3019\) и \(R(100,10000)=29750422\) перед финальным вычислением.
Пусть \(K_{\max}=\lfloor n^2/\ell^2\rfloor\). Внешний перебор рассматривает только нечётные числа до \(K_{\max}\), то есть выполняет \(O(K_{\max})\) основных итераций. В каждой итерации делается константное число точных вычислений границ и не более трёх интервальных сумм \(T_d(a,b)\).
Если \(d\) является квадратом, сумма считается за \(O(1)\). Иначе используется описанная выше мемоизированная рекурсия по квадратическим иррациональностям; глубина рекурсии на практике невелика, а повторяющиеся состояния активно переиспользуются. Поэтому безопасная сводка такова:
$$O\!\bigl(K_{\max}\cdot T_{\mathrm{surd}}\bigr)\ \text{времени},$$
где \(T_{\mathrm{surd}}\) означает амортизированную стоимость одного мемоизированного запроса Beatty, а память пропорциональна числу закешированных состояний. Для реальных параметров \(K_{\max}\approx 2{,}5\cdot 10^5\), поэтому метод остаётся практичным.
تعرف مسألة Project Euler 372 الدالة \(R(M,N)\) بأنها عدد نقاط الشبكة \((x,y)\) التي تحقق \(M \lt x \le N\) و \(M \lt y \le N\) ويكون فيها
$$\left\lfloor \frac{y^2}{x^2} \right\rfloor$$
عددًا فرديًا. ملفات الحل تستعمل الاسمين m وn، لذلك سنكتب أدناه \(R(m,n)\) ونضع
\(\ell=m+1\). إذن المطلوب هو عد جميع الأزواج الصحيحة
$$\ell \le x \le n,\qquad \ell \le y \le n,$$
بحيث يقع مربع الميل \((y/x)^2\) داخل شريط فردي من الصورة \([k,k+1)\).
لكل عدد صحيح فردي \(k \ge 1\) نعرّف
$$A_k=\left\{(x,y)\in \mathbb{Z}^2:\ell \le x,y \le n,\ k \le \frac{y^2}{x^2} \lt k+1\right\}.$$
هذه المجموعات متباينة، وكل نقطة صالحة تنتمي إلى واحدة منها فقط لأن قيمة الـ floor محددة بشكل وحيد. لذلك
$$R(m,n)=\sum_{\substack{k \ge 1\\k \equiv 1 \pmod 2}} |A_k|.$$
وبما أن \(x \ge \ell\) و\(y \le n\)، فإن
$$\frac{y^2}{x^2}\le \frac{n^2}{\ell^2},$$
ومن ثم لا يمكن أن يسهم إلا \(k\) الفردي حتى الحد
$$K_{\max}=\left\lfloor\frac{n^2}{\ell^2}\right\rfloor.$$
وهذا هو حد الحلقة الخارجية في جميع الحلول.
لنثبت عددًا فرديًا \(k\) وقيمة \(x\). الشرط
$$k \le \frac{y^2}{x^2} \lt k+1$$
يكافئ
$$x\sqrt{k}\le y \lt x\sqrt{k+1}.$$
وبما أن \(y\) يجب أن يكون عددًا صحيحًا وأن يحقق أيضًا \(y \le n\)، فإن عدد القيم المسموح بها هو
$$C_k(x)=\max\!\left(0,\ \min\!\bigl(n+1,\lceil x\sqrt{k+1}\rceil\bigr)-\lceil x\sqrt{k}\rceil\right).$$
صيغة \(n+1\) مفيدة لأن عدد الأعداد الصحيحة في المجال شبه المفتوح \([a,b)\cap\mathbb{Z}\) يساوي \(\lceil b\rceil-\lceil a\rceil\)، وشرط \(y \le n\) يعني فقط قطع الحد العلوي عند \(n+1\).
يصبح القيد السفلي \(y \ge \lceil x\sqrt{k}\rceil\) مستحيلاً عندما \(x\sqrt{k} \gt n\). لذلك نعرّف
$$u_k=\left\lfloor\frac{n}{\sqrt{k}}\right\rfloor.$$
إذن لا يمكن أن تسهم إلا القيم \(x \le u_k\).
أما الحد العلوي فيتغير سلوكه عندما يتجاوز \(x\sqrt{k+1}\) القيمة \(n+1\). لهذا نعرّف
$$v_k=\left\lfloor\frac{n+1}{\sqrt{k+1}}\right\rfloor.$$
وعندها تظهر حالتان بالضبط:
$$C_k(x)=\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil \qquad (\ell \le x \le \min(u_k,v_k)),$$
$$C_k(x)=(n+1)-\lceil x\sqrt{k}\rceil \qquad (\max(\ell,v_k+1)\le x \le u_k).$$
ولهذا السبب يقسم الكود مساهمة كل \(k\) فردي إلى المجالين \([\ell,\min(u_k,v_k)]\) و\([\max(\ell,v_k+1),u_k]\).
لنعرّف
$$T_d(a,b)=\sum_{x=a}^{b}\lceil x\sqrt{d}\rceil.$$
ولعدد فردي ثابت \(k\) نضع
$$x_1=\min(u_k,v_k),\qquad x_2=\max(\ell,v_k+1).$$
فتصبح مساهمة هذا \(k\) هي
$$\sum_{x=\ell}^{x_1}\bigl(\lceil x\sqrt{k+1}\rceil-\lceil x\sqrt{k}\rceil\bigr) +\sum_{x=x_2}^{u_k}\bigl((n+1)-\lceil x\sqrt{k}\rceil\bigr),$$
أي
$$T_{k+1}(\ell,x_1)-T_k(\ell,x_1)+(u_k-x_2+1)(n+1)-T_k(x_2,u_k).$$
وهكذا تتحول مسألة العد الهندسي إلى تقييم سريع لمجاميع على فترات من النوع \(T_d(a,b)\).
إذا كان \(d=s^2\) مربعًا كاملًا، فإن \(\lceil x\sqrt{d}\rceil=sx\)، ومن ثم
$$T_d(a,b)=s\sum_{x=a}^{b}x=s\frac{(a+b)(b-a+1)}{2}.$$
الحالة المهمة هي عندما لا يكون \(d\) مربعًا. نعرّف المجموع التراكمي
$$F_d(N)=\sum_{x=1}^{N}\lfloor x\sqrt{d}\rfloor.$$
ولأن \(d\) غير مربع، فإن كل \(x\sqrt{d}\) يكون غير نسبي، وبالتالي \(\lceil x\sqrt{d}\rceil=\lfloor x\sqrt{d}\rfloor+1\). لذا
$$T_d(a,b)=F_d(b)-F_d(a-1)+(b-a+1).$$
الفئة المساعدة تحسب الكمية الأعم
$$F_{d,p,q}(N)=\sum_{i=1}^{N}\left\lfloor i\frac{\sqrt{d}+p}{q}\right\rfloor,$$
والقيمة المطلوبة لدينا هي \(F_d(N)=F_{d,0,1}(N)\). نكتب
$$\alpha=\frac{\sqrt{d}+p}{q}=a+\beta,\qquad a=\lfloor \alpha \rfloor,\qquad 0 \lt \beta \lt 1.$$
أما المعاملات المحولة التي يستعملها الكود فهي
$$p_1=aq-p,\qquad q_1=\frac{d-p_1^2}{q}.$$
وعندها
$$\beta=\frac{\sqrt{d}-p_1}{q},\qquad \frac{1}{\beta}=\frac{\sqrt{d}+p_1}{q_1}.$$
إذا وضعنا
$$r=\left\lfloor N\beta \right\rfloor=\left\lfloor \frac{N\sqrt{d}-Np_1}{q}\right\rfloor,$$
فإن هوية Beatty التكميلية تعطي
$$\sum_{i=1}^{N}\lfloor i\beta\rfloor = Nr - F_{d,p_1,q_1}(r).$$
ومع \(\lfloor i\alpha\rfloor = ai + \lfloor i\beta\rfloor\) نحصل مباشرة على العلاقة المستخدمة حرفيًا في التنفيذ:
$$\boxed{F_{d,p,q}(N)=a\frac{N(N+1)}{2}+Nr-F_{d,p_1,q_1}(r).}$$
التخزين المؤقت للحالات بالمفتاح \((d,p,q,N)\) يجعل الاستدعاءات المتكررة رخيصة، وهذا ضروري لأن الطبقات المتجاورة في \(k\) تعيد استخدام مجاميع الجذور نفسها مرات كثيرة.
لنأخذ \(m=0\) و\(n=5\). عندئذ \(\ell=1\) و\(K_{\max}=\lfloor 25/1\rfloor=25\). عمليًا لا يسهم إلا \(k=1\). بالنسبة إلى \(k=1\) نحصل على
$$u_1=\left\lfloor\frac{5}{1}\right\rfloor=5,\qquad v_1=\left\lfloor\frac{6}{\sqrt{2}}\right\rfloor=4.$$
فتصبح المساهمة
$$\sum_{x=1}^{4}\bigl(\lceil x\sqrt{2}\rceil-\lceil x\rceil\bigr) +\sum_{x=5}^{5}\bigl(6-\lceil x\rceil\bigr),$$
أي
$$[(2-1)+(3-2)+(5-3)+(6-4)] + (6-5)=6+1=7.$$
والنقاط الصالحة هي \((1,1),(2,2),(3,3),(3,4),(4,4),(4,5),(5,5)\)، وهذا يطابق التحليل تمامًا.
نسخ C++ وPython وJava تطبق البنية نفسها.
أولاً، يحدد k_max = (n*n)/(l*l) أقصى الطبقات الفردية. ولكل \(k\) فردي تقوم الدوال
floor_div_by_sqrt وfloor_div_sqrt وisqrt((t*t)/d) بحساب
\(\lfloor t/\sqrt{d}\rfloor\) بدقة باستخدام الجذر التربيعي الصحيح، وبذلك تتجنب أخطاء الحواف المرتبطة بالأعداد
العشرية العائمة.
بعد ذلك تعيد sum_ceil_sqrt القيمة \(T_d(a,b)\). إذا كان \(d\) مربعًا كاملًا فالحساب مباشر بصيغة
متسلسلة حسابية. وإذا لم يكن كذلك، يطلب البرنامج من BeattySumSqrt مجاميع prefix من نوع floor ثم
يحولها إلى مجاميع ceiling. الدالة العودية تبدأ بتقدير عشري تقريبي لـ
\(\left\lfloor\frac{a\sqrt{d}+b}{c}\right\rfloor\)، ثم تصححه بمقارنات صحيحة دقيقة فقط، ولذلك تكون النتيجة النهائية
صحيحة تمامًا.
وفي النهاية تضاف مساهمة كل \(k\) فردي إلى مجمّع كبير. نسخة C++ تتحقق أيضًا من قيمتي الفحص المنشورتين \(R(0,100)=3019\) و\(R(100,10000)=29750422\) قبل حساب الجواب النهائي.
ليكن \(K_{\max}=\lfloor n^2/\ell^2\rfloor\). التعداد الخارجي يمر فقط على الأعداد الفردية حتى \(K_{\max}\)، ولذلك لدينا \(O(K_{\max})\) دورة رئيسية. في كل دورة يوجد عدد ثابت من حسابات الحدود الدقيقة، وما يصل إلى ثلاثة مجاميع على فترات من الشكل \(T_d(a,b)\).
إذا كان \(d\) مربعًا كاملًا فالمجموع يحسب في \(O(1)\). وإلا فإنه يحسب بالعلاقة العودية المذكرة أعلاه على الأعداد غير النسبية التربيعية؛ عمق العودية صغير عمليًا، والحالات المتكررة يعاد استخدامها بكثافة. لذلك يمكن تلخيص الزمن بشكل آمن بأنه
$$O\!\bigl(K_{\max}\cdot T_{\mathrm{surd}}\bigr)$$
حيث تمثل \(T_{\mathrm{surd}}\) الكلفة الوسطية لاستعلام Beatty واحد مع التخزين المؤقت، أما الذاكرة فهي متناسبة مع عدد الحالات المخزنة. في القيم الحقيقية للمسألة يكون \(K_{\max}\approx 2.5\times 10^5\)، لذا تصبح الطريقة عملية.