Problem 882: Removing Bits

Project Euler Problem 882 Solution

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

Problem Summary For each positive integer $x$, consider all integers obtained by deleting exactly one binary digit from the binary expansion of $x$. These smaller integers act as predecessors. A deleted $1$-bit gives a lower constraint, and a deleted $0$-bit gives an upper constraint. Starting from the base value $v_0=0$, the value $v_x$ is defined to be the simplest dyadic rational compatible with all those constraints. Here a dyadic rational means a number of the form $\frac{m}{2^k}$ with $m\in\mathbb{Z}$ and $k\ge 0$. The target quantity is $$S(n)=\left\lceil \sum_{i=1}^{n} i\,v_i \right\rceil.$$ Mathematical Approach The recurrence is local: once the deleted-bit predecessors of $x$ are known, $v_x$ is forced to lie in an open interval, and the simplest dyadic number inside that interval is chosen. Step 1: Deleted-Bit Predecessors Create the Interval Let the deleted-bit predecessors of $x$ be $y_1,\dots,y_\ell$, where $\ell$ is the binary length of $x$. Every predecessor is smaller than $x$, because deleting one bit shortens the representation to at most $\ell-1$ bits, so $$y_j<2^{\ell-1}\le x.$$ This is why the values can be computed in increasing order of $x$. Let $A_x$ be the set of predecessor values coming from deleted $1$-bits, and let $B_x$ be the set coming from deleted $0$-bits....

Detailed mathematical approach

Problem Summary

For each positive integer $x$, consider all integers obtained by deleting exactly one binary digit from the binary expansion of $x$. These smaller integers act as predecessors. A deleted $1$-bit gives a lower constraint, and a deleted $0$-bit gives an upper constraint.

Starting from the base value $v_0=0$, the value $v_x$ is defined to be the simplest dyadic rational compatible with all those constraints. Here a dyadic rational means a number of the form $\frac{m}{2^k}$ with $m\in\mathbb{Z}$ and $k\ge 0$.

The target quantity is

$$S(n)=\left\lceil \sum_{i=1}^{n} i\,v_i \right\rceil.$$

Mathematical Approach

The recurrence is local: once the deleted-bit predecessors of $x$ are known, $v_x$ is forced to lie in an open interval, and the simplest dyadic number inside that interval is chosen.

Step 1: Deleted-Bit Predecessors Create the Interval

Let the deleted-bit predecessors of $x$ be $y_1,\dots,y_\ell$, where $\ell$ is the binary length of $x$. Every predecessor is smaller than $x$, because deleting one bit shortens the representation to at most $\ell-1$ bits, so

$$y_j<2^{\ell-1}\le x.$$

This is why the values can be computed in increasing order of $x$.

Let $A_x$ be the set of predecessor values coming from deleted $1$-bits, and let $B_x$ be the set coming from deleted $0$-bits. Define

$$L_x=\max A_x,$$

and, when $B_x\neq\varnothing$, define

$$R_x=\min B_x.$$

The defining inequalities are

$$v_x>L_x,$$

and, if $B_x\neq\varnothing$, also

$$v_x<R_x.$$

Because the leading binary digit is always $1$, the set $A_x$ is never empty, so every $x\ge 1$ has at least a lower bound.

Step 2: Turn the Open Interval into Integer Inequalities

Suppose first that $B_x\neq\varnothing$ and test a dyadic candidate $\frac{m}{2^k}$. The condition

$$L_x<\frac{m}{2^k}<R_x$$

is equivalent to

$$\left\lfloor 2^kL_x\right\rfloor+1\le m\le \left\lceil 2^kR_x\right\rceil-1.$$

So, for a fixed denominator $2^k$, admissible numerators are exactly the integers in that closed interval. The smallest $k$ for which such an integer exists gives the coarsest denominator and therefore the simplest dyadic level.

If $B_x=\varnothing$, then the binary expansion of $x$ contains no $0$-bit. There is no upper bound, and the simplest admissible dyadic is just the smallest integer strictly larger than $L_x$:

$$v_x=\lfloor L_x\rfloor+1.$$

Step 3: Choose the Simplest Dyadic and Reduce It

Once the minimal denominator $2^k$ is known, there may be several admissible numerators. The implementations select the feasible integer closest to $0$: choose $m=0$ if possible, otherwise choose the admissible integer with smallest absolute value. This reproduces the intended notion of the simplest dyadic inside the interval.

After that, remove all common factors of $2$ from numerator and denominator. Each value is therefore stored in reduced form

$$v_x=\frac{a_x}{2^{e_x}},$$

where $a_x$ is odd unless the value is $0$.

This reduction is important because dyadic numbers can then be compared exactly by lifting them to a common power-of-two denominator and comparing the scaled numerators as integers. No floating-point approximation is needed anywhere.

Step 4: Evaluate the Weighted Ceiling Sum Exactly

After computing $v_1,\dots,v_n$, write each one as $v_i=\frac{a_i}{2^{e_i}}$ in reduced form, and let

$$E=\max_{1\le i\le n} e_i.$$

Then

$$\sum_{i=1}^{n} i\,v_i=\frac{1}{2^E}\sum_{i=1}^{n} i\,a_i\,2^{E-e_i}=\frac{N_n}{2^E},$$

where $N_n$ is an integer.

In this recurrence every $v_i$ is positive, so $N_n\ge 0$. Therefore

$$S(n)=\left\lceil \frac{N_n}{2^E}\right\rceil=\left\lfloor \frac{N_n+2^E-1}{2^E}\right\rfloor.$$

The final ceiling can thus be computed with integer arithmetic only.

Worked Example

The first few values come directly from the interval rule:

$$\begin{aligned} v_1&=1,\\ v_2&=\frac{1}{2},\\ v_3&=2,\\ v_4&=\frac{1}{4},\\ v_5&=\frac{3}{2}. \end{aligned}$$

For example, $x=5$ has binary form $101_2$. Deleting the leftmost $1$ gives $1$, deleting the middle $0$ gives $3$, and deleting the last $1$ gives $2$. So the lower bound is

$$L_5=\max(v_1,v_2)=1,$$

the upper bound is

$$R_5=v_3=2,$$

and the simplest dyadic in $(1,2)$ is

$$v_5=\frac{3}{2}.$$

Using the first five values,

$$S(5)=\left\lceil 1\cdot 1+2\cdot \frac{1}{2}+3\cdot 2+4\cdot \frac{1}{4}+5\cdot \frac{3}{2}\right\rceil=\left\lceil \frac{33}{2}\right\rceil=17.$$

Continuing the same recurrence gives

$$v_6=1,\qquad v_7=3,\qquad v_8=\frac{1}{8},\qquad v_9=\frac{5}{4},\qquad v_{10}=\frac{3}{4},$$

hence

$$S(10)=64,$$

which matches the implementation checkpoint.

How the Code Works

The C++, Python, and Java implementations iterate from $1$ up to $n$. For each $x$, they scan all bit positions, form each deleted-bit predecessor, and update the best lower and upper bounds using exact dyadic comparisons. Because every predecessor is smaller than $x$, all required values are already known when $x$ is processed.

They then test denominators $1,2,4,\dots$ until the interval contains an admissible dyadic rational. The chosen value is reduced immediately so the stored denominator is as small as possible. After the table is complete, the implementations lift every term $i\,v_i$ to a common denominator, add the integer numerators, and apply the ceiling formula exactly.

Complexity Analysis

For a given $x$, there are $\lfloor\log_2 x\rfloor+1$ deleted-bit predecessors to inspect, so building the interval for $v_x$ costs $O(\log x)$. The additional search over dyadic denominators is bounded by a small fixed constant in the implementations, so the total running time up to $n$ is $O(n\log n)$. The table of dyadic values uses $O(n)$ memory.

Footnotes and References

Problem page: https://projecteuler.net/problem=882
Dyadic rational: Wikipedia - Dyadic rational
Binary numeral system: Wikipedia - Binary number
Floor and ceiling functions: Wikipedia - Floor and ceiling functions

Problem 882 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

using i64 = std::int64_t;
using i128 = __int128_t;

struct Dyad {
    i64 num;
    int exp;
};

static Dyad normalize(i128 num, int exp) {
    if (num == 0) return {0, 0};
    while (exp > 0 && (num & 1) == 0) {
        num >>= 1;
        --exp;
    }
    return {static_cast<i64>(num), exp};
}

static i128 floor_div_pow2(i128 x, int sh) {
    if (sh == 0) return x;
    if (x >= 0) return x >> sh;
    i128 a = -x;
    return -((a + ((static_cast<i128>(1) << sh) - 1)) >> sh);
}

static int cmp_dyad(const Dyad& a, const Dyad& b) {
    i128 na = a.num;
    i128 nb = b.num;
    if (a.exp < b.exp) na <<= (b.exp - a.exp);
    else if (b.exp < a.exp) nb <<= (a.exp - b.exp);
    if (na < nb) return -1;
    if (na > nb) return 1;
    return 0;
}

static i128 floor_dyad(const Dyad& a) {
    return floor_div_pow2(static_cast<i128>(a.num), a.exp);
}

static i128 ceil_dyad(const Dyad& a) {
    Dyad neg{-a.num, a.exp};
    return -floor_dyad(neg);
}

static i128 floor_scaled(const Dyad& a, int k) {
    i128 n = static_cast<i128>(a.num);
    if (k >= a.exp) return n << (k - a.exp);
    return floor_div_pow2(n, a.exp - k);
}

static i128 ceil_scaled(const Dyad& a, int k) {
    Dyad neg{-a.num, a.exp};
    return -floor_scaled(neg, k);
}

static Dyad simplest_between(const Dyad& L, bool hasL, const Dyad& R, bool hasR) {
    if (!hasL && !hasR) return {0, 0};
    if (hasL && !hasR) return {static_cast<i64>(floor_dyad(L) + 1), 0};
    if (!hasL && hasR) return {static_cast<i64>(ceil_dyad(R) - 1), 0};

    assert(cmp_dyad(L, R) < 0);

    for (int k = 0; k <= 80; ++k) {
        i128 lo = floor_scaled(L, k) + 1;
        i128 hi = ceil_scaled(R, k) - 1;
        if (lo > hi) continue;

        i128 m;
        if (lo <= 0 && hi >= 0) m = 0;
        else if (hi < 0) m = hi;
        else m = lo;

        return normalize(m, k);
    }

    assert(false);
    return {0, 0};
}

static int bit_length(i64 x) {
    if (x == 0) return 1;
    return 64 - __builtin_clzll(static_cast<unsigned long long>(x));
}

static i64 delete_bit(i64 x, int len, int idx_from_left) {
    int p = len - 1 - idx_from_left;
    i64 hi = x >> (p + 1);
    i64 lo = (p == 0) ? 0 : (x & ((1LL << p) - 1));
    return (hi << p) | lo;
}

static std::vector<Dyad> compute_values(int n) {
    std::vector<Dyad> v(static_cast<std::size_t>(n + 1), Dyad{0, 0});

    for (int x = 1; x <= n; ++x) {
        int len = bit_length(x);
        Dyad bestL{0, 0}, bestR{0, 0};
        bool hasL = false, hasR = false;

        for (int i = 0; i < len; ++i) {
            int bit = (x >> (len - 1 - i)) & 1;
            int y = static_cast<int>(delete_bit(x, len, i));
            const Dyad& vy = v[static_cast<std::size_t>(y)];

            if (bit == 1) {
                if (!hasL || cmp_dyad(vy, bestL) > 0) {
                    bestL = vy;
                    hasL = true;
                }
            } else {
                if (!hasR || cmp_dyad(vy, bestR) < 0) {
                    bestR = vy;
                    hasR = true;
                }
            }
        }

        v[static_cast<std::size_t>(x)] = simplest_between(bestL, hasL, bestR, hasR);
    }

    return v;
}

static i128 S_value(int n) {
    auto vals = compute_values(n);
    int max_exp = 0;
    for (int i = 1; i <= n; ++i) max_exp = std::max(max_exp, vals[static_cast<std::size_t>(i)].exp);

    i128 scaled_sum = 0;
    for (int i = 1; i <= n; ++i) {
        const auto& d = vals[static_cast<std::size_t>(i)];
        int sh = max_exp - d.exp;
        i128 term = static_cast<i128>(i) * static_cast<i128>(d.num);
        scaled_sum += (term << sh);
    }

    i128 den = static_cast<i128>(1) << max_exp;
    return (scaled_sum + den - 1) / den;
}

static std::string to_string_i128(i128 x) {
    if (x == 0) return "0";
    bool neg = x < 0;
    if (neg) x = -x;
    std::string s;
    while (x > 0) {
        int digit = static_cast<int>(x % 10);
        s.push_back(static_cast<char>('0' + digit));
        x /= 10;
    }
    if (neg) s.push_back('-');
    std::reverse(s.begin(), s.end());
    return s;
}

int main() {
    assert(S_value(2) == 2);
    assert(S_value(5) == 17);
    assert(S_value(10) == 64);

    std::cout << to_string_i128(S_value(100'000)) << '\n';
    return 0;
}

Python

def normalize(num, exp):
    if num == 0:
        return 0, 0
    while exp > 0 and (num & 1) == 0:
        num >>= 1
        exp -= 1
    return num, exp

def floor_div_pow2(x, sh):
    if sh == 0:
        return x
    if x >= 0:
        return x >> sh
    a = -x
    return -((a + ((1 << sh) - 1)) >> sh)

def cmp_dyad(num_a, exp_a, num_b, exp_b):
    na = num_a
    nb = num_b
    if exp_a < exp_b:
        na <<= (exp_b - exp_a)
    elif exp_b < exp_a:
        nb <<= (exp_a - exp_b)
    
    if na < nb:
        return -1
    if na > nb:
        return 1
    return 0

def floor_dyad(num, exp):
    return floor_div_pow2(num, exp)

def ceil_dyad(num, exp):
    return -floor_div_pow2(-num, exp)

def floor_scaled(num, exp, k):
    if k >= exp:
        return num << (k - exp)
    return floor_div_pow2(num, exp - k)

def ceil_scaled(num, exp, k):
    return -floor_scaled(-num, exp, k)

def simplest_between(num_L, exp_L, hasL, num_R, exp_R, hasR):
    if not hasL and not hasR:
        return 0, 0
    if hasL and not hasR:
        return floor_dyad(num_L, exp_L) + 1, 0
    if not hasL and hasR:
        return ceil_dyad(num_R, exp_R) - 1, 0
        
    for k in range(81):
        lo = floor_scaled(num_L, exp_L, k) + 1
        hi = ceil_scaled(num_R, exp_R, k) - 1
        if lo > hi:
            continue
            
        if lo <= 0 and hi >= 0:
            m = 0
        elif hi < 0:
            m = hi
        else:
            m = lo
            
        return normalize(m, k)
    
    return 0, 0

def bit_length_custom(x):
    if x == 0:
        return 1
    return x.bit_length()

def delete_bit(x, length, idx_from_left):
    p = length - 1 - idx_from_left
    hi = x >> (p + 1)
    lo = 0 if p == 0 else (x & ((1 << p) - 1))
    return (hi << p) | lo

def compute_values(n):
    v = [(0, 0)] * (n + 1)
    
    for x in range(1, n + 1):
        length = bit_length_custom(x)
        bestL_num, bestL_exp = 0, 0
        bestR_num, bestR_exp = 0, 0
        hasL = False
        hasR = False
        
        for i in range(length):
            bit = (x >> (length - 1 - i)) & 1
            y = delete_bit(x, length, i)
            vy_num, vy_exp = v[y]
            
            if bit == 1:
                if not hasL or cmp_dyad(vy_num, vy_exp, bestL_num, bestL_exp) > 0:
                    bestL_num = vy_num
                    bestL_exp = vy_exp
                    hasL = True
            else:
                if not hasR or cmp_dyad(vy_num, vy_exp, bestR_num, bestR_exp) < 0:
                    bestR_num = vy_num
                    bestR_exp = vy_exp
                    hasR = True
                    
        v[x] = simplest_between(bestL_num, bestL_exp, hasL, bestR_num, bestR_exp, hasR)
        
    return v

def S_value(n):
    vals = compute_values(n)
    max_exp = 0
    for i in range(1, n + 1):
        max_exp = max(max_exp, vals[i][1])
        
    scaled_sum = 0
    for i in range(1, n + 1):
        num, exp = vals[i]
        sh = max_exp - exp
        term = i * num
        scaled_sum += (term << sh)
        
    den = 1 << max_exp
    return (scaled_sum + den - 1) // den

def solve():
    return str(S_value(100000))

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

Java

import java.math.BigInteger;

public class Euler882 {

    static class Dyad {
        BigInteger num;
        int exp;

        Dyad(BigInteger num, int exp) {
            this.num = num;
            this.exp = exp;
        }
    }

    static Dyad normalize(BigInteger num, int exp) {
        if (num.equals(BigInteger.ZERO))
            return new Dyad(BigInteger.ZERO, 0);
        while (exp > 0 && !num.testBit(0)) {
            num = num.shiftRight(1);
            --exp;
        }
        return new Dyad(num, exp);
    }

    static BigInteger floorDivPow2(BigInteger x, int sh) {
        if (sh == 0)
            return x;
        if (x.compareTo(BigInteger.ZERO) >= 0)
            return x.shiftRight(sh);
        BigInteger a = x.negate();
        BigInteger oneMask = BigInteger.ONE.shiftLeft(sh).subtract(BigInteger.ONE);
        return a.add(oneMask).shiftRight(sh).negate();
    }

    static int cmpDyad(Dyad a, Dyad b) {
        BigInteger na = a.num;
        BigInteger nb = b.num;
        if (a.exp < b.exp)
            na = na.shiftLeft(b.exp - a.exp);
        else if (b.exp < a.exp)
            nb = nb.shiftLeft(a.exp - b.exp);
        return na.compareTo(nb);
    }

    static BigInteger floorDyad(Dyad a) {
        return floorDivPow2(a.num, a.exp);
    }

    static BigInteger ceilDyad(Dyad a) {
        Dyad neg = new Dyad(a.num.negate(), a.exp);
        return floorDyad(neg).negate();
    }

    static BigInteger floorScaled(Dyad a, int k) {
        BigInteger n = a.num;
        if (k >= a.exp)
            return n.shiftLeft(k - a.exp);
        return floorDivPow2(n, a.exp - k);
    }

    static BigInteger ceilScaled(Dyad a, int k) {
        Dyad neg = new Dyad(a.num.negate(), a.exp);
        return floorScaled(neg, k).negate();
    }

    static Dyad simplestBetween(Dyad L, boolean hasL, Dyad R, boolean hasR) {
        if (!hasL && !hasR)
            return new Dyad(BigInteger.ZERO, 0);
        if (hasL && !hasR)
            return new Dyad(floorDyad(L).add(BigInteger.ONE), 0);
        if (!hasL && hasR)
            return new Dyad(ceilDyad(R).subtract(BigInteger.ONE), 0);

        for (int k = 0; k <= 80; ++k) {
            BigInteger lo = floorScaled(L, k).add(BigInteger.ONE);
            BigInteger hi = ceilScaled(R, k).subtract(BigInteger.ONE);
            if (lo.compareTo(hi) > 0)
                continue;

            BigInteger m;
            if (lo.compareTo(BigInteger.ZERO) <= 0 && hi.compareTo(BigInteger.ZERO) >= 0)
                m = BigInteger.ZERO;
            else if (hi.compareTo(BigInteger.ZERO) < 0)
                m = hi;
            else
                m = lo;

            return normalize(m, k);
        }

        return new Dyad(BigInteger.ZERO, 0);
    }

    static int bitLength(int x) {
        if (x == 0)
            return 1;
        return 32 - Integer.numberOfLeadingZeros(x);
    }

    static int deleteBit(int x, int len, int idxFromLeft) {
        int p = len - 1 - idxFromLeft;
        int hi = x >>> (p + 1);
        int lo = (p == 0) ? 0 : (x & ((1 << p) - 1));
        return (hi << p) | lo;
    }

    static Dyad[] computeValues(int n) {
        Dyad[] v = new Dyad[n + 1];
        v[0] = new Dyad(BigInteger.ZERO, 0);

        for (int x = 1; x <= n; ++x) {
            int len = bitLength(x);
            Dyad bestL = new Dyad(BigInteger.ZERO, 0);
            Dyad bestR = new Dyad(BigInteger.ZERO, 0);
            boolean hasL = false, hasR = false;

            for (int i = 0; i < len; ++i) {
                int bit = (x >> (len - 1 - i)) & 1;
                int y = deleteBit(x, len, i);
                Dyad vy = v[y];

                if (bit == 1) {
                    if (!hasL || cmpDyad(vy, bestL) > 0) {
                        bestL = vy;
                        hasL = true;
                    }
                } else {
                    if (!hasR || cmpDyad(vy, bestR) < 0) {
                        bestR = vy;
                        hasR = true;
                    }
                }
            }

            v[x] = simplestBetween(bestL, hasL, bestR, hasR);
        }

        return v;
    }

    static BigInteger SValue(int n) {
        Dyad[] vals = computeValues(n);
        int maxExp = 0;
        for (int i = 1; i <= n; ++i) {
            if (vals[i].exp > maxExp)
                maxExp = vals[i].exp;
        }

        BigInteger scaledSum = BigInteger.ZERO;
        for (int i = 1; i <= n; ++i) {
            Dyad d = vals[i];
            int sh = maxExp - d.exp;
            BigInteger term = BigInteger.valueOf(i).multiply(d.num);
            scaledSum = scaledSum.add(term.shiftLeft(sh));
        }

        BigInteger den = BigInteger.ONE.shiftLeft(maxExp);
        return scaledSum.add(den).subtract(BigInteger.ONE).divide(den);
    }

    public static String solve() {
        return SValue(100000).toString();
    }

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