Problem 890: Binary Partitions

Project Euler Problem 890 Solution

EulerSolve provides an optimized solution for Project Euler Problem 890, Binary Partitions, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $p(n)$ denote the number of ways to write $n$ as a sum of powers of $2$, with unlimited repetition allowed. Problem 890 asks for $$p\left(7^{777}\right)\pmod{10^9+7}.$$ The classical generating function is $$\prod_{j\ge 0}\frac{1}{1-x^{2^j}}=\sum_{n\ge 0}p(n)x^n.$$ Because $7^{777}$ is enormous, the implementation cannot iterate up to $n$. Instead, it reads the binary digits of $n$ and counts valid carry patterns. That turns the problem into a digit DP whose size depends on the bit-length of $n$, not on $n$ itself. Mathematical Approach The key observation is that a binary partition can be read bit by bit. At each bit position we only need to know how many lower-power terms have paired up and carried into the next position. Step 1: Encode a Partition as a Carry Process Write $$n=\sum_{k=0}^{L-1} b_k2^k,\qquad b_k\in\{0,1\}.$$ If a partition uses $x_k$ copies of $2^k$, then after combining pairs of $2^k$-terms into $2^{k+1}$-terms we get the balance equation $$x_k+c_k=b_k+2c_{k+1},$$ where $c_k$ is the carry entering bit $k$, and $c_{k+1}$ is the carry leaving it....

Detailed mathematical approach

Problem Summary

Let $p(n)$ denote the number of ways to write $n$ as a sum of powers of $2$, with unlimited repetition allowed. Problem 890 asks for

$$p\left(7^{777}\right)\pmod{10^9+7}.$$

The classical generating function is

$$\prod_{j\ge 0}\frac{1}{1-x^{2^j}}=\sum_{n\ge 0}p(n)x^n.$$

Because $7^{777}$ is enormous, the implementation cannot iterate up to $n$. Instead, it reads the binary digits of $n$ and counts valid carry patterns. That turns the problem into a digit DP whose size depends on the bit-length of $n$, not on $n$ itself.

Mathematical Approach

The key observation is that a binary partition can be read bit by bit. At each bit position we only need to know how many lower-power terms have paired up and carried into the next position.

Step 1: Encode a Partition as a Carry Process

Write

$$n=\sum_{k=0}^{L-1} b_k2^k,\qquad b_k\in\{0,1\}.$$

If a partition uses $x_k$ copies of $2^k$, then after combining pairs of $2^k$-terms into $2^{k+1}$-terms we get the balance equation

$$x_k+c_k=b_k+2c_{k+1},$$

where $c_k$ is the carry entering bit $k$, and $c_{k+1}$ is the carry leaving it. For fixed $c_k$ and $c_{k+1}$, there is exactly one possible value of $x_k$, namely

$$x_k=b_k+2c_{k+1}-c_k.$$

This value is valid if and only if it is nonnegative, so

$$c_k\le 2c_{k+1}+b_k.$$

Now define $F_k(c)$ to be the number of ways to satisfy the lowest $k$ bits of $n$ and end with carry $c$ into bit $k$. Then for $k\ge 1$,

$$F_{k+1}(u)=\sum_{c=0}^{2u+b_k}F_k(c).$$

After the least significant bit there is always exactly one choice once the outgoing carry is fixed, so

$$F_1(c)=1\qquad(c\ge 0).$$

This is also why the least significant bit does not need a special case later: both $p(2m)$ and $p(2m+1)$ start from the same initial carry polynomial.

Step 2: Expand Each Carry Function in the Binomial Basis

The recurrence is a prefix sum evaluated at $2u+b_k$, so starting from the constant function $F_1(c)=1$, each step raises the degree by at most one. Therefore $F_k(c)$ is always a polynomial of degree at most $k-1$.

The implementation stores this polynomial in the basis

$$\binom{c}{0},\binom{c}{1},\binom{c}{2},\dots$$

so that

$$F_k(c)=\sum_{j=0}^{k-1}A_{k,j}\binom{c}{j}.$$

This basis is ideal because of the hockey-stick identity

$$\sum_{c=0}^{M}\binom{c}{j}=\binom{M+1}{j+1}.$$

Substituting the binomial expansion into the recurrence gives

$$F_{k+1}(u)=\sum_{j=0}^{k-1}A_{k,j}\binom{2u+b_k+1}{j+1}.$$

So the whole problem becomes: given the coefficient vector $(A_{k,0},\dots,A_{k,k-1})$, compute the next coefficient vector efficiently.

Step 3: Derive the Transition for a 0 Bit

If the current bit is $b_k=0$, then

$$F_{k+1}(u)=\sum_{j=0}^{k-1}A_{k,j}\binom{2u+1}{j+1}.$$

Using Pascal's identity,

$$\binom{2u+1}{j+1}=\binom{2u}{j+1}+\binom{2u}{j}.$$

Hence if we define an intermediate coefficient list by

$$B_r=A_{k,r}+A_{k,r-1},$$

with the convention that terms outside the valid range are $0$, then

$$F_{k+1}(u)=\sum_{r=0}^{k}B_r\binom{2u}{r}.$$

We now need to re-expand $\binom{2u}{r}$ in the basis $\binom{u}{j}$. Start from

$$\sum_{r\ge 0}\binom{2u}{r}t^r=(1+t)^{2u}=((1+t)^2)^u=(1+2t+t^2)^u.$$

Expanding again,

$$ (1+2t+t^2)^u=\sum_{j\ge 0}\binom{u}{j}(2t+t^2)^j=\sum_{j\ge 0}\binom{u}{j}\sum_{m=0}^{j}\binom{j}{m}2^{j-m}t^{j+m}. $$

Therefore the coefficient of $\binom{u}{j}$ inside $\binom{2u}{j+m}$ is

$$w_j(m)=\binom{j}{m}2^{j-m}.$$

So the next coefficient vector is

$$A_{k+1,j}=\sum_{m\ge 0} B_{j+m}\,w_j(m).$$

In the finite arrays used by the implementation, the sum stops at

$$m_{\max}=\min(j,\;k-j).$$

Step 4: Derive the Transition for a 1 Bit

If the current bit is $b_k=1$, then

$$F_{k+1}(u)=\sum_{j=0}^{k-1}A_{k,j}\binom{2u+2}{j+1}.$$

Applying Pascal once more gives

$$\binom{2u+2}{j+1}=\binom{2u+1}{j+1}+\binom{2u+1}{j},$$

so the same intermediate coefficients $B_r$ appear and

$$F_{k+1}(u)=\sum_{r=0}^{k}B_r\binom{2u+1}{r}.$$

Now use

$$\binom{2u+1}{r}=\binom{2u}{r}+\binom{2u}{r-1}.$$

This means the kernel for bit $1$ is just the sum of two neighboring bit-$0$ kernels:

$$w^{(1)}_j(m)=w_j(m)+w_j(m-1),$$

where $w_j(m)=0$ whenever $m<0$ or $m>j$. Thus

$$A_{k+1,j}=\sum_{m\ge 0} B_{j+m}\,w^{(1)}_j(m).$$

In the finite implementation, the upper limit is

$$m_{\max}=\min(j+1,\;k-j).$$

Step 5: Extract the Final Answer

After all $L$ bits have been processed, there can be no carry beyond the most significant bit. Therefore the desired count is

$$p(n)=F_L(0).$$

In the binomial basis, this is especially simple:

$$F_L(0)=\sum_{j=0}^{L-1}A_{L,j}\binom{0}{j}=A_{L,0},$$

because $\binom{0}{0}=1$ and $\binom{0}{j}=0$ for every $j>0$. That is why the implementation returns the first coefficient of the last state vector.

Worked Example: $n=7$

Take $n=7=111_2$. We begin with

$$F_1(c)=1.$$

The second bit is $1$, so

$$F_2(u)=\sum_{c=0}^{2u+1}1=2u+2.$$

In binomial form,

$$F_2(u)=2\binom{u}{0}+2\binom{u}{1}.$$

The third bit is again $1$, therefore

$$F_3(u)=\sum_{c=0}^{2u+1}(2c+2)=(2u+2)(2u+3)=4u^2+10u+6.$$

Convert this polynomial back to the binomial basis:

$$F_3(u)=6\binom{u}{0}+14\binom{u}{1}+8\binom{u}{2}.$$

Hence

$$p(7)=F_3(0)=6,$$

which matches the known value and the implementation checkpoint.

How the Code Works

The C++, Python, and Java implementations all target the same digit-DP formula. The C++ and Java implementations first convert $n$ to binary, least significant bit first. They then precompute powers of $2$ modulo $10^9+7$, Pascal coefficients modulo $10^9+7$, and two lower-triangular transition tables corresponding to the formulas for a current bit of $0$ and $1$.

The running state is the coefficient list of $F_k(c)$ in the basis $\binom{c}{j}$. For each new bit, the implementation first applies the Pascal update $B_r=A_r+A_{r-1}$, then performs the appropriate convolution against the precomputed kernel. The C++ implementation can split the outer coefficient range across several threads, while the Java implementation performs the same arithmetic serially.

The Python implementation is a thin execution bridge: it compiles and runs the C++ solver when necessary, then parses the numeric output. In every language, the published answer is the first coefficient after the most significant bit has been processed.

Complexity Analysis

Let $L=\lfloor\log_2 n\rfloor+1$. Precomputing powers of $2$, Pascal coefficients, and the two triangular transition tables costs $O(L^2)$ time and $O(L^2)$ memory. At stage $k$, the convolution examines $\Theta(k^2)$ coefficient pairs, so the full dynamic program costs

$$\sum_{k=1}^{L-1}\Theta(k^2)=\Theta(L^3)$$

time. The optional threading in the C++ implementation improves wall-clock time but does not change the asymptotic bound.

Footnotes and References

Project Euler Problem 890: https://projecteuler.net/problem=890
OEIS A000123, binary partitions: https://oeis.org/A000123
Generating function: Wikipedia — Generating function
Binomial coefficient: Wikipedia — Binomial coefficient
Pascal's triangle: Wikipedia — Pascal's triangle

Problem 890 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <thread>
#include <vector>

#include <boost/multiprecision/cpp_int.hpp>

using boost::multiprecision::cpp_int;

namespace {

constexpr int64_t kMod = 1'000'000'007;

std::vector<int> ToBitsLSB(cpp_int n) {
  std::vector<int> bits;
  if (n == 0) {
    return bits;
  }
  while (n > 0) {
    bits.push_back(static_cast<int>(n & 1));
    n >>= 1;
  }
  return bits;
}

int64_t ModAdd(int64_t a, int64_t b) {
  int64_t s = a + b;
  if (s >= kMod) {
    s -= kMod;
  }
  return s;
}

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

int64_t ComputePartitions(const cpp_int &n, int thread_count = 1) {
  if (n == 0) {
    return 1;
  }

  std::vector<int> bits = ToBitsLSB(n);
  if (bits.size() == 1) {
    return 1;
  }

  const int L = static_cast<int>(bits.size());

  std::vector<int64_t> pow2(L + 3, 1);
  for (int i = 1; i < static_cast<int>(pow2.size()); ++i) {
    pow2[i] = ModAdd(pow2[i - 1], pow2[i - 1]);
  }

  std::vector<std::vector<int64_t>> C(L + 1, std::vector<int64_t>(L + 1, 0));
  for (int i = 0; i <= L; ++i) {
    C[i][0] = 1;
    for (int j = 1; j <= i; ++j) {
      C[i][j] = ModAdd(C[i - 1][j - 1], C[i - 1][j]);
    }
  }

  std::vector<std::vector<int64_t>> w(L + 1);
  std::vector<std::vector<int64_t>> w2(L + 1);
  for (int j = 0; j <= L; ++j) {
    w[j].resize(j + 1);
    for (int m = 0; m <= j; ++m) {
      w[j][m] = ModMul(C[j][m], pow2[j - m]);
    }
    w2[j].resize(j + 2);
    w2[j][0] = w[j][0];
    for (int m = 1; m <= j; ++m) {
      int64_t val = ModAdd(w[j][m], w[j][m - 1]);
      w2[j][m] = val;
    }
    w2[j][j + 1] = w[j][j];
  }

  std::vector<std::vector<int64_t>>().swap(C);

  std::vector<int64_t> a(1, 1);

  auto compute_range = [&](const std::vector<int64_t> &b,
                           std::vector<int64_t> &next, int bit, int d, int start,
                           int end) {
    for (int j = start; j < end; ++j) {
      int max_m = 0;
      if (bit == 0) {
        max_m = std::min(j, d + 1 - j);
      } else {
        max_m = std::min(j + 1, d + 1 - j);
      }
      int64_t sum = 0;
      if (bit == 0) {
        const auto &wj = w[j];
        for (int m = 0; m <= max_m; ++m) {
          int i = j + m;
          sum += ModMul(b[i], wj[m]);
          if (sum >= kMod) {
            sum -= kMod;
          }
        }
      } else {
        const auto &w2j = w2[j];
        for (int m = 0; m <= max_m; ++m) {
          int i = j + m;
          sum += ModMul(b[i], w2j[m]);
          if (sum >= kMod) {
            sum -= kMod;
          }
        }
      }
      next[j] = sum;
    }
  };

  const int min_parallel = 128;
  if (thread_count < 1) {
    thread_count = 1;
  }

  for (int idx = 1; idx < L; ++idx) {
    const int bit = bits[idx];
    const int d = static_cast<int>(a.size()) - 1;

    std::vector<int64_t> b(d + 2, 0);
    b[0] = a[0];
    for (int i = 1; i <= d; ++i) {
      b[i] = ModAdd(a[i], a[i - 1]);
    }
    b[d + 1] = a[d];

    std::vector<int64_t> next(d + 2, 0);

    if (thread_count == 1 || d + 2 < min_parallel) {
      compute_range(b, next, bit, d, 0, d + 2);
    } else {
      const int threads = std::min(thread_count, d + 2);
      const int chunk = (d + 2 + threads - 1) / threads;
      std::vector<std::thread> pool;
      pool.reserve(threads);
      for (int t = 0; t < threads; ++t) {
        int start = t * chunk;
        int end = std::min(d + 2, start + chunk);
        if (start >= end) {
          continue;
        }
        pool.emplace_back(compute_range, std::cref(b), std::ref(next), bit, d,
                          start, end);
      }
      for (auto &th : pool) {
        th.join();
      }
    }

    a.swap(next);
  }

  return a[0] % kMod;
}

bool Validate() {
  cpp_int n7 = 7;
  if (ComputePartitions(n7) != 6) {
    std::cerr << "Validation failed: p(7) != 6\n";
    return false;
  }

  cpp_int n7_7 = 1;
  for (int i = 0; i < 7; ++i) {
    n7_7 *= 7;
  }
  if (ComputePartitions(n7_7) != 144548435) {
    std::cerr << "Validation failed: p(7^7) != 144548435\n";
    return false;
  }

  return true;
}

}  // namespace

int main() {
  if (!Validate()) {
    return 1;
  }

  cpp_int n = 1;
  for (int i = 0; i < 777; ++i) {
    n *= 7;
  }

  int threads = static_cast<int>(std::thread::hardware_concurrency());
  if (threads <= 0) {
    threads = 1;
  }

  int64_t result = ComputePartitions(n, threads);
  std::cout << result << "\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

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

public class Euler890 {

    static final long kMod = 1000000007L;

    static List<Integer> toBitsLSB(BigInteger n) {
        List<Integer> bits = new ArrayList<>();
        if (n.equals(BigInteger.ZERO)) {
            return bits;
        }
        while (n.compareTo(BigInteger.ZERO) > 0) {
            bits.add(n.testBit(0) ? 1 : 0);
            n = n.shiftRight(1);
        }
        return bits;
    }

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

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

    static long computePartitions(BigInteger n) {
        if (n.equals(BigInteger.ZERO)) {
            return 1;
        }

        List<Integer> bits = toBitsLSB(n);
        if (bits.size() == 1) {
            return 1;
        }

        int L = bits.size();

        long[] pow2 = new long[L + 3];
        pow2[0] = 1;
        for (int i = 1; i < pow2.length; ++i) {
            pow2[i] = modAdd(pow2[i - 1], pow2[i - 1]);
        }

        long[][] C = new long[L + 1][L + 1];
        for (int i = 0; i <= L; ++i) {
            C[i][0] = 1;
            for (int j = 1; j <= i; ++j) {
                C[i][j] = modAdd(C[i - 1][j - 1], C[i - 1][j]);
            }
        }

        long[][] w = new long[L + 1][];
        long[][] w2 = new long[L + 1][];
        for (int j = 0; j <= L; ++j) {
            w[j] = new long[j + 1];
            for (int m = 0; m <= j; ++m) {
                w[j][m] = modMul(C[j][m], pow2[j - m]);
            }
            w2[j] = new long[j + 2];
            w2[j][0] = w[j][0];
            for (int m = 1; m <= j; ++m) {
                w2[j][m] = modAdd(w[j][m], w[j][m - 1]);
            }
            w2[j][j + 1] = w[j][j];
        }

        long[] a = { 1 };

        for (int idx = 1; idx < L; ++idx) {
            int bit = bits.get(idx);
            int d = a.length - 1;

            long[] b = new long[d + 2];
            b[0] = a[0];
            for (int i = 1; i <= d; ++i) {
                b[i] = modAdd(a[i], a[i - 1]);
            }
            b[d + 1] = a[d];

            long[] nextA = new long[d + 2];

            for (int j = 0; j < d + 2; ++j) {
                int maxM = 0;
                if (bit == 0) {
                    maxM = Math.min(j, d + 1 - j);
                } else {
                    maxM = Math.min(j + 1, d + 1 - j);
                }

                long sum = 0;
                if (bit == 0) {
                    long[] wj = w[j];
                    for (int m = 0; m <= maxM; ++m) {
                        int i = j + m;
                        sum += modMul(b[i], wj[m]);
                        if (sum >= kMod) {
                            sum -= kMod;
                        }
                    }
                } else {
                    long[] w2j = w2[j];
                    for (int m = 0; m <= maxM; ++m) {
                        int i = j + m;
                        sum += modMul(b[i], w2j[m]);
                        if (sum >= kMod) {
                            sum -= kMod;
                        }
                    }
                }
                nextA[j] = sum;
            }

            a = nextA;
        }

        return a[0] % kMod;
    }

    public static String solve() {
        BigInteger n = BigInteger.valueOf(7).pow(777);
        return Long.toString(computePartitions(n));
    }

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