Problem 341: Golomb's Self-describing Sequence

Project Euler Problem 341 Solution

EulerSolve provides an optimized solution for Project Euler Problem 341, Golomb's Self-describing Sequence, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The Golomb sequence $G(n)$ is the unique nondecreasing sequence of positive integers in which the integer $n$ appears exactly $G(n)$ times. Its beginning is $$1,2,2,3,3,4,4,4,5,5,5,6,6,6,6,\dots$$ Project Euler 341 asks for $$\sum_{1\le n \lt 10^6} G(n^3).$$ The official checkpoints $G(10^3)=86$, $G(10^6)=6137$, and $\sum_{1\le n \lt 10^3}G(n^3)=153506976$ show that the queried positions become enormous very quickly. A naive approach that expands the sequence all the way to the largest queried position $(10^6-1)^3$ would be hopelessly too slow and too memory-hungry. Mathematical Approach Self-Describing Structure By definition, value $k$ occurs exactly $G(k)$ times. The implementation generates the sequence with the classical recurrence $$G(1)=1,\qquad G(n)=1+G\bigl(n-G(G(n-1))\bigr)\quad(n\ge 2).$$ This recurrence is equivalent to the usual Project Euler formulation and allows us to precompute $G(1),G(2),\dots$ once. First Prefix Sum: Value Boundaries Define the first prefix sum $$S(k)=\sum_{i=1}^{k}G(i).$$ Since the value $k$ appears $G(k)$ times, the positions holding $k$ are exactly $$S(k-1)+1,\ S(k-1)+2,\ \dots,\ S(k).$$ Therefore $$\boxed{G(n)=k\iff S(k-1)\lt n\le S(k).}$$ So $S(k)$ marks where the block of value $k$ ends. Second Prefix Sum: Blocking the Indices Themselves Now look at the indices $m$ for which $G(m)=k$....

Detailed mathematical approach

Problem Summary

Golomb's self-describing sequence $G$ is the only nondecreasing sequence of natural numbers in which the number $n$ occurs exactly $G(n)$ times. Its first fifteen terms are

$$1,2,2,3,3,4,4,4,5,5,5,6,6,6,6,\ldots$$

and they already describe themselves: $1$ occurs once, $2$ occurs twice and $4$ occurs three times, exactly as $G(1)=1$, $G(2)=2$ and $G(4)=3$ demand. We must evaluate

$$\sum_{1\le n\lt 10^6}G(n^3).$$

The statement gives three checkpoints, ^[1]

$$G(10^3)=86,\qquad G(10^6)=6137,$$

$$\sum_{1\le n\lt 10^3}G(n^3)=153506976.$$

The largest position asked about is

$$(10^6-1)^3=999997000002999999,$$

so writing the sequence out that far is impossible: even four bytes per term would come to four exabytes, and the terms near that position need more than $32$ bits anyway.

The whole solution in one sentence: equal terms of $G$ sit in consecutive blocks, the blocks of length at most $k$ fill exactly the first $P(k)=\sum_{i\le k}i\,G(i)$ positions, and since $P$ passes $10^{18}$ after only about $10^7$ steps, a table of that size puts any requested position into one such group of blocks, where a single division finishes the query.

Mathematical Approach

1. Every term is forced

Nothing has to be guessed: the definition determines the sequence term by term. First, $G(1)=1$. Suppose instead $G(1)=c\ge2$. The sequence is nondecreasing, so every term is at least $c\ge2$ and the number $1$ never occurs; but it must occur $G(1)=c\ge2$ times, a contradiction. So $G(1)=1$ and the value $1$ occurs exactly once.

Next, every natural number occurs: the terms are natural numbers, so $G(n)\ge1$, and $n$ occurs $G(n)\ge1$ times. Since the sequence is nondecreasing, all copies of a value are consecutive, so it splits into blocks with no gaps and no value missing: one copy of $1$, then $G(2)$ copies of $2$, then $G(3)$ copies of $3$, and so on.

The blocks can therefore be laid down from the left, which is already an algorithm. $G(2)\ne1$, which would be a second copy of $1$, so $G(2)\ge2$; the block of $2$ runs from position $2$ to $1+G(2)$ and hence contains position $2$, forcing $G(2)=2$. The start is $1,2,2$. Nothing is circular afterwards: for $k\ge3$ the blocks of $1,\ldots,k-1$ fill at least $1+2+(k-3)=k$ positions, since every block has positive length, so position $k$ is already written when the block of $k$ is appended and its length $G(k)$ is read off. Every step is forced, which proves existence and uniqueness.

2. The first prefix sum locates a value

Let

$$S(k)=\sum_{i=1}^{k}G(i),\qquad S(0)=0.$$

The blocks of $1,\ldots,k$ have lengths $G(1),\ldots,G(k)$, so together they occupy positions $1$ to $S(k)$, and the block of $k$ occupies positions $S(k-1)+1$ to $S(k)$. Hence $G(n)=k$ holds precisely when

$$S(k-1)\lt n\le S(k),$$

or equivalently $G(n)=\min\{m:S(m)\ge n\}$, because $S$ is strictly increasing. The first values of $S$ are

$$1,3,5,8,11,15,19,23,28,33,\ldots$$

and $S(k)$ is exactly the position of the last copy of $k$. ^[4]

3. Neighbouring terms differ by at most one

Positions $n$ and $n+1$ either lie in the same block or in two consecutive blocks, since no block is empty. Therefore

$$G(n+1)-G(n)\in\{0,1\}.$$

A program confirms it for all $n\le3\,000\,000$.

4. Mallows' recurrence builds the table

The block construction of Section 1 is a correct algorithm, and the verification below uses it. The published programs instead use the explicit recurrence found by Colin Mallows, which starts from $G(1)=1$ and then sets ^[2] ^[3]

$$G(n+1)=1+G\bigl(n+1-G(G(n))\bigr).$$

Here is a proof from the block picture. Put $k=G(n)$, so $S(k-1)\lt n\le S(k)$, and note that $G(G(n))=G(k)$ is the length of the block of the value $k$.

If $n\lt S(k)$ then $k\ge2$, since $S(1)=1$ leaves no room below it, and $n+1$ is still in block $k$, so $G(n+1)=k$; we must show $n+1-G(k)$ lies in block $k-1$. From $n\le S(k)-1$, $n+1-G(k)\le S(k)-G(k)$${}=S(k-1)$. From $n\ge S(k-1)+1$, $n+1-G(k)$${}\ge S(k-1)+2-G(k)$, which exceeds $S(k-1)-G(k-1)=S(k-2)$ exactly when $G(k)\lt G(k-1)+2$, true by Section 3. So $G(n+1-G(k))=k-1$.
If $n=S(k)$, then $n+1$ starts block $k+1$, so $G(n+1)=k+1$, and $n+1-G(k)=S(k)+1-G(k)$, which is $S(k-1)+1$, the first position of block $k$, whose value is $k$.

In both cases $G(n+1)=1+G(n+1-G(G(n)))$. All three lookups use indices at most $n$, because $G(n)\le n$ and $G(k)\ge1$, so one left-to-right pass fills a table of $G$ in linear time.

5. The second prefix sum

Inverting $S$ directly is hopeless here: the least $m$ with $S(m)\ge x$ is $G(x)$ itself, which for $x=(10^6-1)^3$ is $160113493649$. Those $1.6\cdot10^{11}$ terms all stay below $2^{32}$, so four bytes each would hold them, but that is $640$ GB.

The way out is to sum one level higher. Define

$$P(k)=\sum_{i=1}^{k}i\,G(i),\qquad P(0)=0,$$

whose first values are $1,5,11,23,38,62,90,122,\ldots$ ^[5] This is $S$ evaluated at the index $S(k)$:

$$P(k)=S\bigl(S(k)\bigr).$$

Indeed, among the positions $1,\ldots,S(k)$ the value $i$ occurs on exactly $G(i)$ of them for each $i\le k$, and on none of them for $i\gt k$, so

$$S(S(k))=\sum_{j\le S(k)}G(j)=\sum_{i\le k}i\,G(i).$$

For $k=5$ this reads $S(S(5))=S(11)=38=P(5)$. Said another way, $P(k)$ ends the last block of length $k$: the step $P(3)-P(2)=11-5=6$ is the combined length of the two blocks of length $3$, those of $4$ and $5$, at positions $6$ to $11$.

6. One division answers a query

Claim. If $P(k-1)\lt x\le P(k)$, then with $d=x-P(k-1)$,

$$G(x)=S(k-1)+\left\lceil\frac{d}{k}\right\rceil.$$

Proof. Write $m=G(x)$, the least index with $S(m)\ge x$. Since $x\le P(k)=S(S(k))$ we get $m\le S(k)$, and since $x\gt P(k-1)=S(S(k-1))$ we get $m\gt S(k-1)$. So $m$ lies in the block of the value $k$, which means $G(j)=k$ for every $j$ with $S(k-1)\lt j\le m$, and therefore

$$S(m)=P(k-1)+\bigl(m-S(k-1)\bigr)k.$$

The condition $S(m)\ge x$ now reads $m-S(k-1)\ge d/k$, and $m$ is the least index satisfying it. Rounding up gives the claim. $\square$

Once the block index $k$ is known, one subtraction, one division and one addition finish the query. For positive integers $a$ and $b$ the ceiling is evaluated with integer arithmetic alone through ^[6]

$$\lceil a/b\rceil=\lfloor(a+b-1)/b\rfloor.$$

7. Worked example: $G(1000)=86$

Continuing both prefix sums to $k=19$, their last three values are

$k=17$: $S(17)=76$, $P(17)=827$;
$k=18$: $S(18)=83$, $P(18)=953$;
$k=19$: $S(19)=90$, $P(19)=1086$.

For $x=1000$ we have $953\lt1000\le1086$, so $k=19$ and $d=1000-953=47$, which gives $G(1000)=83+\lceil47/19\rceil=86$, the value promised by the statement. The same two lines give $G(10)=3+\lceil5/3\rceil=5$ with $k=3$ and $d=5$, and $G(100)=19+\lceil10/8\rceil=21$ with $k=8$ and $d=10$; both agree with the sequence read off directly.

8. Why one extra level is enough

The growth rates explain the saving. $G(n)$ is known to behave like $\varphi^{2-\varphi}n^{\varphi-1}$ with $\varphi=(1+\sqrt5)/2$, ^[2] ^[3] so $S(k)$ grows like $k^{\varphi}\approx k^{1.618}$ and $P(k)$ like $k^{\varphi+1}\approx k^{2.618}$. Summing one level higher therefore reduces the table length from $x^{0.618}$ to $x^{0.382}$. Numerically, $\varphi^{2-\varphi}\approx1.201783$ predicts $G(10^3)\approx85.887$ and $G(10^6)\approx6138.05$ against the exact $86$ and $6137$.

Building $P$ until it passes $(10^6-1)^3$ takes exactly $K=10105256$ terms, with

$$P(K)=999997139415781941$$

and $P(K-1)=999996880327123357$, which is below $(10^6-1)^3$. Writing $C=\varphi^{2-\varphi}$ and summing $i\,G(i)\approx C\,i^{\varphi}$ gives

$$P(k)\approx\frac{C}{\varphi+1}k^{\varphi+1},$$

so the stopping index for a largest query $x$ is about $\bigl((\varphi+1)x/C\bigr)^{1/(\varphi+1)}$, which at $x=(10^6-1)^3$ is $10105299$, within $0.0005\%$ of $K$. The estimate only explains the scale; the programs stop on the exact test $P(k)\ge x$. The table is $K+1=10105257$ four-byte entries, $40.42$ MB instead of $640$ GB.

9. The queries arrive in increasing order

The requested positions $1^3,2^3,\ldots,(10^6-1)^3$ increase and $P$ increases, so the block index $k$ never moves backwards. Keeping $k$ with $S(k-1)$, $S(k)$, $P(k-1)$ and $P(k)$ in a small state and advancing it only when $P(k)\lt x$ costs $K-1=10105255$ advances in total, one for each block index from $2$ to $K$, as a run confirms. No binary search is needed. Writing $L$ for the limit, the sweep costs $O(L+K)$, that is $O(1+K/L)$ per query, and here $K/L\approx10.11$. The ratio is not bounded as $L$ grows, since $K$ scales like $L^{3/(\varphi+1)}\approx L^{1.146}$; the formula itself always costs constant time.

10. How large the numbers become

The table stays small: its largest entry is $G(K)=25639$, which fits in $16$ bits, so $32$-bit entries are generous. The prefix sums are the large quantities, $S(K)=160113507445$ and $P(K)$ just under $10^{18}$, and both fit in $64$ bits.

The requested total fits too. Since $G$ is nondecreasing, every term obeys $G(n^3)\le G((10^6-1)^3)$${}=160113493649$, and there are $999999$ terms, so the total is below

$$999999\cdot160113493649\lt1.61\cdot10^{17},$$

well inside the signed $64$-bit range $2^{63}-1\approx9.22\cdot10^{18}$.

How the Code Works

precompute_golomb_until implements Sections 4 and 5 in one loop. ^[9] It keeps a $1$-indexed table with golomb[0] as an unused placeholder and golomb[1]${}=1$, and a running products that is $P(i)$. For each $i$ it reads g_prev${}=G(i-1)$, then nested${}=G(G(i-1))$, sets idx${}=i-$nested and appends current${}=1+G($idx$)$. All three reads are below $i$, as Section 4 requires. The loop stops at the first $i$ with $P(i)\ge$ max_query_position, which is precisely $K$, so the table is exactly as long as the largest query needs. The C++ version reserves $11\,000\,000$ entries; the real need, $10105257$, fits, and std::vector would grow anyway if it did not.

The caller is main, which sets limit${}=10^6$ and passes max_query_position${}=($limit${}-1)^3$, the largest cube that will be asked for, so the table length follows the hardest query rather than a fixed guess.

QueryState is the small state of Section 9: index is $k$, sum and products are $S(k)$ and $P(k)$, and prev_sum and prev_products are $S(k-1)$ and $P(k-1)$. It starts at $k=1$ with $S(1)=P(1)=1$ and $S(0)=P(0)=0$.

golomb_at(x, golomb, st) is the query of Section 6. Its while loop advances index, shifting the current sums into the prev_ fields, until products${}\ge x$; then it returns prev_sum plus the integer ceiling (x - prev_products + index - 1) / index. For $x=1$ nothing is advanced and the result is $S(0)+1=1$.

sum_golomb_cubes(limit, golomb) creates one QueryState, walks $n=1,\ldots,$limit${}-1$, forms cube${}=n^3$ and adds golomb_at(cube, ...). Sharing the state across the loop is what makes the sweep amortised; a fresh state per query would restart every search at $k=1$.

run_checkpoints verifies the three values from the statement before the real computation, each with a fresh QueryState. That matters because the state only ever moves forward, while the third checkpoint sweeps from $n=1$ again. A failed check is reported on standard error. The flag --skip-checkpoints suppresses them and main rejects any other argument. The Python and Java files have no such checks; they go straight to solve.

The three versions share the recurrence and the sweep but differ in their checks, their names and their arithmetic; Java calls the same functions precomputeGolombUntil, golombAt and sumGolombCubes. The C++ accumulates products, the prefix sums and the result in unsigned __int128, ^[7] which is wider than needed but never has to be re-examined, and prints it with to_string_u128 because the standard streams cannot. Java keeps the table and the precomputation index in int and the positions, the state fields and the total in long, which Section 10 shows is safe; its int[] starts at $11\,000\,000$ entries and would double with Arrays.copyOf if that were too few. The doubling never happens here, but the closing Arrays.copyOf(golomb, size) allocates a second array of $40.42$ MB while the first is still alive. Python relies on unlimited integer precision ^[8] and keeps the table in a list, which costs far more than four bytes per entry.

Complexity and Verification

Building the table costs $O(K)$ time with three array reads per step, and the sweep costs $O($limit${}+K)$ because the block index advances $K-1$ times in total. With $K\approx1.01\cdot10^7$ the whole program is linear in $K$. Memory is $O(K)$ entries: the $K+1$ used ones are $40.42$ MB at four bytes each, the C++ reserve of $11\,000\,000$ is $44$ MB, and Java briefly holds both of its arrays, $84.42$ MB of payload. The three reads have good locality: golomb[i-1] sits next to the write position, golomb[idx] at most $638$ entries or $2552$ bytes behind it, and golomb[g_prev] at an index below $25639$, inside the first $100$ KB of the table.

Measured on one machine, the C++ program spends about $12$ ms building the $10105256$ table entries and $8$ ms on the sweep, $0.02$ s in all, with $42$ MB resident. Java takes $0.05$ s and $130$ MB, Python $2.0$ s and $421$ MB, the price of ten million integer objects. Elapsed times and resident sizes depend on the machine and the toolchain; the operation counts do not.

The following was checked by running code for this article. The sequence was rebuilt from its definition alone, by appending the value $k$ exactly $G(k)$ times, for $10^7$ positions: it is nondecreasing, the value $k$ occurs exactly $G(k)$ times for every $k\le25473$, and it agrees with Mallows' recurrence at every position. The block formula of Section 6 reproduces that independently built sequence at all $10^7$ positions, and $P(k)=S(S(k))$ holds for every one of the $25473$ indices $k$ with $S(k)\le10^7$. The identity $G(P(k))=S(k)$, which follows from Section 6, was confirmed for $k\le30$.

The checkpoints were then reproduced twice over, by two different inversions. Reading $G(x)$ off the first prefix sum by a search in $S$ — affordable for the checkpoint queries $x\le999^3$, where $437931$ terms suffice because $S(437931)=997003923$ — and reading it off the block formula give the same $G(10^3)=86$, $G(10^6)=6137$, the same $\sum_{1\le n\lt10^3}G(n^3)=153506976$ and the same $\sum_{1\le n\lt10^2}G(n^3)=211947$. For the smallest case the nine terms $G(1),G(8),G(27),\ldots,G(729)$ are $1,4,9,16,24,33,44,57,71$, summing to $259$. Finally, the three programs were run on the full problem and their outputs compared byte for byte: all three agree, and the C++ checkpoint routine passes.

Footnotes and References

The references below supply the background facts used above. The block structure, Mallows' recurrence as proved in Section 4, the identity $P(k)=S(S(k))$, the query formula and the counts and bounds are derived explicitly in this article.

Project Euler 341 — Golomb's Self-describing Sequence. The official statement defines G as the only nondecreasing sequence of natural numbers in which n appears exactly G(n) times, lists its first fifteen terms, gives G(10³) = 86, G(10⁶) = 6137 and the sum of G(n³) over 1 ≤ n < 10³ as 153506976, and asks for the same sum over 1 ≤ n < 10⁶. These three values test an implementation without disclosing the requested one.
Wikipedia — Golomb sequence. The article describes the sequence, also called Silverman's sequence, as the increasing sequence in which the n-th term is the number of times n appears — "increasing" there means nondecreasing, as the repeated terms show — attributes the recurrence of Section 4 to Colin Mallows, states that the terms grow asymptotically like φ^(2−φ)·n^(φ−1) with φ the golden ratio, and notes that the sequence is named after Solomon W. Golomb.
OEIS A001462 — Golomb's sequence. The entry lists the first terms, repeats Mallows' recurrence, and gives the asymptotic form a(n) = φ^(2−φ)·n^(φ−1) + E(n) with the error term bounded by Vardi as O(n^(φ−1)/log n), which is the growth rate used in Section 8.
OEIS A001463 — partial sums of Golomb's sequence. The entry defines the sequence as the partial sums of A001462 and states that its n-th term is the position of the last occurrence of n in A001462, which is the block boundary S(k) of Section 2.
OEIS A143125 — weighted sums of Golomb's sequence. The entry is defined as the sum of j·A001462(j) for j = 1 to n, with first terms 1, 5, 11, 23, 38, 62, 90, 122, …, the quantity P(k) of Section 5.
Wikipedia — Floor and ceiling functions. Among the quotient identities the article states ⌈n/m⌉ = ⌊(n+m−1)/m⌋ for positive m, which is how the ceiling in the query formula is evaluated with integer arithmetic alone.
GCC manual — 128-bit integers. The page documents __int128 and unsigned __int128 as a compiler extension available on targets with a wide enough integer mode; the C++ solution uses the unsigned form for its prefix sums and its result.
Java Language Specification — primitive types fixes the range of long as −9223372036854775808 to 9223372036854775807, which the bound of Section 10 stays well inside, and Python — numeric types states that integers have unlimited precision, so the Python version needs no bound at all.
C++ — Euler341.cpp. The linked immutable revision contains precompute_golomb_until, golomb_at, sum_golomb_cubes, run_checkpoints and the --skip-checkpoints flag described above. It is the implementation record for this article; the independent rebuild and the two-inversion comparison are extra validation, not tests packaged in that file.

Problem 341 source code

C++

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

namespace {

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

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

std::vector<u32> precompute_golomb_until(const u64 max_query_position) {
    // 1-indexed: golomb[n] = G(n)
    std::vector<u32> golomb;
    golomb.reserve(11000000);
    golomb.push_back(0U);
    golomb.push_back(1U);

    u128 products = 1;  // sum_{i=1..k} i * G(i)

    for (u64 i = 2; products < static_cast<u128>(max_query_position); ++i) {
        const u64 prev = i - 1;
        const u32 g_prev = golomb[static_cast<std::size_t>(prev)];
        const u32 nested = golomb[static_cast<std::size_t>(g_prev)];
        const u64 idx = i - static_cast<u64>(nested);
        const u32 current = static_cast<u32>(1U + golomb[static_cast<std::size_t>(idx)]);

        golomb.push_back(current);
        products += static_cast<u128>(i) * static_cast<u128>(current);
    }

    return golomb;
}

struct QueryState {
    u64 index = 1;      // current k
    u128 sum = 1;       // S(k) = sum G(i)
    u128 products = 1;  // P(k) = sum i*G(i)

    u128 prev_sum = 0;       // S(k-1)
    u128 prev_products = 0;  // P(k-1)
};

u64 golomb_at(const u64 x, const std::vector<u32>& golomb, QueryState& st) {
    while (st.products < static_cast<u128>(x)) {
        ++st.index;
        st.prev_sum = st.sum;
        st.prev_products = st.products;

        const u32 g = golomb[static_cast<std::size_t>(st.index)];
        st.sum += static_cast<u128>(g);
        st.products += static_cast<u128>(st.index) * static_cast<u128>(g);
    }

    const u128 span_index = static_cast<u128>(st.index);
    const u128 offset = (static_cast<u128>(x) - st.prev_products + span_index - 1U) / span_index;
    return static_cast<u64>(st.prev_sum + offset);
}

u128 sum_golomb_cubes(const u64 limit, const std::vector<u32>& golomb) {
    QueryState st;
    u128 sum = 0;

    for (u64 n = 1; n < limit; ++n) {
        const u64 cube = n * n * n;
        sum += static_cast<u128>(golomb_at(cube, golomb, st));
    }

    return sum;
}

bool run_checkpoints(const std::vector<u32>& golomb) {
    {
        QueryState st;
        if (golomb_at(1000ULL, golomb, st) != 86ULL) {
            std::cerr << "Checkpoint failed: G(10^3)\n";
            return false;
        }
    }

    {
        QueryState st;
        if (golomb_at(1000000ULL, golomb, st) != 6137ULL) {
            std::cerr << "Checkpoint failed: G(10^6)\n";
            return false;
        }
    }

    {
        const u128 sample = sum_golomb_cubes(1000ULL, golomb);
        if (sample != static_cast<u128>(153506976ULL)) {
            std::cerr << "Checkpoint failed: sum_{1<=n<10^3} G(n^3)\n";
            return false;
        }
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    const u64 limit = 1000000ULL;
    const u64 max_query_position = (limit - 1ULL) * (limit - 1ULL) * (limit - 1ULL);
    const std::vector<u32> golomb = precompute_golomb_until(max_query_position);

    if (!skip_checkpoints && !run_checkpoints(golomb)) {
        return 2;
    }

    const u128 answer = sum_golomb_cubes(limit, golomb);
    std::cout << to_string_u128(answer) << '\n';
    return 0;
}

Python

def precompute_golomb_until(max_query_position):
    golomb = [0, 1]
    
    products = 1
    i = 2
    
    while products < max_query_position:
        prev = i - 1
        g_prev = golomb[prev]
        nested = golomb[g_prev]
        idx = i - nested
        current = 1 + golomb[idx]
        
        golomb.append(current)
        products += i * current
        i += 1
        
    return golomb

class QueryState:
    def __init__(self):
        self.index = 1
        self.sum = 1
        self.products = 1
        self.prev_sum = 0
        self.prev_products = 0

def golomb_at(x, golomb, st):
    while st.products < x:
        st.index += 1
        st.prev_sum = st.sum
        st.prev_products = st.products
        
        g = golomb[st.index]
        st.sum += g
        st.products += st.index * g
        
    span_index = st.index
    offset = (x - st.prev_products + span_index - 1) // span_index
    return st.prev_sum + offset

def sum_golomb_cubes(limit, golomb):
    st = QueryState()
    total_sum = 0
    
    for n in range(1, limit):
        cube = n * n * n
        total_sum += golomb_at(cube, golomb, st)
        
    return total_sum

def solve():
    limit = 1000000
    max_query_position = (limit - 1) ** 3
    golomb = precompute_golomb_until(max_query_position)
    ans = sum_golomb_cubes(limit, golomb)
    return str(ans)

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

Java

import java.util.*;

public class Euler341 {

    static int[] precomputeGolombUntil(long maxQueryPosition) {
        int[] golomb = new int[11000000];
        golomb[0] = 0;
        golomb[1] = 1;
        int size = 2;

        long products = 1;
        int i = 2;

        while (products < maxQueryPosition) {
            int prev = i - 1;
            int g_prev = golomb[prev];
            int nested = golomb[g_prev];
            int idx = i - nested;
            int current = 1 + golomb[idx];

            if (size == golomb.length) {
                golomb = Arrays.copyOf(golomb, golomb.length * 2);
            }
            golomb[size++] = current;
            products += (long) i * current;
            i++;
        }

        return Arrays.copyOf(golomb, size);
    }

    static class QueryState {
        long index = 1;
        long sum = 1;
        long products = 1;
        long prev_sum = 0;
        long prev_products = 0;
    }

    static long golombAt(long x, int[] golomb, QueryState st) {
        while (st.products < x) {
            st.index++;
            st.prev_sum = st.sum;
            st.prev_products = st.products;

            int g = golomb[(int) st.index];
            st.sum += g;
            st.products += st.index * g;
        }

        long spanIndex = st.index;
        long offset = (x - st.prev_products + spanIndex - 1) / spanIndex;
        return st.prev_sum + offset;
    }

    static long sumGolombCubes(long limit, int[] golomb) {
        QueryState st = new QueryState();
        long sum = 0;

        for (long n = 1; n < limit; n++) {
            long cube = n * n * n;
            sum += golombAt(cube, golomb, st);
        }

        return sum;
    }

    public static String solve() {
        long limit = 1000000;
        long maxQueryPosition = (limit - 1) * (limit - 1) * (limit - 1);
        int[] golomb = precomputeGolombUntil(maxQueryPosition);
        long ans = sumGolombCubes(limit, golomb);
        return String.valueOf(ans);
    }

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