Problem 508: Integers in Base $i-1$

View on Project Euler

Project Euler Problem 508 Solution

EulerSolve provides an optimized solution for Project Euler Problem 508, Integers in Base $i-1$, with C++, Python, Java, and a step-by-step mathematical explanation.

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

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

Problem 508 source code

C++

#include <algorithm>
#include <cstdint>
#include <cstdlib>
#include <future>
#include <iostream>
#include <thread>
#include <unordered_map>

using namespace std;

namespace {

constexpr int64_t MOD = 1000000007LL;
constexpr int64_t BRUTE_LIMIT = 1024;

struct Key {
    int64_t a = 0;
    int64_t b = 0;
    int64_t c = 0;
    int64_t d = 0;

    bool operator==(const Key& other) const {
        return a == other.a && b == other.b && c == other.c && d == other.d;
    }
};

struct KeyHash {
    size_t operator()(const Key& k) const {
        size_t h = std::hash<int64_t>{}(k.a);
        h ^= std::hash<int64_t>{}(k.b) + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2);
        h ^= std::hash<int64_t>{}(k.c) + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2);
        h ^= std::hash<int64_t>{}(k.d) + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2);
        return h;
    }
};

int64_t floor_div(int64_t a, int64_t b) {
    int64_t q = a / b;
    int64_t r = a % b;
    if (r != 0 && a < 0) --q;
    return q;
}

int64_t ceil_div(int64_t a, int64_t b) {
    return -floor_div(-a, b);
}

void count_even_odd(int64_t l, int64_t r, int64_t& even, int64_t& odd) {
    if (l > r) {
        even = 0;
        odd = 0;
        return;
    }
    int64_t first_even = (l % 2 == 0) ? l : l + 1;
    int64_t last_even = (r % 2 == 0) ? r : r - 1;
    if (first_even > last_even) {
        even = 0;
    } else {
        even = (last_even - first_even) / 2 + 1;
    }
    int64_t total = r - l + 1;
    odd = total - even;
}

int64_t mod_add(int64_t a, int64_t b) {
    a += b;
    if (a >= MOD) a -= MOD;
    return a;
}

int64_t ones_in_base(int64_t a, int64_t b) {
    int64_t count = 0;
    while (a != 0 || b != 0) {
        int64_t r = (a - b) & 1LL;
        count += r;
        int64_t na = (b - a + r) / 2;
        int64_t nb = (r - a - b) / 2;
        a = na;
        b = nb;
    }
    return count;
}

class Solver {
  public:
    int64_t compute_square(int64_t L) { return sum_b(-L, L, -L, L); }
    int64_t compute_a(int64_t umin, int64_t umax, int64_t vmin, int64_t vmax) {
        return sum_a(umin, umax, vmin, vmax);
    }

  private:
    unordered_map<Key, int64_t, KeyHash> memo_b_;
    unordered_map<Key, int64_t, KeyHash> memo_a_;

    int64_t brute_b(int64_t amin, int64_t amax, int64_t bmin, int64_t bmax) {
        int64_t total = 0;
        for (int64_t a = amin; a <= amax; ++a) {
            for (int64_t b = bmin; b <= bmax; ++b) {
                total += ones_in_base(a, b);
            }
        }
        return total % MOD;
    }

    int64_t sum_b(int64_t amin, int64_t amax, int64_t bmin, int64_t bmax) {
        if (amin > amax || bmin > bmax) return 0;
        Key key{amin, amax, bmin, bmax};
        auto it = memo_b_.find(key);
        if (it != memo_b_.end()) return it->second;

        __int128 area = static_cast<__int128>(amax - amin + 1) * static_cast<__int128>(bmax - bmin + 1);
        if (area <= BRUTE_LIMIT) {
            int64_t res = brute_b(amin, amax, bmin, bmax);
            memo_b_[key] = res;
            return res;
        }

        int64_t a_even, a_odd, b_even, b_odd;
        count_even_odd(amin, amax, a_even, a_odd);
        count_even_odd(bmin, bmax, b_even, b_odd);

        __int128 odd_pairs = static_cast<__int128>(a_even) * b_odd + static_cast<__int128>(a_odd) * b_even;
        int64_t res = static_cast<int64_t>(odd_pairs % MOD);

        int64_t u0_min = -amax;
        int64_t u0_max = -amin;
        int64_t u1_min = 1 - amax;
        int64_t u1_max = 1 - amin;

        res = mod_add(res, sum_a(u0_min, u0_max, bmin, bmax));
        res = mod_add(res, sum_a(u1_min, u1_max, bmin, bmax));

        memo_b_[key] = res;
        return res;
    }

    int64_t sum_a(int64_t umin, int64_t umax, int64_t vmin, int64_t vmax) {
        if (umin > umax || vmin > vmax) return 0;
        Key key{umin, umax, vmin, vmax};
        auto it = memo_a_.find(key);
        if (it != memo_a_.end()) return it->second;

        int64_t u_even, u_odd, v_even, v_odd;
        count_even_odd(umin, umax, u_even, u_odd);
        count_even_odd(vmin, vmax, v_even, v_odd);

        __int128 odd_pairs = static_cast<__int128>(u_odd) * v_odd;
        int64_t res = static_cast<int64_t>(odd_pairs % MOD);

        if (u_even && v_even) {
            int64_t x_min = ceil_div(umin, 2);
            int64_t x_max = floor_div(umax, 2);
            int64_t y_min = ceil_div(vmin, 2);
            int64_t y_max = floor_div(vmax, 2);
            if (x_min <= x_max && y_min <= y_max) {
                int64_t amin = -y_max;
                int64_t amax = -y_min;
                int64_t bmin = -x_max;
                int64_t bmax = -x_min;
                res = mod_add(res, sum_b(amin, amax, bmin, bmax));
            }
        }

        if (u_odd && v_odd) {
            int64_t x_min = ceil_div(umin - 1, 2);
            int64_t x_max = floor_div(umax - 1, 2);
            int64_t y_min = ceil_div(vmin - 1, 2);
            int64_t y_max = floor_div(vmax - 1, 2);
            if (x_min <= x_max && y_min <= y_max) {
                int64_t amin = -y_max;
                int64_t amax = -y_min;
                int64_t bmin = -x_max;
                int64_t bmax = -x_min;
                res = mod_add(res, sum_b(amin, amax, bmin, bmax));
            }
        }

        memo_a_[key] = res;
        return res;
    }
};

int64_t brute_square(int64_t L) {
    int64_t total = 0;
    for (int64_t a = -L; a <= L; ++a) {
        for (int64_t b = -L; b <= L; ++b) {
            total += ones_in_base(a, b);
        }
    }
    return total % MOD;
}

int64_t compute_parallel(int64_t L, unsigned threads) {
    if (threads <= 1) {
        Solver solver;
        return solver.compute_square(L);
    }

    int64_t amin = -L;
    int64_t amax = L;
    int64_t bmin = -L;
    int64_t bmax = L;

    int64_t a_even, a_odd, b_even, b_odd;
    count_even_odd(amin, amax, a_even, a_odd);
    count_even_odd(bmin, bmax, b_even, b_odd);

    __int128 odd_pairs = static_cast<__int128>(a_even) * b_odd + static_cast<__int128>(a_odd) * b_even;
    int64_t base = static_cast<int64_t>(odd_pairs % MOD);

    int64_t u0_min = -amax;
    int64_t u0_max = -amin;
    int64_t u1_min = 1 - amax;
    int64_t u1_max = 1 - amin;

    Solver solver0;
    Solver solver1;

    auto fut0 = async(launch::async, [&]() { return solver0.compute_a(u0_min, u0_max, bmin, bmax); });
    auto fut1 = async(launch::async, [&]() { return solver1.compute_a(u1_min, u1_max, bmin, bmax); });

    int64_t res = base;
    res = mod_add(res, fut0.get());
    res = mod_add(res, fut1.get());
    return res;
}

bool validate() {
    if (ones_in_base(11, 24) != 9) {
        cerr << "Validation failed: f(11+24i) mismatch." << '\n';
        return false;
    }
    if (ones_in_base(24, -11) != 7) {
        cerr << "Validation failed: f(24-11i) mismatch." << '\n';
        return false;
    }
    if (ones_in_base(8, 0) != 3) {
        cerr << "Validation failed: f(8+0i) mismatch." << '\n';
        return false;
    }
    if (ones_in_base(-5, 0) != 5) {
        cerr << "Validation failed: f(-5+0i) mismatch." << '\n';
        return false;
    }

    if (brute_square(1) != 20) {
        cerr << "Validation failed: B(1) mismatch." << '\n';
        return false;
    }
    if (brute_square(2) != 75) {
        cerr << "Validation failed: B(2) mismatch." << '\n';
        return false;
    }

    Solver solver;
    int64_t brute20 = brute_square(20);
    int64_t fast20 = solver.compute_square(20);
    if (fast20 != brute20) {
        cerr << "Validation failed: B(20) recursion mismatch." << '\n';
        return false;
    }
    int64_t fast500 = solver.compute_square(500);
    if (fast500 != 10795060) {
        cerr << "Validation failed: B(500) mismatch." << '\n';
        return false;
    }
    return true;
}

} // namespace

int main(int argc, char** argv) {
    if (!validate()) return 1;

    int64_t L = 1000000000000000LL;
    if (argc > 1) L = max<int64_t>(0, atoll(argv[1]));
    unsigned threads = thread::hardware_concurrency();
    if (argc > 2) threads = max(1, atoi(argv[2]));

    unsigned use_threads = min<unsigned>(2, threads);
    int64_t result = compute_parallel(L, use_threads);
    cout << result % MOD << '\n';
    return 0;
}

Python

MOD = 1000000007
BRUTE_LIMIT = 1024

def floor_div(a, b):
    q = a // b
    return q

def ceil_div(a, b):
    return -floor_div(-a, b)

def count_even_odd(l, r):
    if l > r: return 0, 0
    first_even = l if l % 2 == 0 else l + 1
    last_even = r if r % 2 == 0 else r - 1
    if first_even > last_even:
        even = 0
    else:
        even = (last_even - first_even) // 2 + 1
    odd = (r - l + 1) - even
    return even, odd

def mod_add(a, b):
    return (a + b) % MOD

def ones_in_base(a, b):
    count = 0
    while a != 0 or b != 0:
        r = (a - b) & 1
        count += r
        na = (b - a + r) // 2
        nb = (r - a - b) // 2
        a, b = na, nb
    return count

class Solver:
    def __init__(self):
        self.memo_b = {}
        self.memo_a = {}
        
    def brute_b(self, amin, amax, bmin, bmax):
        total = 0
        for a in range(amin, amax + 1):
            for b in range(bmin, bmax + 1):
                total += ones_in_base(a, b)
        return total % MOD

    def sum_b(self, amin, amax, bmin, bmax):
        if amin > amax or bmin > bmax:
            return 0
        key = (amin, amax, bmin, bmax)
        if key in self.memo_b:
            return self.memo_b[key]
            
        area = (amax - amin + 1) * (bmax - bmin + 1)
        if area <= BRUTE_LIMIT:
            res = self.brute_b(amin, amax, bmin, bmax)
            self.memo_b[key] = res
            return res
            
        a_even, a_odd = count_even_odd(amin, amax)
        b_even, b_odd = count_even_odd(bmin, bmax)
        
        odd_pairs = a_even * b_odd + a_odd * b_even
        res = odd_pairs % MOD
        
        u0_min = -amax
        u0_max = -amin
        u1_min = 1 - amax
        u1_max = 1 - amin
        
        res = mod_add(res, self.sum_a(u0_min, u0_max, bmin, bmax))
        res = mod_add(res, self.sum_a(u1_min, u1_max, bmin, bmax))
        
        self.memo_b[key] = res
        return res

    def sum_a(self, umin, umax, vmin, vmax):
        if umin > umax or vmin > vmax:
            return 0
        key = (umin, umax, vmin, vmax)
        if key in self.memo_a:
            return self.memo_a[key]
            
        u_even, u_odd = count_even_odd(umin, umax)
        v_even, v_odd = count_even_odd(vmin, vmax)
        
        odd_pairs = u_odd * v_odd
        res = odd_pairs % MOD
        
        if u_even > 0 and v_even > 0:
            x_min = ceil_div(umin, 2)
            x_max = floor_div(umax, 2)
            y_min = ceil_div(vmin, 2)
            y_max = floor_div(vmax, 2)
            if x_min <= x_max and y_min <= y_max:
                res = mod_add(res, self.sum_b(-y_max, -y_min, -x_max, -x_min))
                
        if u_odd > 0 and v_odd > 0:
            x_min = ceil_div(umin - 1, 2)
            x_max = floor_div(umax - 1, 2)
            y_min = ceil_div(vmin - 1, 2)
            y_max = floor_div(vmax - 1, 2)
            if x_min <= x_max and y_min <= y_max:
                res = mod_add(res, self.sum_b(-y_max, -y_min, -x_max, -x_min))
                
        self.memo_a[key] = res
        return res

    def compute_square(self, L):
        return self.sum_b(-L, L, -L, L)

def compute_parallel(L):
    solver = Solver()
    return solver.compute_square(L) % MOD

def solve():
    L = 1000000000000000
    ans = compute_parallel(L)
    return str(ans)

if __name__ == '__main__':
    print(solve())

Java

import java.util.HashMap;
import java.util.Map;

public class Euler508 {
    private static final long MOD = 1000000007L;
    private static final long BRUTE_LIMIT = 1024L;

    static class Key {
        long a, b, c, d;

        Key(long a, long b, long c, long d) {
            this.a = a;
            this.b = b;
            this.c = c;
            this.d = d;
        }

        @Override
        public boolean equals(Object o) {
            if (this == o)
                return true;
            if (!(o instanceof Key))
                return false;
            Key key = (Key) o;
            return a == key.a && b == key.b && c == key.c && d == key.d;
        }

        @Override
        public int hashCode() {
            int result = (int) (a ^ (a >>> 32));
            result = 31 * result + (int) (b ^ (b >>> 32));
            result = 31 * result + (int) (c ^ (c >>> 32));
            result = 31 * result + (int) (d ^ (d >>> 32));
            return result;
        }
    }

    private static long floorDiv(long a, long b) {
        long q = a / b;
        long r = a % b;
        if (r != 0 && ((a < 0) ^ (b < 0))) {
            q--;
        }
        return q;
    }

    private static long ceilDiv(long a, long b) {
        return -floorDiv(-a, b);
    }

    private static long[] countEvenOdd(long l, long r) {
        if (l > r)
            return new long[] { 0, 0 };
        long firstEven = (l % 2 == 0) ? l : l + 1;
        long lastEven = (r % 2 == 0) ? r : r - 1;
        long even = 0;
        if (firstEven <= lastEven) {
            even = (lastEven - firstEven) / 2 + 1;
        }
        long odd = (r - l + 1) - even;
        return new long[] { even, odd };
    }

    private static long modAdd(long a, long b) {
        a += b;
        if (a >= MOD)
            a -= MOD;
        return a;
    }

    private static long onesInBase(long a, long b) {
        long count = 0;
        while (a != 0 || b != 0) {
            long r = (a - b) & 1L;
            count += r;
            long na = (b - a + r) / 2;
            long nb = (r - a - b) / 2;
            a = na;
            b = nb;
        }
        return count;
    }

    static class Solver {
        Map<Key, Long> memoB = new HashMap<>();
        Map<Key, Long> memoA = new HashMap<>();

        long bruteB(long amin, long amax, long bmin, long bmax) {
            long total = 0;
            for (long a = amin; a <= amax; ++a) {
                for (long b = bmin; b <= bmax; ++b) {
                    total += onesInBase(a, b);
                }
            }
            return total % MOD;
        }

        long sumB(long amin, long amax, long bmin, long bmax) {
            if (amin > amax || bmin > bmax)
                return 0;
            Key key = new Key(amin, amax, bmin, bmax);
            if (memoB.containsKey(key))
                return memoB.get(key);

            long width = amax - amin + 1;
            long height = bmax - bmin + 1;
            if (width <= BRUTE_LIMIT && height <= BRUTE_LIMIT && width * height <= BRUTE_LIMIT) {
                long res = bruteB(amin, amax, bmin, bmax);
                memoB.put(key, res);
                return res;
            }

            long[] aCounts = countEvenOdd(amin, amax);
            long[] bCounts = countEvenOdd(bmin, bmax);

            long aEven = aCounts[0], aOdd = aCounts[1];
            long bEven = bCounts[0], bOdd = bCounts[1];

            long oddPairs1 = (aEven % MOD) * (bOdd % MOD) % MOD;
            long oddPairs2 = (aOdd % MOD) * (bEven % MOD) % MOD;
            long res = (oddPairs1 + oddPairs2) % MOD;

            long u0Min = -amax;
            long u0Max = -amin;
            long u1Min = 1 - amax;
            long u1Max = 1 - amin;

            res = modAdd(res, sumA(u0Min, u0Max, bmin, bmax));
            res = modAdd(res, sumA(u1Min, u1Max, bmin, bmax));

            memoB.put(key, res);
            return res;
        }

        long sumA(long umin, long umax, long vmin, long vmax) {
            if (umin > umax || vmin > vmax)
                return 0;
            Key key = new Key(umin, umax, vmin, vmax);
            if (memoA.containsKey(key))
                return memoA.get(key);

            long[] uCounts = countEvenOdd(umin, umax);
            long[] vCounts = countEvenOdd(vmin, vmax);

            long uEven = uCounts[0], uOdd = uCounts[1];
            long vEven = vCounts[0], vOdd = vCounts[1];

            long res = (uOdd % MOD) * (vOdd % MOD) % MOD;

            if (uEven > 0 && vEven > 0) {
                long xMin = ceilDiv(umin, 2);
                long xMax = floorDiv(umax, 2);
                long yMin = ceilDiv(vmin, 2);
                long yMax = floorDiv(vmax, 2);
                if (xMin <= xMax && yMin <= yMax) {
                    res = modAdd(res, sumB(-yMax, -yMin, -xMax, -xMin));
                }
            }

            if (uOdd > 0 && vOdd > 0) {
                long xMin = ceilDiv(umin - 1, 2);
                long xMax = floorDiv(umax - 1, 2);
                long yMin = ceilDiv(vmin - 1, 2);
                long yMax = floorDiv(vmax - 1, 2);
                if (xMin <= xMax && yMin <= yMax) {
                    res = modAdd(res, sumB(-yMax, -yMin, -xMax, -xMin));
                }
            }

            memoA.put(key, res);
            return res;
        }

        long computeSquare(long L) {
            return sumB(-L, L, -L, L);
        }
    }

    public static void main(String[] args) {
        long L = 1000000000000000L;
        Solver solver = new Solver();
        System.out.println(solver.computeSquare(L));
    }
}