Problem 506: Clock Sequence

View on Project Euler

Project Euler Problem 506 Solution

EulerSolve provides an optimized solution for Project Euler Problem 506, Clock Sequence, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary An infinite digit stream repeats $$1,2,3,4,3,2,1,2,3,4,3,2,\dots$$ The stream is continuous: after \(v_n\) is built, the next term starts exactly where the previous one stopped. For each positive integer \(n\), \(v_n\) is the decimal integer obtained by concatenating successive digits until their digit-sum is exactly \(n\). The first terms are $$v_1=1,\quad v_2=2,\quad v_3=3,\quad v_4=4,\quad v_5=32,\quad v_6=123,\dots$$ The goal is to compute $$S(N)=\sum_{n=1}^{N} v_n \pmod{123454321}$$ for an enormous value of \(N\). A direct simulation would require producing every term one after another, so the solution must exploit the periodicity of the digit stream and the arithmetic structure of the digit sums. Mathematical Approach The decisive observation is that the sequence does not need to be handled term by term forever. Once the terms are grouped by their residue class modulo \(15\), each group follows a simple affine recurrence. Step 1: Track the Starting Position of Each Term After the first \(m\) terms have been completed, the total digit-sum consumed from the stream is $$A_m=1+2+\cdots+m=\frac{m(m+1)}{2}.$$ One complete clock cycle uses the six digits \(1,2,3,4,3,2\), whose sum is $$1+2+3+4+3+2=15.$$ The partial sums inside one cycle are \(0,1,3,6,10,13,15\), and these are all distinct modulo \(15\)....

Detailed mathematical approach

Problem Summary

An infinite digit stream repeats

$$1,2,3,4,3,2,1,2,3,4,3,2,\dots$$

The stream is continuous: after \(v_n\) is built, the next term starts exactly where the previous one stopped. For each positive integer \(n\), \(v_n\) is the decimal integer obtained by concatenating successive digits until their digit-sum is exactly \(n\). The first terms are

$$v_1=1,\quad v_2=2,\quad v_3=3,\quad v_4=4,\quad v_5=32,\quad v_6=123,\dots$$

The goal is to compute

$$S(N)=\sum_{n=1}^{N} v_n \pmod{123454321}$$

for an enormous value of \(N\). A direct simulation would require producing every term one after another, so the solution must exploit the periodicity of the digit stream and the arithmetic structure of the digit sums.

Mathematical Approach

The decisive observation is that the sequence does not need to be handled term by term forever. Once the terms are grouped by their residue class modulo \(15\), each group follows a simple affine recurrence.

Step 1: Track the Starting Position of Each Term

After the first \(m\) terms have been completed, the total digit-sum consumed from the stream is

$$A_m=1+2+\cdots+m=\frac{m(m+1)}{2}.$$

One complete clock cycle uses the six digits \(1,2,3,4,3,2\), whose sum is

$$1+2+3+4+3+2=15.$$

The partial sums inside one cycle are \(0,1,3,6,10,13,15\), and these are all distinct modulo \(15\). Therefore the position inside the six-digit cycle is determined by the consumed digit-sum modulo \(15\).

The term \(v_n\) starts after the first \(n-1\) terms, so its starting position depends on \(A_{n-1}\). Now

$$A_{n+14}-A_{n-1}=\frac{(n+14)(n+15)-(n-1)n}{2}=15(n+7),$$

which is divisible by \(15\). Hence \(v_n\) and \(v_{n+15}\) always start at the same point of the repeating digit stream. This is why the natural classification is by \(n \bmod 15\).

Step 2: Terms in the Same Residue Class Differ by One Six-Digit Block

Fix a residue \(r \in \{1,\dots,15\}\) and write

$$n=r+15k.$$

The terms \(v_r,v_{r+15},v_{r+30},\dots\) all start at the same stream position. Their target digit-sums differ by \(15\), and the next complete cycle from any fixed position contains exactly six digits whose sum is \(15\). Therefore, once \(v_{r+15k}\) has been read, the term \(v_{r+15(k+1)}\) is obtained by appending one more six-digit rotation of the clock sequence.

If \(B=10^6\), appending six decimal digits multiplies the previous number by \(B\). So there exist constants \(b_r\) and \(a_r\), depending only on the residue class \(r\), such that

$$b_r=v_r,$$

$$v_{r+15(k+1)}=B\,v_{r+15k}+a_r.$$

The quantity \(a_r\) is the six-digit block appended in that residue class.

Step 3: Solve the Affine Recurrence

The recurrence

$$x_{k+1}=B\,x_k+a_r,\qquad x_0=b_r$$

has the closed form

$$x_k=b_r B^k+a_r\sum_{j=0}^{k-1}B^j.$$

Therefore

$$v_{r+15k}=b_r B^k+a_r\frac{B^k-1}{B-1}.$$

This formula is exact over the integers and is the core reason the huge input becomes manageable.

Step 4: Sum One Residue Class at a Time

For a fixed residue \(r\), the number of indices of the form \(r+15k\) that do not exceed \(N\) is

$$K_r=\left\lfloor\frac{N-r}{15}\right\rfloor+1 \qquad (N \ge r).$$

Define the two standard sums

$$P(K)=\sum_{k=0}^{K-1} B^k=\frac{B^K-1}{B-1},$$

$$Q(K)=\sum_{k=0}^{K-1}\frac{B^k-1}{B-1}=\frac{P(K)-K}{B-1}.$$

Then the total contribution of residue class \(r\) is

$$\sum_{k=0}^{K_r-1} v_{r+15k}=b_r P(K_r)+a_r Q(K_r).$$

Summing over all possible residues gives

$$S(N)=\sum_{r=1}^{\min(15,N)}\left(b_r P(K_r)+a_r Q(K_r)\right)\pmod{123454321}.$$

Step 5: Modular Form of the Formula

The implementation works modulo

$$M=123454321.$$

Since

$$\gcd(B-1,M)=\gcd(999999,123454321)=1,$$

the factor \(B-1\) has a modular inverse modulo \(M\). Thus the divisions in \(P(K)\) and \(Q(K)\) are carried out safely as multiplications by \((B-1)^{-1} \pmod{M}\), while \(B^K \pmod{M}\) is obtained by fast exponentiation.

Worked Example: The Residue Class \(r=5\)

The exact terms in this class begin as

$$v_5=32,\qquad v_{20}=32123432,\qquad v_{35}=32123432123432.$$

Therefore

$$a_5=v_{20}-10^6 v_5=32123432-32000000=123432.$$

The closed form becomes

$$v_{5+15k}=32 \cdot 10^{6k}+123432\frac{10^{6k}-1}{10^6-1}.$$

For \(k=2\), this gives

$$32 \cdot 10^{12}+123432(10^6+1)=32123432123432,$$

which matches the third exact term. As a small full-sequence checkpoint, the first eleven terms satisfy

$$S(11)=36120.$$

How the Code Works

The C++, Python, and Java implementations first generate the first \(45\) terms exactly. That is enough to obtain three terms from each residue class modulo \(15\): \(v_r\), \(v_{r+15}\), and \(v_{r+30}\).

From these values, the implementation extracts the initial value \(b_r=v_r\) and the appended six-digit block

$$a_r=v_{r+15}-10^6 v_r.$$

It then checks that the same block also satisfies

$$v_{r+30}=10^6 v_{r+15}+a_r,$$

so the affine law is verified before the large-\(N\) computation begins.

After that, the program uses fast modular exponentiation to evaluate \(B^{K_r} \pmod{M}\), computes the modular inverse of \(B-1\), forms \(P(K_r)\) and \(Q(K_r)\) modulo \(M\), and accumulates the \(15\) residue-class contributions. The exact prefix generation is only a fixed-size calibration step; the enormous target \(N\) is handled entirely by closed forms.

Complexity Analysis

Generating the first \(45\) exact terms is constant work. The main computation processes only \(15\) residue classes, and each class requires a constant number of modular operations plus one modular exponentiation \(B^{K_r}\). Therefore the total running time is

$$O(\log N),$$

and the memory usage is

$$O(1).$$

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=506
  2. Geometric series: Wikipedia — Geometric series
  3. Recurrence relation: Wikipedia — Recurrence relation
  4. Modular multiplicative inverse: Wikipedia — Modular multiplicative inverse
  5. Triangular number: Wikipedia — Triangular number

Problem 506 source code

C++

#include <array>
#include <cstdint>
#include <iostream>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

constexpr u64 kMod = 123'454'321ULL;
constexpr u64 kBlockBase = 1'000'000ULL;  // 10^6
constexpr std::array<int, 6> kDigits = {1, 2, 3, 4, 3, 2};

struct DigitStream {
    int pos = 0;

    int next() {
        const int d = kDigits[static_cast<std::size_t>(pos)];
        pos = (pos + 1) % static_cast<int>(kDigits.size());
        return d;
    }
};

u64 mod_pow(u64 base, u64 exp, const u64 mod) {
    u64 result = 1ULL % mod;
    u64 cur = base % mod;
    u64 e = exp;
    while (e > 0ULL) {
        if (e & 1ULL) {
            result = static_cast<u64>((static_cast<u128>(result) * cur) % mod);
        }
        cur = static_cast<u64>((static_cast<u128>(cur) * cur) % mod);
        e >>= 1ULL;
    }
    return result;
}

std::int64_t extended_gcd(const std::int64_t a, const std::int64_t b, std::int64_t& x,
                          std::int64_t& y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    std::int64_t x1 = 0;
    std::int64_t y1 = 0;
    const std::int64_t g = extended_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

u64 mod_inverse(const u64 a, const u64 mod) {
    std::int64_t x = 0;
    std::int64_t y = 0;
    const std::int64_t g = extended_gcd(static_cast<std::int64_t>(a),
                                        static_cast<std::int64_t>(mod), x, y);
    if (g != 1) {
        return 0ULL;
    }
    const std::int64_t m = static_cast<std::int64_t>(mod);
    std::int64_t res = x % m;
    if (res < 0) {
        res += m;
    }
    return static_cast<u64>(res);
}

std::vector<u128> direct_terms_exact(const int count) {
    DigitStream stream;
    std::vector<u128> out(static_cast<std::size_t>(count + 1), 0);
    for (int n = 1; n <= count; ++n) {
        int sum = 0;
        u128 value = 0;
        while (sum < n) {
            const int d = stream.next();
            sum += d;
            value = value * 10 + static_cast<u128>(d);
        }
        if (sum != n) {
            return {};
        }
        out[static_cast<std::size_t>(n)] = value;
    }
    return out;
}

u64 direct_sum_mod(const int n_max) {
    DigitStream stream;
    u64 total = 0ULL;
    for (int n = 1; n <= n_max; ++n) {
        int sum = 0;
        u64 value_mod = 0ULL;
        while (sum < n) {
            const int d = stream.next();
            sum += d;
            value_mod = (value_mod * 10ULL + static_cast<u64>(d)) % kMod;
        }
        if (sum != n) {
            return 0ULL;
        }
        total += value_mod;
        total %= kMod;
    }
    return total;
}

u64 direct_sum_exact_small(const int n_max) {
    DigitStream stream;
    u64 total = 0ULL;
    for (int n = 1; n <= n_max; ++n) {
        int sum = 0;
        u64 value = 0ULL;
        while (sum < n) {
            const int d = stream.next();
            sum += d;
            value = value * 10ULL + static_cast<u64>(d);
        }
        if (sum != n) {
            return 0ULL;
        }
        total += value;
    }
    return total;
}

u64 geometric_sum(const u64 ratio, const u64 terms, const u64 mod, const u64 inv_ratio_minus_one) {
    if (terms == 0ULL) {
        return 0ULL;
    }
    const u64 p = mod_pow(ratio, terms, mod);
    const u64 numerator = (p + mod - 1ULL) % mod;
    return static_cast<u64>((static_cast<u128>(numerator) * inv_ratio_minus_one) % mod);
}

u64 fast_sum_mod(const u64 n_max) {
    const std::vector<u128> terms = direct_terms_exact(45);
    if (terms.empty()) {
        return 0ULL;
    }

    std::array<u64, 15> base{};
    std::array<u64, 15> add{};
    for (int r = 1; r <= 15; ++r) {
        const u128 v0 = terms[static_cast<std::size_t>(r)];
        const u128 v1 = terms[static_cast<std::size_t>(r + 15)];
        const u128 v2 = terms[static_cast<std::size_t>(r + 30)];
        const u128 c = v1 - v0 * static_cast<u128>(kBlockBase);
        if (v2 != v1 * static_cast<u128>(kBlockBase) + c) {
            return 0ULL;
        }
        base[static_cast<std::size_t>(r - 1)] = static_cast<u64>(v0);
        add[static_cast<std::size_t>(r - 1)] = static_cast<u64>(c);
    }

    const u64 inv_b_minus_1 = mod_inverse(kBlockBase - 1ULL, kMod);
    if (inv_b_minus_1 == 0ULL) {
        return 0ULL;
    }

    u64 total = 0ULL;
    for (int r = 1; r <= 15; ++r) {
        if (n_max < static_cast<u64>(r)) {
            break;
        }
        const u64 count = (n_max - static_cast<u64>(r)) / 15ULL + 1ULL;
        const u64 g = geometric_sum(kBlockBase, count, kMod, inv_b_minus_1);
        const u64 count_mod = count % kMod;
        const u64 g_minus_count = (g + kMod - count_mod) % kMod;
        const u64 t = static_cast<u64>((static_cast<u128>(g_minus_count) * inv_b_minus_1) % kMod);

        const u64 part_base =
            static_cast<u64>((static_cast<u128>(base[static_cast<std::size_t>(r - 1)] % kMod) * g) %
                             kMod);
        const u64 part_add =
            static_cast<u64>((static_cast<u128>(add[static_cast<std::size_t>(r - 1)] % kMod) * t) %
                             kMod);
        total += part_base;
        total %= kMod;
        total += part_add;
        total %= kMod;
    }
    return total;
}

bool run_checkpoints() {
    if (direct_sum_exact_small(11) != 36'120ULL) {
        std::cerr << "Checkpoint failed: S(11)\n";
        return false;
    }
    if (direct_sum_mod(1'000) != 18'232'686ULL) {
        std::cerr << "Checkpoint failed: S(1000) mod M\n";
        return false;
    }
    if (fast_sum_mod(5'000ULL) != direct_sum_mod(5'000)) {
        std::cerr << "Checkpoint failed: fast/direct consistency at 5000\n";
        return false;
    }
    return true;
}

}  // namespace

int main() {
    if (!run_checkpoints()) {
        return 1;
    }

    constexpr u64 n = 100'000'000'000'000ULL;
    std::cout << fast_sum_mod(n) << '\n';
    return 0;
}

Python

def mod_pow(base, exp, mod):
    return pow(base, exp, mod)

def extended_gcd(a, b):
    if b == 0:
        return a, 1, 0
    g, x1, y1 = extended_gcd(b, a % b)
    x = y1
    y = x1 - (a // b) * y1
    return g, x, y

def mod_inverse(a, mod):
    g, x, y = extended_gcd(a, mod)
    if g != 1: return 0
    res = x % mod
    if res < 0: res += mod
    return res

def direct_terms_exact(count):
    kDigits = [1, 2, 3, 4, 3, 2]
    out = [0] * (count + 1)
    pos = 0
    for n in range(1, count + 1):
        s = 0
        val = 0
        while s < n:
            d = kDigits[pos]
            pos = (pos + 1) % 6
            s += d
            val = val * 10 + d
        if s != n: return []
        out[n] = val
    return out

def geometric_sum(ratio, terms, mod, inv_ratio_minus_one):
    if terms == 0: return 0
    p = mod_pow(ratio, terms, mod)
    numerator = (p + mod - 1) % mod
    return (numerator * inv_ratio_minus_one) % mod

def fast_sum_mod(n_max):
    kMod = 123454321
    kBlockBase = 1000000
    terms = direct_terms_exact(45)
    if not terms: return 0
    
    base = [0] * 15
    add = [0] * 15
    for r in range(1, 16):
        v0 = terms[r]
        v1 = terms[r + 15]
        v2 = terms[r + 30]
        c = v1 - v0 * kBlockBase
        if v2 != v1 * kBlockBase + c: return 0
        base[r - 1] = v0
        add[r - 1] = c
        
    inv_b_minus_1 = mod_inverse(kBlockBase - 1, kMod)
    if inv_b_minus_1 == 0: return 0
    
    total = 0
    for r in range(1, 16):
        if n_max < r: break
        count = (n_max - r) // 15 + 1
        g = geometric_sum(kBlockBase, count, kMod, inv_b_minus_1)
        count_mod = count % kMod
        g_minus_count = (g + kMod - count_mod) % kMod
        t = (g_minus_count * inv_b_minus_1) % kMod
        
        part_base = ((base[r - 1] % kMod) * g) % kMod
        part_add = ((add[r - 1] % kMod) * t) % kMod
        total = (total + part_base) % kMod
        total = (total + part_add) % kMod
        
    return total

def solve():
    return str(fast_sum_mod(100000000000000))

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

Java

import java.math.BigInteger;

public class Euler506 {
    private static final long kMod = 123454321L;
    private static final BigInteger kBlockBase = BigInteger.valueOf(1000000L);
    private static final int[] kDigits = { 1, 2, 3, 4, 3, 2 };

    private static long modPow(long base, long exp, long mod) {
        long result = 1L % mod;
        long cur = base % mod;
        long e = exp;
        while (e > 0L) {
            if ((e & 1L) != 0L) {
                result = (result * cur) % mod;
            }
            cur = (cur * cur) % mod;
            e >>= 1L;
        }
        return result;
    }

    private static long[] extendedGcd(long a, long b) {
        if (b == 0)
            return new long[] { a, 1, 0 };
        long[] res = extendedGcd(b, a % b);
        long g = res[0];
        long x1 = res[1];
        long y1 = res[2];
        long x = y1;
        long y = x1 - (a / b) * y1;
        return new long[] { g, x, y };
    }

    private static long modInverse(long a, long mod) {
        long[] res = extendedGcd(a, mod);
        if (res[0] != 1)
            return 0L;
        long x = res[1] % mod;
        if (x < 0)
            x += mod;
        return x;
    }

    private static BigInteger[] directTermsExact(int count) {
        BigInteger[] out = new BigInteger[count + 1];
        int pos = 0;
        for (int n = 1; n <= count; ++n) {
            int s = 0;
            BigInteger val = BigInteger.ZERO;
            while (s < n) {
                int d = kDigits[pos];
                pos = (pos + 1) % 6;
                s += d;
                val = val.multiply(BigInteger.TEN).add(BigInteger.valueOf(d));
            }
            if (s != n)
                return new BigInteger[0];
            out[n] = val;
        }
        return out;
    }

    private static long geometricSum(long ratio, long terms, long mod, long invRatioMinusOne) {
        if (terms == 0)
            return 0L;
        long p = modPow(ratio, terms, mod);
        long numerator = (p + mod - 1L) % mod;
        return (numerator * invRatioMinusOne) % mod;
    }

    public static void main(String[] args) {
        long nMax = 100000000000000L;
        BigInteger[] terms = directTermsExact(45);
        if (terms.length == 0)
            return;

        long[] base = new long[15];
        long[] add = new long[15];

        for (int r = 1; r <= 15; ++r) {
            BigInteger v0 = terms[r];
            BigInteger v1 = terms[r + 15];
            BigInteger v2 = terms[r + 30];
            BigInteger c = v1.subtract(v0.multiply(kBlockBase));
            if (!v2.equals(v1.multiply(kBlockBase).add(c)))
                return;
            base[r - 1] = v0.longValue();
            add[r - 1] = c.longValue();
        }

        long invBMinus1 = modInverse(kBlockBase.longValue() - 1L, kMod);
        if (invBMinus1 == 0L)
            return;

        long total = 0L;
        for (int r = 1; r <= 15; ++r) {
            if (nMax < r)
                break;
            long count = (nMax - r) / 15L + 1L;
            long g = geometricSum(kBlockBase.longValue(), count, kMod, invBMinus1);
            long countMod = count % kMod;
            long gMinusCount = (g + kMod - countMod) % kMod;
            long t = (gMinusCount * invBMinus1) % kMod;

            long partBase = ((base[r - 1] % kMod) * g) % kMod;
            long partAdd = ((add[r - 1] % kMod) * t) % kMod;

            total = (total + partBase) % kMod;
            total = (total + partAdd) % kMod;
        }
        System.out.println(total);
    }
}