Problem 760: Sum over Bitwise Operators

Project Euler Problem 760 Solution

EulerSolve provides an optimized solution for Project Euler Problem 760, Sum over Bitwise Operators, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The implementations first collapse the original three-operator expression with the identity $$(u \oplus v) + (u \wedge v) + \operatorname{OR}(u,v)=2\operatorname{OR}(u,v).$$ So the target quantity can be written as $$G(n)=2\sum_{x=0}^{n}\sum_{k=0}^{x}\operatorname{OR}\bigl(k,x-k\bigr),$$ where $\operatorname{OR}$ denotes bitwise OR. The real input is $n=10^{18}$, so a direct double loop is hopeless. The solution therefore derives recurrences that replace $n$ by roughly $n/2$ at every step. Mathematical Approach Introduce the diagonal sum $$T(x)=\sum_{k=0}^{x}\operatorname{OR}\bigl(k,x-k\bigr)$$ and its prefix sum $$Q(n)=\sum_{x=0}^{n}T(x).$$ Then the desired answer is simply $$G(n)=2Q(n).$$ Step 1: Use bitwise parity identities Writing integers as $2a$ or $2a+1$ isolates the least significant bit. The four cases are $$\operatorname{OR}(2a,2b)=2\operatorname{OR}(a,b),$$ $$\operatorname{OR}(2a+1,2b)=2\operatorname{OR}(a,b)+1,$$ $$\operatorname{OR}(2a,2b+1)=2\operatorname{OR}(a,b)+1,$$ $$\operatorname{OR}(2a+1,2b+1)=2\operatorname{OR}(a,b)+1.$$ Each identity strips off one binary digit and replaces the original pair by a smaller pair. That is the key divide-and-conquer observation. Step 2: Derive the recurrence for one diagonal $T(x)$ Take an even argument $x=2m$. The pairs $(k,2m-k)$ split into even-even and odd-odd cases....

Detailed mathematical approach

Problem Summary

The implementations first collapse the original three-operator expression with the identity

$$(u \oplus v) + (u \wedge v) + \operatorname{OR}(u,v)=2\operatorname{OR}(u,v).$$

So the target quantity can be written as

$$G(n)=2\sum_{x=0}^{n}\sum_{k=0}^{x}\operatorname{OR}\bigl(k,x-k\bigr),$$

where $\operatorname{OR}$ denotes bitwise OR. The real input is $n=10^{18}$, so a direct double loop is hopeless. The solution therefore derives recurrences that replace $n$ by roughly $n/2$ at every step.

Mathematical Approach

Introduce the diagonal sum

$$T(x)=\sum_{k=0}^{x}\operatorname{OR}\bigl(k,x-k\bigr)$$

and its prefix sum

$$Q(n)=\sum_{x=0}^{n}T(x).$$

Then the desired answer is simply

$$G(n)=2Q(n).$$

Step 1: Use bitwise parity identities

Writing integers as $2a$ or $2a+1$ isolates the least significant bit. The four cases are

$$\operatorname{OR}(2a,2b)=2\operatorname{OR}(a,b),$$

$$\operatorname{OR}(2a+1,2b)=2\operatorname{OR}(a,b)+1,$$

$$\operatorname{OR}(2a,2b+1)=2\operatorname{OR}(a,b)+1,$$

$$\operatorname{OR}(2a+1,2b+1)=2\operatorname{OR}(a,b)+1.$$

Each identity strips off one binary digit and replaces the original pair by a smaller pair. That is the key divide-and-conquer observation.

Step 2: Derive the recurrence for one diagonal $T(x)$

Take an even argument $x=2m$. The pairs $(k,2m-k)$ split into even-even and odd-odd cases.

If $k=2a$, then $2m-k=2(m-a)$, so the contribution is $2\operatorname{OR}(a,m-a)$.

If $k=2a+1$, then $2m-k=2(m-a-1)+1$, so the contribution is $2\operatorname{OR}(a,m-a-1)+1$.

Summing both families gives

$$\begin{aligned} T(2m) &=\sum_{a=0}^{m}2\operatorname{OR}(a,m-a)+\sum_{a=0}^{m-1}\left(2\operatorname{OR}(a,m-1-a)+1\right)\\ &=2T(m)+2T(m-1)+m. \end{aligned}$$

For an odd argument $x=2m+1$, every pair has opposite parity, and both parity orders reduce to the same smaller problem. Therefore

$$\begin{aligned} T(2m+1) &=\sum_{a=0}^{m}\left(2\operatorname{OR}(a,m-a)+1\right)+\sum_{a=0}^{m}\left(2\operatorname{OR}(a,m-a)+1\right)\\ &=4T(m)+2(m+1). \end{aligned}$$

Step 3: Turn $T(x)$ into a recurrence for the prefix sum $Q(n)$

Since $Q(n)$ is the prefix sum of $T$, an even index can be grouped by parity:

$$Q(2m)=\sum_{j=0}^{m}T(2j)+\sum_{j=0}^{m-1}T(2j+1).$$

Substituting the formulas from Step 2 and collecting the repeated terms yields

$$Q(2m)=2Q(m)+6Q(m-1)+\frac{3m(m+1)}{2}.$$

For odd indices, start from $Q(2m+1)=Q(2m)+T(2m+1)$. Because

$$T(m)=Q(m)-Q(m-1),$$

the odd case simplifies to

$$Q(2m+1)=6Q(m)+2Q(m-1)+\frac{3m^2+7m+4}{2}.$$

Step 4: Final recurrence system and boundary values

The whole computation is determined by the zero boundary rule

$$T(n)=Q(n)=0 \qquad \text{for } n\le 0,$$

together with

$$T(2m)=2T(m)+2T(m-1)+m,$$

$$T(2m+1)=4T(m)+2(m+1),$$

$$Q(2m)=2Q(m)+6Q(m-1)+\frac{3m(m+1)}{2},$$

$$Q(2m+1)=6Q(m)+2Q(m-1)+\frac{3m^2+7m+4}{2}.$$

The polynomial parts are always integers: $m(m+1)$ is even, and $3m^2+7m+4 \equiv m(m+1)\pmod 2$ is even as well. This is why the modular implementation can safely replace division by $2$ with multiplication by the modular inverse of $2$.

Step 5: Worked example for $n=10$

Start from the small diagonal values:

$$T(1)=2,\qquad T(2)=5,\qquad T(4)=16,\qquad T(5)=26.$$

The even recurrence then gives

$$T(10)=2T(5)+2T(4)+5=2\cdot 26+2\cdot 16+5=89.$$

For the prefix sums, one obtains

$$Q(4)=35,\qquad Q(5)=61,$$

$$Q(10)=2Q(5)+6Q(4)+\frac{3\cdot 5\cdot 6}{2}=2\cdot 61+6\cdot 35+45=377.$$

Finally,

$$G(10)=2Q(10)=754,$$

which matches the checkpoint used by the implementations. A second checkpoint is $G(100)=583766$.

How the Code Works

The C++, Python, and Java implementations memoize the pair $(T(n),Q(n))$ for each requested argument. When a new state is needed, they first evaluate the two smaller states at $\lfloor n/2\rfloor$ and $\lfloor n/2\rfloor-1$, then apply the even or odd formula above.

All arithmetic is performed modulo $10^9+7$. The coefficients involving $\frac{1}{2}$ are handled by multiplying with the modular inverse of $2$, so the recurrence remains exact in modular arithmetic.

Once the memoized prefix value $Q(n)$ is known, the final answer is returned as $2Q(n)\bmod 10^9+7$.

Complexity Analysis

Each recurrence step replaces $n$ by $\lfloor n/2\rfloor$ or $\lfloor n/2\rfloor-1$, so only a logarithmic number of distinct arguments are memoized. That gives $O(\log n)$ time and $O(\log n)$ memory. The direct definition, by contrast, needs $\Theta(n^2)$ bitwise evaluations.

Footnotes and References

Problem page: https://projecteuler.net/problem=760
Bitwise operation: Wikipedia — Bitwise operation
Exclusive or: Wikipedia — Exclusive or
Divide-and-conquer algorithm: Wikipedia — Divide-and-conquer algorithm
Memoization: Wikipedia — Memoization
Recurrence relation: Wikipedia — Recurrence relation

Problem 760 source code

C++

#include <cassert>
#include <cstdint>
#include <iostream>
#include <unordered_map>
#include <utility>

namespace {

using i64 = long long;
using i128 = __int128_t;

constexpr i64 MOD = 1'000'000'007LL;
constexpr i64 INV2 = (MOD + 1) / 2;

struct HS {
    i64 h;
    i64 s;
};

std::unordered_map<i64, HS> memo;

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

i64 mul_mod(i64 a, i64 b) {
    return static_cast<i64>((static_cast<i128>(a) * b) % MOD);
}

HS solve(i64 n) {
    if (n <= 0) {
        return HS{0, 0};
    }
    auto it = memo.find(n);
    if (it != memo.end()) {
        return it->second;
    }

    const i64 m = n / 2;
    const HS sm = solve(m);
    const HS sm1 = solve(m - 1);
    const i64 mm = m % MOD;

    i64 h = 0;
    i64 s = 0;

    if ((n & 1LL) == 0LL) {
        h = (mul_mod(2, sm.h) + mul_mod(2, sm1.h) + mm) % MOD;

        i64 poly = mul_mod(3, mul_mod(mm, (mm + 1) % MOD));
        poly = mul_mod(poly, INV2);
        s = (mul_mod(2, sm.s) + mul_mod(6, sm1.s) + poly) % MOD;
    } else {
        h = (mul_mod(4, sm.h) + mul_mod(2, (mm + 1) % MOD)) % MOD;

        i64 poly = 0;
        poly = (poly + mul_mod(3, mul_mod(mm, mm))) % MOD;
        poly = (poly + mul_mod(7, mm)) % MOD;
        poly = (poly + 4) % MOD;
        poly = mul_mod(poly, INV2);
        s = (mul_mod(6, sm.s) + mul_mod(2, sm1.s) + poly) % MOD;
    }

    const HS out{h, s};
    memo.emplace(n, out);
    return out;
}

i64 brute_G(int n) {
    i64 total = 0;
    for (int x = 0; x <= n; ++x) {
        for (int k = 0; k <= x; ++k) {
            total += 2LL * static_cast<i64>(k | (x - k));
        }
    }
    return total;
}

}  // namespace

int main() {
    assert(brute_G(10) == 754);
    assert(brute_G(100) == 583766);

    memo.reserve(256);
    const i64 g10 = mul_mod(2, solve(10).s);
    const i64 g100 = mul_mod(2, solve(100).s);
    assert(g10 == 754);
    assert(g100 == 583766);

    const i64 answer = mul_mod(2, solve(1'000'000'000'000'000'000LL).s);
    std::cout << answer << '\n';
    return 0;
}

Python

import sys

sys.setrecursionlimit(2000)

MOD = 1000000007
INV2 = (MOD + 1) // 2

memo = {}

def add_mod(a, b):
    r = a + b
    if r >= MOD:
        r -= MOD
    return r

def mul_mod(a, b):
    return (a * b) % MOD

def solve_rec(n):
    if n <= 0:
        return (0, 0)
    if n in memo:
        return memo[n]
        
    m = n // 2
    sm_h, sm_s = solve_rec(m)
    sm1_h, sm1_s = solve_rec(m - 1)
    mm = m % MOD
    
    if (n & 1) == 0:
        h = (mul_mod(2, sm_h) + mul_mod(2, sm1_h) + mm) % MOD
        
        poly = mul_mod(3, mul_mod(mm, (mm + 1) % MOD))
        poly = mul_mod(poly, INV2)
        s = (mul_mod(2, sm_s) + mul_mod(6, sm1_s) + poly) % MOD
    else:
        h = (mul_mod(4, sm_h) + mul_mod(2, (mm + 1) % MOD)) % MOD
        
        poly = 0
        poly = (poly + mul_mod(3, mul_mod(mm, mm))) % MOD
        poly = (poly + mul_mod(7, mm)) % MOD
        poly = (poly + 4) % MOD
        poly = mul_mod(poly, INV2)
        s = (mul_mod(6, sm_s) + mul_mod(2, sm1_s) + poly) % MOD
        
    out = (h, s)
    memo[n] = out
    return out

def solve():
    n = 1000000000000000000
    h, s = solve_rec(n)
    ans = mul_mod(2, s)
    return str(ans)

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

Java

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

public class Euler760 {
    static final long MOD = 1000000007L;
    static final long INV2 = (MOD + 1) / 2;

    static class HS {
        long h, s;

        HS(long h, long s) {
            this.h = h;
            this.s = s;
        }
    }

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

    static long addMod(long a, long b) {
        long r = a + b;
        return (r >= MOD) ? (r - MOD) : r;
    }

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

    static HS solveRec(long n) {
        if (n <= 0) {
            return new HS(0, 0);
        }
        HS cached = memo.get(n);
        if (cached != null) {
            return cached;
        }

        long m = n / 2;
        HS sm = solveRec(m);
        HS sm1 = solveRec(m - 1);
        long mm = m % MOD;

        long h = 0;
        long s = 0;

        if ((n & 1L) == 0L) {
            h = (mulMod(2, sm.h) + mulMod(2, sm1.h) + mm) % MOD;

            long poly = mulMod(3, mulMod(mm, (mm + 1) % MOD));
            poly = mulMod(poly, INV2);
            s = (mulMod(2, sm.s) + mulMod(6, sm1.s) + poly) % MOD;
        } else {
            h = (mulMod(4, sm.h) + mulMod(2, (mm + 1) % MOD)) % MOD;

            long poly = 0;
            poly = (poly + mulMod(3, mulMod(mm, mm))) % MOD;
            poly = (poly + mulMod(7, mm)) % MOD;
            poly = (poly + 4) % MOD;
            poly = mulMod(poly, INV2);
            s = (mulMod(6, sm.s) + mulMod(2, sm1.s) + poly) % MOD;
        }

        HS out = new HS(h, s);
        memo.put(n, out);
        return out;
    }

    public static String solve() {
        long n = 1000000000000000000L;
        HS res = solveRec(n);
        long ans = mulMod(2, res.s);
        return Long.toString(ans);
    }

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