Problem 340: Crazy Function

Project Euler Problem 340 Solution

EulerSolve provides an optimized solution for Project Euler Problem 340, Crazy Function, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The function $F(n)$ is defined by a four-layer nested recursion depending on parameters $(a,b,c)$, and the goal is to compute $$S(a,b,c)=\sum_{n=0}^{b}F(n)$$ for extremely large inputs. In the Project Euler instance, $a=21^7$, $b=7^{21}$, and $c=12^7$, so neither direct recursion nor iterating over all $0 \le n \le b$ is realistic. The solution works only because the recursion collapses to a closed form. Mathematical Approach Start from the definition $$F(n)=\begin{cases} n-c, & n \gt b, \\ F\bigl(a+F(a+F(a+F(a+n)))\bigr), & n \le b. \end{cases}$$ Let $$d=a-c,$$ so $d \gt 0$. This quantity is the key simplifier throughout the derivation. Outside the Recursive Region Whenever an argument $x$ is already larger than $b$, the recursion stops immediately and the rule becomes linear: $$F(x)=x-c,\qquad x \gt b.$$ Now look at one outer layer of the nested expression. It always has the shape $F(a+\text{something})$. If the inner value is already above $b$, then that layer contributes $$F(a+y)=a+y-c=y+d.$$ So after the computation has escaped above $b$, each remaining outer layer adds exactly $d$. The whole problem is therefore about determining how many layers are still genuinely recursive before we cross the boundary. Base Strip: $b-a \lt n \le b$ If $n$ lies in the top strip just below $b$, then $n+a \gt b$....

Detailed mathematical approach

Problem Summary

Three integers $a$, $b$ and $c$ are fixed, and the crazy function $F$ is defined by two rules. For every $n\gt b$ it is simply $F(n)=n-c$, and for every $n\le b$ it is

$$F(n)=F\bigl(a+F(a+F(a+F(a+n)))\bigr).$$

We must find the last nine digits of

$$S(a,b,c)=\sum_{n=0}^{b}F(n)$$

for $a=21^7=1801088541$, $c=12^7=35831808$ and $b=7^{21}$, that is $b=558545864083284007$. The statement also supplies three checkpoints for the small parameters $(a,b,c)=(50,2000,40)$: $F(0)=3240$, $F(2000)=2040$ and $S(50,2000,40)=5204240$. ^[1] The sum has $b+1\approx5.6\cdot10^{17}$ terms, each defined by a fourfold nested recursion, so nothing can be evaluated term by term.

The whole solution in one sentence: whatever the innermost call returns always lands above a fixed threshold, so the three calls around it stop recursing and each simply adds $a-c$, which turns $F$ into a simple formula on bands of $a$ consecutive inputs whose sum is a handful of arithmetic.

Mathematical Approach

1. Reading the definition

The first rule is explicit. The second is not: its innermost call starts at the larger argument $a+n$, while the other three arguments depend on values $F$ has yet to return. Reaching an argument above $b$ ends one call, but that a whole nest ends — and that the two rules pin $F$ down at all — has to be proved; Sections 3 to 5 do it. Such nested recursion, where the argument of the outer $F$ is built from $F$, resists the usual reasoning; the classical example is the McCarthy 91 function, invented as a verification test case and likewise owning a one-line closed form. ^[2]

Both can fail. For $a=4$, $b=10$, $c=6$ the evaluation of $F(9)$ obtains $F(13)=7$, then $F(11)=5$, and then asks for $F(9)$ again: it cycles. The equations themselves are satisfied there by many functions, for instance by $F(n)=n-6$ above $10$ together with any constant value $k\le6$ on every $n\le10$, as substitution confirms. We therefore assume what the published programs demand, namely

$$0\le c\lt a\le b,\qquad d=a-c\ge1,$$

which the target parameters satisfy, since $21^7\gt12^7$ and $7^{21}\gt21^7$. The unfolding never asks for an argument below $n$, so for $n\ge0$ negative arguments never arise.

2. The only rule that can finish a layer

A layer of the nest has the shape $F(a+x)$, where $x$ is the value computed one level deeper. If $x\gt b-a$ then $a+x\gt b$, the first rule applies, and

$$F(a+x)=a+x-c=x+d.$$

So a layer whose inner value exceeds the threshold $b-a$ contributes exactly $+d$, with no recursion. The threshold is $b-a$, not $b$: the layer adds $a$ before testing against $b$.

3. The top band, by hand

Let $n$ lie in the top band $b-a\lt n\le b$. By Section 2 the innermost layer gives $F(a+n)=n+d$, which again exceeds $b-a$, so the second layer gives $n+2d$, the third $n+3d$ and the fourth $n+4d$:

$$F(n)=n+4d\qquad (b-a\lt n\le b).$$

A concrete case: take $a=4$, $b=10$, $c=1$, so $d=3$ and the threshold is $b-a=6$. For $n=7$ the four layers are $F(11)=10$, $F(14)=13$, $F(17)=16$, $F(20)=19$: each added $3$, and indeed $F(7)=19=7+4d$. A checkpoint agrees too: for $(50,2000,40)$, $d=10$ and $n=2000$ lies in the top band $1951\le n\le2000$, so $F(2000)=2040$.

4. One step down the ladder

Now let $0\le n\le b-a$, so $m=n+a\le b$ and $F(m)$ may be thought of as already computed. The innermost layer of $F(n)$ is exactly $F(m)$, so if $F(m)\gt b-a$, Section 2 makes each of the three remaining layers add $d$ and the nest collapses to

$$F(n)=F(n+a)+3d.$$

The descent rule is the engine of the solution, and its only premise is the overshoot $F(m)\gt b-a$, which cannot be waved away: $F(m)\ge m$ would not help, since $m$ itself may lie far below $b$. It is proved together with the closed form.

5. The closed form

For $0\le n\le b$ define the band index

$$q=q(n)=\Bigl\lfloor\frac{b-n}{a}\Bigr\rfloor,$$

the number of times $a$ can be added to $n$ without passing $b$. Integers with the same band index form a block of $a$ consecutive values, except the lowest block, which may be shorter.

Proposition. For every $0\le n\le b$,

$$F(n)=n+4d+(a+3d)\,q(n),$$

and moreover $F(n)\gt b-a$.

Proof by induction on $q(n)$; the two statements must be carried together, because the formula in one band uses the overshoot in the band above. ^[3]

Base $q=0$. This is the band $b-a\lt n\le b$, where Section 3 gives $F(n)=n+4d$, matching the formula. Since $n\gt b-a$ and $d\ge1$, also $F(n)\gt b-a$.
Step $q\ge1$. Then $b-n\ge a$, so $m=n+a\le b$ and $q(m)=q-1$. The induction hypothesis gives $F(m)\gt b-a$ and $F(m)=m+4d+(a+3d)(q-1)$, which simplifies to $F(m)=n+d+(a+3d)q$. The overshoot licenses the descent rule of Section 4, so $F(n)=F(m)+3d$${}=n+4d+(a+3d)q$, as claimed. For the overshoot at $n$, rewrite the formula as $F(n)=(n+aq)+4d+3dq$. The definition of $q$ says $n+aq\le b\lt n+a(q+1)$, hence $n+aq\gt b-a$ and therefore $F(n)\gt b-a$. □

The induction answers Section 1 as well: band by band it forces the value of any function obeying the two rules, so under $0\le c\lt a\le b$ the rules determine $F$ completely, and the evaluation terminates — the top band needs the first rule alone, each lower band one call in the band above plus three uses of the first rule. It also explains the two constants: $4d$ is the four layers of the nest, each worth $+d$ once the climb is over, and the coefficient $a+3d$ adds the $3d$ one descent step pays to the $a$ by which $n$ lies below the band above. With $k$ nested calls the same argument gives $F(n)=n+kd+(a+(k-1)d)q$, confirmed by brute force for $k=1,\dots,6$.

6. What the function looks like

Inside a band $F(n+1)=F(n)+1$; at a boundary the index drops, so $F(n)$ exceeds $F(n+1)$ by $a+3d-1$ and the graph is a sawtooth of rising unit steps. For $a=4$, $b=10$, $c=1$ the three bands are

$q=0$, $n=7,8,9,10$: $F=19,20,21,22$;
$q=1$, $n=3,4,5,6$: $F=28,29,30,31$;
$q=2$, $n=0,1,2$: $F=38,39,40$.

The drop from $F(2)=40$ to $F(3)=28$ is $12=a+3d-1$, and the lowest band holds only the $r+1=3$ inputs $0\le n\le r$ with $r=b\bmod a=2$. The checkpoint for $F(0)$ follows at once: for $(50,2000,40)$ the band index of $0$ is $q=\lfloor2000/50\rfloor=40$, so $F(0)=F(2000)+40\cdot 3d$${}=2040+1200=3240$.

7. Summing the closed form

Insert the closed form into the sum and split it into three pieces:

$$S(a,b,c)=\frac{b(b+1)}{2}+4d(b+1)+(a+3d)\,T,$$

where $T=\sum_{n=0}^{b}q(n)$ collects the band indices. The first piece is $0+1+\cdots+b$: added to its own reversal it pairs every term to $b$, so twice the sum is $b(b+1)$. The second is $4d$ taken $b+1$ times. In $T$, substituting $k=b-n$ only reverses the order of summation, so $T=\sum_{k=0}^{b}\lfloor k/a\rfloor$.

8. The floor sum

Write the Euclidean division $b=aq_0+r$, with unique quotient $q_0=\lfloor b/a\rfloor$ and remainder $0\le r\lt a$. ^[4] For each $0\le j\lt q_0$ the value $\lfloor k/a\rfloor=j$ occurs for exactly the $a$ integers $k=aj,\dots,aj+a-1$, and for $j=q_0$ for the $r+1$ integers $k=aq_0,\dots,b$. Hence

$$\begin{aligned}T&=a\sum_{j=0}^{q_0-1}j+(r+1)\,q_0\\&=\frac{a\,q_0(q_0-1)}{2}+(r+1)\,q_0.\end{aligned}$$

For $a=4$ and $b=10$ the floors are $0,0,0,0,1,1,1,1,2,2,2$ with sum $10$, and the formula gives $4\cdot2\cdot1/2+3\cdot2=10$. Both worked examples now come out in three terms. For $(4,10,1)$: $55+132+13\cdot10=317$, and a brute force over the definition sums to $317$ as well. For $(50,2000,40)$: $q_0=40$, $r=0$, $T=39040$, and

$$2001000+80040+80\cdot39040=5204240,$$

the last checkpoint of the statement. Everything is now $O(1)$ arithmetic.

9. Taking the last nine digits

The last nine digits of a non-negative integer are its residue modulo $10^9$, written with leading zeros if the residue is below $10^8$. Congruences are compatible with addition, subtraction and multiplication, so each of the three pieces may be reduced on its own. ^[5]

Halving is the trap. An inverse of $2$ exists only for moduli coprime to $2$, and $\gcd(2,10^9)=2$, so $b(b+1)/2$ must not be obtained by reducing $b(b+1)$ and then dividing. For $b=999999999$ the reduced factors $999999999$ and $0$ have product $0$, which halves to $0$, whereas $b(b+1)/2\equiv500000000$. The way out: one of $b$, $b+1$ is even, so halve that factor exactly before reducing, giving $(b/2)(b+1)$ for even $b$ and $b\,\bigl((b+1)/2\bigr)$ otherwise. The same objection hits the $a\,q_0(q_0-1)/2$ inside $T$, where the programs take the other route and evaluate $T$ exactly in a wide integer type, reducing only the finished value: for the target $T\approx8.7\cdot10^{25}$ needs $87$ bits, past the $64$ available, and the exact $S\approx7.7\cdot10^{35}$ needs $120$. Once every factor is below $10^9$, each product stays below $10^{18}$ and fits in a signed 64-bit integer.

How the Code Works

The C++ program keeps the parameters in an Options struct whose defaults are the target values $21^7$, $7^{21}$ and $12^7$. parse_u64_after_prefix and parse_arguments accept --a=, --b=, --c= and --skip-checkpoints; parse_arguments returns false unless the parsed values satisfy $a\gt c$ and $b\ge a$, the hypotheses of Section 1, and main then exits with status $1$; failed checkpoints exit with status $2$. Neither check looks for overflow, in the parsing or in the arithmetic: the accepted input $a=b=2^{62}$, $c=0$ makes $a+3d$ exactly $2^{64}$, which wraps to zero, so the floor term vanishes and the residue comes out too small by $709551616$. The target stays far inside the safe range.

F_closed is the Proposition of Section 5: it returns n - c above $b$ and otherwise n + 4*d + (a + 3*d)*q with q = (b - n)/a. S_closed_mod_1e9 is Section 9 line by line: it branches on the parity of $b$ to halve the even factor of $b(b+1)$, forms term_linear from $(b+1)$ and $4d$ with both factors reduced first, accumulates sum_floor as $T$ in a signed 128-bit integer, and multiplies its residue by $(a+3d)\bmod 10^9$. The helper mod_mul widens to __uint128_t before the multiplication, a compiler extension of GCC and Clang. ^[6] One detail is worth noting: sum_n is first assigned a reduced product divided by $2$, then overwritten by whichever parity branch applies, so only the branch value is used — fortunately, since the discarded expression is exactly the invalid halving of Section 9.

F_bruteforce implements the definition literally with a std::function and an unordered_map memo, S_bruteforce sums it, and S_closed_exact_small evaluates $S$ exactly in 128 bits and returns it as a 64-bit value, which is safe only for the small parameters it is called with. run_checkpoints verifies $F(0)=3240$, $F(2000)=2040$ and $S(50,2000,40)=5204240$, then compares the closed form against the definition for all $33$ pairs $(a,c)$ with $4\le a\le9$ and $1\le c\lt a$, using $b=60+a$: that is $2245$ comparisons of $F$ and $33$ of $S$.

The Python version inlines the arithmetic of Section 9 in solve, with the same parity branch; its // is floor division and its integers have unlimited precision, so $T$ needs no special type. ^[7] The Java version mirrors the C++ in sClosedMod1e9 but builds $T$ from BigInteger values with multiply, subtract, divide and add, then remainder and longValue. ^[8] Its modMul needs both arguments below $10^9$, which keeps the product below $10^{18}$ where a long suffices; every call site reduces beforehand, including $4(d\bmod 10^9)$. Neither port carries the checkpoints, and both call solve with the target parameters fixed in the source, which the C++ reads from the command line instead. All three print the residue as it stands, so a residue below $10^8$ would need a leading zero to be read as nine digits.

Complexity and Verification

With the closed form, $F(n)$ and $S(a,b,c)$ cost a fixed number of additions, multiplications and divisions in a fixed number of variables. For the target the operands fit in $128$ bits, so this is $O(1)$ machine operations and $O(1)$ memory; for arbitrary parameters the operands have $O(L)$ bits, where $L$ is the bit length of $b$, hence $O(L^2)$ bit operations with schoolbook multiplication. The definition instead would ask for all $b+1$ values, applying the recursive rule $q(n)+1$ times for the input $n$: once in the top band, up to $\lfloor b/a\rfloor+1$, about $3.1\cdot10^8$, times near $n=0$. The C++ checkpoints make $2245$ comparisons of $F$ and $33$ of $S$; since S_bruteforce re-sums every triple, $4490$ brute-force evaluations are started in all, each with its own memo table and all with $b\le69$. The whole program finished in about five milliseconds in the build used here.

Verified independently for this article. A literal brute force over the two rules reproduces the closed form of Section 5 on all $4030$ triples with $1\le a\le12$, $0\le c\lt a$ and $a\le b\le59$, that is $139373$ values of $F$, and the sum formula of Sections 7 and 8 on every one of those triples. In the same family the top-band formula holds in all $32916$ cases, the descent rule in all $106457$, and the overshoot at all $139373$ values; the floor-sum identity was checked for $1\le a\le19$, $0\le b\le399$, and the $k$-fold generalisation for $k=1,\dots,6$ on all triples with $a\le9$, $b\le39$, that is $214470$ values. The modular pipeline of Section 9 reproduces the exact sum modulo $10^9$ on $405$ triples: five fixed, among them $b=10^{13}+7$, and $400$ random ones with $a\lt10^6$, $b\lt10^{12}$. A cycle detector confirmed the circularity at $(a,b,c)=(4,10,6)$. The three published programs print the same line, the C++ checkpoints pass, and with --a=50 --b=2000 --c=40 and --a=4 --b=10 --c=1 the C++ returns $5204240$ and $317$. ^[9]

Footnotes and References

The references below supply the background facts used above; the closed form, the overshoot, the descent rule and the floor sum are derived explicitly in this article.

Project Euler 340 — Crazy Function. The official statement fixes integers a, b, c, defines F(n) = n − c for n > b and F(n) = F(a + F(a + F(a + F(a + n)))) for n ≤ b, sets S(a, b, c) as the sum of F(n) for 0 ≤ n ≤ b, gives F(0) = 3240, F(2000) = 2040 and S(50, 2000, 40) = 5204240 for a = 50, b = 2000, c = 40, and asks for the last nine digits of S(21^7, 7^21, 12^7).
Wikipedia — McCarthy 91 function. M(n) = n − 10 for n > 100 and M(M(n + 11)) for n ≤ 100; Manna, Pnueli and McCarthy introduced it in 1970 as a challenge for automated program verification, chosen because it is nested-recursive rather than defined from f(n − 1). Despite the nesting M(n) = 91 for all n ≤ 100, the collapse that Section 5 establishes for the crazy function.
Wikipedia — Mathematical induction. The method proves a base case and then shows that validity at one index implies validity at the next; the article also describes variants that induct on a derived quantity such as ⌊k/2⌋. Section 5 applies ordinary induction to the band index, carrying the closed form and the overshoot together from one band to the next.
Wikipedia — Euclidean division. For integers a and b with b ≠ 0 there exist unique integers q and r with a = bq + r and 0 ≤ r < |b|. Section 8 uses this uniqueness to split 0, 1, …, b into blocks of a integers sharing one quotient, plus a final block of r + 1 integers.
Wikipedia — Modular arithmetic. Congruence modulo m is compatible with addition, subtraction and multiplication, and a has a multiplicative inverse modulo m if and only if a is coprime to m. This is why the three pieces of S may be reduced separately but the division by 2 may not be performed modulo 10⁹.
GCC manual — 128-bit integers. __int128 and unsigned __int128 are an extension providing 128-bit integer scalar types on targets with a wide enough integer mode. The C++ program uses them for the exact value of T and inside mod_mul.
Python documentation — numeric types. Python integers have unlimited precision, and the // operator is integer division whose result is always rounded towards minus infinity, which is the floor used throughout this article.
Java API — BigInteger. Immutable arbitrary-precision integers whose arithmetic operations mimic Java's primitive integer semantics, including multiply, divide, add and remainder. The Java port uses them for T, which does not fit in a long.
C++ — Euler340.cpp. The linked immutable revision contains F_closed, S_closed_mod_1e9, S_closed_exact_small, F_bruteforce, S_bruteforce and run_checkpoints, with the target parameters as defaults. It is the implementation record for this article; the wider cross-checks above were run separately and are not part of that file.

Problem 340 source code

C++

#include <cstdint>
#include <functional>
#include <iostream>
#include <string>
#include <unordered_map>

namespace {

using u64 = std::uint64_t;
using i128 = __int128_t;

struct Options {
    u64 a = 1'801'088'541ULL;                    // 21^7
    u64 b = 558'545'864'083'284'007ULL;          // 7^21
    u64 c = 35'831'808ULL;                       // 12^7
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    u64 parsed = 0ULL;
    for (char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        parsed = parsed * 10ULL + static_cast<u64>(ch - '0');
    }
    value = parsed;
    return true;
}

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

u64 F_closed(const u64 n, const u64 a, const u64 b, const u64 c) {
    if (n > b) {
        return n - c;
    }
    const u64 d = a - c;
    const u64 q = (b - n) / a;
    return n + 4ULL * d + (a + 3ULL * d) * q;
}

u64 S_closed_mod_1e9(const u64 a, const u64 b, const u64 c) {
    constexpr u64 MOD = 1'000'000'000ULL;
    const u64 d = a - c;
    const u64 q = b / a;
    const u64 r = b % a;

    const auto mod_mul = [](u64 x, u64 y) -> u64 {
        return static_cast<u64>((static_cast<__uint128_t>(x) * y) % 1'000'000'000ULL);
    };

    u64 sum_n = static_cast<u64>((static_cast<__uint128_t>(b % MOD) * ((b + 1ULL) % MOD) / 2ULL) % MOD);
    if ((b & 1ULL) == 0ULL) {
        sum_n = mod_mul((b / 2ULL) % MOD, (b + 1ULL) % MOD);
    } else {
        sum_n = mod_mul(b % MOD, ((b + 1ULL) / 2ULL) % MOD);
    }

    u64 term_linear = mod_mul((b + 1ULL) % MOD, (4ULL * (d % MOD)) % MOD);

    i128 sum_floor =
        static_cast<i128>(a) * static_cast<i128>(q) * static_cast<i128>(q - 1ULL) / 2 +
        static_cast<i128>(r + 1ULL) * static_cast<i128>(q);
    const u64 sum_floor_mod = static_cast<u64>(sum_floor % MOD);
    u64 term_floor = mod_mul((a + 3ULL * d) % MOD, sum_floor_mod);

    return (sum_n + term_linear + term_floor) % MOD;
}

u64 S_closed_exact_small(const u64 a, const u64 b, const u64 c) {
    const u64 d = a - c;
    const u64 q = b / a;
    const u64 r = b % a;
    const i128 sum_floor =
        static_cast<i128>(a) * static_cast<i128>(q) * static_cast<i128>(q - 1ULL) / 2 +
        static_cast<i128>(r + 1ULL) * static_cast<i128>(q);
    const i128 result =
        static_cast<i128>(b) * static_cast<i128>(b + 1ULL) / 2 +
        static_cast<i128>(b + 1ULL) * static_cast<i128>(4ULL * d) +
        static_cast<i128>(a + 3ULL * d) * sum_floor;
    return static_cast<u64>(result);
}

u64 F_bruteforce(const u64 n, const u64 a, const u64 b, const u64 c) {
    std::unordered_map<u64, u64> memo;
    std::function<u64(u64)> rec = [&](u64 x) -> u64 {
        auto it = memo.find(x);
        if (it != memo.end()) {
            return it->second;
        }
        u64 out = 0ULL;
        if (x > b) {
            out = x - c;
        } else {
            out = rec(a + rec(a + rec(a + rec(a + x))));
        }
        memo.emplace(x, out);
        return out;
    };
    return rec(n);
}

u64 S_bruteforce(const u64 a, const u64 b, const u64 c) {
    u64 sum = 0ULL;
    for (u64 n = 0ULL; n <= b; ++n) {
        sum += F_bruteforce(n, a, b, c);
    }
    return sum;
}

bool run_checkpoints() {
    if (F_closed(0, 50, 2000, 40) != 3240ULL) {
        std::cerr << "Checkpoint failed: F(0) for sample parameters" << '\n';
        return false;
    }
    if (F_closed(2000, 50, 2000, 40) != 2040ULL) {
        std::cerr << "Checkpoint failed: F(2000) for sample parameters" << '\n';
        return false;
    }
    if (S_closed_exact_small(50, 2000, 40) != 5'204'240ULL) {
        std::cerr << "Checkpoint failed: S(50,2000,40)" << '\n';
        return false;
    }
    for (u64 a = 4; a <= 9; ++a) {
        for (u64 c = 1; c < a; ++c) {
            const u64 b = 60 + a;
            for (u64 n = 0; n <= b; ++n) {
                if (F_closed(n, a, b, c) != F_bruteforce(n, a, b, c)) {
                    std::cerr << "Checkpoint failed: brute-force cross-check for F(n)" << '\n';
                    return false;
                }
            }
            if (S_closed_exact_small(a, b, c) != S_bruteforce(a, b, c)) {
                std::cerr << "Checkpoint failed: brute-force cross-check for S" << '\n';
                return false;
            }
        }
    }
    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;
    }
    std::cout << S_closed_mod_1e9(options.a, options.b, options.c) << '\n';
    return 0;
}

Python

def solve():
    a = 1801088541  # 21^7
    b = 558545864083284007  # 7^21
    c = 35831808  # 12^7
    MOD = 1_000_000_000

    d = a - c
    q = b // a
    r = b % a

    def mod_mul(x, y):
        return (x * y) % MOD

    if b % 2 == 0:
        sum_n = mod_mul((b // 2) % MOD, (b + 1) % MOD)
    else:
        sum_n = mod_mul(b % MOD, ((b + 1) // 2) % MOD)

    term_linear = mod_mul((b + 1) % MOD, (4 * (d % MOD)) % MOD)

    sum_floor = a * q * (q - 1) // 2 + (r + 1) * q
    sum_floor_mod = sum_floor % MOD
    term_floor = mod_mul((a + 3 * d) % MOD, sum_floor_mod)

    return str((sum_n + term_linear + term_floor) % MOD)

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

Java

import java.util.*;

public class Euler340 {

    static final long MOD = 1000000000L;

    static long modMul(long x, long y) {
        // Can overflow 64-bit if x,y ~ 10^9, so use 128-bit via BigInteger or split.
        // Or since x, y < 10^9, x*y < 10^18, which fits in long perfectly!
        return (x * y) % MOD;
    }

    static long sClosedMod1e9(long a, long b, long c) {
        long d = a - c;
        long q = b / a;
        long r = b % a;

        long sumN = 0;
        if ((b & 1) == 0) {
            sumN = modMul((b / 2) % MOD, (b + 1) % MOD);
        } else {
            sumN = modMul(b % MOD, ((b + 1) / 2) % MOD);
        }

        long termLinear = modMul((b + 1) % MOD, (4 * (d % MOD)) % MOD);

        java.math.BigInteger aBig = java.math.BigInteger.valueOf(a);
        java.math.BigInteger qBig = java.math.BigInteger.valueOf(q);
        java.math.BigInteger sumFloor = aBig.multiply(qBig).multiply(qBig.subtract(java.math.BigInteger.ONE))
                .divide(java.math.BigInteger.valueOf(2))
                .add(java.math.BigInteger.valueOf(r + 1).multiply(qBig));

        long sumFloorMod = sumFloor.remainder(java.math.BigInteger.valueOf(MOD)).longValue();
        long termFloor = modMul((a + 3 * d) % MOD, sumFloorMod);

        return (sumN + termLinear + termFloor) % MOD;
    }

    public static String solve() {
        long a = 1801088541L;
        long b = 558545864083284007L;
        long c = 35831808L;
        long ans = sClosedMod1e9(a, b, c);
        return String.valueOf(ans);
    }

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