Problem 889: Rational Blancmange

Project Euler Problem 889 Solution

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

Problem Summary For fixed integers $k$, $t$, and $r$, define $$\alpha = 2^t + 1,\qquad \beta = 2^k + 1,\qquad M = 10^9 + 62031.$$ For each integer $n$ with $0 \le n < k$, let $R_n$ be the least nonnegative residue of $2^n\alpha^r$ modulo $\beta$, and let $$D_n = \min(R_n,\beta - R_n).$$ The required value is $$F(k,t,r)=\sum_{n=0}^{k-1} 2^{k-n} D_n \pmod{M}.$$ The hard part is that $k$ and $t$ are enormous, so iterating through all $n$ directly is impossible. The solution works by expanding $\alpha^r$, folding powers across the modulus $2^k+1$, and proving that almost all $n$-values inside each interval contribute the same weighted amount. Mathematical Approach The implementations use the target value $r=62$, and the given parameters satisfy $rt<k$. That inequality is crucial: after expansion, every exponent lies below $2k$, so any term can cross the $k$-boundary at most once. Step 1: Expand $(2^t+1)^r$ and track the exponents By the binomial theorem, $$\alpha^r=(1+2^t)^r=\sum_{i=0}^{r}\binom{r}{i}2^{it}.$$ Multiplying by $2^n$ gives $$2^n\alpha^r=\sum_{i=0}^{r}\binom{r}{i}2^{n+it}.$$ So the entire problem is controlled by the exponents $$e_i(n)=n+it.$$ Because $0 \le n < k$ and $rt<k$, every $e_i(n)$ lies in the half-open range $[0,2k)$....

Detailed mathematical approach

Problem Summary

For fixed integers $k$, $t$, and $r$, define

$$\alpha = 2^t + 1,\qquad \beta = 2^k + 1,\qquad M = 10^9 + 62031.$$

For each integer $n$ with $0 \le n < k$, let $R_n$ be the least nonnegative residue of $2^n\alpha^r$ modulo $\beta$, and let

$$D_n = \min(R_n,\beta - R_n).$$

The required value is

$$F(k,t,r)=\sum_{n=0}^{k-1} 2^{k-n} D_n \pmod{M}.$$

The hard part is that $k$ and $t$ are enormous, so iterating through all $n$ directly is impossible. The solution works by expanding $\alpha^r$, folding powers across the modulus $2^k+1$, and proving that almost all $n$-values inside each interval contribute the same weighted amount.

Mathematical Approach

The implementations use the target value $r=62$, and the given parameters satisfy $rt<k$. That inequality is crucial: after expansion, every exponent lies below $2k$, so any term can cross the $k$-boundary at most once.

Step 1: Expand $(2^t+1)^r$ and track the exponents

By the binomial theorem,

$$\alpha^r=(1+2^t)^r=\sum_{i=0}^{r}\binom{r}{i}2^{it}.$$

Multiplying by $2^n$ gives

$$2^n\alpha^r=\sum_{i=0}^{r}\binom{r}{i}2^{n+it}.$$

So the entire problem is controlled by the exponents

$$e_i(n)=n+it.$$

Because $0 \le n < k$ and $rt<k$, every $e_i(n)$ lies in the half-open range $[0,2k)$. Therefore each term is either already below $k$, or it exceeds $k$ once and can be folded back using the congruence $2^k\equiv -1 \pmod{\beta}$.

Step 2: Split the range of $n$ into intervals with one fixed folding pattern

For $s=1,2,\dots,r$, define

$$I_s=\left[k-st,\ k-(s-1)t-1\right]\cap [0,k-1].$$

There is also a last interval

$$I_{r+1}=[0,\ k-rt-1]\cap [0,k-1],$$

which may be empty for some parameters.

If $n\in I_s$, then exactly the terms with $i\ge s$ satisfy $e_i(n)\ge k$, while the terms with $i<s$ stay below $k$. Hence

$$2^{e_i(n)} \equiv \begin{cases} 2^{n+it}, & i<s,\\[4pt] -2^{n+it-k}, & i\ge s, \end{cases}\pmod{\beta}$$

and therefore

$$R_n \equiv E_s(n):=\sum_{i=0}^{s-1}\binom{r}{i}2^{n+it}-\sum_{i=s}^{r}\binom{r}{i}2^{n+it-k}\pmod{\beta}.$$

On the last interval $I_{r+1}$, nothing folds, so only the first sum remains. This interval decomposition is the key structural simplification: instead of treating each $n$ separately, we only need to understand one pattern per interval.

Step 3: After weighting by $2^{k-n}$, the bulk contribution of an interval is constant

Suppose that for some $n\in I_s$ the centered representative is simply $D_n=E_s(n)$. Then multiplying by the outer weight gives

$$2^{k-n}D_n=\sum_{i=0}^{s-1}\binom{r}{i}2^{k+it}-\sum_{i=s}^{r}\binom{r}{i}2^{it}.$$

The right-hand side no longer depends on $n$. So every non-exceptional $n$ in the same interval contributes the same weighted amount.

For the last interval, the constant becomes

$$2^{k-n}D_n=\sum_{i=0}^{r}\binom{r}{i}2^{k+it},$$

again independent of $n$, whenever the centered residue is taken without switching to $\beta-R_n$.

This is why the program can add most of the contribution of a huge interval in one bulk step instead of visiting every index inside it.

Step 4: Only the last 62 positions of an interval can be ambiguous

Fix an interval $I_s$, and let the leading positive exponent be

$$e_{\max}=n+(s-1)t$$

for $s\le r$, or $e_{\max}=n+rt$ on the final interval. Write

$$\delta = k-e_{\max}.$$

When $\delta \ge 63$, every other term in the signed sum is at least $2^{63}$ times smaller than the boundary scale $2^k$, while the total coefficient mass is bounded by

$$\sum_{i=0}^{r}\binom{r}{i}=2^r=2^{62}.$$

So the leading term dominates and the entire signed value sits strictly between $0$ and $\beta/2$. In that regime there is no need to compare $R_n$ with $\beta-R_n$: the nearer value is automatically the unreversed one.

Consequently, only the values with $\delta\le 62$ need exact centered-remainder logic. Those are exactly the last $62$ positions of each nonempty interval, which is why the implementations peel off that short suffix and treat the earlier part as bulk.

Step 5: Recover the exact boundary cases from the leading binomial term

For an exceptional $n$, let $j=s-1$ for $s\le r$, and $j=r$ on the last interval. The dominant term is

$$\binom{r}{j}2^{k-\delta}.$$

Write the coefficient in base $2^\delta$:

$$\binom{r}{j}=q\,2^\delta+\ell,\qquad 0\le \ell < 2^\delta.$$

Then

$$\binom{r}{j}2^{k-\delta}=q2^k+\ell 2^{k-\delta}=q(2^k+1)+\ell 2^{k-\delta}-q.$$

This identity shows that the exact quotient by $\beta=2^k+1$ can be read off from the bit shift

$$q=\left\lfloor \frac{\binom{r}{j}}{2^\delta}\right\rfloor.$$

After subtracting $q\beta$, the low bits $\ell$ determine on which side of $\beta/2$ the remainder lies. If $\ell<2^{\delta-1}$, the remainder is below half the modulus. If $\ell>2^{\delta-1}$, the complementary distance $\beta-R_n$ is smaller. In the tie case $\ell=2^{\delta-1}$, the lower-order terms decide the outcome. This is exactly what the implementations evaluate for the short exceptional suffix of each interval.

Worked Example: $F(3,1,1)=42$

Here

$$\alpha=2^1+1=3,\qquad \beta=2^3+1=9.$$

Since $r=1$, we only need the values $n=0,1,2$:

$$\begin{aligned} n=0&:&&2^0\alpha^1=3\equiv 3 \pmod 9,\quad D_0=\min(3,6)=3,\quad 2^{3-0}D_0=24,\\ n=1&:&&2^1\alpha^1=6\equiv 6 \pmod 9,\quad D_1=\min(6,3)=3,\quad 2^{3-1}D_1=12,\\ n=2&:&&2^2\alpha^1=12\equiv 3 \pmod 9,\quad D_2=\min(3,6)=3,\quad 2^{3-2}D_2=6. \end{aligned}$$

Therefore

$$F(3,1,1)=24+12+6=42,$$

which matches the validation value used by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical plan. They first precompute the binomial coefficients $\binom{r}{i}$, their values modulo $M$, and the modular powers $2^{it}\pmod M$. They also compute $2^k\pmod M$ and the modular inverse of $2^k$ modulo $M$, because folded terms naturally introduce factors of $2^{-k}$ when everything is tracked directly modulo $M$.

Next, the implementation forms prefix sums for the unfolded part and suffix sums for the folded part. That makes the constant bulk contribution of an interval available in $O(1)$ modular time once the interval length is known. Only the last $62$ values of each interval are handled one by one.

For those exceptional positions, the implementation rebuilds the signed folded expression, extracts the quotient $q$ from the dominant binomial coefficient by a bit shift, subtracts the corresponding multiple of $\beta$, decides whether the remainder or its complement is closer to zero, and finally multiplies by the outer factor $2^{k-n}$. The C++ version parallelizes these independent exceptional evaluations; the Python and Java versions evaluate the same cases sequentially.

Complexity Analysis

There are only $r+1$ interval patterns. Each interval contributes one bulk term and at most $62$ exceptional positions. Each exceptional evaluation scans the $r+1$ binomial terms, so the arithmetic core uses $O(r^2)$ modular operations and $O(r)$ memory.

For the target problem $r=62$ is fixed, so the runtime is effectively constant with respect to the sizes of $k$ and $t$; it does not grow linearly with either parameter. The only remaining dependence on their magnitudes comes from modular exponentiation, which costs logarithmic time in the exponent values.

Footnotes and References

Problem page: Project Euler 889
Takagi curve (Blancmange curve): Wikipedia — Takagi curve
Binomial theorem: Wikipedia — Binomial theorem
Modular arithmetic: Wikipedia — Modular arithmetic
Exponentiation by squaring: Wikipedia — Exponentiation by squaring

Problem 889 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <thread>
#include <vector>
#include <boost/multiprecision/cpp_int.hpp>

using namespace std;
using boost::multiprecision::cpp_int;

namespace {

constexpr int64_t MOD = 1'000'062'031LL;

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 base, uint64_t exp) {
    int64_t result = 1 % MOD;
    int64_t cur = base % MOD;
    while (exp > 0) {
        if (exp & 1) result = mod_mul(result, cur);
        cur = mod_mul(cur, cur);
        exp >>= 1;
    }
    return result;
}

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

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

int64_t brute_mod(int64_t k, int64_t t, int r) {
    cpp_int A = (cpp_int(1) << t) + 1;
    cpp_int Apow = 1;
    for (int i = 0; i < r; ++i) {
        Apow *= A;
    }
    cpp_int B = (cpp_int(1) << k) + 1;
    cpp_int total = 0;
    for (int64_t n = 0; n < k; ++n) {
        cpp_int rn = (Apow << n) % B;
        cpp_int dn = rn;
        cpp_int comp = B - rn;
        if (comp < dn) dn = comp;
        total += dn << (k - n);
    }
    int64_t result = (total % MOD).convert_to<long long>();
    if (result < 0) result += MOD;
    return result;
}

struct FastContext {
    int r = 0;
    int64_t k = 0;
    int64_t t = 0;
    vector<uint64_t> binom;
    vector<int64_t> binom_mod;
    vector<int64_t> pow2_ti;
    int64_t pow2_k = 0;
    int64_t pow2_k_inv = 0;
    int64_t b_mod = 0;
    vector<int64_t> prefix_k;
    vector<int64_t> suffix;
};

struct Task {
    int64_t n = 0;
    int s = 0;
};

int64_t exception_contrib(const FastContext& ctx, int64_t n, int s) {
    int64_t pow2_n = mod_pow(2, static_cast<uint64_t>(n));
    int64_t S_mod = 0;
    for (int i = 0; i <= ctx.r; ++i) {
        int64_t base = mod_mul(ctx.binom_mod[i], mod_mul(pow2_n, ctx.pow2_ti[i]));
        if (i < s) {
            S_mod = mod_add(S_mod, base);
        } else {
            int64_t term = mod_mul(base, ctx.pow2_k_inv);
            S_mod = mod_sub(S_mod, term);
        }
    }

    int idx = (s <= ctx.r) ? (s - 1) : ctx.r;
    int64_t e_max = n + ctx.t * static_cast<int64_t>(idx);
    int64_t shift = ctx.k - e_max;
    uint64_t C_top = ctx.binom[idx];

    uint64_t Q = 0;
    if (shift >= 0 && shift < 63) {
        Q = C_top >> shift;
    }
    int64_t r_mod = mod_sub(S_mod, mod_mul(static_cast<int64_t>(Q % MOD), ctx.b_mod));

    bool greater = false;
    if (shift < 63) {
        uint64_t low_mask = (static_cast<uint64_t>(1) << shift) - 1;
        uint64_t low = C_top & low_mask;
        uint64_t half = static_cast<uint64_t>(1) << (shift - 1);
        if (low < half) {
            greater = false;
        } else if (low > half) {
            greater = true;
        } else {
            greater = (s >= 2);
        }
    }

    int64_t d_mod = greater ? mod_sub(ctx.b_mod, r_mod) : r_mod;
    int64_t pow2_kn = mod_pow(2, static_cast<uint64_t>(ctx.k - n));
    return mod_mul(d_mod, pow2_kn);
}

int64_t compute_fast(int64_t k, int64_t t, int r) {
    FastContext ctx;
    ctx.r = r;
    ctx.k = k;
    ctx.t = t;
    ctx.binom.assign(r + 1, 0);
    ctx.binom_mod.assign(r + 1, 0);
    ctx.pow2_ti.assign(r + 1, 0);
    ctx.prefix_k.assign(r + 2, 0);
    ctx.suffix.assign(r + 2, 0);

    ctx.binom[0] = 1;
    for (int i = 1; i <= r; ++i) {
        __int128 num = static_cast<__int128>(ctx.binom[i - 1]) * (r - i + 1);
        ctx.binom[i] = static_cast<uint64_t>(num / i);
    }
    for (int i = 0; i <= r; ++i) {
        ctx.binom_mod[i] = static_cast<int64_t>(ctx.binom[i] % MOD);
    }

    int64_t pow2_t = mod_pow(2, static_cast<uint64_t>(t));
    ctx.pow2_ti[0] = 1;
    for (int i = 1; i <= r; ++i) {
        ctx.pow2_ti[i] = mod_mul(ctx.pow2_ti[i - 1], pow2_t);
    }

    ctx.pow2_k = mod_pow(2, static_cast<uint64_t>(k));
    ctx.pow2_k_inv = mod_pow(ctx.pow2_k, MOD - 2);
    ctx.b_mod = mod_add(ctx.pow2_k, 1);

    vector<int64_t> term(r + 1, 0);
    vector<int64_t> term_k(r + 1, 0);
    for (int i = 0; i <= r; ++i) {
        term[i] = mod_mul(ctx.binom_mod[i], ctx.pow2_ti[i]);
        term_k[i] = mod_mul(term[i], ctx.pow2_k);
    }

    for (int i = 0; i <= r; ++i) {
        ctx.prefix_k[i + 1] = mod_add(ctx.prefix_k[i], term_k[i]);
    }
    for (int i = r; i >= 0; --i) {
        ctx.suffix[i] = mod_add(ctx.suffix[i + 1], term[i]);
    }

    vector<Task> tasks;
    tasks.reserve(static_cast<size_t>((r + 1) * 62));
    int64_t total = 0;

    auto add_bulk = [&](int64_t bulk, int64_t constant) {
        if (bulk <= 0) return;
        int64_t bulk_mod = static_cast<int64_t>(bulk % MOD);
        total = mod_add(total, mod_mul(bulk_mod, constant));
    };

    for (int s = 1; s <= r; ++s) {
        __int128 start_raw = static_cast<__int128>(k) - static_cast<__int128>(s) * t;
        __int128 end_raw = static_cast<__int128>(k) - static_cast<__int128>(s - 1) * t - 1;
        if (end_raw < 0) continue;
        int64_t start = start_raw > 0 ? static_cast<int64_t>(start_raw) : 0;
        int64_t end = static_cast<int64_t>(end_raw);
        if (end > k - 1) end = k - 1;
        if (start > end) continue;

        int64_t length = end - start + 1;
        int64_t exc_start = max(start, end - 61);
        int64_t exc_count = (exc_start <= end) ? (end - exc_start + 1) : 0;
        int64_t bulk = length - exc_count;

        int64_t constant = mod_sub(ctx.prefix_k[s], ctx.suffix[s]);
        add_bulk(bulk, constant);

        for (int64_t n = exc_start; n <= end; ++n) {
            tasks.push_back({n, s});
        }
    }

    {
        int s = r + 1;
        __int128 end_raw = static_cast<__int128>(k) - static_cast<__int128>(r) * t - 1;
        if (end_raw >= 0) {
            int64_t start = 0;
            int64_t end = static_cast<int64_t>(end_raw);
            if (end > k - 1) end = k - 1;
            if (start <= end) {
                int64_t length = end - start + 1;
                int64_t exc_start = max(start, end - 61);
                int64_t exc_count = (exc_start <= end) ? (end - exc_start + 1) : 0;
                int64_t bulk = length - exc_count;

                int64_t constant = ctx.prefix_k[r + 1];
                add_bulk(bulk, constant);

                for (int64_t n = exc_start; n <= end; ++n) {
                    tasks.push_back({n, s});
                }
            }
        }
    }

    size_t task_count = tasks.size();
    if (task_count > 0) {
        int thread_count = static_cast<int>(thread::hardware_concurrency());
        if (thread_count <= 0) thread_count = 1;
        if (thread_count > static_cast<int>(task_count)) {
            thread_count = static_cast<int>(task_count);
        }
        vector<int64_t> partial(thread_count, 0);
        vector<thread> workers;
        workers.reserve(thread_count);

        for (int t_idx = 0; t_idx < thread_count; ++t_idx) {
            workers.emplace_back([&, t_idx]() {
                int64_t acc = 0;
                for (size_t i = t_idx; i < task_count; i += thread_count) {
                    acc += exception_contrib(ctx, tasks[i].n, tasks[i].s);
                    if (acc >= MOD) acc -= MOD;
                }
                partial[t_idx] = acc;
            });
        }
        for (auto& th : workers) th.join();
        for (int t_idx = 0; t_idx < thread_count; ++t_idx) {
            total = mod_add(total, partial[t_idx]);
        }
    }

    return total % MOD;
}

}  // namespace

int main() {
    if (brute_mod(3, 1, 1) != 42) {
        cerr << "Validation failure: F(3,1,1)\n";
        return 1;
    }
    if (brute_mod(13, 3, 3) != 23093880) {
        cerr << "Validation failure: F(13,3,3)\n";
        return 1;
    }
    if (brute_mod(103, 13, 6) != 878922518) {
        cerr << "Validation failure: F(103,13,6)\n";
        return 1;
    }

    const int64_t k = 1'000'000'000'000'000'000LL + 31;
    const int64_t t = 100'000'000'000'000LL + 31;
    const int r = 62;

    cout << compute_fast(k, t, r) << '\n';
    return 0;
}

Python

def solve():
    MOD = 1000062031
    k = 1000000000000000031; t = 100000000000031; r = 62

    def mm(a, b): return a*b%MOD
    def mp(base, exp):
        r = 1; base %= MOD
        while exp > 0:
            if exp&1: r=r*base%MOD
            base=base*base%MOD; exp >>= 1
        return r
    def ma(a, b):
        v = a+b; return v-MOD if v>=MOD else v
    def ms(a, b):
        v = a-b; return v+MOD if v<0 else v

    binom = [0]*(r+1); binom[0] = 1
    for i in range(1, r+1): binom[i] = binom[i-1]*(r-i+1)//i
    binom_mod = [b%MOD for b in binom]

    pow2_t = mp(2, t)
    pow2_ti = [0]*(r+1); pow2_ti[0] = 1
    for i in range(1, r+1): pow2_ti[i] = mm(pow2_ti[i-1], pow2_t)
    pow2_k = mp(2, k); pow2_k_inv = mp(pow2_k, MOD-2); b_mod = ma(pow2_k, 1)

    term = [mm(binom_mod[i], pow2_ti[i]) for i in range(r+1)]
    term_k = [mm(term[i], pow2_k) for i in range(r+1)]
    prefix_k = [0]*(r+2); suffix = [0]*(r+2)
    for i in range(r+1): prefix_k[i+1] = ma(prefix_k[i], term_k[i])
    for i in range(r, -1, -1): suffix[i] = ma(suffix[i+1], term[i])

    def exc_contrib(n, s):
        p2n = mp(2, n); S = 0
        for i in range(r+1):
            base = mm(binom_mod[i], mm(p2n, pow2_ti[i]))
            if i < s: S = ma(S, base)
            else: S = ms(S, mm(base, pow2_k_inv))
        idx = min(s-1, r) if s <= r+1 else r
        e_max = n + t*idx; shift = k - e_max
        Q = 0
        if 0 <= shift < 63: Q = binom[idx] >> shift
        r_mod = ms(S, mm(Q%MOD, b_mod))
        greater = False
        if shift < 63:
            low_mask = (1<<shift)-1; low = binom[idx] & low_mask; half = 1<<(shift-1)
            if low < half: greater = False
            elif low > half: greater = True
            else: greater = s >= 2
        d_mod = ms(b_mod, r_mod) if greater else r_mod
        return mm(d_mod, mp(2, k-n))

    total = 0
    for s in range(1, r+1):
        start_raw = k - s*t; end_raw = k - (s-1)*t - 1
        if end_raw < 0: continue
        start = max(0, start_raw); end = min(end_raw, k-1)
        if start > end: continue
        length = end - start + 1
        exc_start = max(start, end-61); exc_count = end-exc_start+1 if exc_start<=end else 0
        bulk = length - exc_count
        if bulk > 0:
            constant = ms(prefix_k[s], suffix[s])
            total = ma(total, mm(bulk%MOD, constant))
        for n in range(exc_start, end+1): total = ma(total, exc_contrib(n, s))

    # s = r+1
    end_raw = k - r*t - 1
    if end_raw >= 0:
        start = 0; end = min(end_raw, k-1)
        if start <= end:
            length = end-start+1
            exc_start = max(start, end-61); exc_count = end-exc_start+1 if exc_start<=end else 0
            bulk = length - exc_count
            if bulk > 0: total = ma(total, mm(bulk%MOD, prefix_k[r+1]))
            for n in range(exc_start, end+1): total = ma(total, exc_contrib(n, r+1))

    return str(total%MOD)

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

Java

import java.math.BigInteger;
import java.util.ArrayList;

public class Euler889 {

    static final long MOD = 1000062031L;

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

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

    static long modMul(long a, long b) {
        return (a * b) % MOD;
    }

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

    static class FastContext {
        int r;
        long k;
        long t;
        long[] binom;
        long[] binomMod;
        long[] pow2Ti;
        long pow2K;
        long pow2KInv;
        long bMod;
        long[] prefixK;
        long[] suffix;
    }

    static class Task {
        long n;
        int s;

        Task(long n, int s) {
            this.n = n;
            this.s = s;
        }
    }

    static long exceptionContrib(FastContext ctx, long n, int s) {
        long pow2N = modPow(2, n);
        long SMod = 0;

        for (int i = 0; i <= ctx.r; ++i) {
            long base = modMul(ctx.binomMod[i], modMul(pow2N, ctx.pow2Ti[i]));
            if (i < s) {
                SMod = modAdd(SMod, base);
            } else {
                long term = modMul(base, ctx.pow2KInv);
                SMod = modSub(SMod, term);
            }
        }

        int idx = (s <= ctx.r) ? (s - 1) : ctx.r;
        long eMax = n + ctx.t * idx;
        long shift = ctx.k - eMax;
        long CTop = ctx.binom[idx];

        long Q = 0;
        if (shift >= 0 && shift < 63) {
            Q = CTop >>> shift;
        }
        long rMod = modSub(SMod, modMul(Q % MOD, ctx.bMod));

        boolean greater = false;
        if (shift >= 0 && shift < 63) {
            long lowMask = (1L << shift) - 1;
            long low = CTop & lowMask;
            long half = (shift == 0) ? 0 : (1L << (shift - 1));

            if (low < half) {
                greater = false;
            } else if (low > half) {
                greater = true;
            } else {
                greater = (s >= 2);
            }
        }

        long dMod = greater ? modSub(ctx.bMod, rMod) : rMod;
        long pow2KN = modPow(2, ctx.k - n);
        return modMul(dMod, pow2KN);
    }

    static long computeFast(long k, long t, int r) {
        FastContext ctx = new FastContext();
        ctx.r = r;
        ctx.k = k;
        ctx.t = t;
        ctx.binom = new long[r + 1];
        ctx.binomMod = new long[r + 1];
        ctx.pow2Ti = new long[r + 1];
        ctx.prefixK = new long[r + 2];
        ctx.suffix = new long[r + 2];

        ctx.binom[0] = 1;
        BigInteger biNum;
        for (int i = 1; i <= r; ++i) {
            biNum = BigInteger.valueOf(ctx.binom[i - 1]).multiply(BigInteger.valueOf(r - i + 1));
            ctx.binom[i] = biNum.divide(BigInteger.valueOf(i)).longValue();
        }
        for (int i = 0; i <= r; ++i) {
            ctx.binomMod[i] = ctx.binom[i] % MOD;
        }

        long pow2T = modPow(2, t);
        ctx.pow2Ti[0] = 1;
        for (int i = 1; i <= r; ++i) {
            ctx.pow2Ti[i] = modMul(ctx.pow2Ti[i - 1], pow2T);
        }

        ctx.pow2K = modPow(2, k);
        ctx.pow2KInv = modPow(ctx.pow2K, MOD - 2);
        ctx.bMod = modAdd(ctx.pow2K, 1);

        long[] term = new long[r + 1];
        long[] termK = new long[r + 1];
        for (int i = 0; i <= r; ++i) {
            term[i] = modMul(ctx.binomMod[i], ctx.pow2Ti[i]);
            termK[i] = modMul(term[i], ctx.pow2K);
        }

        for (int i = 0; i <= r; ++i) {
            ctx.prefixK[i + 1] = modAdd(ctx.prefixK[i], termK[i]);
        }
        for (int i = r; i >= 0; --i) {
            ctx.suffix[i] = modAdd(ctx.suffix[i + 1], term[i]);
        }

        ArrayList<Task> tasks = new ArrayList<>();
        long[] total = { 0 };

        java.util.function.BiConsumer<Long, Long> addBulk = (bulk, constant) -> {
            if (bulk <= 0)
                return;
            long bulkMod = bulk % MOD;
            total[0] = modAdd(total[0], modMul(bulkMod, constant));
        };

        for (int s = 1; s <= r; ++s) {
            BigInteger startRaw = BigInteger.valueOf(k).subtract(BigInteger.valueOf(s).multiply(BigInteger.valueOf(t)));
            BigInteger endRaw = BigInteger.valueOf(k)
                    .subtract(BigInteger.valueOf(s - 1).multiply(BigInteger.valueOf(t))).subtract(BigInteger.ONE);

            if (endRaw.compareTo(BigInteger.ZERO) < 0)
                continue;

            long start = startRaw.compareTo(BigInteger.ZERO) > 0 ? startRaw.longValue() : 0;
            long end = endRaw.longValue();
            if (end > k - 1)
                end = k - 1;
            if (start > end)
                continue;

            long length = end - start + 1;
            long excStart = Math.max(start, end - 61);
            long excCount = (excStart <= end) ? (end - excStart + 1) : 0;
            long bulk = length - excCount;

            long constant = modSub(ctx.prefixK[s], ctx.suffix[s]);
            addBulk.accept(bulk, constant);

            for (long n = excStart; n <= end; ++n) {
                tasks.add(new Task(n, s));
            }
        }

        int s = r + 1;
        BigInteger endRaw = BigInteger.valueOf(k).subtract(BigInteger.valueOf(r).multiply(BigInteger.valueOf(t)))
                .subtract(BigInteger.ONE);
        if (endRaw.compareTo(BigInteger.ZERO) >= 0) {
            long start = 0;
            long end = endRaw.longValue();
            if (end > k - 1)
                end = k - 1;

            if (start <= end) {
                long length = end - start + 1;
                long excStart = Math.max(start, end - 61);
                long excCount = (excStart <= end) ? (end - excStart + 1) : 0;
                long bulk = length - excCount;

                long constant = ctx.prefixK[r + 1];
                addBulk.accept(bulk, constant);

                for (long n = excStart; n <= end; ++n) {
                    tasks.add(new Task(n, s));
                }
            }
        }

        for (Task task : tasks) {
            total[0] = modAdd(total[0], exceptionContrib(ctx, task.n, task.s));
        }

        return total[0] % MOD;
    }

    public static String solve() {
        long k = 1000000000000000000L + 31L;
        long t = 100000000000000L + 31L;
        int r = 62;
        return Long.toString(computeFast(k, t, r));
    }

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