Problem 284: Steady Squares

Project Euler Problem 284 Solution

EulerSolve provides an optimized solution for Project Euler Problem 284, Steady Squares, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary In base $14$, a positive integer $x$ is a steady square of length $n$ if its square ends in the same $n$ base-$14$ digits, i.e. $$x^2 \equiv x \pmod{14^n}.$$ The code sums the base-$14$ digit sums of all positive steady squares with at most $n_{\max}$ digits and prints the final total in base $14$. The Project Euler final value is intentionally omitted here. Mathematical Approach 1) Idempotents modulo $14^n$. The congruence $x^2 \equiv x \pmod{14^n}$ is equivalent to $$x(x-1)\equiv 0 \pmod{2^n7^n}.$$ Because $\gcd(x,x-1)=1$, for the factor $2^n$ one of $x$ and $x-1$ must be divisible by $2^n$, so $$x \equiv 0 \text{ or } 1 \pmod{2^n}.$$ The same argument modulo $7^n$ gives $$x \equiv 0 \text{ or } 1 \pmod{7^n}.$$ By the Chinese Remainder Theorem there are exactly four idempotent residue classes modulo $14^n$: $$ (0,0),\quad (1,1),\quad (0,1),\quad (1,0) $$ with respect to $(\bmod\,2^n,\bmod\,7^n)$. For $n=1$ these are exactly $$0,\ 1,\ 7,\ 8 \pmod{14}.$$ The branches starting at $0$ and $1$ are trivial; the two nontrivial infinite branches start at $7$ and $8$. That is why the implementation only tracks those two sequences, while adding the one-digit steady square $1$ separately. 2) One-step lifting in base $14$....

Detailed mathematical approach

Problem Summary

In base $14$ the digits are $0,\dots,9$ and $\mathtt a,\mathtt b,\mathtt c,\mathtt d$ for $10,\dots,13$; base-14 numerals are set in typewriter type, such as $\mathtt{c37}$, and decimal numbers in ordinary type. A number with $n$ digits, the first of which is not $0$, is a steady square if its square ends in the same $n$ digits. In base $10$, $376$ is one, because $376^2=141376$. In base $14$, $\mathtt{c37}$ is one, because $\mathtt{c37}^2=\mathtt{aa0c37}$, and its digit sum is $\mathtt c+3+7=\mathtt{18}$, that is $22$. For $1\le n\le9$, the digit sums of all $n$-digit steady squares in base $14$ add up to $\mathtt{2d8}=582$. We must find the same total for $1\le n\le10000$ and write it in base $14$, with lower-case letters for the digits above $9$. ^[1] Steady squares are better known as automorphic numbers. ^[2]

The whole solution in one sentence: apart from $1$, every steady square is an ending of one of two infinite digit strings that end in $7$ and in $8$, each new digit of these strings is found from a simple congruence modulo $14$, and so the program builds both strings one digit at a time up to length $10000$, adding a string's running digit sum whenever its new front digit is not $0$.

Mathematical Approach

We turn the digit condition into a congruence and then build its solutions one digit at a time.

1. Steady squares solve a congruence

Digit positions are numbered from the right, starting with $0$: position $i$ is worth $14^i$. A number $x$ with $n$ digits satisfies $14^{n-1}\le x\lt14^n$, and its last $n$ digits form the remainder of $x$ modulo $14^n$, which is $x$ itself, while the last $n$ digits of $x^2$ form the remainder of $x^2$. So $x$ is a steady square exactly when $x^2\equiv x\pmod{14^n}$, that is, when

$$14^n\mid x^2-x=x(x-1).$$

For one digit we can simply try everything. The squares of $0,1,2,\dots,13$ leave the remainders

$$0,1,4,9,2,11,8,7,8,11,2,9,4,1$$

after division by $14$, so only for $0$, $1$, $7$ and $8$ is the remainder of the square the number itself. The one-digit steady squares are $1$, $7$ and $8$, and the total for $n=1$ is $1+7+8=16$; whether $0$ is counted makes no difference, as its digit sum is $0$. In base $10$ the same test gives $1$, $5$, $6$ for one digit, then $25$, $76$, $376$, $625$, $9376$, …, the automorphic numbers of OEIS A003226. ^[3]

2. Four solutions for every length

Since $x-(x-1)=1$, the numbers $x$ and $x-1$ have no common prime factor. So if $2^n$ divides $x(x-1)$, all the factors $2$ sit in one of them, and $2^n$ divides $x$ or $x-1$; the same holds for $7^n$. Hence $14^n=2^n7^n$ divides $x(x-1)$ exactly when both remainders are $0$ or $1$:

$$\begin{aligned}x\bmod 2^n&\in\{0,1\},\\ x\bmod 7^n&\in\{0,1\}.\end{aligned}$$

Because $2^n$ and $7^n$ are coprime, the Chinese remainder theorem says that each of the four combinations is met by exactly one $x$ with $0\le x\lt14^n$. ^[4] Two of them are $0$ and $1$. The other two we call

$a_n$, with $a_n\equiv0\pmod{2^n}$ and $a_n\equiv1\pmod{7^n}$, and
$b_n$, with $b_n\equiv1\pmod{2^n}$ and $b_n\equiv0\pmod{7^n}$.

For $n=1$ these are $a_1=8$ and $b_1=7$. In the same way a base with $k$ different prime factors has $2^k$ solutions for every length. ^[2] The decimal example of the statement is of the first kind, with $5$ in place of $7$: $376=8\cdot47$ is divisible by $2^3$, and $376=3\cdot125+1$ leaves the remainder $1$ modulo $5^3$. For $n\ge2$ the numbers $0$ and $1$ have fewer than $n$ digits, so the $n$-digit steady squares are those of $a_n$ and $b_n$ whose digit in position $n-1$ is not $0$, together with $1$ when $n=1$.

3. Each family is one infinite string of digits

The remainder of $a_{n+1}$ modulo $14^n$ is again a solution modulo $14^n$, and it is still $\equiv0\pmod{2^n}$ and $\equiv1\pmod{7^n}$. By the uniqueness in Section 2 it equals $a_n$. So $a_{n+1}=a_n+t\cdot14^n$ for a digit $t$: passing from $a_n$ to $a_{n+1}$ only writes one new digit in front, and the same holds for $b_n$. Each family is therefore the list of endings of one infinite string of digits, read from the right: ^[2]

$$\begin{aligned}\cdots\mathtt{6aa633d1a8}&\quad(a_n),\\ \cdots\mathtt{7337aa0c37}&\quad(b_n).\end{aligned}$$

We write $\alpha_i$ and $\beta_i$ for the digits in position $i$ of these strings, so that $a_n=\sum_{i\lt n}\alpha_i14^i$ and $b_n=\sum_{i\lt n}\beta_i14^i$.

4. One new digit from a linear congruence

Suppose that $0\le x\lt M=14^k$ with $k\ge1$ solves $x^2\equiv x\pmod M$, and write $x^2-x=fM$ with an integer $f$. Every number below $14M$ that ends in the same $k$ digits as $x$ is $x'=x+tM$ with a digit $t$ from $0$ to $13$, and multiplying out gives

$$\begin{aligned}&x'^2-x'\\&\quad=fM+(2x-1)tM+t^2M^2\\&\quad=M\bigl(f+(2x-1)t+t^2M\bigr).\end{aligned}$$

This is divisible by $14M$ exactly when the bracket is divisible by $14$. The term $t^2M$ always is, because $14\mid M$, so the condition is the linear congruence

$$f+(2x-1)\,t\equiv0\pmod{14}.$$

Its coefficient is its own inverse modulo $14$, because

$$\begin{aligned}(2x-1)^2&=4(x^2-x)+1\\&\equiv1\pmod{14}.\end{aligned}$$

Multiplying the congruence by $2x-1$ therefore leaves exactly one solution digit, $t\equiv-(2x-1)f\pmod{14}$. This is the lifting step of Hensel's lemma for the polynomial $x^2-x$, whose derivative $2x-1$ is invertible; the lemma is usually stated for powers of a prime, but the computation above works modulo $14$ directly. ^[5] Every member of the family of $7$ ends in the digit $7$, so there $2x-1\equiv13\equiv-1$ and $t\equiv f$; in the family of $8$, $2x-1\equiv15\equiv1$ and $t\equiv-f\pmod{14}$.

5. Carrying the quotient instead of squaring

For this $t$ the bracket is divisible by $14$, and the new quotient

$$f'=\frac{f+(2x-1)\,t+t^2M}{14}$$

satisfies $x'^2-x'=f'\cdot14M$. So the state $(x,f,M)$ passes from one length to the next with additions, multiplications by small numbers and one exact division by $14$; the long numbers are never squared. Both $x$ and $f$ stay below $M$, because $f=x(x-1)/M\le x$.

6. Worked example: five digits of each family

Family of $7$: we start with $x=7$, $M=14$ and $f=(49-7)/14=3$, so $t\equiv f$ gives $t=3$, $x'=7+3\cdot14=49=\mathtt{37}$ and $f'=(3+13\cdot3+9\cdot14)/14=12$. Indeed $49^2=2401=\mathtt{c37}$ ends in $\mathtt{37}$. Next, $t\equiv f=12$ is the digit $\mathtt c$, so $x'=49+12\cdot196=2401=\mathtt{c37}$, and the bracket $12+97\cdot12+144\cdot196=29400$ gives $f'=29400/14=2100$. Now $2100=150\cdot14$, so $t=0$: the ending $\mathtt{0c37}$ still solves $x^2\equiv x\pmod{14^4}$, but it is not a four-digit steady square. The state goes on with $f'=2100/14=150$, and $t\equiv150\equiv10$ gives $\mathtt{a0c37}=386561$.

Family of $8$: we start with $x=8$ and $f=(64-8)/14=4$, so $t\equiv-4\equiv10$, $x'=148=\mathtt{a8}$ and $f'=(4+15\cdot10+100\cdot14)/14=111$. Then $t\equiv-111\equiv1$ gives $\mathtt{1a8}=344$ with $f'=43$, next $t\equiv-43\equiv13$ gives $\mathtt{d1a8}=36016$ with $f'=33765$, and $t\equiv-33765\equiv3$ gives $\mathtt{3d1a8}=151264$.

7. Digit sums and leading zeros

Adding $t\cdot14^k$ to a number below $14^k$ writes the digit $t$ in front of its digits without any carry, so the digit sum grows by exactly $t$, and the new ending of length $k+1$ counts exactly when $t\ne0$. The endings that count are:

$n=1$: $1$, $7$, $8$; running total $16$;
$n=2$: $\mathtt{37}$ and $\mathtt{a8}$, digit sums $10$ and $18$; total $44$;
$n=3$: $\mathtt{c37}$ and $\mathtt{1a8}$, digit sums $22$ and $19$; total $85$;
$n=4$: only $\mathtt{d1a8}$, digit sum $32$; total $117$;
$n=5$: $\mathtt{a0c37}$ and $\mathtt{3d1a8}$, digit sums $32$ and $35$; total $184$.

Going on in the same way, the totals for the lengths up to $6$, $7$, $8$ and $9$ are $264$, $357$, $463$ and $582=\mathtt{2d8}$, the checkpoint of the statement.

8. The two families mirror each other

The sum $a_n+b_n$ is $\equiv1$ modulo $2^n$ and modulo $7^n$, hence $\equiv1\pmod{14^n}$. Both numbers are positive and below $14^n$, and the only number between $2$ and $2\cdot14^n$ of this kind is $14^n+1$, so

$$b_n=\bigl(14^n-1-a_n\bigr)+2.$$

The number $14^n-1$ is written with $n$ digits $\mathtt d$, so subtracting $a_n$ needs no borrowing and leaves the digits $13-\alpha_i$; adding $2$ then changes only the last digit, from $13-8=5$ to $7$. Therefore

$$\beta_0=7,\qquad\beta_i=13-\alpha_i\quad(i\ge1).$$

For instance $\mathtt c+1=3+\mathtt a=13$ for $\mathtt{c37}$ and $\mathtt{1a8}$; in base $10$, $376+625=1001$ is the same relation. ^[3] Consequently, for $n\ge2$ the leading digits of $a_n$ and $b_n$ add up to $13$, so they are never both $0$ and every length has one or two steady squares; and the digit sums of $a_n$ and $b_n$ add up to $15+13(n-1)=13n+2$. Up to $n=9$ the only leading zero is that of $\mathtt{0c37}$, with digit sum $22$, so the checkpoint is also

$$\begin{aligned}&1+\sum_{n=1}^{9}(13n+2)-22\\&\quad=1+603-22=582.\end{aligned}$$

The programs do not use this relation; it serves below as an independent check.

How the Code Works

The C++ program ^[9] starts in main. parse_arguments fills an Options record with the largest length $N$, which is $10000$ unless --n-max= gives another positive value, and with the flag --skip-checkpoints. Unless that flag is set, run_checkpoints runs first; then solve_sum_of_digit_sums(N) computes the total, and to_base14 writes it in base $14$ by repeated division, taking the digits from a character array holding 0123456789abcd.

One family is held in a SequenceState: the numbers $x$, $f$ and $M=14^k$ of Sections 4 and 5 as x, f and mod, of the arbitrary-precision type cpp_int from Boost.Multiprecision, ^[7] together with the newest digit leading_digit and the running digit_sum. make_initial_state(root) sets $x$ to the root, $M=14$, $f=(x^2-x)/14$, which is $3$ for the root $7$ and $4$ for the root $8$, and the digit sum to the root. The method step performs one extension of Section 4. It computes $a=(2x-1)\bmod14$ with the helper mod14, which also turns a negative remainder into one between $0$ and $13$, looks up the inverse of $a$ modulo $14$, and sets $t$ to the product of $(14-f\bmod14)\bmod14$ and that inverse, reduced modulo $14$. It then replaces $x$, $f$ and $M$ by $x+tM$, $(f+(2x-1)t+t^2M)/14$ and $14M$, stores $t$ in leading_digit and adds it to digit_sum, exactly as in Sections 5 and 7.

solve_sum_of_digit_sums first fills the table inv_mod14 by trying $b=1,2,\dots$ for each $a$ until $ab\equiv1\pmod{14}$: the six numbers $1,3,5,9,11,13$ coprime to $14$ receive the inverses $1,5,3,11,9,13$, ^[6] and the other entries stay $-1$. They are never read, because $a$ is always $13$ in the family of $7$ and $1$ in the family of $8$. The function then creates seq7 and seq8 and sets total to $1$ for the steady square $1$. For $n=1,\dots,N$ the loop adds the digit_sum of each family whose leading_digit is not $0$, and then extends both families by one digit, except after the last length. run_checkpoints requires the total $582$ for $N=9$ and checks that $x^2-x$ is divisible by mod for the first $20$ states of both families.

The Python and Java versions contain the same SequenceState class and the same loop, named solve_sum_of_digit_sums in Python and solveSumOfDigitSums in Java, but no options and no checkpoints: Python's solve and Java's main always compute the total for $N=10000$ and print it in base $14$. Python's integers have unlimited precision, and its mod14 is simply %. Java uses BigInteger; since its remainder takes the sign of the dividend, mod14 adds $14$ to a negative result as in C++, although no value here is negative. ^[8] The total adds at most $2N+1$ digit sums, each at most $13N$, so for $N=10000$ it stays below $2^{32}$; C++ keeps the sums in std::uint64_t and Java in long.

Complexity and Verification

The loop makes $N-1=9999$ extensions per family. Extension $k$ works with numbers of about $k$ base-14 digits, or $k\log_2 14\approx3.81k$ bits; at the end $14^{10000}$ has $38074$ bits. Each extension makes a fixed number of passes over such numbers: forming $2x-1$, two remainders modulo $14$, a few additions, multiplications by $t$, $t^2$ or $14$, and one exact division by $14$. Each pass costs time proportional to the length, so extension $k$ costs $O(k)$ and the whole run $O(N^2)$ operations on machine words, not $O(N)$ as a count of big-number operations would suggest. The memory is $O(N)$ bits for a handful of numbers of at most about $5$ kB each. On an Apple M4 Max the C++ program (clang++ -O2) takes about $0.2$ seconds, and the Java (OpenJDK 17) and Python 3.11 versions about $0.3$ seconds each. With the option --n-max= the C++ program needs about $0.9$, $3.7$ and $14$ seconds for $N=20000$, $40000$ and $80000$: doubling $N$ multiplies the time by four, as a quadratic cost predicts.

The C++ program checks the total $582$ of the statement and the congruence $x^2\equiv x$ for the first $20$ states of both families. Beyond these checks, a brute-force search over all numbers below $14^8\approx1.48\cdot10^9$ reproduces the steady squares of Section 7 and the running totals $16$, $44$, $85$, $117$, $184$, $264$, $357$ and $463$ for the lengths $1$ to $8$, which the C++ program also prints, in base $14$, with --n-max=1 to --n-max=8. The strings can also be built without lifting: let $u$ be the inverse of $2^N$ modulo $7^N$, taken between $1$ and $7^N-1$, which exists because the two powers are coprime. ^[6] Then $2^Nu$ is divisible by $2^N$, leaves the remainder $1$ modulo $7^N$ and lies below $14^N$, so by Section 2 it is $a_N$; for $N=2$, $u=37$ because $4\cdot37=3\cdot49+1$, and $a_2=148=\mathtt{a8}$. Computed in this way for $N=10000$, the endings of $a_N$ agree with the same construction carried out for every length $n\le300$ and, together with $b_n=14^n+1-a_n$, with the states $x$ of both families after each of the $9999$ extensions, where every $f$ equals $(x^2-x)/14^k$. The mirror relation of Section 8 holds at all $9999$ positions $i\ge1$, the first ten digits of both strings agree with the Wikipedia table, ^[2] and the total recomputed from these digits equals the common output of the three programs. The new leading digit is $0$ at $731$ of the $10000$ lengths in the family of $8$ and at $702$ in the family of $7$, about one in $14$. These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts used above. The lifting step, the mirror relation of the two families and all counts and timings are derived or computed in this article.

Project Euler 284 — Steady Squares. The official statement defines steady squares, with the examples 376² = 141376 in base 10 and c37² = aa0c37 in base 14, where c37 has the digit sum 18 in base 14, forbids leading zeros, gives the total 2d8 (582) of the digit sums for 1 ≤ n ≤ 9 and asks for the total for 1 ≤ n ≤ 10000, written in base 14 with lower-case letters.
Wikipedia — Automorphic number. A k-digit number in base b is automorphic if its square ends in the same k digits. A base b has 2^ω(b) fixed points, where ω(b) is the number of distinct prime factors of b, and by Hensel's lemma they lift to every power of b and to the b-adic integers. The table for base 14 lists the fixed points 0, 1, 7, 8 modulo 14, the automorphic numbers 0, 1, 7, 8, 37, A8, 1A8, C37, D1A8, … and the 14-adic fixed points …7337AA0C37 and …6AA633D1A8.
OEIS A003226 — Automorphic numbers: m² ends with m. The entry lists 0, 1, 5, 6, 25, 76, 376, 625, 9376, 90625, …, that is, 0 and 1 together with the families ending in 5 and in 6. A comment points out the complement relation, for example 10³ + 1 − 376 = 625: for every d-digit term n other than 0 and 1, the number 10ᵈ + 1 − n is again a term, which is the decimal form of the mirror relation of Section 8.
Wikipedia — Chinese remainder theorem. If the remainders of an integer after division by pairwise coprime moduli are prescribed, there is exactly one such integer between 0 and the product of the moduli minus 1; Section 2 applies this to the moduli 2ⁿ and 7ⁿ.
Wikipedia — Hensel's lemma. If f(r) ≡ 0 modulo pᵏ for a prime p and f′(r) is not divisible by p, then r lifts to a unique root s modulo pᵏ⁺¹ with s ≡ r modulo pᵏ, found by writing s = r + t·pᵏ and solving a linear congruence for t. Section 4 carries out the same computation for f(x) = x² − x modulo powers of 14.
Wikipedia — Modular multiplicative inverse. An inverse of a modulo m is an integer x with ax ≡ 1 modulo m, and it exists if and only if gcd(a, m) = 1. The programs tabulate the inverses modulo 14, and the independent check uses the inverse of 2ᴺ modulo 7ᴺ.
Boost.Multiprecision — cpp_int. The type cpp_int is a signed integer of arbitrary precision without a fixed bit limit; it holds x, f and 14ᵏ.
Python — numeric types states that integers have unlimited precision. Java — BigInteger describes immutable arbitrary-precision integers whose remainder follows the sign of the dividend, unlike mod, which is never negative.
C++ — Euler284.cpp. The linked immutable revision contains SequenceState, solve_sum_of_digit_sums, run_checkpoints and the options --n-max= and --skip-checkpoints. It is the implementation record for this article.

Problem 284 source code

C++

#include <array>
#include <cstdint>
#include <iostream>
#include <string>

#include <boost/multiprecision/cpp_int.hpp>

namespace {

using u64 = std::uint64_t;
using boost::multiprecision::cpp_int;

constexpr int kBase = 14;
constexpr int kDefaultN = 10000;

struct Options {
    int n_max = kDefaultN;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--n-max=", options.n_max)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n_max >= 1;
}

int mod14(const cpp_int& v) {
    cpp_int r = v % kBase;
    if (r < 0) {
        r += kBase;
    }
    return static_cast<int>(r);
}

std::string to_base14(u64 value) {
    constexpr std::array<char, 14> digits = {
        '0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 'a', 'b', 'c', 'd'};
    if (value == 0) {
        return "0";
    }
    std::string out;
    while (value > 0) {
        out.push_back(digits[static_cast<std::size_t>(value % 14)]);
        value /= 14;
    }
    std::reverse(out.begin(), out.end());
    return out;
}

struct SequenceState {
    cpp_int x;
    cpp_int f;
    cpp_int mod;
    int leading_digit;
    u64 digit_sum;

    void step(const std::array<int, kBase>& inv_mod14) {
        const cpp_int x_old = x;
        const cpp_int mod_old = mod;

        const int a = mod14(2 * x_old - 1);
        const int inv = inv_mod14[static_cast<std::size_t>(a)];
        const int f_mod = mod14(f);
        const int t = (kBase - f_mod) % kBase * inv % kBase;

        x = x_old + static_cast<cpp_int>(t) * mod_old;
        f = (f + (2 * x_old - 1) * t + static_cast<cpp_int>(t) * t * mod_old) / kBase;
        mod = mod_old * kBase;

        leading_digit = t;
        digit_sum += static_cast<u64>(t);
    }
};

SequenceState make_initial_state(int root) {
    SequenceState s;
    s.x = root;
    s.mod = kBase;
    s.f = (s.x * s.x - s.x) / kBase;
    s.leading_digit = root;
    s.digit_sum = static_cast<u64>(root);
    return s;
}

u64 solve_sum_of_digit_sums(const int n_max) {
    std::array<int, kBase> inv_mod14{};
    inv_mod14.fill(-1);
    for (int a = 1; a < kBase; ++a) {
        for (int b = 1; b < kBase; ++b) {
            if ((a * b) % kBase == 1) {
                inv_mod14[static_cast<std::size_t>(a)] = b;
                break;
            }
        }
    }

    SequenceState seq7 = make_initial_state(7);
    SequenceState seq8 = make_initial_state(8);

    u64 total = 1;  // steady square "1" for n = 1

    for (int n = 1; n <= n_max; ++n) {
        if (seq7.leading_digit != 0) {
            total += seq7.digit_sum;
        }
        if (seq8.leading_digit != 0) {
            total += seq8.digit_sum;
        }

        if (n == n_max) {
            break;
        }

        seq7.step(inv_mod14);
        seq8.step(inv_mod14);
    }

    return total;
}

bool run_checkpoints() {
    if (solve_sum_of_digit_sums(9) != 582ULL) {
        std::cerr << "Checkpoint failed for n<=9" << '\n';
        return false;
    }

    // Verify steady-square condition during the first few lifts.
    SequenceState seq7 = make_initial_state(7);
    SequenceState seq8 = make_initial_state(8);
    std::array<int, kBase> inv_mod14{};
    inv_mod14.fill(-1);
    for (int a = 1; a < kBase; ++a) {
        for (int b = 1; b < kBase; ++b) {
            if ((a * b) % kBase == 1) {
                inv_mod14[static_cast<std::size_t>(a)] = b;
                break;
            }
        }
    }

    for (int n = 1; n <= 20; ++n) {
        if ((seq7.x * seq7.x - seq7.x) % seq7.mod != 0) {
            std::cerr << "Idempotence checkpoint failed for root 7 at n=" << n << '\n';
            return false;
        }
        if ((seq8.x * seq8.x - seq8.x) % seq8.mod != 0) {
            std::cerr << "Idempotence checkpoint failed for root 8 at n=" << n << '\n';
            return false;
        }
        seq7.step(inv_mod14);
        seq8.step(inv_mod14);
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    const u64 answer = solve_sum_of_digit_sums(options.n_max);
    std::cout << to_base14(answer) << '\n';
    return 0;
}

Python

def mod14(v):
    return v % 14

def to_base14(value):
    digits = "0123456789abcd"
    if value == 0:
        return "0"
    out = []
    while value > 0:
        out.append(digits[value % 14])
        value //= 14
    return "".join(reversed(out))

class SequenceState:
    def __init__(self, root):
        self.x = root
        self.mod = 14
        self.f = (self.x * self.x - self.x) // 14
        self.leading_digit = root
        self.digit_sum = root

    def step(self, inv_mod14):
        x_old = self.x
        mod_old = self.mod
        
        a = mod14(2 * x_old - 1)
        inv = inv_mod14[a]
        f_mod = mod14(self.f)
        t = ((14 - f_mod) % 14 * inv) % 14
        
        self.x = x_old + t * mod_old
        self.f = (self.f + (2 * x_old - 1) * t + t * t * mod_old) // 14
        self.mod = mod_old * 14
        
        self.leading_digit = t
        self.digit_sum += t

def solve_sum_of_digit_sums(n_max):
    inv_mod14 = [-1] * 14
    for a in range(1, 14):
        for b in range(1, 14):
            if (a * b) % 14 == 1:
                inv_mod14[a] = b
                break
                
    seq7 = SequenceState(7)
    seq8 = SequenceState(8)
    
    total = 1
    
    for n in range(1, n_max + 1):
        if seq7.leading_digit != 0:
            total += seq7.digit_sum
        if seq8.leading_digit != 0:
            total += seq8.digit_sum
            
        if n == n_max:
            break
            
        seq7.step(inv_mod14)
        seq8.step(inv_mod14)
        
    return total

def solve():
    ans = solve_sum_of_digit_sums(10000)
    return to_base14(ans)

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

Java

import java.math.BigInteger;

public class Euler284 {

    static final int kBase = 14;

    static String toBase14(long value) {
        char[] digits = "0123456789abcd".toCharArray();
        if (value == 0)
            return "0";
        StringBuilder sb = new StringBuilder();
        while (value > 0) {
            sb.append(digits[(int) (value % 14)]);
            value /= 14;
        }
        return sb.reverse().toString();
    }

    static int mod14(BigInteger v) {
        BigInteger r = v.remainder(BigInteger.valueOf(kBase));
        if (r.signum() < 0) {
            r = r.add(BigInteger.valueOf(kBase));
        }
        return r.intValue();
    }

    static class SequenceState {
        BigInteger x;
        BigInteger f;
        BigInteger mod;
        int leadingDigit;
        long digitSum;

        SequenceState(int root) {
            this.x = BigInteger.valueOf(root);
            this.mod = BigInteger.valueOf(kBase);
            this.f = this.x.multiply(this.x).subtract(this.x).divide(BigInteger.valueOf(kBase));
            this.leadingDigit = root;
            this.digitSum = root;
        }

        void step(int[] invMod14) {
            BigInteger xOld = x;
            BigInteger modOld = mod;

            int a = mod14(xOld.multiply(BigInteger.TWO).subtract(BigInteger.ONE));
            int inv = invMod14[a];
            int fMod = mod14(f);
            int t = ((kBase - fMod) % kBase * inv) % kBase;

            x = xOld.add(modOld.multiply(BigInteger.valueOf(t)));
            BigInteger term2 = xOld.multiply(BigInteger.TWO).subtract(BigInteger.ONE).multiply(BigInteger.valueOf(t));
            BigInteger term3 = modOld.multiply(BigInteger.valueOf((long) t * t));
            f = f.add(term2).add(term3).divide(BigInteger.valueOf(kBase));
            mod = modOld.multiply(BigInteger.valueOf(kBase));

            leadingDigit = t;
            digitSum += t;
        }
    }

    static long solveSumOfDigitSums(int nMax) {
        int[] invMod14 = new int[kBase];
        for (int i = 0; i < kBase; i++)
            invMod14[i] = -1;
        for (int a = 1; a < kBase; a++) {
            for (int b = 1; b < kBase; b++) {
                if ((a * b) % kBase == 1) {
                    invMod14[a] = b;
                    break;
                }
            }
        }

        SequenceState seq7 = new SequenceState(7);
        SequenceState seq8 = new SequenceState(8);

        long total = 1;

        for (int n = 1; n <= nMax; n++) {
            if (seq7.leadingDigit != 0) {
                total += seq7.digitSum;
            }
            if (seq8.leadingDigit != 0) {
                total += seq8.digitSum;
            }

            if (n == nMax)
                break;

            seq7.step(invMod14);
            seq8.step(invMod14);
        }

        return total;
    }

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