Every Gaussian integer \(z=a+bi\) has a representation in base \(i-1\) using only the digits \(0\) and \(1\). Let \(f(a,b)\) be the number of digits equal to \(1\) in that expansion, and define
$$B(L)=\sum_{a=-L}^{L}\sum_{b=-L}^{L} f(a,b).$$
The goal is to compute \(B(10^{15}) \bmod (10^9+7)\). A direct scan of the full square would involve about \(4\times 10^{30}\) lattice points, so the only workable method is to exploit the arithmetic of division by \(i-1\).
The implementations do not enumerate Gaussian integers one by one. Instead, they turn the digit-extraction rule for a single point into mutually recursive formulas for sums over entire rectangles.
Write one digit step as
$$a+bi=r_0+(i-1)(a_1+b_1 i),\qquad r_0\in\{0,1\}.$$
Expanding \((i-1)(a_1+b_1 i)\) and matching real and imaginary parts gives
$$r_0\equiv a-b \pmod{2},$$
$$a_1=\frac{b-a+r_0}{2},\qquad b_1=\frac{r_0-a-b}{2}.$$
Therefore one digit contributes \(r_0\) ones, and the remaining tail is the representation of the smaller Gaussian integer \(a_1+b_1 i\):
$$f(a,b)=r_0+f(a_1,b_1).$$
This is the basic recurrence used everywhere in the solution.
For an axis-aligned rectangle define
$$S([a_{\min},a_{\max}]\times[b_{\min},b_{\max}])=\sum_{a=a_{\min}}^{a_{\max}}\sum_{b=b_{\min}}^{b_{\max}} f(a,b).$$
Then the required answer is simply
$$B(L)=S([-L,L]\times[-L,L]).$$
It is also convenient to write \(E([x_1,x_2])\) for the number of even integers in an interval and \(O([x_1,x_2])\) for the number of odd integers in that interval. Those counts are available by closed formulas, so the algorithm never loops just to count parity classes.
The first digit is \(r_0=(a-b)\bmod 2\). Hence all lattice points with \(a-b\) odd contribute an immediate \(+1\). The number of such points in a rectangle is
$$E([a_{\min},a_{\max}])\,O([b_{\min},b_{\max}])+O([a_{\min},a_{\max}])\,E([b_{\min},b_{\max}]).$$
What remains is the tail \(f(a_1,b_1)\). To keep the image of a rectangle axis-aligned, introduce auxiliary coordinates
$$u=a_1+b_1=r_0-a,\qquad v=a_1-b_1=b.$$
For fixed \(r_0\), the map \((a,b)\mapsto (u,v)\) sends a rectangle to another rectangle. This motivates an auxiliary sum \(T\) over \((u,v)\)-rectangles. The first recursion becomes
$$\begin{aligned} S([a_{\min},a_{\max}]\times[b_{\min},b_{\max}])&=N_{\mathrm{odd}}\\ &\quad +T([-a_{\max},-a_{\min}]\times[b_{\min},b_{\max}])\\ &\quad +T([1-a_{\max},1-a_{\min}]\times[b_{\min},b_{\max}]), \end{aligned}$$
where \(N_{\mathrm{odd}}\) is the parity-count term above. The two \(T\)-rectangles correspond to the two possible values \(r_0=0\) and \(r_0=1\).
Since \(a_1=(u+v)/2\) and \(b_1=(u-v)/2\), only pairs with the same parity can occur. Let their common parity be \(r_1\in\{0,1\}\). Then the second digit equals
$$r_1\equiv a_1-b_1\equiv v \pmod{2}.$$
Removing that digit once more gives
$$a_2=\frac{b_1-a_1+r_1}{2}=-\frac{v-r_1}{2},\qquad b_2=\frac{r_1-a_1-b_1}{2}=-\frac{u-r_1}{2}.$$
So the odd-odd points in a \((u,v)\)-rectangle contribute one more immediate \(+1\), while the even-even and odd-odd classes map back to smaller rectangles in the original \((a,b)\)-coordinates. If \(U=[u_{\min},u_{\max}]\) and \(V=[v_{\min},v_{\max}]\), then
$$\begin{aligned} T(U\times V)&=O(U)\,O(V)\\ &\quad +S\!\left(\left[-\left\lfloor\frac{v_{\max}}{2}\right\rfloor,-\left\lceil\frac{v_{\min}}{2}\right\rceil\right]\times\left[-\left\lfloor\frac{u_{\max}}{2}\right\rfloor,-\left\lceil\frac{u_{\min}}{2}\right\rceil\right]\right)\\ &\quad +S\!\left(\left[-\left\lfloor\frac{v_{\max}-1}{2}\right\rfloor,-\left\lceil\frac{v_{\min}-1}{2}\right\rceil\right]\times\left[-\left\lfloor\frac{u_{\max}-1}{2}\right\rfloor,-\left\lceil\frac{u_{\min}-1}{2}\right\rceil\right]\right). \end{aligned}$$
Any interval with lower bound greater than upper bound is simply empty and contributes \(0\).
The point \(8+0i\) is a small but informative example. Repeated digit extraction yields
$$8\to -4-4i\to 4i\to 2-2i\to -2\to 1+i\to -i\to i\to 1\to 0.$$
The corresponding digit sequence is
$$0,0,0,0,0,0,1,1,1,$$
so
$$f(8,0)=3.$$
Two larger checkpoints used by the implementations are
$$B(1)=20,\qquad B(2)=75.$$
They confirm that the rectangle recurrences reproduce the direct brute-force totals on small squares.
The desired value is
$$B(L)=S([-L,L]\times[-L,L]).$$
Each pair of recursion layers roughly halves the coordinate scale. That is why a huge square can be handled through repeated parity counting, affine transforms, and memoized subrectangles instead of explicit enumeration.
The C++, Python, and Java implementations all follow the same plan. First they provide a direct evaluator for a single Gaussian integer by repeatedly extracting the next base-\(i-1\) digit with the recurrence from Step 1. That direct evaluator is only used on very small rectangles.
For larger rectangles the implementation switches to the two recursive summation routines described above: one routine works in the original \((a,b)\)-coordinates, and the other works in the auxiliary \((u,v)\)-coordinates. Each call counts the immediate parity contribution by closed formulas, transforms the remaining work into one or two smaller rectangles, and stores the result in a memoization table keyed by the four rectangle boundaries.
The brute-force cutoff is an area of at most \(1024\) lattice points. Below that threshold, direct evaluation is cheaper than further recursive splitting. The C++ implementation also evaluates the two top-level transformed subproblems independently, because they are mathematically disjoint; the Python and Java implementations use the same recurrence without that extra parallel split. All reported values are reduced modulo \(10^9+7\).
For a single Gaussian integer, repeated division by \(i-1\) takes \(O(\log (|a|+|b|+1))\) steps because the coordinate scale shrinks rapidly. For the square sum, every two recursive layers reduce interval sizes by about a factor of \(2\), so the recursion depth is \(O(\log L)\). The total running time is proportional to the number of distinct memoized rectangles that actually appear, and the memory usage is of the same order. The exact state count is messy to write in closed form, but it is dramatically smaller than the \((2L+1)^2\) points of a direct scan, which is why the method remains practical even for \(L=10^{15}\).
Jede gaußsche ganze Zahl \(z=a+bi\) besitzt eine Darstellung zur Basis \(i-1\), in der nur die Ziffern \(0\) und \(1\) vorkommen. Sei \(f(a,b)\) die Anzahl der Ziffern \(1\) in dieser Darstellung, und definiere
$$B(L)=\sum_{a=-L}^{L}\sum_{b=-L}^{L} f(a,b).$$
Gesucht ist \(B(10^{15}) \bmod (10^9+7)\). Ein direkter Durchlauf über das ganze Quadrat würde etwa \(4\times 10^{30}\) Gitterpunkte berühren und ist damit ausgeschlossen. Entscheidend ist also, wie die Division durch \(i-1\) das Gitter transformiert.
Die Implementierungen zählen gaußsche ganze Zahlen nicht einzeln. Stattdessen wird die Ziffernrekursion eines einzelnen Punkts in zwei gegenseitig rekursive Formeln für ganze Rechteckssummen umgeschrieben.
Ein Ziffernschritt hat die Form
$$a+bi=r_0+(i-1)(a_1+b_1 i),\qquad r_0\in\{0,1\}.$$
Durch Ausmultiplizieren und Vergleich von Real- und Imaginärteil erhält man
$$r_0\equiv a-b \pmod{2},$$
$$a_1=\frac{b-a+r_0}{2},\qquad b_1=\frac{r_0-a-b}{2}.$$
Also liefert diese erste Ziffer genau \(r_0\) Einsen, und der Rest ist die Darstellung der kleineren gaußschen Zahl \(a_1+b_1 i\):
$$f(a,b)=r_0+f(a_1,b_1).$$
Diese Rekursion ist der Kern der gesamten Lösung.
Für ein achsenparalleles Rechteck definieren wir
$$S([a_{\min},a_{\max}]\times[b_{\min},b_{\max}])=\sum_{a=a_{\min}}^{a_{\max}}\sum_{b=b_{\min}}^{b_{\max}} f(a,b).$$
Dann ist die gesuchte Größe einfach
$$B(L)=S([-L,L]\times[-L,L]).$$
Außerdem schreiben wir \(E([x_1,x_2])\) für die Anzahl gerader und \(O([x_1,x_2])\) für die Anzahl ungerader Zahlen in einem Intervall. Diese Paritätszahlen lassen sich geschlossen berechnen; dafür braucht der Algorithmus keine Schleifen.
Die erste Ziffer ist \(r_0=(a-b)\bmod 2\). Damit liefern alle Gitterpunkte mit ungeradem \(a-b\) sofort einen Beitrag \(+1\). Ihre Anzahl in einem Rechteck ist
$$E([a_{\min},a_{\max}])\,O([b_{\min},b_{\max}])+O([a_{\min},a_{\max}])\,E([b_{\min},b_{\max}]).$$
Übrig bleibt der Restterm \(f(a_1,b_1)\). Damit das Bild eines Rechtecks wieder rechteckig bleibt, führen wir Hilfskoordinaten ein:
$$u=a_1+b_1=r_0-a,\qquad v=a_1-b_1=b.$$
Für festes \(r_0\) wird ein Rechteck unter \((a,b)\mapsto (u,v)\) wieder auf ein Rechteck abgebildet. Deshalb definieren wir eine Hilfssumme \(T\) über \((u,v)\)-Rechtecke. Die erste Rekursion lautet dann
$$\begin{aligned} S([a_{\min},a_{\max}]\times[b_{\min},b_{\max}])&=N_{\mathrm{odd}}\\ &\quad +T([-a_{\max},-a_{\min}]\times[b_{\min},b_{\max}])\\ &\quad +T([1-a_{\max},1-a_{\min}]\times[b_{\min},b_{\max}]), \end{aligned}$$
wobei \(N_{\mathrm{odd}}\) der Paritätsbeitrag aus der vorherigen Formel ist. Die beiden \(T\)-Rechtecke entsprechen den Fällen \(r_0=0\) und \(r_0=1\).
Aus \(a_1=(u+v)/2\) und \(b_1=(u-v)/2\) folgt, dass nur Paare gleicher Parität vorkommen können. Sei \(r_1\in\{0,1\}\) diese gemeinsame Parität. Dann ist die zweite Ziffer
$$r_1\equiv a_1-b_1\equiv v \pmod{2}.$$
Nach dem Abspalten dieser Ziffer erhält man
$$a_2=\frac{b_1-a_1+r_1}{2}=-\frac{v-r_1}{2},\qquad b_2=\frac{r_1-a_1-b_1}{2}=-\frac{u-r_1}{2}.$$
Also geben die ungerade-ungerade Punkte in einem \((u,v)\)-Rechteck noch ein unmittelbares \(+1\), während die Klassen gerade-gerade und ungerade-ungerade wieder auf kleinere Rechtecke in den ursprünglichen \((a,b)\)-Koordinaten zurückfallen. Für \(U=[u_{\min},u_{\max}]\) und \(V=[v_{\min},v_{\max}]\) gilt
$$\begin{aligned} T(U\times V)&=O(U)\,O(V)\\ &\quad +S\!\left(\left[-\left\lfloor\frac{v_{\max}}{2}\right\rfloor,-\left\lceil\frac{v_{\min}}{2}\right\rceil\right]\times\left[-\left\lfloor\frac{u_{\max}}{2}\right\rfloor,-\left\lceil\frac{u_{\min}}{2}\right\rceil\right]\right)\\ &\quad +S\!\left(\left[-\left\lfloor\frac{v_{\max}-1}{2}\right\rfloor,-\left\lceil\frac{v_{\min}-1}{2}\right\rceil\right]\times\left[-\left\lfloor\frac{u_{\max}-1}{2}\right\rfloor,-\left\lceil\frac{u_{\min}-1}{2}\right\rceil\right]\right). \end{aligned}$$
Ist ein Intervall leer, also untere Grenze größer als obere, trägt es einfach \(0\) bei.
Der Punkt \(8+0i\) ist klein, zeigt die Rekursion aber gut. Wiederholtes Ziffernabspalten ergibt
$$8\to -4-4i\to 4i\to 2-2i\to -2\to 1+i\to -i\to i\to 1\to 0.$$
Die zugehörige Ziffernfolge lautet
$$0,0,0,0,0,0,1,1,1,$$
also
$$f(8,0)=3.$$
Zwei weitere Kontrollwerte aus den Implementierungen sind
$$B(1)=20,\qquad B(2)=75.$$
Damit lässt sich prüfen, dass die Rechteckrekursionen auf kleinen Quadraten exakt mit der direkten Auswertung übereinstimmen.
Gesucht ist
$$B(L)=S([-L,L]\times[-L,L]).$$
Jeweils zwei Rekursionsebenen halbieren die Koordinatenskala ungefähr. Dadurch kann ein riesiges Quadrat über Paritätszählung, affine Transformationen und Memoisierung von Teilrechtecken behandelt werden, statt Punkte einzeln aufzuzählen.
Die C++-, Python- und Java-Implementierungen folgen demselben Bauplan. Zuerst gibt es einen direkten Auswerter für einen einzelnen gaußschen Integer, der Schritt für Schritt die nächste Basis-\(i-1\)-Ziffer bestimmt. Dieser direkte Auswerter wird nur auf sehr kleinen Rechtecken verwendet.
Für größere Rechtecke schaltet die Implementierung auf die beiden oben beschriebenen rekursiven Summen um: eine in den ursprünglichen \((a,b)\)-Koordinaten und eine in den Hilfskoordinaten \((u,v)\). Jeder Aufruf zählt den unmittelbaren Paritätsbeitrag per geschlossener Formel, transformiert die Restarbeit in ein oder zwei kleinere Rechtecke und speichert das Ergebnis in einer Memoisierungstabelle, deren Schlüssel genau aus den vier Rechteckgrenzen besteht.
Die Brute-Force-Grenze liegt bei einer Fläche von höchstens \(1024\) Gitterpunkten. Unterhalb dieser Schwelle ist direkte Auswertung billiger als weiteres Rekursionssplitten. Die C++-Implementierung bearbeitet zudem die beiden obersten transformierten Teilprobleme unabhängig voneinander parallel; Python und Java verwenden dieselbe Mathematik ohne diese zusätzliche Aufteilung. Alle Ergebnisse werden modulo \(10^9+7\) reduziert.
Für einen einzelnen gaußschen Integer benötigt die wiederholte Division durch \(i-1\) \(O(\log (|a|+|b|+1))\) Schritte, weil die Koordinatenskala schnell schrumpft. Für die Quadratsumme halbieren zwei Rekursionsebenen die Intervallgrößen ungefähr, also ist die Rekursionstiefe \(O(\log L)\). Die Gesamtzeit ist proportional zur Anzahl tatsächlich auftretender memoisierten Rechtecke, und der Speicherverbrauch ist von derselben Größenordnung. Eine saubere geschlossene Formel für diese Zustandszahl ist unhandlich, aber sie ist drastisch kleiner als die \((2L+1)^2\) Punkte eines direkten Gitterscans. Genau deshalb bleibt die Methode auch für \(L=10^{15}\) praktikabel.
Her Gaussian tamsayı \(z=a+bi\), yalnızca \(0\) ve \(1\) rakamları kullanılarak \(i-1\) tabanında yazılabilir. Bu yazımdaki \(1\) rakamlarının sayısını \(f(a,b)\) ile gösterelim ve
$$B(L)=\sum_{a=-L}^{L}\sum_{b=-L}^{L} f(a,b)$$
tanımını yapalım. Amaç \(B(10^{15}) \bmod (10^9+7)\) değerini bulmaktır. Tüm kareyi doğrudan taramak yaklaşık \(4\times 10^{30}\) kafes noktası gerektirir; dolayısıyla çözüm, \(i-1\) ile bölmenin sayıları nasıl dönüştürdüğünü kullanmak zorundadır.
Uygulamalar Gaussian tamsayıları tek tek saymaz. Bunun yerine, tek bir nokta için geçerli olan basamak çıkarma kuralı tüm dikdörtgenler üzerindeki toplamlar için iki karşılıklı yinelemeli formüle dönüştürülür.
Bir basamak adımı
$$a+bi=r_0+(i-1)(a_1+b_1 i),\qquad r_0\in\{0,1\}$$
şeklinde yazılır. \((i-1)(a_1+b_1 i)\) açılıp reel ve imajiner kısımlar eşleştirilince
$$r_0\equiv a-b \pmod{2},$$
$$a_1=\frac{b-a+r_0}{2},\qquad b_1=\frac{r_0-a-b}{2}$$
elde edilir. Böylece ilk basamak tam olarak \(r_0\) adet \(1\) katkısı yapar, geri kalan kuyruk ise daha küçük Gaussian tamsayı \(a_1+b_1 i\)'nin yazımıdır:
$$f(a,b)=r_0+f(a_1,b_1).$$
Tüm çözüm bu temel bağıntıya dayanır.
Eksenlere paralel bir dikdörtgen için
$$S([a_{\min},a_{\max}]\times[b_{\min},b_{\max}])=\sum_{a=a_{\min}}^{a_{\max}}\sum_{b=b_{\min}}^{b_{\max}} f(a,b)$$
tanımını yapalım. O zaman aranan değer
$$B(L)=S([-L,L]\times[-L,L])$$
olur. Ayrıca \(E([x_1,x_2])\) bir aralıktaki çift, \(O([x_1,x_2])\) tek sayı adedi olsun. Bu sayılar kapalı formülle bulunabildiği için algoritma parite saymak adına döngü kurmaz.
İlk basamak \(r_0=(a-b)\bmod 2\) olduğundan, \(a-b\) tek olan tüm noktalar hemen \(+1\) katkı verir. Bir dikdörtgendeki böyle noktaların sayısı
$$E([a_{\min},a_{\max}])\,O([b_{\min},b_{\max}])+O([a_{\min},a_{\max}])\,E([b_{\min},b_{\max}])$$
olur. Geriye \(f(a_1,b_1)\) terimi kalır. Bu kalanı işlerken dikdörtgen yapısını korumak için yardımcı koordinatlar tanıtılır:
$$u=a_1+b_1=r_0-a,\qquad v=a_1-b_1=b.$$
Sabit \(r_0\) için \((a,b)\mapsto (u,v)\) dönüşümü bir dikdörtgeni yine bir dikdörtgene götürür. Bu yüzden \((u,v)\)-dikdörtgenleri üzerinde tanımlı yardımcı bir \(T\) toplamı kullanılır. İlk yineleme böylece
$$\begin{aligned} S([a_{\min},a_{\max}]\times[b_{\min},b_{\max}])&=N_{\mathrm{odd}}\\ &\quad +T([-a_{\max},-a_{\min}]\times[b_{\min},b_{\max}])\\ &\quad +T([1-a_{\max},1-a_{\min}]\times[b_{\min},b_{\max}]) \end{aligned}$$
biçimini alır. Buradaki \(N_{\mathrm{odd}}\), yukarıdaki tek fark sayımıdır. İki \(T\) dikdörtgeni sırasıyla \(r_0=0\) ve \(r_0=1\) durumlarını temsil eder.
\(a_1=(u+v)/2\) ve \(b_1=(u-v)/2\) olduğundan yalnızca aynı pariteli \((u,v)\) çiftleri geçerlidir. Bu ortak pariteyi \(r_1\in\{0,1\}\) diyelim. İkinci basamak
$$r_1\equiv a_1-b_1\equiv v \pmod{2}$$
olur. Bu basamak da çıkarıldığında
$$a_2=\frac{b_1-a_1+r_1}{2}=-\frac{v-r_1}{2},\qquad b_2=\frac{r_1-a_1-b_1}{2}=-\frac{u-r_1}{2}$$
elde edilir. Dolayısıyla \((u,v)\)-dikdörtgenindeki tek-tek noktalar bir ek \(+1\) daha verir; çift-çift ve tek-tek sınıfları ise tekrar daha küçük \((a,b)\)-dikdörtgenlerine döner. \(U=[u_{\min},u_{\max}]\) ve \(V=[v_{\min},v_{\max}]\) için
$$\begin{aligned} T(U\times V)&=O(U)\,O(V)\\ &\quad +S\!\left(\left[-\left\lfloor\frac{v_{\max}}{2}\right\rfloor,-\left\lceil\frac{v_{\min}}{2}\right\rceil\right]\times\left[-\left\lfloor\frac{u_{\max}}{2}\right\rfloor,-\left\lceil\frac{u_{\min}}{2}\right\rceil\right]\right)\\ &\quad +S\!\left(\left[-\left\lfloor\frac{v_{\max}-1}{2}\right\rfloor,-\left\lceil\frac{v_{\min}-1}{2}\right\rceil\right]\times\left[-\left\lfloor\frac{u_{\max}-1}{2}\right\rfloor,-\left\lceil\frac{u_{\min}-1}{2}\right\rceil\right]\right). \end{aligned}$$
Alt sınırı üst sınırdan büyük olan her aralık boş kabul edilir ve katkısı \(0\)'dır.
\(8+0i\) noktası küçük ama öğretici bir örnektir. Basamak çıkarma işlemi tekrarlanınca
$$8\to -4-4i\to 4i\to 2-2i\to -2\to 1+i\to -i\to i\to 1\to 0$$
zinciri elde edilir. Karşılık gelen basamak dizisi
$$0,0,0,0,0,0,1,1,1$$
olduğu için
$$f(8,0)=3$$
olur. Uygulamaların kullandığı iki daha büyük kontrol değeri de
$$B(1)=20,\qquad B(2)=75$$
şeklindedir. Böylece dikdörtgen yinelemelerinin küçük karelerde brute-force toplamlarla tam uyuştuğu görülür.
İstenen değer
$$B(L)=S([-L,L]\times[-L,L])$$
şeklindedir. Her iki yineleme katmanı koordinat ölçeğini yaklaşık yarıya indirdiği için, devasa bir kare nokta nokta taranmak yerine parite sayımı, affine dönüşümler ve önbelleğe alınmış alt dikdörtgenler üzerinden çözülebilir.
C++, Python ve Java uygulamaları aynı planı izler. Önce tek bir Gaussian tamsayı için \(i-1\) tabanındaki bir sonraki basamağı tekrar tekrar çıkaran doğrudan bir değerlendirici vardır. Bu doğrudan yöntem yalnızca çok küçük dikdörtgenlerde kullanılır.
Daha büyük dikdörtgenlerde uygulama, yukarıdaki iki yinelemeli toplam yordamına geçer: biri özgün \((a,b)\) koordinatlarında, diğeri yardımcı \((u,v)\) koordinatlarında çalışır. Her çağrı, anlık parite katkısını kapalı formülle hesaplar, kalan işi bir veya iki küçük dikdörtgene dönüştürür ve sonucu dört sınır değeriyle anahtarlanan bir önbellekte saklar.
Brute-force eşiği en fazla \(1024\) kafes noktasıdır. Bu eşik altında doğrudan değerlendirme, tekrar bölmekten daha ucuzdur. C++ uygulaması ayrıca üst düzeyde ortaya çıkan iki bağımsız dönüştürülmüş alt problemi paralel işler; Python ve Java aynı matematiği ek paralelleştirme olmadan kullanır. Tüm sonuçlar \(10^9+7\) modunda tutulur.
Tek bir Gaussian tamsayı için \(i-1\) ile art arda bölme işlemi, koordinat ölçeği hızla küçüldüğü için \(O(\log (|a|+|b|+1))\) adım sürer. Kare toplamı için iki yineleme katmanı aralık boyutlarını yaklaşık ikiye böler; bu yüzden derinlik \(O(\log L)\)'dir. Toplam süre, gerçekten oluşan farklı önbellekli dikdörtgen sayısı ile orantılıdır ve bellek kullanımı da aynı mertebededir. Bu durum sayısının tam kapalı formu karmaşıktır; ancak \((2L+1)^2\) adet doğrudan kafes noktasından çok daha küçüktür. Yöntemin \(L=10^{15}\) için de uygulanabilir kalmasının nedeni budur.
Todo entero gaussiano \(z=a+bi\) tiene una representación en base \(i-1\) que usa solo los dígitos \(0\) y \(1\). Sea \(f(a,b)\) el número de dígitos iguales a \(1\) en esa expansión, y definamos
$$B(L)=\sum_{a=-L}^{L}\sum_{b=-L}^{L} f(a,b).$$
La meta es calcular \(B(10^{15}) \bmod (10^9+7)\). Un barrido directo del cuadrado completo tocaría aproximadamente \(4\times 10^{30}\) puntos de la red, así que la solución debe aprovechar la estructura aritmética de dividir por \(i-1\).
Las implementaciones no recorren los enteros gaussianos uno por uno. En su lugar, convierten la regla de extracción de dígitos para un solo punto en dos fórmulas recursivas mutuas para sumas sobre rectángulos enteros.
Un paso de dígitos se escribe como
$$a+bi=r_0+(i-1)(a_1+b_1 i),\qquad r_0\in\{0,1\}.$$
Al expandir y comparar partes real e imaginaria se obtiene
$$r_0\equiv a-b \pmod{2},$$
$$a_1=\frac{b-a+r_0}{2},\qquad b_1=\frac{r_0-a-b}{2}.$$
Por tanto, ese dígito aporta exactamente \(r_0\) unos y el resto de la representación es la del entero gaussiano más pequeño \(a_1+b_1 i\):
$$f(a,b)=r_0+f(a_1,b_1).$$
Esa es la recurrencia fundamental de toda la solución.
Para un rectángulo alineado con los ejes definimos
$$S([a_{\min},a_{\max}]\times[b_{\min},b_{\max}])=\sum_{a=a_{\min}}^{a_{\max}}\sum_{b=b_{\min}}^{b_{\max}} f(a,b).$$
Entonces la cantidad pedida es
$$B(L)=S([-L,L]\times[-L,L]).$$
También escribimos \(E([x_1,x_2])\) para el número de enteros pares en un intervalo y \(O([x_1,x_2])\) para el número de impares. Esas cantidades se calculan mediante fórmulas cerradas, de modo que el algoritmo no itera solo para contar paridades.
El primer dígito es \(r_0=(a-b)\bmod 2\). Por ello, todos los puntos con \(a-b\) impar aportan inmediatamente \(+1\). Su número dentro de un rectángulo es
$$E([a_{\min},a_{\max}])\,O([b_{\min},b_{\max}])+O([a_{\min},a_{\max}])\,E([b_{\min},b_{\max}]).$$
Queda el término residual \(f(a_1,b_1)\). Para que la imagen de un rectángulo siga siendo otro rectángulo, introducimos coordenadas auxiliares:
$$u=a_1+b_1=r_0-a,\qquad v=a_1-b_1=b.$$
Con \(r_0\) fijo, la transformación \((a,b)\mapsto(u,v)\) envía un rectángulo a otro rectángulo. Por eso se define una suma auxiliar \(T\) sobre rectángulos en \((u,v)\). La primera recurrencia queda
$$\begin{aligned} S([a_{\min},a_{\max}]\times[b_{\min},b_{\max}])&=N_{\mathrm{odd}}\\ &\quad +T([-a_{\max},-a_{\min}]\times[b_{\min},b_{\max}])\\ &\quad +T([1-a_{\max},1-a_{\min}]\times[b_{\min},b_{\max}]), \end{aligned}$$
donde \(N_{\mathrm{odd}}\) es el conteo de paridad anterior. Los dos rectángulos de \(T\) corresponden a los casos \(r_0=0\) y \(r_0=1\).
Como \(a_1=(u+v)/2\) y \(b_1=(u-v)/2\), solo pueden aparecer pares con la misma paridad. Llamemos \(r_1\in\{0,1\}\) a esa paridad común. Entonces el segundo dígito es
$$r_1\equiv a_1-b_1\equiv v \pmod{2}.$$
Al quitar también ese dígito se obtiene
$$a_2=\frac{b_1-a_1+r_1}{2}=-\frac{v-r_1}{2},\qquad b_2=\frac{r_1-a_1-b_1}{2}=-\frac{u-r_1}{2}.$$
Así, los puntos impar-impar de un rectángulo en \((u,v)\) aportan un \(+1\) adicional, mientras que las clases par-par e impar-impar regresan a rectángulos más pequeños en las coordenadas originales \((a,b)\). Si \(U=[u_{\min},u_{\max}]\) y \(V=[v_{\min},v_{\max}]\), entonces
$$\begin{aligned} T(U\times V)&=O(U)\,O(V)\\ &\quad +S\!\left(\left[-\left\lfloor\frac{v_{\max}}{2}\right\rfloor,-\left\lceil\frac{v_{\min}}{2}\right\rceil\right]\times\left[-\left\lfloor\frac{u_{\max}}{2}\right\rfloor,-\left\lceil\frac{u_{\min}}{2}\right\rceil\right]\right)\\ &\quad +S\!\left(\left[-\left\lfloor\frac{v_{\max}-1}{2}\right\rfloor,-\left\lceil\frac{v_{\min}-1}{2}\right\rceil\right]\times\left[-\left\lfloor\frac{u_{\max}-1}{2}\right\rfloor,-\left\lceil\frac{u_{\min}-1}{2}\right\rceil\right]\right). \end{aligned}$$
Cualquier intervalo cuyo extremo inferior supere al superior se considera vacío y aporta \(0\).
El punto \(8+0i\) es pequeño pero ilustrativo. La extracción repetida de dígitos produce
$$8\to -4-4i\to 4i\to 2-2i\to -2\to 1+i\to -i\to i\to 1\to 0.$$
La secuencia de dígitos correspondiente es
$$0,0,0,0,0,0,1,1,1,$$
de modo que
$$f(8,0)=3.$$
Otros dos puntos de control usados por las implementaciones son
$$B(1)=20,\qquad B(2)=75.$$
Eso verifica que las recurrencias rectangulares reproducen exactamente las sumas por fuerza bruta en cuadrados pequeños.
La cantidad buscada es
$$B(L)=S([-L,L]\times[-L,L]).$$
Cada par de niveles recursivos reduce aproximadamente a la mitad la escala de las coordenadas. Por eso un cuadrado gigantesco puede resolverse mediante conteos de paridad, transformaciones afines y memoización de subrectángulos, en lugar de enumerar punto por punto.
Las implementaciones en C++, Python y Java siguen exactamente el mismo esquema. Primero disponen de un evaluador directo para un solo entero gaussiano, que extrae repetidamente el siguiente dígito en base \(i-1\). Ese evaluador directo solo se usa en rectángulos muy pequeños.
Para rectángulos grandes, la implementación cambia a las dos rutinas recursivas descritas arriba: una opera en las coordenadas originales \((a,b)\) y la otra en las coordenadas auxiliares \((u,v)\). Cada llamada cuenta la contribución inmediata de la paridad mediante fórmulas cerradas, transforma el trabajo restante en uno o dos rectángulos más pequeños y guarda el resultado en una tabla de memoización indexada por las cuatro fronteras del rectángulo.
El umbral de fuerza bruta es un área de como máximo \(1024\) puntos de la red. Por debajo de ese tamaño, evaluar directamente es más barato que seguir dividiendo. La implementación en C++ además resuelve de manera independiente los dos subproblemas transformados del nivel superior, porque son disjuntos desde el punto de vista matemático; las versiones en Python y Java usan la misma recurrencia sin esa división paralela adicional. Todos los resultados se reducen módulo \(10^9+7\).
Para un solo entero gaussiano, la división repetida por \(i-1\) requiere \(O(\log (|a|+|b|+1))\) pasos, porque la escala de las coordenadas disminuye rápidamente. Para la suma sobre el cuadrado, cada dos niveles recursivos reducen aproximadamente a la mitad el tamaño de los intervalos, así que la profundidad es \(O(\log L)\). El tiempo total es proporcional al número de rectángulos distintos que realmente se memorizan, y la memoria tiene el mismo orden. Es incómodo escribir una fórmula cerrada exacta para ese número de estados, pero es muchísimo menor que los \((2L+1)^2\) puntos de un barrido directo. Esa es la razón de que el método siga siendo viable incluso para \(L=10^{15}\).
每个高斯整数 \(z=a+bi\) 都可以用仅含数字 \(0\) 与 \(1\) 的 \(i-1\) 进制表示。记 \(f(a,b)\) 为该表示中数字 \(1\) 的个数，并定义
$$B(L)=\sum_{a=-L}^{L}\sum_{b=-L}^{L} f(a,b).$$
题目要求计算 \(B(10^{15}) \bmod (10^9+7)\)。如果直接遍历整个正方形区域，大约需要处理 \(4\times 10^{30}\) 个格点，这是完全不可行的。因此必须利用“除以 \(i-1\)”时高斯整数坐标所满足的递推结构。
三份实现的核心都不是逐点枚举，而是先把单个高斯整数的“取下一位”公式写出来，再把它推广成对整块矩形区域求和的递推。这样才能把超大正方形分解成大量更小的重复子问题。
把一个数字提取步骤写成
$$a+bi=r_0+(i-1)(a_1+b_1 i),\qquad r_0\in\{0,1\}.$$
将 \((i-1)(a_1+b_1 i)\) 展开，并分别比较实部与虚部，可以得到
$$r_0\equiv a-b \pmod{2},$$
$$a_1=\frac{b-a+r_0}{2},\qquad b_1=\frac{r_0-a-b}{2}.$$
这说明第一位数字完全由 \(a-b\) 的奇偶性决定。若 \(a-b\) 为偶数，则第一位是 \(0\)；若为奇数，则第一位是 \(1\)。因此
$$f(a,b)=r_0+f(a_1,b_1).$$
也就是说，当前这一步只贡献 \(r_0\) 个 \(1\)，剩余部分就是较小高斯整数 \(a_1+b_1 i\) 的同类问题。
对任意轴对齐矩形，定义
$$S([a_{\min},a_{\max}]\times[b_{\min},b_{\max}])=\sum_{a=a_{\min}}^{a_{\max}}\sum_{b=b_{\min}}^{b_{\max}} f(a,b).$$
那么题目要求的答案就是
$$B(L)=S([-L,L]\times[-L,L]).$$
为了后面的奇偶分类，我们再记 \(E([x_1,x_2])\) 为区间中的偶整数个数，\(O([x_1,x_2])\) 为区间中的奇整数个数。这两个数量都可以用整除公式直接算出，因此程序无需为了数奇偶而遍历区间中的每个整数。
第一位数字是 \(r_0=(a-b)\bmod 2\)，所以凡是满足 \(a-b\) 为奇数的点，都会立刻贡献一个 \(1\)。这类点在矩形中的总数为
$$E([a_{\min},a_{\max}])\,O([b_{\min},b_{\max}])+O([a_{\min},a_{\max}])\,E([b_{\min},b_{\max}]).$$
剩下还没算进去的是尾部项 \(f(a_1,b_1)\)。如果直接把矩形在 \((a_1,b_1)\) 平面里追踪，其边界会变得不方便处理。实现里的关键技巧是引入辅助坐标
$$u=a_1+b_1=r_0-a,\qquad v=a_1-b_1=b.$$
这样做的好处是：当 \(r_0\) 固定时，原来的一个轴对齐矩形在 \((u,v)\) 平面中仍然是轴对齐矩形。因此可以定义一个辅助求和函数 \(T\)，专门处理 \((u,v)\)-矩形上的剩余贡献。于是第一次递推变成
$$\begin{aligned} S([a_{\min},a_{\max}]\times[b_{\min},b_{\max}])&=N_{\mathrm{odd}}\\ &\quad +T([-a_{\max},-a_{\min}]\times[b_{\min},b_{\max}])\\ &\quad +T([1-a_{\max},1-a_{\min}]\times[b_{\min},b_{\max}]), \end{aligned}$$
其中 \(N_{\mathrm{odd}}\) 就是上面的奇差点数。两个 \(T\) 矩形分别对应第一位数字 \(r_0=0\) 与 \(r_0=1\) 两种情况。
由
$$a_1=\frac{u+v}{2},\qquad b_1=\frac{u-v}{2}$$
可知，只有当 \(u\) 与 \(v\) 同奇同偶时，\(a_1,b_1\) 才是整数。所以在辅助矩形里，异奇偶的点根本不会产生有效状态。设这共同的奇偶性为 \(r_1\in\{0,1\}\)，那么第二位数字满足
$$r_1\equiv a_1-b_1\equiv v \pmod{2}.$$
再去掉这一位之后，新的原坐标变为
$$a_2=\frac{b_1-a_1+r_1}{2}=-\frac{v-r_1}{2},\qquad b_2=\frac{r_1-a_1-b_1}{2}=-\frac{u-r_1}{2}.$$
这就带来一个非常整齐的结构：
奇奇点会额外贡献一个 \(1\)，因为这时第二位数字 \(r_1=1\)。
偶偶点和奇奇点在去掉第二位之后，又都回到了更小的原始 \((a,b)\)-矩形上。
若 \(U=[u_{\min},u_{\max}]\)、\(V=[v_{\min},v_{\max}]\)，则辅助递推可写为
$$\begin{aligned} T(U\times V)&=O(U)\,O(V)\\ &\quad +S\!\left(\left[-\left\lfloor\frac{v_{\max}}{2}\right\rfloor,-\left\lceil\frac{v_{\min}}{2}\right\rceil\right]\times\left[-\left\lfloor\frac{u_{\max}}{2}\right\rfloor,-\left\lceil\frac{u_{\min}}{2}\right\rceil\right]\right)\\ &\quad +S\!\left(\left[-\left\lfloor\frac{v_{\max}-1}{2}\right\rfloor,-\left\lceil\frac{v_{\min}-1}{2}\right\rceil\right]\times\left[-\left\lfloor\frac{u_{\max}-1}{2}\right\rfloor,-\left\lceil\frac{u_{\min}-1}{2}\right\rceil\right]\right). \end{aligned}$$
如果某个区间的左端点大于右端点，就说明该子矩形为空，贡献自然记为 \(0\)。
\(8+0i\) 是一个很适合手算的例子。不断应用步骤 1 的公式，会得到
$$8\to -4-4i\to 4i\to 2-2i\to -2\to 1+i\to -i\to i\to 1\to 0.$$
对应的数字序列为
$$0,0,0,0,0,0,1,1,1,$$
因此
$$f(8,0)=3.$$
实现中还使用了两个更大的校验值：
$$B(1)=20,\qquad B(2)=75.$$
这两个值可以用小范围暴力直接验证，因此也说明矩形递推在小规模上与直接求和完全一致。
最后要求的就是
$$B(L)=S([-L,L]\times[-L,L]).$$
观察步骤 3 与步骤 4 可以发现：每经过两层递推，坐标范围大致缩小一半。于是，大正方形的求和问题会很快分裂成许多更小的矩形，而且这些小矩形大量重复。正因如此，记忆化缓存才会非常有效。
C++、Python 和 Java 三份实现使用的是同一套思路。首先，它们都保留了一个针对单个高斯整数的直接计算器，反复执行“确定当前位并更新到下一状态”的过程。这个直接计算器只在矩形非常小时使用。
一旦矩形变大，程序就切换到上面推导出的两个递归求和器：一个处理原坐标 \((a,b)\)，另一个处理辅助坐标 \((u,v)\)。每次调用都先用闭式公式计算本层的奇偶贡献，再把剩余部分转换为一个或两个更小的矩形，最后把结果按四个边界值缓存起来，这样同一子问题以后就不必重复计算。
暴力阈值是面积不超过 \(1024\) 个格点的矩形。低于这个阈值时，直接枚举反而更划算。C++ 实现还把最顶层产生的两个互不依赖的辅助子问题拆开并行计算；Python 与 Java 则保持同样的数学递推，但不额外做这一步并行化。所有输出最终都对 \(10^9+7\) 取模。
对于单个高斯整数，连续除以 \(i-1\) 的步数是 \(O(\log (|a|+|b|+1))\)，因为坐标尺度会很快下降。对于整个正方形求和，每两层递推大约把区间长度减半，因此递归深度是 \(O(\log L)\)。总时间与实际出现的不同缓存矩形数量成正比，内存复杂度也是同一量级。这个状态数很难写成漂亮的精确闭式，但它远小于直接扫描需要处理的 \((2L+1)^2\) 个格点，这正是算法能够应付 \(L=10^{15}\) 的根本原因。
Каждое гауссово целое \(z=a+bi\) можно записать в системе счисления с основанием \(i-1\), используя только цифры \(0\) и \(1\). Пусть \(f(a,b)\) обозначает количество цифр \(1\) в такой записи, а
$$B(L)=\sum_{a=-L}^{L}\sum_{b=-L}^{L} f(a,b).$$
Требуется вычислить \(B(10^{15}) \bmod (10^9+7)\). Прямой перебор всего квадрата означал бы примерно \(4\times 10^{30}\) точек решётки, поэтому решение должно использовать специальную структуру деления на \(i-1\).
Реальные программы не перебирают гауссовы целые по одному. Вместо этого правило извлечения очередной цифры для одной точки превращается в две взаимные рекурсии для сумм по целым прямоугольникам.
Один шаг записи имеет вид
$$a+bi=r_0+(i-1)(a_1+b_1 i),\qquad r_0\in\{0,1\}.$$
После раскрытия скобок и сравнения действительной и мнимой частей получаем
$$r_0\equiv a-b \pmod{2},$$
$$a_1=\frac{b-a+r_0}{2},\qquad b_1=\frac{r_0-a-b}{2}.$$
Следовательно, текущая цифра вносит ровно \(r_0\) единиц, а хвост записи совпадает с записью меньшего гауссова целого \(a_1+b_1 i\):
$$f(a,b)=r_0+f(a_1,b_1).$$
Это базовая рекуррентная формула всей конструкции.
Для прямоугольника, параллельного осям, введём
$$S([a_{\min},a_{\max}]\times[b_{\min},b_{\max}])=\sum_{a=a_{\min}}^{a_{\max}}\sum_{b=b_{\min}}^{b_{\max}} f(a,b).$$
Тогда искомая величина равна
$$B(L)=S([-L,L]\times[-L,L]).$$
Также удобно обозначить через \(E([x_1,x_2])\) число чётных целых в интервале, а через \(O([x_1,x_2])\) число нечётных. Эти количества находятся по явным формулам, поэтому код не тратит время на отдельный перебор ради подсчёта паритетов.
Первая цифра равна \(r_0=(a-b)\bmod 2\). Значит, все точки с нечётным \(a-b\) сразу дают вклад \(+1\). Их число в прямоугольнике равно
$$E([a_{\min},a_{\max}])\,O([b_{\min},b_{\max}])+O([a_{\min},a_{\max}])\,E([b_{\min},b_{\max}]).$$
Остаётся хвост \(f(a_1,b_1)\). Чтобы образ прямоугольника после преобразования снова был прямоугольником, вводятся вспомогательные координаты
$$u=a_1+b_1=r_0-a,\qquad v=a_1-b_1=b.$$
При фиксированном \(r_0\) отображение \((a,b)\mapsto(u,v)\) переводит осевой прямоугольник в осевой прямоугольник. Поэтому вводится вспомогательная сумма \(T\) по прямоугольникам в координатах \((u,v)\). Первая рекурсия принимает вид
$$\begin{aligned} S([a_{\min},a_{\max}]\times[b_{\min},b_{\max}])&=N_{\mathrm{odd}}\\ &\quad +T([-a_{\max},-a_{\min}]\times[b_{\min},b_{\max}])\\ &\quad +T([1-a_{\max},1-a_{\min}]\times[b_{\min},b_{\max}]), \end{aligned}$$
где \(N_{\mathrm{odd}}\) и есть только что подсчитанный вклад нечётной разности. Два прямоугольника для \(T\) соответствуют случаям \(r_0=0\) и \(r_0=1\).
Поскольку
$$a_1=\frac{u+v}{2},\qquad b_1=\frac{u-v}{2},$$
возможны только пары \((u,v)\) одинаковой чётности. Обозначим эту общую чётность через \(r_1\in\{0,1\}\). Тогда вторая цифра равна
$$r_1\equiv a_1-b_1\equiv v \pmod{2}.$$
После удаления и этой цифры новые исходные координаты становятся такими:
$$a_2=\frac{b_1-a_1+r_1}{2}=-\frac{v-r_1}{2},\qquad b_2=\frac{r_1-a_1-b_1}{2}=-\frac{u-r_1}{2}.$$
Отсюда следует удобная структура: точки нечётность-нечётность дают ещё один непосредственный вклад \(+1\), а классы чётность-чётность и нечётность-нечётность переводятся обратно в меньшие прямоугольники в исходных координатах \((a,b)\). Если \(U=[u_{\min},u_{\max}]\) и \(V=[v_{\min},v_{\max}]\), то
$$\begin{aligned} T(U\times V)&=O(U)\,O(V)\\ &\quad +S\!\left(\left[-\left\lfloor\frac{v_{\max}}{2}\right\rfloor,-\left\lceil\frac{v_{\min}}{2}\right\rceil\right]\times\left[-\left\lfloor\frac{u_{\max}}{2}\right\rfloor,-\left\lceil\frac{u_{\min}}{2}\right\rceil\right]\right)\\ &\quad +S\!\left(\left[-\left\lfloor\frac{v_{\max}-1}{2}\right\rfloor,-\left\lceil\frac{v_{\min}-1}{2}\right\rceil\right]\times\left[-\left\lfloor\frac{u_{\max}-1}{2}\right\rfloor,-\left\lceil\frac{u_{\min}-1}{2}\right\rceil\right]\right). \end{aligned}$$
Если нижняя граница интервала оказывается больше верхней, это просто пустой прямоугольник с нулевым вкладом.
Точка \(8+0i\) мала, но хорошо показывает механику метода. Повторяя шаг извлечения цифры, получаем
$$8\to -4-4i\to 4i\to 2-2i\to -2\to 1+i\to -i\to i\to 1\to 0.$$
Соответствующая последовательность цифр такова:
$$0,0,0,0,0,0,1,1,1,$$
поэтому
$$f(8,0)=3.$$
Ещё две контрольные величины, которые используют реализации:
$$B(1)=20,\qquad B(2)=75.$$
Они подтверждают, что рекурсии по прямоугольникам точно совпадают с прямым перебором на малых квадратах.
Нужно вычислить
$$B(L)=S([-L,L]\times[-L,L]).$$
Каждые два уровня рекурсии уменьшают масштаб координат примерно вдвое. Поэтому огромный квадрат удаётся обработать через подсчёт паритетов, аффинные преобразования и мемоизацию подпрямоугольников, а не через явный просмотр каждой точки.
Реализации на C++, Python и Java устроены одинаково. Сначала в них есть прямой вычислитель для одного гауссова целого, который многократно извлекает следующую цифру в базе \(i-1\). Этот прямой путь используется только тогда, когда прямоугольник совсем мал.
Для больших прямоугольников программа переключается на две рекурсивные суммы, выведенные выше: одна работает в исходных координатах \((a,b)\), другая во вспомогательных \((u,v)\). Каждый вызов по явным формулам считает немедленный вклад чётности, преобразует остаток задачи в один или два меньших прямоугольника и сохраняет результат в таблице мемоизации, ключом в которой служат четыре граничных значения.
Порог прямого перебора равен площади не более \(1024\) узлов решётки. Ниже этого порога прямой расчёт дешевле, чем дальнейшее рекурсивное разбиение. Реализация на C++ дополнительно запускает два верхних независимых подзадачи отдельно, поскольку они математически не пересекаются; версии на Python и Java используют ту же рекурсию без этого дополнительного распараллеливания. Все ответы считаются по модулю \(10^9+7\).
Для одного гауссова целого повторное деление на \(i-1\) требует \(O(\log (|a|+|b|+1))\) шагов, потому что масштаб координат быстро уменьшается. Для суммы по квадрату каждые два рекурсивных уровня примерно вдвое уменьшают размеры интервалов, так что глубина рекурсии равна \(O(\log L)\). Полное время работы пропорционально числу различных прямоугольников, реально попавших в мемоизацию, а память имеет тот же порядок. Точную красивую формулу для этого числа состояний написать трудно, но оно на порядки меньше, чем \((2L+1)^2\) точек полного перебора. Именно поэтому метод остаётся практичным даже при \(L=10^{15}\).
كل عدد غاوسي \(z=a+bi\) يملك تمثيلًا بالأساس \(i-1\) يستخدم فقط الرقمين \(0\) و\(1\). لنعرّف \(f(a,b)\) بأنه عدد الخانات التي تساوي \(1\) في هذا التمثيل، ثم نضع
$$B(L)=\sum_{a=-L}^{L}\sum_{b=-L}^{L} f(a,b).$$
المطلوب هو حساب \(B(10^{15}) \bmod (10^9+7)\). المسح المباشر للمربع كله يعني التعامل مع نحو \(4\times 10^{30}\) نقطة شبكية، وهذا غير ممكن عمليًا، لذلك يجب استغلال البنية الخاصة للقسمة على \(i-1\).
التنفيذات لا تعدّ الأعداد الغاوسية نقطةً نقطة. بدلًا من ذلك، تأخذ قاعدة استخراج الخانة التالية لعدد واحد وتحولها إلى علاقتين عوديتين متبادلتين لمجاميع على مستطيلات كاملة.
نكتب خطوة واحدة على الصورة
$$a+bi=r_0+(i-1)(a_1+b_1 i),\qquad r_0\in\{0,1\}.$$
وبعد التوسيع ومطابقة الجزأين الحقيقي والتخيلي نحصل على
$$r_0\equiv a-b \pmod{2},$$
$$a_1=\frac{b-a+r_0}{2},\qquad b_1=\frac{r_0-a-b}{2}.$$
إذن الخانة الحالية تضيف بالضبط \(r_0\) إلى عدد الآحاد، وما يبقى هو تمثيل العدد الغاوسي الأصغر \(a_1+b_1 i\):
$$f(a,b)=r_0+f(a_1,b_1).$$
هذه هي العلاقة الأساسية التي تعتمد عليها الحلول كلها.
لأي مستطيل موازي للمحورين نعرّف
$$S([a_{\min},a_{\max}]\times[b_{\min},b_{\max}])=\sum_{a=a_{\min}}^{a_{\max}}\sum_{b=b_{\min}}^{b_{\max}} f(a,b).$$
وعندئذ تصبح القيمة المطلوبة
$$B(L)=S([-L,L]\times[-L,L]).$$
كما نرمز بـ \(E([x_1,x_2])\) إلى عدد الأعداد الزوجية في مجال صحيح، وبـ \(O([x_1,x_2])\) إلى عدد الأعداد الفردية فيه. هذه الأعداد تُحسب بصيغ مغلقة، لذلك لا تحتاج الخوارزمية إلى المرور على المجال عنصرًا عنصرًا لمجرد معرفة الزوجية.
الخانة الأولى هي \(r_0=(a-b)\bmod 2\). ولذلك فإن كل نقطة تحقق أن \(a-b\) فردي تعطي مساهمة مباشرة مقدارها \(+1\). وعدد هذه النقاط داخل مستطيل يساوي
$$E([a_{\min},a_{\max}])\,O([b_{\min},b_{\max}])+O([a_{\min},a_{\max}])\,E([b_{\min},b_{\max}]).$$
يبقى بعد ذلك الحد \(f(a_1,b_1)\). وللحفاظ على كون صورة المستطيل مستطيلًا محورياً، ندخل إحداثيين مساعدين:
$$u=a_1+b_1=r_0-a,\qquad v=a_1-b_1=b.$$
عند تثبيت \(r_0\)، فإن التحويل \((a,b)\mapsto(u,v)\) يرسل مستطيلًا إلى مستطيل. لذلك نعرّف مجموعًا مساعدًا \(T\) على مستطيلات \((u,v)\). وتصبح العودية الأولى
$$\begin{aligned} S([a_{\min},a_{\max}]\times[b_{\min},b_{\max}])&=N_{\mathrm{odd}}\\ &\quad +T([-a_{\max},-a_{\min}]\times[b_{\min},b_{\max}])\\ &\quad +T([1-a_{\max},1-a_{\min}]\times[b_{\min},b_{\max}]), \end{aligned}$$
حيث إن \(N_{\mathrm{odd}}\) هو حدّ العدّ الفردي السابق. والمستطيلان في \(T\) يمثلان الحالتين \(r_0=0\) و\(r_0=1\).
من
$$a_1=\frac{u+v}{2},\qquad b_1=\frac{u-v}{2}$$
نستنتج أن الأزواج الممكنة فقط هي التي لها الزوجية نفسها. ولنسم هذه الزوجية المشتركة \(r_1\in\{0,1\}\). عندئذ تكون الخانة الثانية
$$r_1\equiv a_1-b_1\equiv v \pmod{2}.$$
وبعد حذف هذه الخانة أيضًا نحصل على
$$a_2=\frac{b_1-a_1+r_1}{2}=-\frac{v-r_1}{2},\qquad b_2=\frac{r_1-a_1-b_1}{2}=-\frac{u-r_1}{2}.$$
وهذا يعني أن نقاط فردي-فردي داخل مستطيل \((u,v)\) تضيف \(+1\) آخر فورًا، بينما تعود حالتا زوجي-زوجي وفردي-فردي بعد التقسيم إلى مستطيلات أصغر في الإحداثيات الأصلية \((a,b)\). فإذا كان \(U=[u_{\min},u_{\max}]\) و\(V=[v_{\min},v_{\max}]\)، فالعلاقة هي
$$\begin{aligned} T(U\times V)&=O(U)\,O(V)\\ &\quad +S\!\left(\left[-\left\lfloor\frac{v_{\max}}{2}\right\rfloor,-\left\lceil\frac{v_{\min}}{2}\right\rceil\right]\times\left[-\left\lfloor\frac{u_{\max}}{2}\right\rfloor,-\left\lceil\frac{u_{\min}}{2}\right\rceil\right]\right)\\ &\quad +S\!\left(\left[-\left\lfloor\frac{v_{\max}-1}{2}\right\rfloor,-\left\lceil\frac{v_{\min}-1}{2}\right\rceil\right]\times\left[-\left\lfloor\frac{u_{\max}-1}{2}\right\rfloor,-\left\lceil\frac{u_{\min}-1}{2}\right\rceil\right]\right). \end{aligned}$$
وأي مجال يكون فيه الحد الأدنى أكبر من الحد الأعلى يُعدّ فارغًا ومساهمته صفر.
النقطة \(8+0i\) صغيرة لكنها توضّح الفكرة جيدًا. إذا كررنا استخراج الخانة نحصل على السلسلة
$$8\to -4-4i\to 4i\to 2-2i\to -2\to 1+i\to -i\to i\to 1\to 0.$$
وتسلسل الخانات الموافق هو
$$0,0,0,0,0,0,1,1,1,$$
ومن ثم
$$f(8,0)=3.$$
وهناك قيمتان إضافيتان للتحقق تستخدمهما التنفيذات:
$$B(1)=20,\qquad B(2)=75.$$
وهذا يؤكد أن عوديات المستطيلات تعيد بالضبط القيم نفسها التي يعطيها العدّ المباشر على المربعات الصغيرة.
القيمة النهائية التي نريدها هي
$$B(L)=S([-L,L]\times[-L,L]).$$
كل مستويين عوديين يخفضان مقياس الإحداثيات تقريبًا إلى النصف. ولهذا يمكن معالجة مربع هائل جدًا عبر عدّ الزوجيات، والتحويلات affine، وتخزين المستطيلات الفرعية، بدلًا من تعداد كل نقطة على حدة.
تنفيذات C++ وPython وJava تتبع الخطة نفسها. أولًا يوجد مقيم مباشر لعدد غاوسي واحد، يكرر عملية استخراج الخانة التالية في الأساس \(i-1\). هذا المسار المباشر يُستخدم فقط عندما يكون المستطيل صغيرًا جدًا.
أما المستطيلات الأكبر، فتُعالج عبر الروتينين العوديين المشتقين أعلاه: روتين في الإحداثيات الأصلية \((a,b)\)، وروتين في الإحداثيات المساعدة \((u,v)\). كل استدعاء يحسب المساهمة الفورية للزوجية بصيغ مغلقة، ثم يحول ما تبقى إلى مستطيل واحد أو مستطيلين أصغر، ثم يخزن النتيجة باستخدام حدود المستطيل الأربعة كمفتاح في جدول تذكّر.
حدّ العدّ المباشر هو مساحة لا تتجاوز \(1024\) نقطة شبكية. تحت هذا الحد يكون التقييم المباشر أرخص من المزيد من التقسيم. تنفيذ C++ يحسب أيضًا المشكلتين العلويتين المستقلتين بصورة منفصلة لأنهما غير متداخلتين رياضيًا، بينما يستخدم Python وJava العودية نفسها من دون هذا التقسيم المتوازي الإضافي. كل القيم تُختزل بترديد \(10^9+7\).
بالنسبة إلى عدد غاوسي واحد، فإن القسمة المتكررة على \(i-1\) تحتاج إلى \(O(\log (|a|+|b|+1))\) خطوة لأن مقياس الإحداثيات يتناقص بسرعة. أما مجموع المربع، فكل مستويين من العودية يقلصان أطوال المجالات بنحو النصف، لذا يكون عمق العودية \(O(\log L)\). وزمن التنفيذ الكلي يتناسب مع عدد المستطيلات المختلفة التي تظهر فعليًا في التذكّر، واستهلاك الذاكرة من المرتبة نفسها. من الصعب كتابة صيغة مغلقة أنيقة لهذا العدد من الحالات، لكنه أصغر بكثير من \((2L+1)^2\) نقطة التي يتطلبها المسح المباشر. ولهذا تبقى الطريقة عملية حتى عند \(L=10^{15}\).