Problem 551: Sum of Digits Sequence
View on Project EulerProject Euler Problem 551 Solution
EulerSolve provides an optimized solution for Project Euler Problem 551, Sum of Digits Sequence, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary The sequence is defined by $$a_1=1,\qquad a_{n+1}=a_n+\sigma(a_n),$$ where \(\sigma(x)\) is the decimal digit sum of \(x\). The target index is so large that simulating every update one by one is infeasible, so the implementations accelerate the process by grouping many ordinary updates into exact carry-to-carry transitions. Mathematical Approach The key observation is that decimal digit sums behave well across powers of \(10\). By isolating the lowest six digits, the recurrence splits into a rapidly changing low part and a slowly changing high part. That split makes it possible to precompute exact jumps to the next carry and then compose those jumps over long intervals. Step 1: Split Each Term into a High Part and a Six-Digit Low Part Set $$B=10^6,\qquad a_n=Bq_n+r_n,\qquad 0\le r_n \lt B.$$ Let \(\sigma_6(r)\) denote the digit sum of the six-digit decimal block of \(r\), allowing leading zeros....
Detailed mathematical approach
Problem Summary
The sequence is defined by
$$a_1=1,\qquad a_{n+1}=a_n+\sigma(a_n),$$
where \(\sigma(x)\) is the decimal digit sum of \(x\). The target index is so large that simulating every update one by one is infeasible, so the implementations accelerate the process by grouping many ordinary updates into exact carry-to-carry transitions.
Mathematical Approach
The key observation is that decimal digit sums behave well across powers of \(10\). By isolating the lowest six digits, the recurrence splits into a rapidly changing low part and a slowly changing high part. That split makes it possible to precompute exact jumps to the next carry and then compose those jumps over long intervals.
Step 1: Split Each Term into a High Part and a Six-Digit Low Part
Set
$$B=10^6,\qquad a_n=Bq_n+r_n,\qquad 0\le r_n \lt B.$$
Let \(\sigma_6(r)\) denote the digit sum of the six-digit decimal block of \(r\), allowing leading zeros. Because \(B\) is a power of \(10\), the decimal digits of \(q_n\) and \(r_n\) do not interfere, so
$$\sigma(a_n)=\sigma(q_n)+\sigma_6(r_n).$$
Therefore one elementary update becomes
$$a_{n+1}=Bq_n+r_n+\sigma(q_n)+\sigma_6(r_n).$$
If the low part stays below \(B\), then
$$q_{n+1}=q_n,\qquad r_{n+1}=r_n+\sigma_6(r_n)+\sigma(q_n).$$
If the low part reaches or exceeds \(B\), one carry occurs:
$$q_{n+1}=q_n+1,\qquad r_{n+1}=r_n+\sigma_6(r_n)+\sigma(q_n)-B.$$
For the target computation, the high part never needs more than \(12\) decimal digits, so \(\sigma(q_n)\le 108\). Since \(\sigma_6(r_n)\le 54\), every elementary update increases the low part by at most \(162\). In particular, a single update can create at most one carry.
Step 2: Jump Directly to the Next Carry When the High Part Is Fixed
While \(q_n\) is fixed, its digit sum is a constant. For a fixed value \(c\), define
$$F_c(x)=x+\sigma_6(x)+c.$$
Starting from a low part \(x\), the next carry happens after
$$\tau_c(x)=\min\{t\ge 1:F_c^{(t)}(x)\ge B\}$$
elementary updates, and the low part immediately after that carry is
$$\rho_c(x)=F_c^{(\tau_c(x))}(x)-B.$$
These quantities are exact: \(\tau_c(x)\) counts how many ordinary recurrence steps are consumed before the boundary \(B\) is crossed, and \(\rho_c(x)\) is the true remainder after subtracting \(B\).
Step 3: After a Carry, the Low Part Returns to a Tiny State Space
At the moment of crossing we have
$$B\le F_c^{(\tau_c(x))}(x)\le B+54+108,$$
so the post-carry remainder always satisfies
$$0\le \rho_c(x)\le 162.$$
This is the main compression step. Although the low part can wander through all values below \(10^6\) before the next carry, the state immediately after a carry lies in the tiny set
$$S=\{0,1,2,\dots,162\}.$$
So a carry event can be viewed as a macro-step that starts from some remainder in \(S\), spends a known number of elementary updates, and returns to another remainder in \(S\).
Step 4: Represent One Carry Event as a Transition Depending Only on the Current High Part
For each current high part \(h\), let \(c=\sigma(h)\). One carry event then induces a transition \(M_h\) on the small remainder set \(S\): for each \(x\in S\), the transition returns the new remainder \(\rho_c(x)\) and the exact cost \(\tau_c(x)\).
If we want to process several carries in a row, starting from high part \(h\), we must compose the transitions
$$M_h,\ M_{h+1},\ M_{h+2},\ \dots.$$
For adjacent intervals of high-part values, composition is exact:
$$\mathcal{T}_{[u,w)}=\mathcal{T}_{[v,w)}\circ \mathcal{T}_{[u,v)}\qquad (u\le v\le w).$$
This means that long runs of carry events can be handled by composing reusable transition blocks instead of replaying each event separately.
Step 5: Build Decimal Blocks of Consecutive High Parts
Consider an aligned interval of length \(10^d\), such as all high-part values from \(2300\) to \(2399\). Every number in that interval shares the same prefix and only the last \(d\) digits vary. Because digit sums split over decimal concatenation, the digit sum of each high part is
$$\sigma(\text{prefix}\cdot 10^d+\text{suffix})=\sigma(\text{prefix})+\sigma(\text{suffix}).$$
Therefore the full transition of such a block depends only on two pieces of information: the block length \(10^d\) and the digit sum of the fixed prefix. The implementations precompute all these aligned decimal blocks up to \(d=12\). An arbitrary interval \([u,v)\) is then decomposed greedily into aligned blocks, so a long interval can be evaluated with only a small number of transition compositions.
Step 6: Use Doubling and Binary Search to Maximize the Number of Carry Events
Suppose \(R\) elementary updates are still available. Let \(C(m)\) be the number of elementary updates consumed by the next \(m\) carry events from the current state. Since each carry event costs at least one ordinary update, \(C(m)\) is strictly increasing.
So the algorithm finds
$$m^*=\max\{m\ge 0:C(m)\le R\}$$
by first doubling \(m\) until the budget is exceeded and then applying binary search. After those \(m^*\) carry events are executed in one shot, the remaining budget is too small to reach the next carry. The last few elementary updates are therefore performed directly with the high part fixed.
Worked Example: One Carry Event and a Short Tail
Take
$$a=17\cdot 10^6+999995.$$
Then \(\sigma(17)=8\) and \(\sigma_6(999995)=50\). One elementary update gives
$$999995+50+8=1000053,$$
so the next term is
$$a'=18\cdot 10^6+53.$$
This single carry event has consumed one ordinary recurrence step and moved the low part back into the small set \(S\). Now the high part is \(18\), so its digit sum is \(9\). If two more elementary updates are taken before the next carry, the low part evolves as
$$53\to 70\to 86,$$
because \(53+\sigma_6(53)+9=53+8+9=70\) and \(70+\sigma_6(70)+9=70+7+9=86\). This is exactly the kind of short leftover tail that the algorithm finishes directly after the carry-jump phase.
How the Code Works
The C++, Python, and Java implementations begin by precomputing the digit sum of every six-digit value from \(0\) to \(999999\). This makes the low-part update \(x\mapsto x+\sigma_6(x)+c\) constant-time for every possible state.
Next, for each possible digit sum \(c\) of the high part, the implementation scans the full low-part range backward and determines the exact number of ordinary updates needed to hit the next carry, together with the post-carry remainder. Only the carry-to-carry states \(0\) through \(162\) need to be retained for later composition, because every carry lands inside that range.
After that, the implementation builds transition blocks for aligned decimal intervals of lengths \(1,10,100,\dots,10^{12}\). Each block stores, for every small remainder, both the new remainder and the number of ordinary updates consumed across that whole interval of consecutive high-part values.
During the actual solve phase, the current value starts at \(1\). The implementation repeatedly asks how many carry events can be afforded under the remaining step budget, answers that question with doubling and binary search, applies the resulting interval transition, and finally performs the tiny leftover suffix directly. The method is exact throughout; nothing is approximated.
Complexity Analysis
Let \(B=10^6\), let \(C=108\) be the maximum possible digit sum of the high part used by the implementations, let \(R=162\) be the post-carry remainder bound, and let \(D=12\) be the maximum decimal block depth. Precomputing all six-digit digit sums costs \(O(B)\) time. Building the next-carry information for every \(c\in[0,C]\) costs \(O(BC)\) time and dominates the preprocessing.
The aligned block tables operate only on the tiny state space of size \(R+1=163\), so their cost is much smaller than the carry-table construction. For the single final query, doubling and binary search use a logarithmic number of interval evaluations, and each interval evaluation touches only a small number of decimal blocks. Memory usage is dominated by the carry tables and is \(O(BC)\).
Footnotes and References
- Problem page: https://projecteuler.net/problem=551
- Digit sum: Wikipedia — Digit sum
- Decimal numeral system: Wikipedia — Decimal
- Function composition and iteration: Wikipedia — Iterated function
- Binary search: Wikipedia — Binary search algorithm
Problem 551 source code
C++
#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <vector>
namespace {
using u8 = std::uint8_t;
using u16 = std::uint16_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;
static constexpr int DMAX = 12; // H fits in <= 12 decimal digits for this problem.
static constexpr int CMAX = 9 * DMAX; // max possible digit sum of H.
static constexpr u32 P = 1'000'000U; // 10^6, keeps carry per step <= 1.
static constexpr int DS6_MAX = 54; // max digit sum for 6 digits.
static constexpr int RMAX = DS6_MAX + CMAX; // max remainder after crossing P.
static int digit_sum_u64(u64 x) {
int s = 0;
while (x) {
s += static_cast<int>(x % 10);
x /= 10;
}
return s;
}
struct Tables {
std::vector<u8> ds6; // digit sum for 0..P-1
// For fixed c in [0..CMAX] and starting x in [0..RMAX], jump to the next carry:
// after steps[c][x] single-step updates, the low part overflows and becomes rem[c][x],
// and the high part increases by 1.
std::array<std::array<u16, RMAX + 1>, CMAX + 1> rem{};
std::array<std::array<u32, RMAX + 1>, CMAX + 1> steps{};
Tables() { build(); }
void build() {
ds6.assign(P, 0);
for (u32 i = 1; i < P; ++i) {
ds6[i] = static_cast<u8>(ds6[i / 10] + (i % 10));
}
std::vector<u32> st(P);
std::vector<u16> rm(P);
for (int c = 0; c <= CMAX; ++c) {
for (int i = static_cast<int>(P) - 1; i >= 0; --i) {
const u32 x = static_cast<u32>(i);
const u32 y = x + static_cast<u32>(ds6[x]) + static_cast<u32>(c);
if (y >= P) {
st[x] = 1;
rm[x] = static_cast<u16>(y - P);
} else {
st[x] = st[y] + 1;
rm[x] = rm[y];
}
}
for (int x = 0; x <= RMAX; ++x) {
steps[c][x] = st[static_cast<u32>(x)];
rem[c][x] = rm[static_cast<u32>(x)];
}
}
}
};
struct Trans {
std::array<u16, RMAX + 1> next{};
std::array<u64, RMAX + 1> cost{}; // number of single-step updates
};
static Trans identity_trans() {
Trans t;
for (int r = 0; r <= RMAX; ++r) {
t.next[r] = static_cast<u16>(r);
t.cost[r] = 0;
}
return t;
}
static Trans compose(const Trans& a, const Trans& b) {
// Apply a then b.
Trans out;
for (int r = 0; r <= RMAX; ++r) {
const u16 mid = a.next[r];
out.next[r] = b.next[mid];
out.cost[r] = a.cost[r] + b.cost[mid];
}
return out;
}
struct BlockBuilder {
const Tables& tab;
std::array<u64, DMAX + 1> pow10{};
std::array<std::array<Trans, CMAX + 1>, DMAX + 1> block{};
explicit BlockBuilder(const Tables& t) : tab(t) { build(); }
Trans single(int digit_sum_h) const {
Trans tr;
for (int r = 0; r <= RMAX; ++r) {
tr.next[r] = tab.rem[digit_sum_h][r];
tr.cost[r] = tab.steps[digit_sum_h][r];
}
return tr;
}
void build() {
pow10[0] = 1;
for (int d = 1; d <= DMAX; ++d) pow10[d] = pow10[d - 1] * 10ULL;
for (int ps = 0; ps <= CMAX; ++ps) block[0][ps] = single(ps);
for (int d = 1; d <= DMAX; ++d) {
for (int ps = 0; ps + 9 * d <= CMAX; ++ps) {
Trans tr = block[d - 1][ps];
for (int dig = 1; dig <= 9; ++dig) tr = compose(tr, block[d - 1][ps + dig]);
block[d][ps] = tr;
}
}
}
Trans range_trans(u64 l, u64 r) const {
// Compose macro-steps for H in [l, r), in increasing order.
Trans res = identity_trans();
u64 cur = l;
while (cur < r) {
int best_d = 0;
for (int d = DMAX; d >= 0; --d) {
const u64 p = pow10[d];
if (p == 0) continue;
if (cur % p == 0 && cur + p <= r) {
best_d = d;
break;
}
}
const u64 p = pow10[best_d];
const u64 hi = cur / p;
const int ps = digit_sum_u64(hi);
res = compose(res, block[best_d][ps]);
cur += p;
}
return res;
}
};
static u64 solve_a(u64 n) {
// a_0=1, a_1=1, and for n>=1: a_{n+1}=a_n+sumdigits(a_n).
if (n == 0) return 1;
if (n == 1) return 1;
static const Tables tables;
static const BlockBuilder builder(tables);
u64 remaining_steps = n - 1; // number of single-step updates from a_1.
u64 H = 0;
u32 low = 1;
// Find maximum number of low-part carry events we can perform within remaining_steps.
// Each carry event processes steps[digit_sum(H)][low] single steps and increases H by 1.
auto cost_after_carries = [&](u64 carry_cnt, u16 start_low, u16& out_low) -> u64 {
const Trans tr = builder.range_trans(H, H + carry_cnt);
out_low = tr.next[start_low];
return tr.cost[start_low];
};
u64 lo = 0, hi = 1;
while (true) {
u16 dummy = 0;
const u64 cost = cost_after_carries(hi, static_cast<u16>(low), dummy);
if (cost > remaining_steps) break;
lo = hi;
hi <<= 1ULL;
}
while (lo + 1 < hi) {
const u64 mid = lo + (hi - lo) / 2;
u16 dummy = 0;
const u64 cost = cost_after_carries(mid, static_cast<u16>(low), dummy);
if (cost <= remaining_steps) lo = mid;
else hi = mid;
}
u16 low_after = 0;
const u64 used = cost_after_carries(lo, static_cast<u16>(low), low_after);
H += lo;
remaining_steps -= used;
low = static_cast<u32>(low_after);
// Finish remaining steps without crossing the next carry boundary.
const int c = digit_sum_u64(H);
for (u64 i = 0; i < remaining_steps; ++i) {
low += static_cast<u32>(tables.ds6[low]) + static_cast<u32>(c);
// By construction we should not overflow.
assert(low < P);
}
return H * static_cast<u64>(P) + static_cast<u64>(low);
}
} // namespace
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
// Checkpoint from the problem statement.
assert(solve_a(1'000'000ULL) == 31054319ULL);
std::cout << solve_a(1'000'000'000'000'000ULL) << '\n';
return 0;
}
Python
def solve():
N = 10**15
DMAX = 12
CMAX = 9 * DMAX
P = 1_000_000
DS6_MAX = 54
RMAX = DS6_MAX + CMAX
def digit_sum(x):
s = 0
while x:
s += x % 10
x //= 10
return s
# Build ds6: digit sum for 0..P-1
ds6 = [0] * P
for i in range(1, P):
ds6[i] = ds6[i // 10] + i % 10
# Build tables: for each carry c, steps/rem to next overflow
steps_tab = [[0]*(RMAX+1) for _ in range(CMAX+1)]
rem_tab = [[0]*(RMAX+1) for _ in range(CMAX+1)]
for c in range(CMAX + 1):
st = [0] * P
rm = [0] * P
for i in range(P - 1, -1, -1):
y = i + ds6[i] + c
if y >= P:
st[i] = 1
rm[i] = y - P
else:
st[i] = st[y] + 1
rm[i] = rm[y]
for x in range(RMAX + 1):
steps_tab[c][x] = st[x]
rem_tab[c][x] = rm[x]
# Transpose compose logic
def make_single(c):
nxt = list(rem_tab[c][:RMAX+1])
cost = list(steps_tab[c][:RMAX+1])
return (nxt, cost)
def compose(a, b):
nxt_a, cost_a = a
nxt_b, cost_b = b
nxt = [0]*(RMAX+1)
cost = [0]*(RMAX+1)
for r in range(RMAX+1):
mid = nxt_a[r]
nxt[r] = nxt_b[mid]
cost[r] = cost_a[r] + cost_b[mid]
return (nxt, cost)
def identity():
return (list(range(RMAX+1)), [0]*(RMAX+1))
pow10 = [1]
for d in range(1, DMAX+1):
pow10.append(pow10[-1] * 10)
block = [[None]*(CMAX+1) for _ in range(DMAX+1)]
for ps in range(CMAX+1):
block[0][ps] = make_single(ps)
for d in range(1, DMAX+1):
for ps in range(CMAX - 9*d + 1):
tr = block[d-1][ps]
for dig in range(1, 10):
tr = compose(tr, block[d-1][ps+dig])
block[d][ps] = tr
def range_trans(l, r):
res = identity()
cur = l
while cur < r:
best_d = 0
for d in range(DMAX, -1, -1):
p = pow10[d]
if cur % p == 0 and cur + p <= r:
best_d = d
break
p = pow10[best_d]
hi = cur // p
ps = digit_sum(hi)
res = compose(res, block[best_d][ps])
cur += p
return res
if N <= 1: return str(1)
remaining = N - 1
H = 0
low = 1
lo_c, hi_c = 0, 1
while True:
tr = range_trans(H, H + hi_c)
cost = tr[1][low]
if cost > remaining: break
lo_c = hi_c
hi_c <<= 1
while lo_c + 1 < hi_c:
mid = lo_c + (hi_c - lo_c) // 2
tr = range_trans(H, H + mid)
cost = tr[1][low]
if cost <= remaining: lo_c = mid
else: hi_c = mid
tr = range_trans(H, H + lo_c)
used = tr[1][low]
low_after = tr[0][low]
H += lo_c
remaining -= used
low = low_after
c = digit_sum(H)
for _ in range(remaining):
low += ds6[low] + c
return str(H * P + low)
if __name__ == '__main__':
print(solve())
Java
public class Euler551 {
static final int DMAX = 12;
static final int CMAX = 9 * DMAX;
static final int P = 1000000;
static final int DS6_MAX = 54;
static final int RMAX = DS6_MAX + CMAX;
static int digitSum(long x) {
int s = 0;
while (x > 0) {
s += (int) (x % 10);
x /= 10;
}
return s;
}
static class Tables {
byte[] ds6 = new byte[P];
short[][] rem = new short[CMAX + 1][RMAX + 1];
long[][] steps = new long[CMAX + 1][RMAX + 1];
Tables() {
for (int i = 1; i < P; i++) {
ds6[i] = (byte) (ds6[i / 10] + (i % 10));
}
int[] st = new int[P];
short[] rm = new short[P];
int[] dsPlusI = new int[P];
for (int i = 0; i < P; i++)
dsPlusI[i] = i + ds6[i];
for (int c = 0; c <= CMAX; c++) {
for (int i = P - 1; i >= 0; i--) {
int y = dsPlusI[i] + c;
if (y >= P) {
st[i] = 1;
rm[i] = (short) (y - P);
} else {
st[i] = st[y] + 1;
rm[i] = rm[y];
}
}
for (int x = 0; x <= RMAX; x++) {
steps[c][x] = st[x];
rem[c][x] = rm[x];
}
}
}
}
static class Trans {
short[] next = new short[RMAX + 1];
long[] cost = new long[RMAX + 1];
}
static Trans identityTrans() {
Trans t = new Trans();
for (int r = 0; r <= RMAX; r++) {
t.next[r] = (short) r;
t.cost[r] = 0;
}
return t;
}
static Trans compose(Trans a, Trans b) {
Trans out = new Trans();
for (int r = 0; r <= RMAX; r++) {
int mid = a.next[r] & 0xFFFF;
out.next[r] = b.next[mid];
out.cost[r] = a.cost[r] + b.cost[mid];
}
return out;
}
static class BlockBuilder {
Tables tab;
long[] pow10 = new long[DMAX + 1];
Trans[][] block = new Trans[DMAX + 1][CMAX + 1];
BlockBuilder(Tables tab) {
this.tab = tab;
pow10[0] = 1;
for (int d = 1; d <= DMAX; d++)
pow10[d] = pow10[d - 1] * 10L;
for (int ps = 0; ps <= CMAX; ps++) {
Trans t = new Trans();
for (int r = 0; r <= RMAX; r++) {
t.next[r] = tab.rem[ps][r];
t.cost[r] = tab.steps[ps][r];
}
block[0][ps] = t;
}
for (int d = 1; d <= DMAX; d++) {
for (int ps = 0; ps + 9 * d <= CMAX; ps++) {
Trans tr = block[d - 1][ps];
for (int dig = 1; dig <= 9; dig++) {
tr = compose(tr, block[d - 1][ps + dig]);
}
block[d][ps] = tr;
}
}
}
Trans rangeTrans(long l, long r) {
Trans res = identityTrans();
long cur = l;
while (cur < r) {
int bestD = 0;
for (int d = DMAX; d >= 0; d--) {
long p = pow10[d];
if (p != 0 && cur % p == 0 && cur + p <= r) {
bestD = d;
break;
}
}
long p = pow10[bestD];
long hi = cur / p;
int ps = digitSum(hi);
res = compose(res, block[bestD][ps]);
cur += p;
}
return res;
}
}
static long solveA(long n) {
if (n <= 1)
return 1;
Tables tables = new Tables();
BlockBuilder builder = new BlockBuilder(tables);
long remainingSteps = n - 1;
long H = 0;
int low = 1;
long lo = 0, hi = 1;
while (true) {
Trans tr = builder.rangeTrans(H, H + hi);
long cost = tr.cost[low];
if (cost > remainingSteps)
break;
lo = hi;
hi <<= 1;
}
while (lo + 1 < hi) {
long mid = lo + (hi - lo) / 2;
Trans tr = builder.rangeTrans(H, H + mid);
long cost = tr.cost[low];
if (cost <= remainingSteps) {
lo = mid;
} else {
hi = mid;
}
}
Trans finalTr = builder.rangeTrans(H, H + lo);
long used = finalTr.cost[low];
int lowAfter = finalTr.next[low] & 0xFFFF;
H += lo;
remainingSteps -= used;
low = lowAfter;
int c = digitSum(H);
for (long i = 0; i < remainingSteps; i++) {
low += (tables.ds6[low] & 0xFF) + c;
}
return H * P + low;
}
public static String solve() {
return Long.toString(solveA(1000000000000000L));
}
public static void main(String[] args) {
System.out.println(solve());
}
}