Problem 884: Removing Cubes

View on Project Euler

Project Euler Problem 884 Solution

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

Problem Summary Let \(d(n)\) be the number of greedy steps needed to reduce \(n\) to \(0\) by repeatedly subtracting the largest perfect cube not exceeding the current value. The quantity required by the problem is $$T(N)=\sum_{n=1}^{N-1} d(n),$$ with the final target \(N=10^{17}\). A direct simulation for every \(n<N\) is hopeless, so the solution groups numbers by the first cube removed and turns the summation into a much smaller recursive dynamic program. Mathematical Approach For bookkeeping it is convenient to set \(d(0)=0\) and to define the partial-sum function $$T(x)=\sum_{n=1}^{x-1} d(n)\qquad (x\ge 1).$$ The answer we want is therefore \(T(10^{17})\). Step 1: Partition by the first removed cube Take any \(x>1\), and let $$m=\left\lfloor\sqrt[3]{x-1}\right\rfloor,\qquad b=m^3,\qquad r=x-b.$$ Then \(b\le x-1 < (m+1)^3\), so every integer in the interval $$b,\ b+1,\ \dots,\ x-1$$ has the same first greedy move: subtract \(b=m^3\)....

Detailed mathematical approach

Problem Summary

Let \(d(n)\) be the number of greedy steps needed to reduce \(n\) to \(0\) by repeatedly subtracting the largest perfect cube not exceeding the current value. The quantity required by the problem is

$$T(N)=\sum_{n=1}^{N-1} d(n),$$

with the final target \(N=10^{17}\). A direct simulation for every \(n<N\) is hopeless, so the solution groups numbers by the first cube removed and turns the summation into a much smaller recursive dynamic program.

Mathematical Approach

For bookkeeping it is convenient to set \(d(0)=0\) and to define the partial-sum function

$$T(x)=\sum_{n=1}^{x-1} d(n)\qquad (x\ge 1).$$

The answer we want is therefore \(T(10^{17})\).

Step 1: Partition by the first removed cube

Take any \(x>1\), and let

$$m=\left\lfloor\sqrt[3]{x-1}\right\rfloor,\qquad b=m^3,\qquad r=x-b.$$

Then \(b\le x-1 < (m+1)^3\), so every integer in the interval

$$b,\ b+1,\ \dots,\ x-1$$

has the same first greedy move: subtract \(b=m^3\).

Any number in that block can be written as \(b+u\) with \(0\le u<r\), and its greedy decomposition begins with

$$b+u \longrightarrow u.$$

Therefore

$$d(b+u)=1+d(u).$$

Step 2: Convert that block identity into a summation recurrence

Summing the previous identity over \(u=0,1,\dots,r-1\) gives

$$T(x)-T(b)=\sum_{u=0}^{r-1}\bigl(1+d(u)\bigr)=r+T(r).$$

Now define the cube-boundary values

$$B_m=T(m^3).$$

Then the main recurrence becomes

$$\boxed{T(x)=B_m+r+T(r),\qquad m=\left\lfloor\sqrt[3]{x-1}\right\rfloor,\ r=x-m^3.}$$

This is the central formula used by the implementations.

Step 3: Derive the boundary recurrence

At an exact cube boundary, take \(x=(m+1)^3\). Then

$$r=(m+1)^3-m^3=3m^2+3m+1.$$

Substituting into the main recurrence yields

$$B_{m+1}=B_m+\bigl((m+1)^3-m^3\bigr)+T\bigl((m+1)^3-m^3\bigr).$$

So once the boundary values are computed in increasing order, each new boundary depends only on an already smaller argument. This turns the original problem into a forward sweep over consecutive cubes.

Step 4: Why the recursive argument shrinks quickly

The remainder satisfies

$$1\le r\le (m+1)^3-m^3=3m^2+3m+1.$$

Since \(m\approx x^{1/3}\), we get the rough bound

$$r\le 3x^{2/3}+3x^{1/3}+1.$$

That means one recursive step reduces an argument of size about \(x\) to a remainder of size about \(x^{2/3}\). Repeating the process shrinks the problem very fast, which is why memoization is enough to make the method practical even for \(10^{17}\).

Step 5: Turn the recurrence into a memoized dynamic program

Let

$$q=\left\lfloor\sqrt[3]{N-1}\right\rfloor.$$

The algorithm first computes all boundary values

$$B_1,\ B_2,\ \dots,\ B_q$$

from left to right. Whenever \(T(x)\) is needed for a non-boundary argument, the same recurrence

$$T(x)=B_m+r+T(r)$$

is applied and the result is stored. Because many boundary gaps and remainders recur, caching turns a potentially branching recursion tree into a directed acyclic graph of already solved states.

Worked Example

Take \(x=10\). Then

$$m=\left\lfloor\sqrt[3]{9}\right\rfloor=2,\qquad b=2^3=8,\qquad r=10-8=2.$$

The recurrence says

$$T(10)=T(8)+2+T(2).$$

Now \(T(8)=d(1)+\cdots+d(7)=1+2+3+4+5+6+7=28\), and \(T(2)=d(1)=1\). Hence

$$T(10)=28+2+1=31.$$

This matches the direct greedy counts: \(d(8)=1\) and \(d(9)=2\), so the last two numbers contribute \(3\) on top of \(T(8)\).

How the Code Works

The C++, Python, and Java implementations all follow the same structure. First they compute an exact floor cube root for a given integer by starting from a floating estimate and then correcting upward or downward until the surrounding cubes are in the right order. That avoids boundary mistakes near perfect cubes.

Next they build an array of cube-boundary totals \(B_m=T(m^3)\) up to \(m=\lfloor\sqrt[3]{N-1}\rfloor\). Each new entry uses the boundary recurrence, so the program only needs the gap between consecutive cubes and one recursive partial sum for that gap.

A hash table or dictionary stores every already computed value of \(T(x)\). When the same remainder appears again, the implementation returns the cached value immediately instead of recomputing the whole chain.

After the boundary sweep is complete, the program applies the same recurrence once more to the target \(N=10^{17}\). The C++ version also includes small brute-force checks for tiny inputs to confirm that the optimized recurrence matches direct simulation.

Complexity Analysis

Let \(q=\lfloor\sqrt[3]{N-1}\rfloor\), and let \(M\) be the number of distinct arguments whose values are actually memoized. The deterministic preprocessing is the forward sweep over the \(q\) cube boundaries, so that part is \(O(q)\).

For any uncached argument \(x\), one recurrence step replaces \(x\) by a remainder \(r\) with size \(O(x^{2/3})\), so an individual recursion chain has depth \(O(\log\log N)\). Since each distinct state is solved only once, the total running time is \(O(M)\) cube-root evaluations plus \(O(q)\) boundary updates, and the memory usage is \(O(M+q)\). In practice \(M\) stays close to the cube-root scale and is tiny compared with \(N\).

Footnotes and References

  1. Project Euler problem page: https://projecteuler.net/problem=884
  2. Cube number: Wikipedia — Cube number
  3. Greedy algorithm: Wikipedia — Greedy algorithm
  4. Memoization: Wikipedia — Memoization
  5. Dynamic programming: Wikipedia — Dynamic programming

Problem 884 source code

C++

#include <cassert>
#include <cmath>
#include <cstdint>
#include <functional>
#include <iostream>
#include <string>
#include <unordered_map>
#include <vector>
#include <algorithm>

using u64 = std::uint64_t;
using u128 = unsigned __int128;

static u64 cube_u64(u64 k) {
    return static_cast<u64>(static_cast<u128>(k) * k * k);
}

static u64 floor_cuberoot(u64 x) {
    long double y = std::cbrt(static_cast<long double>(x));
    u64 k = static_cast<u64>(y);
    while (cube_u64(k + 1) <= x) ++k;
    while (cube_u64(k) > x) --k;
    return k;
}

static u64 D_bruteforce(u64 n) {
    u64 steps = 0;
    while (n > 0) {
        u64 k = floor_cuberoot(n);
        n -= cube_u64(k);
        ++steps;
    }
    return steps;
}

static u128 S_bruteforce(u64 N) {
    u128 s = 0;
    for (u64 n = 1; n < N; ++n) s += D_bruteforce(n);
    return s;
}

static u128 S_fast(u64 N) {
    if (N <= 1) return 0;

    u64 K = floor_cuberoot(N - 1);
    std::vector<u128> A(static_cast<std::size_t>(K + 1), 0);

    std::unordered_map<u64, u128> memo;
    memo.reserve(static_cast<std::size_t>(K * 3));
    memo[0] = 0;
    memo[1] = 0;
    if (K >= 1) A[1] = 0;

    std::function<u128(u64, u64)> solve = [&](u64 x, u64 ready_k) -> u128 {
        if (x <= 1) return 0;
        auto it = memo.find(x);
        if (it != memo.end()) return it->second;

        u64 k = floor_cuberoot(x - 1);
        assert(k <= ready_k);
        u64 c = cube_u64(k);
        u64 t = x - c;

        u128 ans = A[static_cast<std::size_t>(k)] + static_cast<u128>(t) + solve(t, ready_k);
        memo[x] = ans;
        return ans;
    };

    for (u64 k = 1; k < K; ++k) {
        u64 delta = 3 * k * k + 3 * k + 1;
        u128 sdelta = solve(delta, k);
        A[static_cast<std::size_t>(k + 1)] = A[static_cast<std::size_t>(k)] + static_cast<u128>(delta) + sdelta;
        memo[cube_u64(k + 1)] = A[static_cast<std::size_t>(k + 1)];
    }

    return solve(N, K);
}

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

int main() {
    assert(S_fast(100) == 512);
    for (u64 n : {2ULL, 10ULL, 100ULL, 500ULL, 1000ULL}) {
        assert(S_fast(n) == S_bruteforce(n));
    }

    std::cout << to_string_u128(S_fast(100'000'000'000'000'000ULL)) << '\n';
    return 0;
}

Python

import sys

sys.setrecursionlimit(2000000)

def floor_cuberoot(x):
    y = int(x ** (1.0/3.0))
    k = y
    while (k + 1) ** 3 <= x:
        k += 1
    while k ** 3 > x:
        k -= 1
    return k

def S_fast(N):
    if N <= 1:
        return 0

    K = floor_cuberoot(N - 1)
    A = [0] * (K + 1)
    
    memo = {}
    memo[0] = 0
    memo[1] = 0

    def solve(x, ready_k):
        if x <= 1:
            return 0
        if x in memo:
            return memo[x]

        k = floor_cuberoot(x - 1)
        c = k ** 3
        t = x - c

        ans = A[k] + t + solve(t, ready_k)
        memo[x] = ans
        return ans

    for k in range(1, K):
        delta = 3 * k * k + 3 * k + 1
        sdelta = solve(delta, k)
        A[k + 1] = A[k] + delta + sdelta
        memo[(k + 1) ** 3] = A[k + 1]

    return solve(N, K)

def solve():
    return str(S_fast(10**17))

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

Java

import java.util.HashMap;

public class Euler884 {

    static long floorCuberoot(long x) {
        long k = (long) Math.cbrt(x);
        while (cube(k + 1) <= x && cube(k + 1) > 0)
            ++k;
        while (cube(k) > x || cube(k) < 0)
            --k;
        return k;
    }

    static long cube(long k) {
        return k * k * k; // Assumes it won't overflow for the values we need
    }

    static class Solver {
        long[] A;
        HashMap<Long, Long> memo = new HashMap<>(30000);

        long solveRec(long x, long readyK) {
            if (x <= 1)
                return 0;
            Long cached = memo.get(x);
            if (cached != null)
                return cached;

            long k = floorCuberoot(x - 1);
            long c = cube(k);
            long t = x - c;

            long ans = A[(int) k] + t + solveRec(t, readyK);
            memo.put(x, ans);
            return ans;
        }

        long SFast(long N) {
            if (N <= 1)
                return 0;

            long K = floorCuberoot(N - 1);
            A = new long[(int) K + 1];

            memo.put(0L, 0L);
            memo.put(1L, 0L);

            for (long k = 1; k < K; ++k) {
                long delta = 3 * k * k + 3 * k + 1;
                long sdelta = solveRec(delta, k);
                A[(int) (k + 1)] = A[(int) k] + delta + sdelta;
                memo.put(cube(k + 1), A[(int) (k + 1)]);
            }

            return solveRec(N, K);
        }
    }

    public static String solve() {
        Solver solver = new Solver();
        return Long.toString(solver.SFast(100000000000000000L));
    }

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