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

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$. From the previous identity, those indices are

$$m=S(k-1)+1,\dots,S(k),$$

and there are exactly $G(k)$ of them. Each such index $m$ contributes the value $k$, so the total length contributed by the entire $k$-block is $k\cdot G(k)$. This motivates the weighted prefix sum

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

After finishing all blocks up to $k$, the occupied sequence positions are precisely $1,\dots,P(k)$.

Inverting a $P$-Block

Suppose $x$ satisfies

$$P(k-1)\lt x\le P(k).$$

Then $G(x)$ must lie inside the block made of the values

$$S(k-1)+1,\dots,S(k),$$

where each value is repeated exactly $k$ times. Let

$$d=x-P(k-1).$$

Every consecutive group of $k$ positions advances the answer by 1, hence the offset inside the block is

$$\left\lceil\frac{d}{k}\right\rceil.$$

So we get the fast query formula

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

This is the core identity used by the function golomb_at.

Worked Example: $x=10$

The first boundaries are

$$S(1)=1,\qquad S(2)=3,\qquad S(3)=5,$$

and the weighted boundaries are

$$P(1)=1,\qquad P(2)=5,\qquad P(3)=11.$$

Since $P(2)=5\lt 10\le 11=P(3)$, the query lies in the $k=3$ block. Therefore

$$G(10)=S(2)+\left\lceil\frac{10-P(2)}{3}\right\rceil =3+\left\lceil\frac{5}{3}\right\rceil =5.$$

Indeed, the tenth term of the Golomb sequence is 5.

Why the Monotone Sweep Works

The problem only asks for $x=n^3$, and these queries arrive in strictly increasing order. Therefore the correct block index $k$ never moves backward. The code keeps the current $k$, $S(k)$, and $P(k)$ in a reusable state object and only advances them when the next cube crosses a new block boundary. Precomputation stops as soon as $P(k)$ reaches the largest required cube $(10^6-1)^3$.

How the Code Works

precompute_golomb_until(max_query_position) builds the table golomb[n]=G(n) using the recurrence above. While doing so, it also tracks products, which is exactly the weighted prefix sum $P(k)$. Thus it stops as soon as the largest needed query position is covered.

QueryState stores the current block index $k$, the current boundaries $S(k)$ and $P(k)$, and the previous boundaries $S(k-1)$ and $P(k-1)$. In golomb_at(x, golomb, st), a simple while loop advances to the unique block satisfying $P(k-1)\lt x\le P(k)$, then evaluates

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

The ceiling is implemented with pure integer arithmetic:

$$\left\lceil\frac{a}{b}\right\rceil=\left\lfloor\frac{a+b-1}{b}\right\rfloor.$$

Finally, sum_golomb_cubes(limit, golomb) loops over $n=1,2,\dots,\text{limit}-1$, computes $x=n^3$, queries $G(x)$, and accumulates the total.

Complexity Analysis

Let $K$ be the largest block index needed so that $P(K)\ge (10^6-1)^3$. Precomputing $G(1),\dots,G(K)$ costs $O(K)$ time and $O(K)$ memory. The cube sweep costs $O(10^6+K)$ amortized, because the pointer over $k$-blocks only moves forward. There is no per-query binary search and no attempt to materialize the sequence up to index $10^{18}$.

Footnotes and References

Problem page: https://projecteuler.net/problem=341
Golomb sequence: Wikipedia - Golomb sequence
Prefix sums: Wikipedia - Prefix sum

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());
    }
}