Detailed mathematical approach

Problem Summary

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\).

Mathematical Approach

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.

Step 1: Extract One Base-\(i-1\) Digit

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.

Step 2: Turn Pointwise Recurrence into Rectangle Sums

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.

Step 3: First Parity Split in the \((a,b)\)-Plane

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\).

Step 4: Second Parity Split in the Auxiliary \((u,v)\)-Plane

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\).

Step 5: Worked Example

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.

Step 6: Final Recursive Formula for the Target Sum

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.

How the Code Works

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\).

Complexity Analysis

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}\).

Footnotes and References

1. Problem page: https://projecteuler.net/problem=508
2. Gaussian integers: Wikipedia — Gaussian integer
3. Complex-base numeral systems: Wikipedia — Complex-base system
4. Floor and ceiling functions: Wikipedia — Floor and ceiling functions
5. Memoization: Wikipedia — Memoization