Problem 792: Too Many Twos

View on Project Euler

Project Euler Problem 792 Solution

EulerSolve provides an optimized solution for Project Euler Problem 792, Too Many Twos, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The problem asks for $$U(10^4)=\sum_{n=1}^{10^4}u(n^3),$$ where the arithmetic function \(u(m)\) is determined by a 2-adic valuation. If we set \(k=m+1\) and define the 2-adically convergent series $$S_k=\sum_{j\ge 0}(-2)^j\binom{2(k+j)}{k+j},$$ then the required quantity is $$u(m)=k+v_2(-3S_k).$$ Here \(v_2(x)\) denotes the exponent of 2 dividing \(x\). The difficulty is not the outer sum itself but recovering the power of 2 in this series for many large cubic inputs without expanding enormous exact integers. Mathematical Approach The key idea is to stay inside \(\mathbb{Z}/2^{64}\mathbb{Z}\). At that precision the relevant 2-adic valuation is still visible, while every arithmetic step fits into machine-word computations. Step 1: Truncate the 2-adic series at fixed precision Let \(C_r=\binom{2r}{r}\). A standard fact for central binomial coefficients is $$v_2(C_r)=s_2(r),$$ where \(s_2(r)\) is the number of 1-bits in the binary expansion of \(r\). Therefore $$v_2\!\left((-2)^j C_{k+j}\right)=j+s_2(k+j)\ge j+1.$$ So modulo \(2^{64}\), all sufficiently late terms vanish automatically. A 64-step accumulation is enough at the chosen precision, which converts an infinite 2-adic series into a fixed-size computation....

Detailed mathematical approach

Problem Summary

The problem asks for

$$U(10^4)=\sum_{n=1}^{10^4}u(n^3),$$

where the arithmetic function \(u(m)\) is determined by a 2-adic valuation. If we set \(k=m+1\) and define the 2-adically convergent series

$$S_k=\sum_{j\ge 0}(-2)^j\binom{2(k+j)}{k+j},$$

then the required quantity is

$$u(m)=k+v_2(-3S_k).$$

Here \(v_2(x)\) denotes the exponent of 2 dividing \(x\). The difficulty is not the outer sum itself but recovering the power of 2 in this series for many large cubic inputs without expanding enormous exact integers.

Mathematical Approach

The key idea is to stay inside \(\mathbb{Z}/2^{64}\mathbb{Z}\). At that precision the relevant 2-adic valuation is still visible, while every arithmetic step fits into machine-word computations.

Step 1: Truncate the 2-adic series at fixed precision

Let \(C_r=\binom{2r}{r}\). A standard fact for central binomial coefficients is

$$v_2(C_r)=s_2(r),$$

where \(s_2(r)\) is the number of 1-bits in the binary expansion of \(r\). Therefore

$$v_2\!\left((-2)^j C_{k+j}\right)=j+s_2(k+j)\ge j+1.$$

So modulo \(2^{64}\), all sufficiently late terms vanish automatically. A 64-step accumulation is enough at the chosen precision, which converts an infinite 2-adic series into a fixed-size computation.

Step 2: Separate powers of 2 from odd parts

For any nonzero integer define its odd part by

$$\operatorname{Odd}(x)=\frac{x}{2^{v_2(x)}}.$$

Then every central binomial coefficient splits as

$$C_r=\operatorname{Odd}(C_r)\,2^{s_2(r)}.$$

The power of 2 is therefore known immediately from the binary digit count of \(r\). The remaining task is to evaluate the odd part efficiently. Using factorials,

$$\operatorname{Odd}(C_r)=\frac{\operatorname{Odd}((2r)!)}{\operatorname{Odd}(r!)^2}.$$

Because the denominator is odd, it is invertible modulo \(2^{64}\), so division can be replaced by multiplication with odd modular inverses.

Step 3: Build the odd part of a factorial recursively

Introduce the odd prefix product

$$F(u)=\prod_{t=0}^{u-1}(2t+1)=1\cdot 3\cdot 5\cdots (2u-1).$$

Every factor in \(n!\) is either odd or twice a smaller integer. After stripping away powers of 2, this gives

$$\operatorname{Odd}(n!)=\operatorname{Odd}\!\left(\left\lfloor\frac{n}{2}\right\rfloor!\right)\,F\!\left(\left\lceil\frac{n}{2}\right\rceil\right).$$

Iterating the recursion yields

$$\operatorname{Odd}(n!)=\prod_{a\ge 0}F\!\left(\left\lfloor\frac{n+2^a}{2^{a+1}}\right\rfloor\right),$$

with only \(O(\log n)\) nontrivial factors. This is why the implementation can evaluate the odd part of very large factorials without ever materializing the factorial itself.

Step 4: Evaluate \(F(u)\) modulo \(2^{64}\) with a Newton series

In the working 2-adic setting, the implementation evaluates \(F(u)\) from a forward-difference expansion

$$F(u)\equiv \sum_{t=0}^{126}\Delta^tF(0)\binom{u}{t}\pmod{2^{64}}.$$

Only the first 127 coefficients are needed at this precision, so they are precomputed once from \(F(0),F(1),\dots,F(126)\). The generalized binomial coefficients \(\binom{u}{t}\) are also formed modulo \(2^{64}\) by separating powers of 2 from odd numerator and denominator parts.

Whenever an odd denominator must be inverted, Newton's method in the 2-adic world is used:

$$x_{r+1}=x_r(2-ax_r)\pmod{2^{64}},\qquad a\text{ odd}.$$

Each iteration doubles the number of correct low bits, so six rounds already reach full 64-bit precision.

Step 5: Update consecutive central binomial coefficients cheaply

Once the odd part of \(C_k\) is known, the next one follows from

$$\frac{C_{r+1}}{C_r}=\frac{2(2r+1)}{r+1}.$$

After removing powers of 2, the odd parts satisfy

$$\operatorname{Odd}(C_{r+1})=\operatorname{Odd}(C_r)\,\frac{2r+1}{\operatorname{Odd}(r+1)}.$$

So the implementation computes the first central binomial term at \(r=k\), then advances through \(r=k,k+1,\dots\) with a short recurrence. At each step it restores the full coefficient by shifting the odd part by \(s_2(r)\), multiplies by the appropriate power of \(-2\), and accumulates modulo \(2^{64}\).

Step 6: Worked example for \(m=4\)

Here \(k=m+1=5\). The first central binomial coefficients are

$$\binom{10}{5}=252=2^2\cdot 63,\qquad \binom{12}{6}=924=2^2\cdot 231,\qquad \binom{14}{7}=3432=2^3\cdot 429.$$

Hence the first terms of the series have valuations

$$v_2\!\left(\binom{10}{5}\right)=2,\qquad v_2\!\left(-2\binom{12}{6}\right)=3,\qquad v_2\!\left(4\binom{14}{7}\right)=5.$$

Modulo \(8\), only the first term survives, so

$$S_5\equiv \binom{10}{5}\equiv 4\pmod{8}.$$

Multiplication by \(-3\) does not change the 2-adic valuation because 3 is odd. Therefore

$$v_2(-3S_5)=2,$$

and thus

$$u(4)=5+2=7.$$

This matches the small checkpoint used by the implementation and shows why the low 2-adic valuation, not the exact integer magnitude, is the critical object.

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. First they precompute the odd prefix products \(F(0),\dots,F(126)\), the associated forward differences, and the inverses of the small odd numbers needed in the Newton expansion. Then, for each cube input \(m=n^3\), the implementation constructs the odd part of the first central binomial coefficient at \(k=m+1\), restores its full power of 2 from the binary digit count, and evaluates the fixed-length alternating series modulo \(2^{64}\).

After the series has been accumulated, the implementation multiplies by \(-3\), counts the exponent of 2 in the resulting residue, adds \(k\), and obtains \(u(m)\). The outer loop sums these values for \(n=1,2,\dots,10^4\). The direct compiled implementations also split that outer range across available processor cores, while the Python version delegates to the same optimized core computation.

Complexity Analysis

The precomputation phase is constant-sized: it stores only 127 forward-difference entries and a matching set of small odd inverses, so its cost is \(O(1)\) for this problem. For one value \(u(m)\), the odd-factorial formula needs \(O(\log m)\) evaluations of \(F\), each performed in constant time with the fixed Newton expansion, and the series summation always uses 64 iterations. Thus one query costs \(O(\log m)\) word operations, and summing \(u(n^3)\) for \(1\le n\le N\) costs

$$O\!\left(\sum_{n=1}^{N}\log(n^3)\right)=O(N\log N)$$

total time with \(O(1)\) extra memory beyond the fixed tables. In practice the constants are small because all arithmetic is carried out modulo \(2^{64}\).

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=792
  2. 2-adic valuation: Wikipedia - P-adic valuation
  3. Central binomial coefficient: Wikipedia - Central binomial coefficient
  4. Legendre's formula: Wikipedia - Legendre's formula
  5. Hensel's lemma: Wikipedia - Hensel's lemma
  6. Newton series: Wikipedia - Newton series

Problem 792 source code

C++

#include <algorithm>
#include <array>
#include <atomic>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <stdexcept>
#include <string>
#include <thread>
#include <vector>

using namespace std;

struct Mod2Pow64 {
    static constexpr int E = 64;
    static constexpr int D = 2 * E - 2;

    array<uint64_t, D + 1> F_small{};

    array<uint64_t, D + 1> dF{};

    array<uint64_t, D + 1> inv_small{};

    static inline uint64_t inv_odd(uint64_t a) {

        if ((a & 1ULL) == 0) {
            throw runtime_error("inv_odd called with even input");
        }
        uint64_t x = 1;
        for (int i = 0; i < 6; ++i) {
            x = x * (uint64_t)(2 - a * x);
        }
        return x;
    }

    Mod2Pow64() {

        F_small[0] = 1;
        for (int u = 1; u <= D; ++u) {
            uint64_t factor = (uint64_t)(2 * (u - 1) + 1);
            F_small[u] = F_small[u - 1] * factor;
        }

        array<uint64_t, D + 1> tmp = F_small;
        for (int k = 0; k <= D; ++k) {
            dF[k] = tmp[0];

            for (int i = 0; i <= D - k - 1; ++i) {
                tmp[i] = tmp[i + 1] - tmp[i];
            }
        }

        inv_small.fill(0);
        for (int x = 1; x <= D; x += 2) {
            inv_small[x] = inv_odd((uint64_t)x);
        }
    }

    inline uint64_t evalF(uint64_t u) const {

        if (u <= (uint64_t)D) return F_small[(size_t)u];

        uint64_t res = dF[0];

        uint64_t binom_odd = 1;
        int exp2 = 0;

        for (int k = 1; k <= D; ++k) {
            uint64_t num = u - (uint64_t)k + 1;
            uint64_t den = (uint64_t)k;

            int tz_num = __builtin_ctzll(num);
            int tz_den = __builtin_ctzll(den);

            exp2 += tz_num - tz_den;

            binom_odd *= (num >> tz_num);
            uint64_t den_odd = (den >> tz_den);
            binom_odd *= inv_small[(size_t)den_odd];

            uint64_t binom = (exp2 >= 64) ? 0ULL : (binom_odd << exp2);
            res += dF[(size_t)k] * binom;
        }
        return res;
    }

    inline uint64_t odd_fact(uint64_t n) const {

        uint64_t res = 1;
        while (n > 0) {
            uint64_t m = (n + 1) >> 1;
            res *= evalF(m);
            n >>= 1;
        }
        return res;
    }

    inline uint64_t odd_central_binom(uint64_t k) const {

        uint64_t a = odd_fact(2 * k);
        uint64_t b = odd_fact(k);
        uint64_t invb = inv_odd(b);
        return a * invb * invb;
    }
};

static inline uint64_t u_value(uint64_t n, const Mod2Pow64& mod) {

    const uint64_t base = n + 1;

    uint64_t O = mod.odd_central_binom(base);
    uint64_t k = base;

    int pc = __builtin_popcountll(k);
    uint64_t B = (pc >= 64) ? 0ULL : (O << pc);

    uint64_t sum = 0;
    for (int j = 0; j < 64; ++j) {
        uint64_t term = B << j;
        if (j & 1) term = (uint64_t)(0) - term;
        sum += term;

        if (j != 63) {

            uint64_t denom = k + 1;
            int tz = __builtin_ctzll(denom);
            uint64_t denom_odd = denom >> tz;
            uint64_t inv = Mod2Pow64::inv_odd(denom_odd);

            O *= (2 * k + 1);
            O *= inv;

            ++k;
            pc = __builtin_popcountll(k);
            B = (pc >= 64) ? 0ULL : (O << pc);
        }
    }

    uint64_t c = (uint64_t)(-3) * sum;

    if (c == 0) {

        throw runtime_error("Need higher 2-adic precision (c==0 mod 2^64)");
    }

    int tz = __builtin_ctzll(c);
    return base + (uint64_t)tz;
}

static inline unsigned long long U_value(int N, const Mod2Pow64& mod) {
    unsigned long long total = 0;
    for (int n = 1; n <= N; ++n) {
        uint64_t x = (uint64_t)n;
        uint64_t cube = x * x * x;
        total += (unsigned long long)u_value(cube, mod);
    }
    return total;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    Mod2Pow64 mod;

    auto check = [&](const string& name, unsigned long long got, unsigned long long expected) {
        if (got != expected) {
            cerr << "Validation failed: " << name << " got " << got
                 << " expected " << expected << "\n";
            exit(1);
        }
    };

    check("u(4)",  (unsigned long long)u_value(4, mod), 7ULL);
    check("u(20)", (unsigned long long)u_value(20, mod), 24ULL);
    check("U(5)",  U_value(5, mod), 241ULL);

    const int N = 10000;
    unsigned int T = thread::hardware_concurrency();
    if (T == 0) T = 4;
    if (T > (unsigned int)N) T = N;

    atomic<unsigned long long> global_sum{0};

    int chunk = (N + (int)T - 1) / (int)T;
    vector<thread> threads;
    threads.reserve(T);

    for (unsigned int t = 0; t < T; ++t) {
        int L = (int)t * chunk + 1;
        int R = min(N, (int)(t + 1) * chunk);
        if (L > R) break;

        threads.emplace_back([&, L, R]() {
            unsigned long long local = 0;
            for (int n = L; n <= R; ++n) {
                uint64_t x = (uint64_t)n;
                uint64_t cube = x * x * x;
                local += (unsigned long long)u_value(cube, mod);
            }
            global_sum.fetch_add(local, memory_order_relaxed);
        });
    }

    for (auto& th : threads) th.join();

    unsigned long long ans = global_sum.load(memory_order_relaxed);

    cout << ans << "\n";
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

public class Euler792 {

    static class Mod2Pow64 {
        static final int E = 64;
        static final int D = 2 * E - 2;

        long[] FSmall = new long[D + 1];
        long[] dF = new long[D + 1];
        long[] invSmall = new long[D + 1];

        static long invOdd(long a) {
            if ((a & 1L) == 0) {
                throw new RuntimeException("invOdd called with even input");
            }
            long x = 1;
            for (int i = 0; i < 6; ++i) {
                x = x * (2L - a * x);
            }
            return x;
        }

        Mod2Pow64() {
            FSmall[0] = 1;
            for (int u = 1; u <= D; ++u) {
                long factor = (long) (2 * (u - 1) + 1);
                FSmall[u] = FSmall[u - 1] * factor;
            }

            long[] tmp = FSmall.clone();
            for (int k = 0; k <= D; ++k) {
                dF[k] = tmp[0];
                for (int i = 0; i <= D - k - 1; ++i) {
                    tmp[i] = tmp[i + 1] - tmp[i];
                }
            }

            for (int x = 1; x <= D; x += 2) {
                invSmall[x] = invOdd((long) x);
            }
        }

        long evalF(long u) {
            if (Long.compareUnsigned(u, D) <= 0)
                return FSmall[(int) u];

            long res = dF[0];
            long binomOdd = 1;
            int exp2 = 0;

            for (int k = 1; k <= D; ++k) {
                long num = u - k + 1;
                long den = k;

                int tzNum = Long.numberOfTrailingZeros(num);
                int tzDen = Long.numberOfTrailingZeros(den);

                exp2 += tzNum - tzDen;

                binomOdd *= (num >>> tzNum);
                long denOdd = (den >>> tzDen);
                binomOdd *= invSmall[(int) denOdd];

                long binom = (exp2 >= 64) ? 0L : (binomOdd << exp2);
                res += dF[k] * binom;
            }
            return res;
        }

        long oddFact(long n) {
            long res = 1;
            while (n > 0) {
                long m = (n + 1) >>> 1;
                res *= evalF(m);
                n >>>= 1;
            }
            return res;
        }

        long oddCentralBinom(long k) {
            long a = oddFact(2 * k);
            long b = oddFact(k);
            long invb = invOdd(b);
            return a * invb * invb;
        }
    }

    static long uValue(long n, Mod2Pow64 mod) {
        long base = n + 1;
        long O = mod.oddCentralBinom(base);
        long k = base;

        int pc = Long.bitCount(k);
        long B = (pc >= 64) ? 0L : (O << pc);

        long sum = 0;
        for (int j = 0; j < 64; ++j) {
            long term = B << j;
            if ((j & 1) == 1)
                term = -term;
            sum += term;

            if (j != 63) {
                long denom = k + 1;
                int tz = Long.numberOfTrailingZeros(denom);
                long denomOdd = denom >>> tz;
                long inv = Mod2Pow64.invOdd(denomOdd);

                O *= (2 * k + 1);
                O *= inv;

                k++;
                pc = Long.bitCount(k);
                B = (pc >= 64) ? 0L : (O << pc);
            }
        }

        long c = -3L * sum;
        if (c == 0) {
            throw new RuntimeException("Need higher 2-adic precision (c==0 mod 2^64)");
        }

        int tz = Long.numberOfTrailingZeros(c);
        return base + tz;
    }

    static long UValue(int N) {
        Mod2Pow64 mod = new Mod2Pow64();
        long[] partials = new long[16];
        Thread[] threads = new Thread[16];

        int numThreads = Runtime.getRuntime().availableProcessors();
        if (numThreads > 16)
            numThreads = 16;
        if (numThreads <= 0)
            numThreads = 1;

        int chunk = (N + numThreads - 1) / numThreads;

        for (int t = 0; t < numThreads; t++) {
            final int tid = t;
            int L = t * chunk + 1;
            int R = Math.min(N, (t + 1) * chunk);

            if (L > R)
                break;

            threads[t] = new Thread(() -> {
                long local = 0;
                for (int n = L; n <= R; n++) {
                    long x = n;
                    long cube = x * x * x;
                    local += uValue(cube, mod);
                }
                partials[tid] = local;
            });
            threads[t].start();
        }

        long total = 0;
        for (int t = 0; t < numThreads; t++) {
            if (threads[t] != null) {
                try {
                    threads[t].join();
                    total += partials[t];
                } catch (InterruptedException e) {
                    e.printStackTrace();
                }
            }
        }

        return total;
    }

    public static String solve() {
        return Long.toString(UValue(10000));
    }

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