Project Euler Problem 341
The Golomb's self-describing sequence (G(n)) is the only nondecreasing sequence of natural numbers such that n appears e
Solution
Answer: 56098610614277014
Let
$G(n)=1+G\bigl(n-G(G(n-1))\bigr)$
be the standard recurrence for the Golomb sequence (valid for $n>1$), with $G(1)=1$.
We want to compute
$$\sum_{1 \le n < 10^6} G(n^3).$$
The challenge is that $n^3$ reaches nearly $10^{18}$, so directly generating the sequence up to that index is impossible.
1. Mathematical analysis
Step 1: Understanding Golomb's sequence
The defining property says:
the integer $k$ appears exactly $G(k)$ times.
So if $G(k)=4$, then the value $k$ occurs four times in the sequence.
Define the cumulative endpoint function
$$S(k)=\sum_{i=1}^{k} G(i).$$
Then:
- $k$ occupies positions
$$S(k-1)+1,\dots,S(k),$$
- therefore
$$G(n)=k \quad\Longleftrightarrow\quad S(k-1)<n\le S(k).$$
So computing $G(n^3)$ reduces to finding the first $k$ with
$$S(k)\ge n^3.$$
Step 2: Avoid generating up to $10^{18}$
A naïve approach would compute all $S(k)$ until $10^{18}$, but that would require about $10^{11}$ values — impossible.
The key observation is:
If $G(r)=L$, then the value $r$ occurs exactly $L$ times among the indices $k$.
That means there is a whole block of $L$ consecutive $k$-values for which
$$G(k)=r.$$
Inside that block:
- $S(k)$ increases by exactly $r$ each step,
- so the cumulative values form an arithmetic progression.
Hence we can process entire blocks at once rather than one $k$ at a time.
Since $G(r)$ only grows to around a few million, generating $G(r)$ up to the necessary size is completely feasible.
Step 3: Efficient block processing
Suppose we are in a block where:
- $G(k)=r$,
- block length is $L=G(r)$.
Then:
$$S(k)$$
advances by increments of $r$.
If we already know the cumulative sum before the block, say $S_0$, then the end of the block is
$$S_0 + Lr.$$
For a cube $c=n^3$ inside this range, the correct index is
$$k = k_0 + \left\lceil \frac{c-S_0}{r} \right\rceil,$$
where $k_0$ is the index before the block starts.
Since there are only $999{,}999$ cubes, we can sweep through them once.
2. Python implementation
def solve(limit=10**6):
"""
Compute:
sum(G(n^3)) for 1 <= n < limit
where G is Golomb's self-describing sequence.
"""
# Golomb values, 1-indexed
G = [0, 1] # G(1) = 1
total = 0
# Current cube pointer: n^3
n = 1
cube = 1
# We have already processed k = 1
# S(1) = 1
current_k = 1
current_S = 1
# Handle cubes <= S(1)
while n < limit and cube <= 1:
total += 1
n += 1
cube = n**3
r = 2
while n < limit:
# Generate G(r)
g_r = 1 + G[r - G[G[r - 1]]]
G.append(g_r)
# This block has length G(r)
block_length = g_r
# End of cumulative range for this block
block_end_S = current_S + block_length * r
# Process all cubes that fall in this block
while n < limit and cube <= block_end_S:
# Ceiling division:
# j = ceil((cube - current_S) / r)
j = (cube - current_S + r - 1) // r
# Corresponding k-value
total += current_k + j
n += 1
cube = n**3
# Move to next block
current_k += block_length
current_S = block_end_S
r += 1
return total
# Verification checks
print(solve(1000)) # should be 153506976
print(solve(10**6))
3. Code walkthrough
Building Golomb values
g_r = 1 + G[r - G[G[r - 1]]]
This is the classic Golomb recurrence:
$$G(r) = 1 + G(r - G(G(r-1))).$$
We only generate as many terms as needed.
Processing blocks
If $G(r)=L$, then the value $r$ repeats $L$ times among consecutive indices $k$.
Instead of iterating over every $k$, we process the whole block:
block_end_S = current_S + block_length * r
because each step increases $S(k)$ by exactly $r$.
Finding where a cube lands
If $c=n^3$ lies in the current block:
j = (cube - current_S + r - 1) // r
computes
$$\left\lceil \frac{c-S_0}{r} \right\rceil,$$
the offset inside the block.
Then:
total += current_k + j
adds the correct Golomb value.
4. Verification
The problem provides:
$$\sum_{1\le n<10^3} G(n^3)=153506976.$$
Running the algorithm for limit=1000 gives exactly:
$$153506976,$$
so the implementation matches the official check.
The final computation for $10^6$ completes quickly and produces a value of the expected magnitude ($\sim 10^{16}$).
Answer: 56098610614277014