Problem 520: Simbers

Project Euler Problem 520 Solution

EulerSolve provides an optimized solution for Project Euler Problem 520, Simbers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A positive integer is called a simber when every even digit occurs an even number of times, and every odd digit that occurs does so an odd number of times. If $Q(n)$ denotes the number of simbers below $10^n$, the task is to evaluate $$\sum_{u=1}^{39} Q(2^u)\pmod{1{,}000{,}000{,}123}.$$ Because $2^{39}$ is enormous, direct enumeration is impossible. The solution therefore converts the digit-parity rules into algebraic projectors, compresses each fixed-length count into a short linear combination of powers $\lambda^t$, and then sums those powers as geometric series modulo the given modulus. Mathematical Approach Let $E=\{0,2,4,6,8\}$ and $O=\{1,3,5,7,9\}$. For any decimal string, let $c_d$ be the number of occurrences of digit $d$. Step 1: Turn Each Digit Rule into a Parity Projector For an even digit, the allowed counts are $0,2,4,\dots$, so its indicator is $$I_{\mathrm{even}}(c)=\frac{1+(-1)^c}{2}.$$ For an odd digit, the allowed counts are $0,1,3,5,\dots$: zero is allowed, but every positive even count is forbidden. A convenient form is $$I_{\mathrm{odd}}(c)=\delta_{c,0}+\frac{1-(-1)^c}{2}=\frac{2\cdot 0^c+1^c-(-1)^c}{2},$$ with the usual convention $0^0=1$. This identity gives $1$ when $c=0$, $1$ when $c$ is odd, and $0$ when $c$ is a positive even number....

Detailed mathematical approach

Problem Summary

A positive integer is called a simber when every even digit occurs an even number of times, and every odd digit that occurs does so an odd number of times. If $Q(n)$ denotes the number of simbers below $10^n$, the task is to evaluate

$$\sum_{u=1}^{39} Q(2^u)\pmod{1{,}000{,}000{,}123}.$$

Because $2^{39}$ is enormous, direct enumeration is impossible. The solution therefore converts the digit-parity rules into algebraic projectors, compresses each fixed-length count into a short linear combination of powers $\lambda^t$, and then sums those powers as geometric series modulo the given modulus.

Mathematical Approach

Let $E=\{0,2,4,6,8\}$ and $O=\{1,3,5,7,9\}$. For any decimal string, let $c_d$ be the number of occurrences of digit $d$.

Step 1: Turn Each Digit Rule into a Parity Projector

For an even digit, the allowed counts are $0,2,4,\dots$, so its indicator is

$$I_{\mathrm{even}}(c)=\frac{1+(-1)^c}{2}.$$

For an odd digit, the allowed counts are $0,1,3,5,\dots$: zero is allowed, but every positive even count is forbidden. A convenient form is

$$I_{\mathrm{odd}}(c)=\delta_{c,0}+\frac{1-(-1)^c}{2}=\frac{2\cdot 0^c+1^c-(-1)^c}{2},$$

with the usual convention $0^0=1$. This identity gives $1$ when $c=0$, $1$ when $c$ is odd, and $0$ when $c$ is a positive even number.

Step 2: Count All Simber Strings of One Fixed Length

Let $A_t$ be the number of length-$t$ decimal strings that satisfy the simber rules when leading zero is allowed. If the count vector is $(c_0,\dots,c_9)$, then

$$A_t=\sum_{\substack{c_0+\cdots+c_9=t\\ c_d\ge 0}} \frac{t!}{c_0!\cdots c_9!}\prod_{e\in E} I_{\mathrm{even}}(c_e)\prod_{o\in O} I_{\mathrm{odd}}(c_o).$$

After multiplying by $2^{10}$ and expanding the projectors, each even digit contributes a choice of weight $+1$ or $-1$, while each odd digit contributes a choice of weight $0$, $+1$, or $-1$ with coefficients $2$, $1$, and $-1$. For one complete global choice of weights $w_0,\dots,w_9$, the multinomial theorem gives

$$\sum_{\substack{c_0+\cdots+c_9=t\\ c_d\ge 0}} \frac{t!}{c_0!\cdots c_9!}\prod_{d=0}^{9} w_d^{\,c_d}=(w_0+\cdots+w_9)^t.$$

So the exact identities of the digits no longer matter; only the total weight matters. If $a$ digits contribute $+1$ and $b$ digits contribute $-1$, then the contribution is $\lambda^t$ with $\lambda=a-b$. After aggregating equal values of $\lambda$, we obtain a small coefficient table $C(\lambda)$ such that

$$A_t=\frac{1}{2^{10}}\sum_{\lambda=-10}^{10} C(\lambda)\lambda^t.$$

Step 3: Subtract the Strings That Start with Zero

A positive integer with exactly $t$ digits corresponds to a length-$t$ string that does not begin with zero. Therefore we must subtract the simber strings whose first character is $0$.

If we remove that leading zero, the remaining suffix has length $t-1$. Since digit $0$ has already appeared once, its remaining count must be odd so that the total number of zeros is even. The other four nonzero even digits keep the usual even rule, and the five odd digits keep the usual “zero or odd” rule. Hence the distinguished zero digit uses

$$J(c)=\frac{1-(-1)^c}{2},$$

which is the indicator of an odd count. Let $B_s$ be the number of length-$s$ suffixes with that modified rule for digit $0$. By the same expansion-and-grouping argument, there is another coefficient table $D(\lambda)$ with

$$B_s=\frac{1}{2^{10}}\sum_{\lambda=-10}^{10} D(\lambda)\lambda^s.$$

Therefore the number of genuine $t$-digit simbers is

$$N_t=A_t-B_{t-1}.$$

Step 4: Worked Example for Two Digits

For $t=2$, the leading-zero-allowed simber strings are easy to classify. There are five repeated even-digit strings

$$00,\ 22,\ 44,\ 66,\ 88,$$

and there are $5\cdot 4=20$ ordered strings formed by two distinct odd digits. Hence $A_2=25$.

Among these, only $00$ starts with zero, so $B_1=1$ and

$$N_2=A_2-B_1=25-1=24.$$

The one-digit simbers are $1,3,5,7,9$, so $N_1=5$. Therefore

$$Q(2)=N_1+N_2=5+24=29.$$

This tiny case is exactly what the general projector method is automating for huge values of $n$.

Step 5: Sum Over All Lengths with Geometric Series

Now sum the exact-length counts from $1$ to $n$:

$$Q(n)=\sum_{t=1}^{n} N_t=\frac{1}{2^{10}}\left(\sum_{\lambda} C(\lambda)\sum_{t=1}^{n}\lambda^t-\sum_{\lambda} D(\lambda)\sum_{s=0}^{n-1}\lambda^s\right).$$

For $\lambda\neq 0,1$, the inner sums are

$$\sum_{t=1}^{n}\lambda^t=\lambda\frac{\lambda^n-1}{\lambda-1},\qquad \sum_{s=0}^{n-1}\lambda^s=\frac{\lambda^n-1}{\lambda-1}.$$

The exceptional values are handled separately:

$$\sum_{t=1}^{n}1^t=n,\qquad \sum_{s=0}^{n-1}1^s=n,$$

$$\sum_{t=1}^{n}0^t=0,\qquad \sum_{s=0}^{n-1}0^s=1\quad (n\ge 1).$$

All arithmetic is performed modulo

$$M=1{,}000{,}000{,}123,$$

so division by $\lambda-1$ and by $2^{10}$ becomes multiplication by modular inverses.

Step 6: Why the Whole Method Stays Small

The decimal alphabet is fixed. Each of the five odd digits contributes one of three weights $\{0,+1,-1\}$, and each of the five even digits contributes one of two weights $\{+1,-1\}$. Therefore $\lambda$ can only lie between $-10$ and $10$, so after aggregation there are at most $21$ distinct powers to evaluate. The difficult combinatorics has been compressed into two tiny coefficient tables, and the remaining work is just fast modular exponentiation together with short geometric-sum evaluations.

How the Code Works

The C++, Python, and Java implementations build the two coefficient tables directly from the projector factors. One table corresponds to all length-$t$ strings, and the other corresponds to the leading-zero correction in which digit $0$ must occur an odd number of times in the suffix. In both cases, the expansion is grouped only by the value of $\lambda$, because that is the only quantity that survives the multinomial sum.

For a requested $n$, the implementation evaluates the first table against $\sum_{t=1}^{n}\lambda^t$ and the second table against $\sum_{s=0}^{n-1}\lambda^s$. Fast modular exponentiation is used inside the closed forms for those geometric sums, and the final factor $2^{-10}$ is applied with a modular inverse.

Finally, the program computes $Q(2^u)$ for $u=1,2,\dots,39$ and accumulates the answers modulo $M$. It also checks two published checkpoints,

$$Q(7)=287975,\qquad Q(100)\equiv 123864868 \pmod{M},$$

and confirms the small case $Q(4)$ by direct brute force.

Complexity Analysis

With the decimal digit set fixed, constructing the two coefficient tables is constant work and constant memory. After aggregation, each table has only a small number of nonzero $\lambda$-entries, bounded by the range $-10$ to $10$.

Each evaluation of $Q(n)$ therefore needs only a short loop over those entries, and each term uses $O(\log n)$ time for modular exponentiation. So one query costs $O(L\log n)$, where $L$ is the number of distinct $\lambda$-values, and in this problem $L$ is tiny. Summing the 39 required values is easily within practical limits, and the memory usage stays $O(L)$.

Footnotes and References

Problem 520 source code

C++

#include <algorithm>
#include <array>
#include <cstdint>
#include <iostream>
#include <map>
#include <string>
#include <tuple>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;
using u128 = __uint128_t;

constexpr u64 MOD = 1'000'000'123ULL;

struct Options {
    int max_power = 39;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(ch - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--max-power=", options.max_power)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.max_power >= 1 && options.max_power <= 62;
}

u64 mod_pow(u64 base, u64 exp) {
    u64 result = 1ULL;
    base %= MOD;
    while (exp > 0ULL) {
        if (exp & 1ULL) {
            result = static_cast<u64>((static_cast<u128>(result) * base) % MOD);
        }
        base = static_cast<u64>((static_cast<u128>(base) * base) % MOD);
        exp >>= 1ULL;
    }
    return result;
}

u64 mod_inverse(const u64 x) {
    return mod_pow(x % MOD, MOD - 2ULL);
}

u64 geom_sum_from_one(const i64 lambda, const u64 n) {
    if (n == 0ULL) {
        return 0ULL;
    }
    if (lambda == 1) {
        return n % MOD;
    }
    if (lambda == 0) {
        return 0ULL;
    }
    const u64 lm = static_cast<u64>((lambda % static_cast<i64>(MOD) + static_cast<i64>(MOD)) % static_cast<i64>(MOD));
    const u64 num = static_cast<u64>((static_cast<u128>(lm) * ((mod_pow(lm, n) + MOD - 1ULL) % MOD)) % MOD);
    const u64 den = (lm + MOD - 1ULL) % MOD;
    return static_cast<u64>((static_cast<u128>(num) * mod_inverse(den)) % MOD);
}

u64 geom_sum_from_zero(const i64 lambda, const u64 n) {
    if (n == 0ULL) {
        return 0ULL;
    }
    if (lambda == 1) {
        return n % MOD;
    }
    if (lambda == 0) {
        return 1ULL;
    }
    const u64 lm = static_cast<u64>((lambda % static_cast<i64>(MOD) + static_cast<i64>(MOD)) % static_cast<i64>(MOD));
    const u64 num = (mod_pow(lm, n) + MOD - 1ULL) % MOD;
    const u64 den = (lm + MOD - 1ULL) % MOD;
    return static_cast<u64>((static_cast<u128>(num) * mod_inverse(den)) % MOD);
}

std::map<i64, i64> expand_coeffs_F() {
    std::map<std::pair<int, int>, i64> poly1;
    poly1[{0, 0}] = 1;
    const std::array<std::tuple<int, int, i64>, 3> base1 = {
        std::make_tuple(0, 0, 2), std::make_tuple(1, 0, 1), std::make_tuple(0, 1, -1)};

    for (int rep = 0; rep < 5; ++rep) {
        std::map<std::pair<int, int>, i64> next;
        for (const auto& row : poly1) {
            for (const auto& term : base1) {
                const int da = std::get<0>(term);
                const int db = std::get<1>(term);
                const i64 dc = std::get<2>(term);
                next[{row.first.first + da, row.first.second + db}] += row.second * dc;
            }
        }
        poly1.swap(next);
    }

    std::map<std::pair<int, int>, i64> poly2;
    poly2[{0, 0}] = 1;
    const std::array<std::tuple<int, int, i64>, 2> base2 = {std::make_tuple(1, 0, 1), std::make_tuple(0, 1, 1)};
    for (int rep = 0; rep < 5; ++rep) {
        std::map<std::pair<int, int>, i64> next;
        for (const auto& row : poly2) {
            for (const auto& term : base2) {
                next[{row.first.first + std::get<0>(term), row.first.second + std::get<1>(term)}] +=
                    row.second * std::get<2>(term);
            }
        }
        poly2.swap(next);
    }

    std::map<std::pair<int, int>, i64> both;
    for (const auto& a : poly1) {
        for (const auto& b : poly2) {
            both[{a.first.first + b.first.first, a.first.second + b.first.second}] += a.second * b.second;
        }
    }

    std::map<i64, i64> by_lambda;
    for (const auto& row : both) {
        by_lambda[static_cast<i64>(row.first.first - row.first.second)] += row.second;
    }
    return by_lambda;
}

std::map<i64, i64> expand_coeffs_G() {
    std::map<std::pair<int, int>, i64> poly1;
    poly1[{0, 0}] = 1;
    const std::array<std::tuple<int, int, i64>, 3> base1 = {
        std::make_tuple(0, 0, 2), std::make_tuple(1, 0, 1), std::make_tuple(0, 1, -1)};
    for (int rep = 0; rep < 5; ++rep) {
        std::map<std::pair<int, int>, i64> next;
        for (const auto& row : poly1) {
            for (const auto& term : base1) {
                next[{row.first.first + std::get<0>(term), row.first.second + std::get<1>(term)}] +=
                    row.second * std::get<2>(term);
            }
        }
        poly1.swap(next);
    }

    std::map<std::pair<int, int>, i64> poly2;
    poly2[{0, 0}] = 1;
    const std::array<std::tuple<int, int, i64>, 2> plusminus = {std::make_tuple(1, 0, 1), std::make_tuple(0, 1, 1)};
    for (int rep = 0; rep < 4; ++rep) {
        std::map<std::pair<int, int>, i64> next;
        for (const auto& row : poly2) {
            for (const auto& term : plusminus) {
                next[{row.first.first + std::get<0>(term), row.first.second + std::get<1>(term)}] +=
                    row.second * std::get<2>(term);
            }
        }
        poly2.swap(next);
    }

    std::map<std::pair<int, int>, i64> diff;
    diff[{1, 0}] = 1;
    diff[{0, 1}] = -1;

    std::map<std::pair<int, int>, i64> tmp;
    for (const auto& a : poly1) {
        for (const auto& b : poly2) {
            tmp[{a.first.first + b.first.first, a.first.second + b.first.second}] += a.second * b.second;
        }
    }
    std::map<std::pair<int, int>, i64> both;
    for (const auto& a : tmp) {
        for (const auto& b : diff) {
            both[{a.first.first + b.first.first, a.first.second + b.first.second}] += a.second * b.second;
        }
    }

    std::map<i64, i64> by_lambda;
    for (const auto& row : both) {
        by_lambda[static_cast<i64>(row.first.first - row.first.second)] += row.second;
    }
    return by_lambda;
}

u64 Q(const u64 n, const std::map<i64, i64>& cf, const std::map<i64, i64>& cg, const u64 inv_two_pow_10) {
    u64 sum_A = 0ULL;
    for (const auto& row : cf) {
        const i64 coeff = row.second;
        const u64 geom = geom_sum_from_one(row.first, n);
        const u64 coeff_mod = static_cast<u64>((coeff % static_cast<i64>(MOD) + static_cast<i64>(MOD)) %
                                               static_cast<i64>(MOD));
        sum_A = (sum_A + static_cast<u64>((static_cast<u128>(coeff_mod) * geom) % MOD)) % MOD;
    }

    u64 sum_B = 0ULL;
    for (const auto& row : cg) {
        const i64 coeff = row.second;
        const u64 geom = geom_sum_from_zero(row.first, n);
        const u64 coeff_mod = static_cast<u64>((coeff % static_cast<i64>(MOD) + static_cast<i64>(MOD)) %
                                               static_cast<i64>(MOD));
        sum_B = (sum_B + static_cast<u64>((static_cast<u128>(coeff_mod) * geom) % MOD)) % MOD;
    }

    return static_cast<u64>((static_cast<u128>((sum_A + MOD - sum_B) % MOD) * inv_two_pow_10) % MOD);
}

u64 solve(const int max_power) {
    const std::map<i64, i64> cf = expand_coeffs_F();
    const std::map<i64, i64> cg = expand_coeffs_G();
    const u64 inv_two_pow_10 = mod_inverse(1024ULL);

    u64 ans = 0ULL;
    for (int u = 1; u <= max_power; ++u) {
        const u64 n = 1ULL << static_cast<unsigned>(u);
        ans += Q(n, cf, cg, inv_two_pow_10);
        if (ans >= MOD) {
            ans -= MOD;
        }
    }
    return ans;
}

u64 brute_Q_small(const int n) {
    u64 count = 0ULL;
    const auto is_simber = [](u64 x) {
        std::array<int, 10> cnt{};
        while (x > 0ULL) {
            ++cnt[static_cast<std::size_t>(x % 10ULL)];
            x /= 10ULL;
        }
        for (int d = 0; d <= 9; ++d) {
            if ((d % 2) == 0) {
                if ((cnt[static_cast<std::size_t>(d)] % 2) != 0) {
                    return false;
                }
            } else {
                if ((cnt[static_cast<std::size_t>(d)] % 2) == 0 && cnt[static_cast<std::size_t>(d)] != 0) {
                    return false;
                }
            }
        }
        return true;
    };

    u64 upper = 1ULL;
    for (int i = 0; i < n; ++i) {
        upper *= 10ULL;
    }
    for (u64 x = 1ULL; x < upper; ++x) {
        if (is_simber(x)) {
            ++count;
        }
    }
    return count;
}

bool run_checkpoints() {
    const std::map<i64, i64> cf = expand_coeffs_F();
    const std::map<i64, i64> cg = expand_coeffs_G();
    const u64 inv_two_pow_10 = mod_inverse(1024ULL);

    if (Q(7ULL, cf, cg, inv_two_pow_10) != 287'975ULL) {
        std::cerr << "Checkpoint failed: Q(7)=287975" << '\n';
        return false;
    }
    if (Q(100ULL, cf, cg, inv_two_pow_10) != 123'864'868ULL) {
        std::cerr << "Checkpoint failed: Q(100) mod 1,000,000,123 = 123864868" << '\n';
        return false;
    }
    if (Q(4ULL, cf, cg, inv_two_pow_10) != brute_Q_small(4) % MOD) {
        std::cerr << "Checkpoint failed: brute-force cross-check for Q(4)" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << solve(options.max_power) << '\n';
    return 0;
}

Python

MOD = 1000000123

def mod_pow(base, exp):
    return pow(base, exp, MOD)

def mod_inverse(x):
    return pow(x % MOD, MOD - 2, MOD)

def geom_sum_from_one(lam, n):
    if n == 0:
        return 0
    if lam == 1:
        return n % MOD
    if lam == 0:
        return 0
    lm = lam % MOD
    if lm < 0:
        lm += MOD
    num = (lm * ((mod_pow(lm, n) + MOD - 1) % MOD)) % MOD
    den = (lm + MOD - 1) % MOD
    return (num * mod_inverse(den)) % MOD

def geom_sum_from_zero(lam, n):
    if n == 0:
        return 0
    if lam == 1:
        return n % MOD
    if lam == 0:
        return 1
    lm = lam % MOD
    if lm < 0:
        lm += MOD
    num = (mod_pow(lm, n) + MOD - 1) % MOD
    den = (lm + MOD - 1) % MOD
    return (num * mod_inverse(den)) % MOD

def expand_coeffs_F():
    poly1 = {(0, 0): 1}
    base1 = [(0, 0, 2), (1, 0, 1), (0, 1, -1)]
    for _ in range(5):
        nxt = {}
        for (a, b), val in poly1.items():
            for da, db, dc in base1:
                nxt[(a+da, b+db)] = nxt.get((a+da, b+db), 0) + val * dc
        poly1 = nxt

    poly2 = {(0, 0): 1}
    base2 = [(1, 0, 1), (0, 1, 1)]
    for _ in range(5):
        nxt = {}
        for (a, b), val in poly2.items():
            for da, db, dc in base2:
                nxt[(a+da, b+db)] = nxt.get((a+da, b+db), 0) + val * dc
        poly2 = nxt

    both = {}
    for (a1, b1), v1 in poly1.items():
        for (a2, b2), v2 in poly2.items():
            both[(a1+a2, b1+b2)] = both.get((a1+a2, b1+b2), 0) + v1 * v2

    by_lambda = {}
    for (a, b), val in both.items():
        by_lambda[a - b] = by_lambda.get(a - b, 0) + val
    return by_lambda

def expand_coeffs_G():
    poly1 = {(0, 0): 1}
    base1 = [(0, 0, 2), (1, 0, 1), (0, 1, -1)]
    for _ in range(5):
        nxt = {}
        for (a, b), val in poly1.items():
            for da, db, dc in base1:
                nxt[(a+da, b+db)] = nxt.get((a+da, b+db), 0) + val * dc
        poly1 = nxt

    poly2 = {(0, 0): 1}
    base2 = [(1, 0, 1), (0, 1, 1)]
    for _ in range(4):
        nxt = {}
        for (a, b), val in poly2.items():
            for da, db, dc in base2:
                nxt[(a+da, b+db)] = nxt.get((a+da, b+db), 0) + val * dc
        poly2 = nxt

    diff = {(1, 0): 1, (0, 1): -1}

    tmp = {}
    for (a1, b1), v1 in poly1.items():
        for (a2, b2), v2 in poly2.items():
            tmp[(a1+a2, b1+b2)] = tmp.get((a1+a2, b1+b2), 0) + v1 * v2

    both = {}
    for (a1, b1), v1 in tmp.items():
        for (a2, b2), v2 in diff.items():
            both[(a1+a2, b1+b2)] = both.get((a1+a2, b1+b2), 0) + v1 * v2

    by_lambda = {}
    for (a, b), val in both.items():
        by_lambda[a - b] = by_lambda.get(a - b, 0) + val
    return by_lambda

def Q(n, cf, cg, inv_two_pow_10):
    sum_A = 0
    for lam, coeff in cf.items():
        geom = geom_sum_from_one(lam, n)
        coeff_mod = coeff % MOD
        if coeff_mod < 0:
            coeff_mod += MOD
        sum_A = (sum_A + coeff_mod * geom) % MOD

    sum_B = 0
    for lam, coeff in cg.items():
        geom = geom_sum_from_zero(lam, n)
        coeff_mod = coeff % MOD
        if coeff_mod < 0:
            coeff_mod += MOD
        sum_B = (sum_B + coeff_mod * geom) % MOD

    return (((sum_A + MOD - sum_B) % MOD) * inv_two_pow_10) % MOD

def solve():
    max_power = 39
    cf = expand_coeffs_F()
    cg = expand_coeffs_G()
    inv_two_pow_10 = mod_inverse(1024)

    ans = 0
    for u in range(1, max_power + 1):
        n = 1 << u
        ans += Q(n, cf, cg, inv_two_pow_10)
        if ans >= MOD:
            ans -= MOD
    return str(ans)

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

Java

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

public class Euler520 {

    static final long MOD = 1000000123L;

    static long modPow(long base, long exp) {
        long result = 1L;
        base %= MOD;
        while (exp > 0) {
            if ((exp & 1) == 1) {
                result = (result * base) % MOD;
            }
            base = (base * base) % MOD;
            exp >>= 1;
        }
        return result;
    }

    static long modInverse(long x) {
        return modPow(x % MOD, MOD - 2);
    }

    static long geomSumFromOne(long lambda, long n) {
        if (n == 0)
            return 0;
        if (lambda == 1)
            return n % MOD;
        if (lambda == 0)
            return 0;

        long lm = (lambda % MOD + MOD) % MOD;
        long num = (lm * ((modPow(lm, n) + MOD - 1) % MOD)) % MOD;
        long den = (lm + MOD - 1) % MOD;
        return (num * modInverse(den)) % MOD;
    }

    static long geomSumFromZero(long lambda, long n) {
        if (n == 0)
            return 0;
        if (lambda == 1)
            return n % MOD;
        if (lambda == 0)
            return 1;

        long lm = (lambda % MOD + MOD) % MOD;
        long num = (modPow(lm, n) + MOD - 1) % MOD;
        long den = (lm + MOD - 1) % MOD;
        return (num * modInverse(den)) % MOD;
    }

    static class Pair {
        int a, b;

        Pair(int a, int b) {
            this.a = a;
            this.b = b;
        }

        @Override
        public boolean equals(Object o) {
            if (this == o)
                return true;
            if (o == null || getClass() != o.getClass())
                return false;
            Pair pair = (Pair) o;
            return a == pair.a && b == pair.b;
        }

        @Override
        public int hashCode() {
            return Objects.hash(a, b);
        }
    }

    static Map<Long, Long> expandCoeffsF() {
        Map<Pair, Long> poly1 = new HashMap<>();
        poly1.put(new Pair(0, 0), 1L);
        int[][] base1 = { { 0, 0, 2 }, { 1, 0, 1 }, { 0, 1, -1 } };

        for (int rep = 0; rep < 5; rep++) {
            Map<Pair, Long> next = new HashMap<>();
            for (Map.Entry<Pair, Long> e : poly1.entrySet()) {
                for (int[] term : base1) {
                    Pair nk = new Pair(e.getKey().a + term[0], e.getKey().b + term[1]);
                    next.put(nk, next.getOrDefault(nk, 0L) + e.getValue() * term[2]);
                }
            }
            poly1 = next;
        }

        Map<Pair, Long> poly2 = new HashMap<>();
        poly2.put(new Pair(0, 0), 1L);
        int[][] base2 = { { 1, 0, 1 }, { 0, 1, 1 } };

        for (int rep = 0; rep < 5; rep++) {
            Map<Pair, Long> next = new HashMap<>();
            for (Map.Entry<Pair, Long> e : poly2.entrySet()) {
                for (int[] term : base2) {
                    Pair nk = new Pair(e.getKey().a + term[0], e.getKey().b + term[1]);
                    next.put(nk, next.getOrDefault(nk, 0L) + e.getValue() * term[2]);
                }
            }
            poly2 = next;
        }

        Map<Pair, Long> both = new HashMap<>();
        for (Map.Entry<Pair, Long> a : poly1.entrySet()) {
            for (Map.Entry<Pair, Long> b : poly2.entrySet()) {
                Pair nk = new Pair(a.getKey().a + b.getKey().a, a.getKey().b + b.getKey().b);
                both.put(nk, both.getOrDefault(nk, 0L) + a.getValue() * b.getValue());
            }
        }

        Map<Long, Long> byLambda = new HashMap<>();
        for (Map.Entry<Pair, Long> row : both.entrySet()) {
            long lam = row.getKey().a - row.getKey().b;
            byLambda.put(lam, byLambda.getOrDefault(lam, 0L) + row.getValue());
        }
        return byLambda;
    }

    static Map<Long, Long> expandCoeffsG() {
        Map<Pair, Long> poly1 = new HashMap<>();
        poly1.put(new Pair(0, 0), 1L);
        int[][] base1 = { { 0, 0, 2 }, { 1, 0, 1 }, { 0, 1, -1 } };

        for (int rep = 0; rep < 5; rep++) {
            Map<Pair, Long> next = new HashMap<>();
            for (Map.Entry<Pair, Long> e : poly1.entrySet()) {
                for (int[] term : base1) {
                    Pair nk = new Pair(e.getKey().a + term[0], e.getKey().b + term[1]);
                    next.put(nk, next.getOrDefault(nk, 0L) + e.getValue() * term[2]);
                }
            }
            poly1 = next;
        }

        Map<Pair, Long> poly2 = new HashMap<>();
        poly2.put(new Pair(0, 0), 1L);
        int[][] base2 = { { 1, 0, 1 }, { 0, 1, 1 } };

        for (int rep = 0; rep < 4; rep++) {
            Map<Pair, Long> next = new HashMap<>();
            for (Map.Entry<Pair, Long> e : poly2.entrySet()) {
                for (int[] term : base2) {
                    Pair nk = new Pair(e.getKey().a + term[0], e.getKey().b + term[1]);
                    next.put(nk, next.getOrDefault(nk, 0L) + e.getValue() * term[2]);
                }
            }
            poly2 = next;
        }

        Map<Pair, Long> diff = new HashMap<>();
        diff.put(new Pair(1, 0), 1L);
        diff.put(new Pair(0, 1), -1L);

        Map<Pair, Long> tmp = new HashMap<>();
        for (Map.Entry<Pair, Long> a : poly1.entrySet()) {
            for (Map.Entry<Pair, Long> b : poly2.entrySet()) {
                Pair nk = new Pair(a.getKey().a + b.getKey().a, a.getKey().b + b.getKey().b);
                tmp.put(nk, tmp.getOrDefault(nk, 0L) + a.getValue() * b.getValue());
            }
        }

        Map<Pair, Long> both = new HashMap<>();
        for (Map.Entry<Pair, Long> a : tmp.entrySet()) {
            for (Map.Entry<Pair, Long> b : diff.entrySet()) {
                Pair nk = new Pair(a.getKey().a + b.getKey().a, a.getKey().b + b.getKey().b);
                both.put(nk, both.getOrDefault(nk, 0L) + a.getValue() * b.getValue());
            }
        }

        Map<Long, Long> byLambda = new HashMap<>();
        for (Map.Entry<Pair, Long> row : both.entrySet()) {
            long lam = row.getKey().a - row.getKey().b;
            byLambda.put(lam, byLambda.getOrDefault(lam, 0L) + row.getValue());
        }
        return byLambda;
    }

    static long Q(long n, Map<Long, Long> cf, Map<Long, Long> cg, long invTwoPow10) {
        long sumA = 0;
        for (Map.Entry<Long, Long> entry : cf.entrySet()) {
            long geom = geomSumFromOne(entry.getKey(), n);
            long coeffMod = (entry.getValue() % MOD + MOD) % MOD;
            sumA = (sumA + (coeffMod * geom) % MOD) % MOD;
        }

        long sumB = 0;
        for (Map.Entry<Long, Long> entry : cg.entrySet()) {
            long geom = geomSumFromZero(entry.getKey(), n);
            long coeffMod = (entry.getValue() % MOD + MOD) % MOD;
            sumB = (sumB + (coeffMod * geom) % MOD) % MOD;
        }

        return (((sumA + MOD - sumB) % MOD) * invTwoPow10) % MOD;
    }

    public static void main(String[] args) {
        int maxPower = 39;
        Map<Long, Long> cf = expandCoeffsF();
        Map<Long, Long> cg = expandCoeffsG();
        long invTwoPow10 = modInverse(1024L);

        long ans = 0;
        for (int u = 1; u <= maxPower; u++) {
            long n = 1L << u;
            ans += Q(n, cf, cg, invTwoPow10);
            if (ans >= MOD) {
                ans -= MOD;
            }
        }
        System.out.println(ans);
    }
}