Problem 847: Jack's Bean

Project Euler Problem 847 Solution

EulerSolve provides an optimized solution for Project Euler Problem 847, Jack's Bean, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary There are three plates containing $a$, $b$, and $c$ beans, and exactly one bean is magical. A question may inspect only a subset taken from a single plate, and the answer is just yes or no. Let $h(a,b,c)$ be the minimum number of such questions that always identifies the magical bean. The task is to evaluate $$H(N)=\sum_{\substack{a,b,c\in\mathbb{Z}_{\ge 0}\\1\le a+b+c\le N}} h(a,b,c) \pmod{10^9+7}.$$ A brute-force search over triples or over decision trees is hopeless for the true input size, so the solution turns the problem into a counting argument over dyadic ranges. Mathematical Approach Write $s=a+b+c$ for the total number of candidate beans. The whole strategy revolves around the threshold $$D=2^{k-1}\qquad\text{whenever}\qquad 2^{k-1}\lt s\le 2^k.$$ At that scale, $k$ is the information-theoretic lower bound, and the only question is whether a given triple attains that bound or needs one extra step. Step 1: Baseline Lower Bound from Binary Questions With $s$ possible locations for the magical bean, any yes/no decision tree of depth $t$ has at most $2^t$ leaves....

Detailed mathematical approach

Problem Summary

There are three plates containing $a$, $b$, and $c$ beans, and exactly one bean is magical. A question may inspect only a subset taken from a single plate, and the answer is just yes or no. Let $h(a,b,c)$ be the minimum number of such questions that always identifies the magical bean. The task is to evaluate

$$H(N)=\sum_{\substack{a,b,c\in\mathbb{Z}_{\ge 0}\\1\le a+b+c\le N}} h(a,b,c) \pmod{10^9+7}.$$

A brute-force search over triples or over decision trees is hopeless for the true input size, so the solution turns the problem into a counting argument over dyadic ranges.

Mathematical Approach

Write $s=a+b+c$ for the total number of candidate beans. The whole strategy revolves around the threshold

$$D=2^{k-1}\qquad\text{whenever}\qquad 2^{k-1}\lt s\le 2^k.$$

At that scale, $k$ is the information-theoretic lower bound, and the only question is whether a given triple attains that bound or needs one extra step.

Step 1: Baseline Lower Bound from Binary Questions

With $s$ possible locations for the magical bean, any yes/no decision tree of depth $t$ has at most $2^t$ leaves. Therefore every strategy satisfies

$$h(a,b,c)\ge \left\lceil \log_2 s \right\rceil.$$

For fixed $s$, the number of nonnegative triples $(a,b,c)$ with total $s$ is the stars-and-bars count

$$\binom{s+2}{2}.$$

So the universal baseline contribution is

$$B(N)=\sum_{s=1}^{N}\left\lceil \log_2 s \right\rceil \binom{s+2}{2}.$$

Since $\lceil\log_2 1\rceil=0$, the $s=1$ term vanishes automatically. The implementations group this sum by dyadic blocks $(2^{k-1},2^k]$ and expand

$$\binom{s+2}{2}=\frac{s^2+3s+2}{2},$$

so each block can be evaluated from closed forms for $\sum s$ and $\sum s^2$.

Step 2: When Can the Lower Bound Be Achieved?

Fix a triple with $2^{k-1}\lt s\le 2^k$, and set $D=2^{k-1}$. If we hope to finish in exactly $k$ questions, then after the first question both branches must leave at most $D$ candidates, because only $k-1$ further yes/no answers remain.

Suppose the first question selects a subset of size $q$ from one plate. Then the two branch sizes are $q$ and $s-q$, so we must have

$$s-D\le q\le D.$$

Inside the cube $0\le a,b,c\le D$, the chosen subset can only come from a plate that contains at least $s-D$ beans. Thus a triple is easy at level $k$ if at least one coordinate is at least $s-D$. It is hard if no coordinate reaches that threshold.

Because $s=a+b+c$, the condition $a\lt s-D$ is equivalent to $b+c\gt D$, and similarly for the other coordinates. Therefore, inside the cube, hardness is exactly the system

$$a+b\gt D,\qquad a+c\gt D,\qquad b+c\gt D.$$

Every hard triple needs $k+1$ questions; every non-hard triple attains the lower bound $k$.

Step 3: Count the Hard Region Inside $[0,D]^3$

Define

$$W_D(X)=\#\left\{(a,b,c)\in\mathbb{Z}_{\ge 0}^3:0\le a,b,c\le D,\ a+b+c\le X,\ a+b\gt D,\ a+c\gt D,\ b+c\gt D\right\}.$$

This is the hard region at level $k$ restricted to triples with total at most $X$. To count it, use inclusion-exclusion on the three “easy” events

$$E_{ab}:\ a+b\le D,\qquad E_{ac}:\ a+c\le D,\qquad E_{bc}:\ b+c\le D.$$

By symmetry,

$$W_D(X)=U_D(X)-3A_D(X)+3B_D(X)-C_D(X),$$

where $U_D(X)$ counts all triples in the box with total at most $X$, $A_D(X)$ counts one easy event such as $E_{ab}$, $B_D(X)$ counts an intersection of two easy events, and $C_D(X)$ counts the intersection of all three.

The simplest term is the universe count:

$$U_D(X)= \begin{cases} \binom{X+3}{3}, & X\le D,\\ \binom{X+3}{3}-3\binom{X-D+2}{3}, & D\lt X\le 2D. \end{cases}$$

The remaining terms reduce to one-dimensional polynomial sums after fixing either a pair sum or one coordinate. That is why the implementation only needs closed forms for arithmetic and quadratic progressions rather than any explicit enumeration of triples.

Step 4: Recursive Reduction When One Plate Exceeds $D$

Now consider a triple with total $s\le 2D$ but with one coordinate larger than $D$. Because $s\le 2D$, at most one coordinate can exceed $D$.

Assume $a\gt D$. To keep both branches within size $D$, the first question must isolate exactly $D$ beans from that oversized plate. The “yes” branch then has exactly $D$ candidates, while the “no” branch leaves the reduced triple

$$\left(a-D,\ b,\ c\right),$$

whose total is $s-D$.

So the original triple is hard at level $k$ precisely when the reduced triple is hard at level $k-1$. If we let $E_k(X)$ denote the number of hard triples with total at most $X$ on level $k$, then

$$E_k(X)= \begin{cases} W_{2^{k-1}}(X), & X\le 2^{k-1},\\ W_{2^{k-1}}(X)+3E_{k-1}(X-2^{k-1}), & 2^{k-1}\lt X\le 2^k. \end{cases}$$

The factor $3$ comes from the choice of which plate is the oversized one.

Step 5: Final Summation Formula

Let $K$ be the largest integer with $2^{K-1}\lt N\le 2^K$. Every triple in the dyadic block $2^{k-1}\lt s\le 2^k$ contributes the baseline $k$, and every hard triple contributes one extra question. Hence

$$H(N)=\sum_{k=1}^{K} k\sum_{s=2^{k-1}+1}^{\min(N,2^k)}\binom{s+2}{2}+\sum_{k=1}^{K}E_k\!\left(\min(N,2^k)\right)\pmod{10^9+7}.$$

This is the exact formula used by the implementation.

Worked Example

Take $s=8$. Then $k=3$ and $D=4$.

For the triple $(3,3,2)$, all pair sums exceed $4$:

$$3+3\gt 4,\qquad 3+2\gt 4,\qquad 3+2\gt 4.$$

So it lies in the hard region. Any three-question strategy would need its first question to split the eight candidates into two branches of size at most four, which forces a queried subset of size exactly $4$. But no plate contains $4$ beans, so such a first move is impossible. Therefore

$$h(3,3,2)=4.$$

By contrast, $(5,2,1)$ is easy at the same level: ask about a four-bean subset of the plate containing five beans. The “yes” branch leaves $4$ candidates, and the “no” branch leaves $(1,2,1)$, whose total is $4$, so two more questions suffice. Thus $h(5,2,1)=3$.

The recursive idea appears one level higher with $(11,3,2)$, where $s=16$ and $D=8$. Asking about an eight-bean subset of the large plate leaves either a block of size $8$ or the reduced triple $(3,3,2)$, already known to be hard. Hence $(11,3,2)$ is hard at level $4$ and needs $5$ questions.

How the Code Works

The C++, Python, and Java implementations work entirely modulo $10^9+7$. They first accumulate the baseline term by grouping totals $s$ into dyadic intervals and evaluating

$$\sum \binom{s+2}{2}=\sum \frac{s^2+3s+2}{2}$$

with closed forms for $\sum s$ and $\sum s^2$. Next they evaluate the hard-region counts $W_D(X)$ by inclusion-exclusion. Each component of that inclusion-exclusion is reduced to a short list of polynomial range sums, so no loop over triples is ever required.

Finally, the implementation memoizes the recursive states $E_k(X)$. Each state performs one within-cube count and, when $X\gt 2^{k-1}$, one recursive call to the previous level. The answer is the baseline term plus the total number of hard triples across all relevant levels.

Complexity Analysis

Let $K=\lceil\log_2 N\rceil$. The number of dyadic levels is $O(K)$, and each recursive state is evaluated once. Every state uses only constant-time modular arithmetic together with closed-form sums, so the total running time is $O(\log N)$. The memoized state table also has size $O(\log N)$.

Footnotes and References

Problem page: https://projecteuler.net/problem=847
Decision tree: Wikipedia — Decision tree
Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
Stars and bars: Wikipedia — Stars and bars
Binomial coefficient: Wikipedia — Binomial coefficient

Problem 847 source code

C++

#include <cstdint>
#include <iostream>
#include <unordered_map>
#include <functional>
#include <vector>

using namespace std;

namespace {

static constexpr int64_t MOD = 1'000'000'007LL;

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

int64_t mod_sub(int64_t a, int64_t b) {
    int64_t r = a - b;
    if (r < 0) r += MOD;
    return r;
}

int64_t mod_mul(int64_t a, int64_t b) {
    return static_cast<int64_t>((__int128)a * b % MOD);
}

int64_t mod_pow(int64_t a, int64_t e) {
    int64_t r = 1 % MOD;
    int64_t x = (a % MOD + MOD) % MOD;
    while (e > 0) {
        if (e & 1) r = mod_mul(r, x);
        x = mod_mul(x, x);
        e >>= 1;
    }
    return r;
}

const int64_t INV2 = mod_pow(2, MOD - 2);
const int64_t INV6 = mod_pow(6, MOD - 2);

int64_t sum1_mod(long long n) {
    if (n < 0) return 0;
    int64_t n_mod = n % MOD;
    return mod_mul(mod_mul(n_mod, (n_mod + 1) % MOD), INV2);
}

int64_t sum2_mod(long long n) {
    if (n < 0) return 0;
    int64_t n_mod = n % MOD;
    int64_t n1 = (n_mod + 1) % MOD;
    int64_t two_n1 = (mod_mul(2 % MOD, n_mod) + 1) % MOD;
    return mod_mul(mod_mul(mod_mul(n_mod, n1), two_n1), INV6);
}

int64_t sum1_range_mod(long long l, long long r) {
    if (l > r) return 0;
    return mod_sub(sum1_mod(r), sum1_mod(l - 1));
}

int64_t sum2_range_mod(long long l, long long r) {
    if (l > r) return 0;
    return mod_sub(sum2_mod(r), sum2_mod(l - 1));
}

// Sum_{s=l..r} (s+1)(N-s+1).
int64_t sum_pair_mod(long long l, long long r, long long N) {
    if (l > r) return 0;
    long long cnt = r - l + 1;
    int64_t S1 = sum1_range_mod(l, r);
    int64_t S2 = sum2_range_mod(l, r);
    int64_t sum_s_plus = mod_add(S1, cnt % MOD);
    int64_t sum_s2_plus = mod_add(S2, S1);
    return mod_sub(mod_mul((N + 1) % MOD, sum_s_plus), sum_s2_plus);
}

// Sum_{s=l..r} (s+1)(s+2)/2.
int64_t sum_tri_mod(long long l, long long r) {
    if (l > r) return 0;
    long long cnt = r - l + 1;
    int64_t S1 = sum1_range_mod(l, r);
    int64_t S2 = sum2_range_mod(l, r);
    int64_t total = mod_add(S2, mod_add(mod_mul(3, S1), mod_mul(2, cnt % MOD)));
    return mod_mul(total, INV2);
}

// Sum_{a=l..r} (D-a+1)^2.
int64_t sum_M1_sq_mod(long long l, long long r, long long D) {
    if (l > r) return 0;
    long long u_low = D - r + 1;
    long long u_high = D - l + 1;
    return sum2_range_mod(u_low, u_high);
}

// Sum of t(t+1)/2 for arithmetic progression: t = t0 + step*i, i=0..n-1.
int64_t sum_t_t1_mod(long long t0, long long step, long long n) {
    if (n <= 0) return 0;
    int64_t n_mod = n % MOD;
    int64_t n1_mod = (n - 1) % MOD;
    if (n1_mod < 0) n1_mod += MOD;
    int64_t sum_i = mod_mul(mod_mul(n_mod, n1_mod), INV2);
    int64_t two_n1 = (mod_mul(2 % MOD, n1_mod) + 1) % MOD;
    int64_t sum_i2 = mod_mul(mod_mul(mod_mul(n_mod, n1_mod), two_n1), INV6);

    int64_t t0_mod = t0 % MOD;
    if (t0_mod < 0) t0_mod += MOD;
    int64_t step_mod = step % MOD;
    if (step_mod < 0) step_mod += MOD;

    int64_t sum_t = mod_add(mod_mul(n_mod, t0_mod), mod_mul(step_mod, sum_i));
    int64_t t0_sq = mod_mul(t0_mod, t0_mod);
    int64_t term1 = mod_mul(n_mod, t0_sq);
    int64_t term2 = mod_mul(mod_mul(mod_mul(2, t0_mod), step_mod), sum_i);
    int64_t term3 = mod_mul(mod_mul(step_mod, step_mod), sum_i2);
    int64_t sum_t2 = mod_add(term1, mod_add(term2, term3));

    return mod_mul(mod_add(sum_t2, sum_t), INV2);
}

int64_t T_mod(long long x) {
    if (x < 0) return 0;
    int64_t xm = x % MOD;
    return mod_mul(mod_mul(mod_mul((xm + 1) % MOD, (xm + 2) % MOD), (xm + 3) % MOD), INV6);
}

// Count triples 0<=a,b,c<=D, a+b+c<=N, and a+b>D, a+c>D, b+c>D.
int64_t count_within_cube_mod(long long N, long long D) {
    if (N < 0) return 0;
    if (N > 2 * D) N = 2 * D;

    int64_t U;
    if (N <= D) {
        U = T_mod(N);
    } else {
        U = mod_sub(T_mod(N), mod_mul(3, T_mod(N - D - 1)));
    }

    int64_t A;
    if (N <= D) {
        A = sum_pair_mod(0, N, N);
    } else {
        long long s1 = N - D;
        int64_t part1 = mod_mul((D + 1) % MOD,
                                mod_mul((s1 + 1) % MOD, (s1 + 2) % MOD) * INV2 % MOD);
        int64_t part2 = sum_pair_mod(s1 + 1, D, N);
        A = mod_add(part1, part2);
    }

    int64_t AB;
    if (N <= D) {
        AB = sum_tri_mod(0, N);
    } else {
        long long a0 = 2 * D - N;
        int64_t mid1 = sum_M1_sq_mod(0, a0 - 1, D);
        int64_t mid2 = sum_t_t1_mod(a0, -1, a0);
        int64_t mid = mod_sub(mid1, mid2);
        int64_t full = sum_M1_sq_mod(a0, D, D);
        AB = mod_add(mid, full);
    }

    int64_t ABC;
    if (N <= D) {
        ABC = sum_tri_mod(0, N);
    } else {
        long long A1 = N - D;
        long long a_half = D / 2;
        int64_t total = 0;

        long long r1_end = min(A1, a_half - 1);
        if (r1_end >= 0) {
            int64_t mid1 = sum_M1_sq_mod(0, r1_end, D);
            long long n = r1_end + 1;
            int64_t mid2 = sum_t_t1_mod(D, -2, n);
            total = mod_add(total, mod_sub(mid1, mid2));
        }

        long long r1_full_start = max(0LL, a_half);
        if (r1_full_start <= A1) {
            total = mod_add(total, sum_M1_sq_mod(r1_full_start, A1, D));
        }

        long long l2 = A1 + 1;
        if (l2 <= D) {
            long long a0 = 2 * D - N;
            long long r2_end = min(a0 - 1, D);
            if (l2 <= r2_end) {
                int64_t mid1 = sum_M1_sq_mod(l2, r2_end, D);
                long long n = r2_end - l2 + 1;
                int64_t mid2 = sum_t_t1_mod(a0 - l2, -1, n);
                total = mod_add(total, mod_sub(mid1, mid2));
            }

            long long l_full = max(l2, a0);
            if (l_full <= D) {
                total = mod_add(total, sum_M1_sq_mod(l_full, D, D));
            }
        }

        ABC = total;
    }

    int64_t res = mod_sub(mod_add(mod_sub(U, mod_mul(3, A)), mod_mul(3, AB)), ABC);
    return res;
}

vector<unordered_map<long long, int64_t>> memo;

int64_t hard_count_mod(int k, long long N) {
    if (N < 0) return 0;
    if (k == 1) {
        if (N <= 2) return 0;
        return mod_sub(T_mod(N), T_mod(2));
    }
    if (N > (1LL << k)) N = 1LL << k;
    auto& mp = memo[k];
    auto it = mp.find(N);
    if (it != mp.end()) return it->second;

    long long D = 1LL << (k - 1);
    int64_t res;
    if (N <= D) {
        res = count_within_cube_mod(N, D);
    } else {
        res = mod_add(count_within_cube_mod(N, D),
                      mod_mul(3, hard_count_mod(k - 1, N - D)));
    }
    mp.emplace(N, res);
    return res;
}

int64_t sum_C_range(long long L, long long R) {
    if (L > R) return 0;
    long long cnt = R - L + 1;
    int64_t S1 = sum1_range_mod(L, R);
    int64_t S2 = sum2_range_mod(L, R);
    int64_t total = mod_add(S2, mod_add(mod_mul(3, S1), mod_mul(2, cnt % MOD)));
    return mod_mul(total, INV2);
}

int64_t compute_H(long long N) {
    if (N <= 1) return 0;
    int max_k = 64 - __builtin_clzll(static_cast<unsigned long long>(N - 1));

    int64_t base = 0;
    for (int k = 1; k <= max_k; ++k) {
        long long L = (1LL << (k - 1)) + 1;
        long long R = min(N, 1LL << k);
        if (L > R) continue;
        int64_t sumC = sum_C_range(L, R);
        base = (base + mod_mul(k % MOD, sumC)) % MOD;
    }

    memo.assign(max_k + 1, {});
    int64_t extra = 0;
    for (int k = 1; k <= max_k; ++k) {
        long long X = min(N, 1LL << k);
        extra = mod_add(extra, hard_count_mod(k, X));
    }

    return mod_add(base, extra);
}

void run_validation() {
    struct Case { long long N; long long expected; };
    const Case cases[] = {
        {6, 203},
        {20, 7718},
        {111, 1634144},
    };
    for (const auto& c : cases) {
        int64_t got = compute_H(c.N);
        if (got != (c.expected % MOD)) {
            cerr << "Validation failed for N=" << c.N
                 << ": got " << got << ", expected " << c.expected << "\n";
            exit(1);
        }
    }
}

} // namespace

int main(int argc, char** argv) {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    long long N = 0;
    for (int i = 0; i < 19; ++i) N = N * 10 + 1; // R_19
    bool validate = true;

    // Optional CLI: ./a.out [N] [validate(0/1)]
    if (argc >= 2) N = stoll(argv[1]);
    if (argc >= 3) validate = (stoi(argv[2]) != 0);

    if (validate) run_validation();

    cout << compute_H(N) << "\n";
    return 0;
}

Python

kMod = 1000000007

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

def mod_sub(a, b):
    return (a - b + kMod) % kMod

def mod_mul(a, b):
    return (a * b) % kMod

def mod_pow(a, e):
    return pow(a, e, kMod)

INV2 = mod_pow(2, kMod - 2)
INV6 = mod_pow(6, kMod - 2)

def sum1_mod(n):
    if n < 0: return 0
    nm = n % kMod
    return mod_mul(mod_mul(nm, (nm + 1) % kMod), INV2)

def sum2_mod(n):
    if n < 0: return 0
    nm = n % kMod
    n1 = (nm + 1) % kMod
    two_n1 = (mod_mul(2, nm) + 1) % kMod
    return mod_mul(mod_mul(mod_mul(nm, n1), two_n1), INV6)

def sum1_range_mod(l, r):
    if l > r: return 0
    return mod_sub(sum1_mod(r), sum1_mod(l - 1))

def sum2_range_mod(l, r):
    if l > r: return 0
    return mod_sub(sum2_mod(r), sum2_mod(l - 1))

def sum_pair_mod(l, r, N):
    if l > r: return 0
    cnt = r - l + 1
    S1 = sum1_range_mod(l, r)
    S2 = sum2_range_mod(l, r)
    sum_s_plus = mod_add(S1, cnt % kMod)
    sum_s2_plus = mod_add(S2, S1)
    return mod_sub(mod_mul((N + 1) % kMod, sum_s_plus), sum_s2_plus)

def sum_tri_mod(l, r):
    if l > r: return 0
    cnt = r - l + 1
    S1 = sum1_range_mod(l, r)
    S2 = sum2_range_mod(l, r)
    total = mod_add(S2, mod_add(mod_mul(3, S1), mod_mul(2, cnt % kMod)))
    return mod_mul(total, INV2)

def sum_M1_sq_mod(l, r, D):
    if l > r: return 0
    u_low = D - r + 1
    u_high = D - l + 1
    return sum2_range_mod(u_low, u_high)

def sum_t_t1_mod(t0, step, n):
    if n <= 0: return 0
    n_mod = n % kMod
    n1_mod = (n - 1) % kMod
    sum_i = mod_mul(mod_mul(n_mod, n1_mod), INV2)
    two_n1 = (mod_mul(2, n1_mod) + 1) % kMod
    sum_i2 = mod_mul(mod_mul(mod_mul(n_mod, n1_mod), two_n1), INV6)
    
    t0_mod = t0 % kMod
    step_mod = step % kMod
    
    sum_t = mod_add(mod_mul(n_mod, t0_mod), mod_mul(step_mod, sum_i))
    t0_sq = mod_mul(t0_mod, t0_mod)
    term1 = mod_mul(n_mod, t0_sq)
    term2 = mod_mul(mod_mul(mod_mul(2, t0_mod), step_mod), sum_i)
    term3 = mod_mul(mod_mul(step_mod, step_mod), sum_i2)
    sum_t2 = mod_add(term1, mod_add(term2, term3))
    
    return mod_mul(mod_add(sum_t2, sum_t), INV2)

def T_mod(x):
    if x < 0: return 0
    xm = x % kMod
    return mod_mul(mod_mul(mod_mul((xm + 1) % kMod, (xm + 2) % kMod), (xm + 3) % kMod), INV6)

def count_within_cube_mod(N, D):
    if N < 0: return 0
    if N > 2 * D: N = 2 * D
    
    if N <= D:
        U = T_mod(N)
    else:
        U = mod_sub(T_mod(N), mod_mul(3, T_mod(N - D - 1)))
        
    if N <= D:
        A = sum_pair_mod(0, N, N)
    else:
        s1 = N - D
        part1 = mod_mul((D + 1) % kMod, mod_mul(mod_mul((s1 + 1) % kMod, (s1 + 2) % kMod), INV2))
        part2 = sum_pair_mod(s1 + 1, D, N)
        A = mod_add(part1, part2)
        
    if N <= D:
        AB = sum_tri_mod(0, N)
    else:
        a0 = 2 * D - N
        mid1 = sum_M1_sq_mod(0, a0 - 1, D)
        mid2 = sum_t_t1_mod(a0, -1, a0)
        mid = mod_sub(mid1, mid2)
        full = sum_M1_sq_mod(a0, D, D)
        AB = mod_add(mid, full)
        
    if N <= D:
        ABC = sum_tri_mod(0, N)
    else:
        A1 = N - D
        a_half = D // 2
        total = 0
        
        r1_end = min(A1, a_half - 1)
        if r1_end >= 0:
            mid1 = sum_M1_sq_mod(0, r1_end, D)
            n = r1_end + 1
            mid2 = sum_t_t1_mod(D, -2, n)
            total = mod_add(total, mod_sub(mid1, mid2))
            
        r1_full_start = max(0, a_half)
        if r1_full_start <= A1:
            total = mod_add(total, sum_M1_sq_mod(r1_full_start, A1, D))
            
        l2 = A1 + 1
        if l2 <= D:
            a0 = 2 * D - N
            r2_end = min(a0 - 1, D)
            if l2 <= r2_end:
                mid1 = sum_M1_sq_mod(l2, r2_end, D)
                n = r2_end - l2 + 1
                mid2 = sum_t_t1_mod(a0 - l2, -1, n)
                total = mod_add(total, mod_sub(mid1, mid2))
                
            l_full = max(l2, a0)
            if l_full <= D:
                total = mod_add(total, sum_M1_sq_mod(l_full, D, D))
                
        ABC = total
        
    return mod_sub(mod_add(mod_sub(U, mod_mul(3, A)), mod_mul(3, AB)), ABC)

memo = {}

def hard_count_mod(k, N):
    if N < 0: return 0
    if k == 1:
        if N <= 2: return 0
        return mod_sub(T_mod(N), T_mod(2))
    if N > (1 << k): N = 1 << k
    
    key = (k, N)
    if key in memo: return memo[key]
    
    D = 1 << (k - 1)
    if N <= D:
        res = count_within_cube_mod(N, D)
    else:
        res = mod_add(count_within_cube_mod(N, D), mod_mul(3, hard_count_mod(k - 1, N - D)))
        
    memo[key] = res
    return res

def sum_C_range(L, R):
    if L > R: return 0
    cnt = R - L + 1
    S1 = sum1_range_mod(L, R)
    S2 = sum2_range_mod(L, R)
    total = mod_add(S2, mod_add(mod_mul(3, S1), mod_mul(2, cnt % kMod)))
    return mod_mul(total, INV2)

def compute_H(N):
    if N <= 1: return 0
    max_k = (N - 1).bit_length()
    
    base = 0
    for k in range(1, max_k + 1):
        L = (1 << (k - 1)) + 1
        R = min(N, 1 << k)
        if L > R: continue
        sumC = sum_C_range(L, R)
        base = mod_add(base, mod_mul(k % kMod, sumC))
        
    global memo
    memo.clear()
    
    extra = 0
    for k in range(1, max_k + 1):
        X = min(N, 1 << k)
        extra = mod_add(extra, hard_count_mod(k, X))
        
    return mod_add(base, extra)

def solve():
    N = int("1" * 19)
    ans = compute_H(N)
    return str(ans)

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

Java

import java.util.HashMap;

public class Euler847 {
    static final long MOD = 1000000007L;

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

    static long modSub(long a, long b) {
        long r = a - b;
        if (r < 0)
            r += MOD;
        return r;
    }

    static long modMul(long a, long b) {
        long aMod = (a % MOD + MOD) % MOD;
        long bMod = (b % MOD + MOD) % MOD;
        return (aMod * bMod) % MOD;
    }

    static long modPow(long a, long e) {
        long r = 1 % MOD;
        long x = (a % MOD + MOD) % MOD;
        while (e > 0) {
            if ((e & 1) == 1)
                r = modMul(r, x);
            x = modMul(x, x);
            e >>= 1;
        }
        return r;
    }

    static final long INV2 = modPow(2, MOD - 2);
    static final long INV6 = modPow(6, MOD - 2);

    static long sum1Mod(long n) {
        if (n < 0)
            return 0;
        long nMod = n % MOD;
        return modMul(modMul(nMod, (nMod + 1) % MOD), INV2);
    }

    static long sum2Mod(long n) {
        if (n < 0)
            return 0;
        long nMod = n % MOD;
        long n1 = (nMod + 1) % MOD;
        long twoN1 = (modMul(2, nMod) + 1) % MOD;
        return modMul(modMul(modMul(nMod, n1), twoN1), INV6);
    }

    static long sum1RangeMod(long l, long r) {
        if (l > r)
            return 0;
        return modSub(sum1Mod(r), sum1Mod(l - 1));
    }

    static long sum2RangeMod(long l, long r) {
        if (l > r)
            return 0;
        return modSub(sum2Mod(r), sum2Mod(l - 1));
    }

    static long sumPairMod(long l, long r, long N) {
        if (l > r)
            return 0;
        long cnt = r - l + 1;
        long S1 = sum1RangeMod(l, r);
        long S2 = sum2RangeMod(l, r);
        long sumSPlus = modAdd(S1, cnt % MOD);
        long sumS2Plus = modAdd(S2, S1);
        return modSub(modMul((N + 1) % MOD, sumSPlus), sumS2Plus);
    }

    static long sumTriMod(long l, long r) {
        if (l > r)
            return 0;
        long cnt = r - l + 1;
        long S1 = sum1RangeMod(l, r);
        long S2 = sum2RangeMod(l, r);
        long total = modAdd(S2, modAdd(modMul(3, S1), modMul(2, cnt % MOD)));
        return modMul(total, INV2);
    }

    static long sumM1SqMod(long l, long r, long D) {
        if (l > r)
            return 0;
        long uLow = D - r + 1;
        long uHigh = D - l + 1;
        return sum2RangeMod(uLow, uHigh);
    }

    static long sumTT1Mod(long t0, long step, long n) {
        if (n <= 0)
            return 0;
        long nMod = n % MOD;
        long n1Mod = (n - 1) % MOD;
        if (n1Mod < 0)
            n1Mod += MOD;
        long sumI = modMul(modMul(nMod, n1Mod), INV2);
        long twoN1 = (modMul(2, n1Mod) + 1) % MOD;
        long sumI2 = modMul(modMul(modMul(nMod, n1Mod), twoN1), INV6);

        long t0Mod = (t0 % MOD + MOD) % MOD;
        long stepMod = (step % MOD + MOD) % MOD;

        long sumT = modAdd(modMul(nMod, t0Mod), modMul(stepMod, sumI));
        long t0Sq = modMul(t0Mod, t0Mod);
        long term1 = modMul(nMod, t0Sq);
        long term2 = modMul(modMul(modMul(2, t0Mod), stepMod), sumI);
        long term3 = modMul(modMul(stepMod, stepMod), sumI2);
        long sumT2 = modAdd(term1, modAdd(term2, term3));

        return modMul(modAdd(sumT2, sumT), INV2);
    }

    static long tMod(long x) {
        if (x < 0)
            return 0;
        long xm = x % MOD;
        return modMul(modMul(modMul((xm + 1) % MOD, (xm + 2) % MOD), (xm + 3) % MOD), INV6);
    }

    static long countWithinCubeMod(long N, long D) {
        if (N < 0)
            return 0;
        if (N > 2 * D)
            N = 2 * D;

        long U;
        if (N <= D) {
            U = tMod(N);
        } else {
            U = modSub(tMod(N), modMul(3, tMod(N - D - 1)));
        }

        long A;
        if (N <= D) {
            A = sumPairMod(0, N, N);
        } else {
            long s1 = N - D;
            long part1 = modMul((D + 1) % MOD, modMul(modMul((s1 + 1) % MOD, (s1 + 2) % MOD), INV2));
            long part2 = sumPairMod(s1 + 1, D, N);
            A = modAdd(part1, part2);
        }

        long AB;
        if (N <= D) {
            AB = sumTriMod(0, N);
        } else {
            long a0 = 2 * D - N;
            long mid1 = sumM1SqMod(0, a0 - 1, D);
            long mid2 = sumTT1Mod(a0, -1, a0);
            long mid = modSub(mid1, mid2);
            long full = sumM1SqMod(a0, D, D);
            AB = modAdd(mid, full);
        }

        long ABC;
        if (N <= D) {
            ABC = sumTriMod(0, N);
        } else {
            long A1 = N - D;
            long aHalf = D / 2;
            long total = 0;

            long r1End = Math.min(A1, aHalf - 1);
            if (r1End >= 0) {
                long mid1 = sumM1SqMod(0, r1End, D);
                long n = r1End + 1;
                long mid2 = sumTT1Mod(D, -2, n);
                total = modAdd(total, modSub(mid1, mid2));
            }

            long r1FullStart = Math.max(0L, aHalf);
            if (r1FullStart <= A1) {
                total = modAdd(total, sumM1SqMod(r1FullStart, A1, D));
            }

            long l2 = A1 + 1;
            if (l2 <= D) {
                long a0 = 2 * D - N;
                long r2End = Math.min(a0 - 1, D);
                if (l2 <= r2End) {
                    long mid1 = sumM1SqMod(l2, r2End, D);
                    long n = r2End - l2 + 1;
                    long mid2 = sumTT1Mod(a0 - l2, -1, n);
                    total = modAdd(total, modSub(mid1, mid2));
                }

                long lFull = Math.max(l2, a0);
                if (lFull <= D) {
                    total = modAdd(total, sumM1SqMod(lFull, D, D));
                }
            }
            ABC = total;
        }

        return modSub(modAdd(modSub(U, modMul(3, A)), modMul(3, AB)), ABC);
    }

    static HashMap<String, Long> memo = new HashMap<>();

    static long hardCountMod(int k, long N) {
        if (N < 0)
            return 0;
        if (k == 1) {
            if (N <= 2)
                return 0;
            return modSub(tMod(N), tMod(2));
        }
        if (N > (1L << k))
            N = 1L << k;

        String key = k + "_" + N;
        if (memo.containsKey(key))
            return memo.get(key);

        long D = 1L << (k - 1);
        long res;
        if (N <= D) {
            res = countWithinCubeMod(N, D);
        } else {
            res = modAdd(countWithinCubeMod(N, D), modMul(3, hardCountMod(k - 1, N - D)));
        }

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

    static long sumCRange(long L, long R) {
        if (L > R)
            return 0;
        long cnt = R - L + 1;
        long S1 = sum1RangeMod(L, R);
        long S2 = sum2RangeMod(L, R);
        long total = modAdd(S2, modAdd(modMul(3, S1), modMul(2, cnt % MOD)));
        return modMul(total, INV2);
    }

    static long computeH(long N) {
        if (N <= 1)
            return 0;
        int maxK = 64 - Long.numberOfLeadingZeros(N - 1);

        long base = 0;
        for (int k = 1; k <= maxK; ++k) {
            long L = (1L << (k - 1)) + 1;
            long R = Math.min(N, 1L << k);
            if (L > R)
                continue;
            long sumC = sumCRange(L, R);
            base = modAdd(base, modMul(k % MOD, sumC));
        }

        memo.clear();
        long extra = 0;
        for (int k = 1; k <= maxK; ++k) {
            long X = Math.min(N, 1L << k);
            extra = modAdd(extra, hardCountMod(k, X));
        }

        return modAdd(base, extra);
    }

    public static String solve() {
        long N = 0;
        for (int i = 0; i < 19; ++i) {
            N = N * 10 + 1;
        }
        return Long.toString(computeH(N));
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}