Problem 335: Gathering the Beans

Project Euler Problem 335 Solution

EulerSolve provides an optimized solution for Project Euler Problem 335, Gathering the Beans, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Place $n$ bowls in a circle, with one bean in each bowl. Starting from a chosen bowl, Peter removes all beans from the current bowl and drops them one by one into successive bowls clockwise. The next move starts from the bowl that received the last bean. Let $M(n)$ be the number of moves required until the initial all-ones configuration appears again. The problem asks for $$\sum_{k=0}^{10^{18}} M(2^k+1)\pmod{7^9}.$$ A direct simulation is easy for small $n$, but the upper limit $10^{18}$ makes brute force completely infeasible. The solution therefore isolates a closed form for the special arguments $2^k+1$, then sums that formula as three geometric series. Mathematical Approach State Model and Direct Definition of $M(n)$ Number the bowls $0,1,\dots,n-1$, let $b_i$ be the bean count in bowl $i$, and let $p$ be the current bowl. If the current bowl contains $t=b_p$ beans, one move performs $$b_p\leftarrow 0,\qquad b_{p+j\pmod n}\leftarrow b_{p+j\pmod n}+1\quad(1\le j\le t),\qquad p\leftarrow p+t\pmod n.$$ The initial state is $(1,1,\dots,1)$. The helper routine simulate_M in the C++ source implements exactly this rule and repeats it until the all-ones state reappears. This gives a literal definition of $M(n)$....

Detailed mathematical approach

Problem Summary

Peter puts $x$ bowls in a circle with one bean in each. He takes all the beans out of one bowl and drops them one by one into the bowls that follow, clockwise; the next move starts at the bowl the last bean landed in, and he carries on until one bean sits in every bowl again. Let $M(x)$ be the number of moves. The statement's animation shows $M(5)=15$ and adds $M(100)=10920$. We must find ^[1]

$$\sum_{k=0}^{10^{18}} M(2^k+1)\pmod{7^9}.$$

The whole solution in one sentence: with $n=2^k+1$ bowls the moves split into $2^k$ blocks that each start at one fixed bowl, block $l$ skips a bowl exactly when the binary digits of $l$ all occur in that bowl's number minus one, counting the skips gives $M(2^k+1)=2^{k+1}+4^k-3^k$, and summing that is three geometric series modulo $7^9$.

Mathematical Approach

1. The move rule

Number the bowls $0,1,\dots,n-1$ clockwise, let $b_j$ be the beans in bowl $j$ and let $p$ be the bowl the next move starts from. With $t=b_p$, one move does

$$b_p\leftarrow 0,\qquad b_{p+r}\leftarrow b_{p+r}+1$$

for $r=1,\dots,t$, all indices modulo $n$, and then $p\leftarrow p+t$. The start is $b_j=1$ for every $j$. The total $\sum_j b_j=n$ never changes, so no move sows more than $n$ beans and the sowing wraps round the circle at most once. This is the relay sowing of mancala games, played on a ring. ^[2] Rotating the circle commutes with the rule, so the bowl Peter starts from only relabels the bowls and cannot change $M(n)$; from here on it is bowl $0$.

2. The starting position always comes back

Every move can be undone. Say one has just finished with the pointer on bowl $p$. If every bowl holds a bean then, the beans totalling $n$, each holds one, and the only move with that result empties $n$ beans out of bowl $p$. Otherwise walk backwards from $p$ while the bowls are non-empty: a move empties the bowl it starts from and raises every bowl from $q+1$ to $p$ by one, so the first empty bowl $q$ is the only possible source, and it sowed $p-q$ modulo $n$ beans.

Take a state to be the bean counts with the pointer, whose bowl is never empty since the pointer follows the last bean. The reconstruction makes the move map injective on this finite set of states, hence a permutation, so repeated moves do return to the start. The all-ones arrangement may reappear earlier, with the pointer elsewhere, and $M(n)$ counts that first moment; Section 8 locates it.

3. Laps

Group the moves into laps: a lap begins with a move from bowl $0$ and ends just before the next such move. Section 8 proves that for $n=2^k+1$ there are exactly $N=2^k$ laps, each beginning with one bean in bowl $0$. Every lap but the last advances the pointer exactly $n$ bowls, so the arcs its moves sweep tile the $n$ bowls once and each bowl gains one bean; bowl $0$, emptied by the first move, takes the last bean and is ready again. The last lap differs: its final move sows $n$ beans right round the circle.

4. How far a move jumps

Let $c_l(j)$ be the beans in bowl $j$ when lap $l$ begins, so $c_l(0)=1$. When the pointer arrives at a bowl $j\ge1$ during lap $l$, that bowl has taken its one bean for the lap and nothing else, so it holds $c_l(j)+1$ beans and the move jumps $c_l(j)+1$ bowls. Before the last lap, a bowl $j\ge1$ passed over gains a bean while one visited is emptied, so $c_{l+1}(j)$ is $0$ if lap $l$ visits bowl $j$ and $c_l(j)+1$ otherwise, with $c_{l+1}(0)=1$. Hence $c_l(j)+1$ counts the laps since bowl $j$ was last visited: a move jumps one bowl further for every lap in which its bowl was skipped.

5. Fifteen moves for five bowls

Take $n=5$. Lap 1 starts from $(1,1,1,1,1)$; its moves start at bowls $0,1,3$ and sow $1,2,2$ beans, passing bowls $2$ and $4$ over, and it ends at $(1,0,2,0,2)$. Lap 2 moves at $0,1,2$ sowing $1,1,3$, ending at $(1,0,0,1,3)$; lap 3 at $0,1,2,3$ sowing $1,1,1,2$, ending at $(1,0,0,0,4)$; lap 4 at $0,1,2,3,4$ sowing $1,1,1,1,5$, the last move spreading all five beans one per bowl. The lap lengths are $3,3,4,5$, so

$$M(5)=3+3+4+5=15,$$

as the statement says. Bowl $4$ is the one to watch: skipped on laps 1, 2 and 3, it grows to $4$ beans, and the move that starts there at last jumps the full circle of $5$.

6. The binary rule

For integers $u,v\ge0$ write $u\preceq v$ when every binary digit $1$ of $u$ is also a $1$ of $v$, that is when $u$ is a submask of $v$. Fix $n=2^k+1$ and $N=2^k$, and write a bowl other than bowl $0$ as $i+1$ with $0\le i\le N-1$.

Claim. The process has exactly $N$ laps, and the moves of lap $l$ start at bowl $0$ and at the bowls $i+1$ with $l\not\preceq i$, in increasing order of $i$, and nowhere else.

So bowl $i+1$ is passed over on lap $l$ exactly when $l\preceq i$. For $n=5$: lap $1$ skips $i\in\{1,3\}$, bowls $2$ and $4$; lap $2$ skips $i\in\{2,3\}$; lap $3$ only $i=3$; lap $4$ nothing, as no $i\le3$ carries digit $4$ — the pattern of Section 5.

With Section 4 the claim fixes every bean count. Let $H(i,l)$ be the length of the longest run $l-1,l-2,\dots$ of integers all submasks of $i$, stopping below $0$. Since $0\preceq i$ always, $H(i,1)=1$, the bean each bowl starts with; for larger $l$ the value $0$ counts only if the run reaches it. As lap $m$ skips bowl $i+1$ exactly when $m\preceq i$, the laps since its last visit come to

$$c_l(i+1)=H(i,l).$$

7. The gap lemma

A small case first. For $l=2$ and $i=1$ the indices $2$ and $3$ both carry digit $2$ and are skipped, the next without it is $4$, and $1,0$ are both submasks of $1$: $H(1,2)=2$ and the gap is $1+2=3$. Note too that adding $1$ to $i$ clears the digits below the lowest digit $0$ of $i$, at position $z$ say, sets digit $z$ and leaves the digits above alone.

Lemma. Let $l\ge1$, let $2^a$ be its lowest digit $1$ and write $l=2^a+L$, so every digit $1$ of $L$ lies above position $a$. Let $i\ge0$ with $l\not\preceq i$. Then $H(i,l)=2^a$ if $L\preceq i$ and $i\equiv 2^a-1\pmod{2^{a+1}}$, and $H(i,l)=0$ otherwise; and in all cases the smallest $i'\gt i$ with $l\not\preceq i'$ is

$$i'=i+1+H(i,l).$$

Proof. Note $l-1=L+(2^a-1)$ carries every digit of $L$ and all digits below $a$.

Case A: $L\not\preceq i$. Then $l-1\not\preceq i$, so $H=0$. Also $l\not\preceq i+1$: a digit $d$ of $L$ is missing from $i$ and the only digit $i+1$ gains is $z$, so $l\preceq i+1$ would force $z=d\gt a$ — but then digit $a$ of $i+1$ was cleared to $0$, where $l$ has a $1$. So $i'=i+1$.
Case B: $L\preceq i$, and some digit below $a$ is $0$ in $i$. Here digit $a$ of $i$ is $0$, again $l-1\not\preceq i$ gives $H=0$, and $z\lt a$ leaves digit $a$ of $i+1$ at $0$, so $i'=i+1$.
Case C: $L\preceq i$, and every digit below $a$ is $1$ in $i$. Digit $a$ of $i$ is again $0$, so $i=2^{a+1}u+2^a-1$. For $1\le r\le 2^a$, $l-r=L+(2^a-r)$ with $0\le 2^a-r\le 2^a-1$, whose digits are $L$'s plus digits below $a$, all present in $i$; so $l-r\preceq i$ and $H\ge 2^a$. The next candidate $l-2^a-1=L-1$ is negative when $L=0$; for $L\gt0$, subtracting one sets every digit below the lowest $1$ of $L$, which lies above position $a$, so $L-1$ has digit $a$ equal to $1$, which $i$ lacks. Either way $H=2^a$. For the gap, $i+r=2^{a+1}u+2^a+(r-1)$ with $1\le r\le 2^a$ has digit $a$ equal to $1$ and the same digits above $a$ as $i$, so $l\preceq i+r$; and $i+2^a+1=2^{a+1}(u+1)$ has digit $a$ equal to $0$, so $i'=i+2^a+1$.

8. Why the binary rule holds

Induct on $l$, carrying $c_l(i+1)=H(i,l)$ along; for $l=1$ the run is the single value $0$, so $H(i,1)=1$ everywhere, one bean per bowl. Assume lap $l$ begins with one bean in bowl $0$ and $H(i,l)$ beans in bowl $i+1$. Its first move carries that bean to bowl $1$, so the second starts at $i=0$, which belongs on the list since $l\not\preceq 0$ for every $l\ge1$. From a bowl $i+1$ with $l\not\preceq i$ the move jumps $H(i,l)+1$ bowls and lands, by the lemma, on bowl $i'+1$ with $i'$ the next integer above $i$ that is not a superset of $l$. Chaining from $i=0$ gives the claim's bowls in increasing order.

For $1\le l\lt N$ every digit of $l$ lies below position $k$, so $l\preceq N-1$ and bowl $N$ is skipped, while $l\not\preceq N$ because $N=2^k$ has no digit below $k$. Every index between the last visited $i$ and $N$ is a superset of $l$ while $N$ is not, so the lemma sends the final move to bowl $N+1$, that is bowl $0$. No earlier move crosses that boundary, so the lap is exactly one circuit and Section 3 applies; the new counts are $0$ where $l\not\preceq i$ and $H(i,l)+1$ elsewhere — in both cases $H(i,l+1)$.

Lap $N$ is the last. Its number $2^k$ has its only digit at position $k$, which no $i\le N-1$ carries, so nothing is skipped, and the lemma gives $H(i,N)=2^k$ for $i=N-1$ and $0$ otherwise: one bean in bowl $0$, empty bowls $1$ to $N-1$, $N$ beans in bowl $N$. That lone bean walks forward one bowl per move to bowl $N$, which then holds $N+1=n$ beans; the next move sows one into each bowl and the start is back, the pointer left on bowl $n-1$ — only the bowls have to look the same again.

This return is the first: restoring the all-ones arrangement needs a move sowing all $n$ beans, since a shorter sowing leaves its own bowl empty, and the lemma keeps every jump on a lap $l\lt N$ at or below $1+2^{k-1}$, under $n$, while on lap $N$ every move but the last sows a single bean.

9. Counting the moves

Lap $l$ makes one move at bowl $0$ plus one for each $i$ in $\{0,\dots,N-1\}$ with $l\not\preceq i$. Let $w(l)$ be the Hamming weight of $l$, its number of binary digits $1$. ^[3] For $l\le N-1$ a superset of $l$ comes from filling the other $k-w(l)$ digits freely, so there are $2^{\,k-w(l)}$ of them and the lap has $1+N-2^{\,k-w(l)}$ moves; for $l=N$ there is none, and the lap has $1+N=n$ moves. Summing,

$$M(n)=N(1+N)-\sum_{l=1}^{N-1}2^{\,k-w(l)}.$$

That sum is a product in disguise: running $l$ over $0,\dots,2^k-1$ chooses each of the $k$ digits independently, a digit $0$ contributing a factor $2$ and a digit $1$ a factor $1$ — the expansion of a product of sums that also proves the binomial theorem: ^[4]

$$\sum_{l=0}^{2^k-1}2^{\,k-w(l)}=\prod_{d=0}^{k-1}(2+1)=3^k.$$

Dropping the term $l=0$, which is $2^k$, leaves $3^k-2^k$, so

$$\begin{aligned}M(2^k+1)&=2^k(2^k+1)-(3^k-2^k)\\&=2^{k+1}+4^k-3^k.\end{aligned}$$

Counting bowl by bowl agrees: bowl $i+1$ is skipped on exactly the $2^{w(i)}-1$ laps with $l\preceq i$, and those skips also total $3^k-2^k$. For $n=9$ they are $0,1,1,3,1,3,3,7$, together $19=3^3-2^3$, while the eight laps take $5,5,7,5,7,7,8,9$ moves: $M(9)=53$.

10. From $M$ to the sum

Put $E=10^{18}$ and $S(E)=\sum_{k=0}^{E}M(2^k+1)$. Splitting the closed form termwise leaves three geometric series, each summed by $\sum_{k=0}^{E}r^k=(r^{E+1}-1)/(r-1)$ for $r\ne1$: ^[5]

$$\begin{aligned}S(E)&=\sum_{k=0}^{E}\left(2^{k+1}+4^k-3^k\right)\\&=2\left(2^{E+1}-1\right)+\frac{4^{E+1}-1}{3}\\&\quad-\frac{3^{E+1}-1}{2}.\end{aligned}$$

A small case checks the bookkeeping: for $E=5$ the right side is $2\cdot63+4095/3-728/2$${}=1127$, and the six terms $2+5+15+53+207+845$ do add to $1127$.

11. Working modulo $7^9$

The exact $S(10^{18})$ has about $6\cdot10^{17}$ decimal digits, hence the modulus $7^9=40353607$. Dividing modulo $m$ means multiplying by an inverse, a number whose product with the divisor leaves remainder $1$: since $3\cdot5\equiv1\pmod 7$, dividing by $3$ modulo $7$ is multiplying by $5$. An inverse exists exactly when divisor and modulus are coprime, ^[6] and $7^9$ is odd and not a multiple of $3$, so both exist: $2^{-1}\equiv(7^9+1)/2=20176804$ and $3^{-1}\equiv(2\cdot7^9+1)/3=26902405$. With $p_r=r^{E+1}\bmod 7^9$,

$$\begin{aligned}S(E)&\equiv 2(p_2-1)+(p_4-1)\,3^{-1}\\&\quad-(p_3-1)\,2^{-1}\pmod{7^9}.\end{aligned}$$

The residues $p_2,p_3,p_4$ stay below $7^9$; the work sits in the exponent. Exponentiation by squaring reads its binary digits, so that $r^{13}=r^{8}r^{4}r$, reducing after every product. ^[7] Here $10^{18}+1$ has $60$ binary digits, $25$ of them $1$, and the C++ and Java loops do $60$ squarings and $25$ accumulator multiplications per base, one of each more than the tightest form of the method needs; Python calls its own built-in. Two shortcuts go unused: $p_4\equiv p_2^{2}$, and Euler's theorem would reduce the exponent modulo $\varphi(7^9)=34588806$ first. ^[8]

How the Code Works

The C++ file ^[9] fixes kExponentLimit $=10^{18}$ and kModulus $=40353607$ and keeps checking apart from computing. simulate_M runs Section 1 on a vector of ones: it empties bowls[pos], advances pos once per bean while filling the bowl reached, and returns the move counter once a scan finds the all-ones position. closed_form_M_for_power_of_two_plus_one returns pow_u64(2, k + 1) + pow_u64(4, k) - pow_u64(3, k), called only for $k\le8$, where even $2^9+4^8=66048$ leaves room to spare in 64 bits and $4^k\ge3^k$ stops the unsigned subtraction underflowing.

run_validations lists twelve ValidationTask records, each a bowl count with the value expected: the statement's $M(5)=15$ and $M(100)=10920$, a check of $M(2)=2$, which the statement does not supply, and nine comparisons of simulate_M against the closed form for $k=0,\dots,8$, that is $n$ up to $257$. Mismatches go to standard error, main exits with code $1$, and --skip-checks jumps to the answer.

solve is Section 11, built from mod_inverse (extended Euclid on signed 64-bit integers, returning $0$ when no inverse exists), pow_mod (binary exponentiation) and mul_mod (factors widened to unsigned __int128), and combining its parts as (sum_2 + sum_4 + kModulus - sum_3) % kModulus, where the extra modulus keeps the subtraction from wrapping around. Python uses pow(base, exp, mod) and pow(a, -1, mod), which needs Python 3.8 or newer, and Java its own powMod with BigInteger.modInverse; Java's mulMod is the bare (a * b) % mod, safe since a product of residues stays below $1.63\cdot10^{15}$, inside a signed 64-bit long, so the C++ widening is a precaution. Neither port runs the checks.

Complexity and Verification

solve does three modular exponentiations with a 60-bit exponent and two extended Euclidean runs, one per inverse, then a constant number of modular operations: with the modulus fixed and every value in a machine word, $O(\log E)$ time and $O(1)$ memory. The checks dominate and are still trivial — the twelve tasks simulate $89499$ moves, the longest $M(257)=59487$, placing $101358$ beans and inspecting $95775$ bowl entries, each scan stopping at the first that differs from one; a simulation with $n$ bowls holds $n$ counts, so $O(n)$ memory and $O(nM(n))$ time. Simulating the sum is hopeless: $10^{18}+1$ terms, the last close to $4^{10^{18}}$ and just below it, since $3^E\gt2^{E+1}$ for $E\ge2$.

Beyond the program's own checks, these were verified for this article. Simulation agrees with $2^{k+1}+4^k-3^k$ for every $k\le12$, the last being $M(4097)=16253967$. For every $k\le7$ the lap picture of Sections 6 and 8 holds item by item — lap count, the bean in bowl $0$ at each start, bowls visited, the counts $H(i,l)$, lap lengths, and that only the final move sows $n$ beans — and the gap lemma holds for every $l\le129$ and $i\lt2100$ with $l\not\preceq i$. Section 2 recovered the unique predecessor of every reachable position for $n\le7$, $M(n)$ ignores the starting bowl for $n\le39$, and Section 11 matches the exact sum modulo $7^9$ for every $E\le200$. The three programs agree with each other and with an independent recomputation.

One gap is worth naming: the program tests the closed form only for $k\le8$ yet applies it up to $k=10^{18}$, and Sections 6 to 9 prove it for every $k$. Three cautions on the code as it stands: mod_inverse returns $0$ when no inverse exists and solve never inspects it, so an even kModulus, or a multiple of $3$, would silently give a wrong result; simulate_M would read past its vector for $0$ bowls, which no task asks for; and the task label for $M(2)$ credits the statement, which gives only $M(5)$ and $M(100)$.

Footnotes and References

The references supply the background facts used above; the lap decomposition, submask rule, gap lemma and move count are derived here.

Project Euler 335 — Gathering the Beans. The statement describes the circle of bowls with one bean each and the move that empties a bowl into those that follow clockwise, gives M(5) = 15 and M(100) = 10920, and asks for the sum of M(2^k + 1) over k up to 10^18, modulo 7^9.
Wikipedia — Mancala. Sowing drops a hole's seeds "one by one into subsequent holes in a motion wrapping around the board"; relay sowing, the last seed starting the next sowing from its hole, is a frequent feature.
Wikipedia — Hamming weight. The Hamming weight, or population count, of a set of bits is "the number of bits set to 1" — the function w of Section 9.
Wikipedia — Binomial theorem. It expands (x + y)^n as a product of n copies of (x + y), one term per choice in each factor — the argument behind Σ 2^(k − w(l)) = 3^k.
Wikipedia — Geometric series. The partial sum with first term a and ratio r is a(1 − r^(n+1))/(1 − r) for r ≠ 1; Section 10 uses ratios 2, 4 and 3.
Wikipedia — Modular multiplicative inverse. An inverse of a modulo m exists "if and only if gcd(a, m) = 1" and follows from the extended Euclidean algorithm's ax + my = 1, as mod_inverse computes it.
Wikipedia — Exponentiation by squaring. The method needs about log₂ n squarings and at most as many further multiplications, one fewer than the exponent's number of 1 bits.
Wikipedia — Euler's totient function. It gives φ(p^k) = p^k − p^(k−1) and Euler's theorem, a^φ(n) ≡ 1 (mod n) for coprime a, so φ(7^9) = 34588806.
C++ — Euler335.cpp. The linked immutable revision holds simulate_M, run_validations and solve, with the comment M(2^k + 1) = 2^(k+1) + 4^k − 3^k.

Problem 335 source code

C++

#include <algorithm>
#include <atomic>
#include <cstdint>
#include <iostream>
#include <limits>
#include <mutex>
#include <sstream>
#include <string>
#include <thread>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr u64 kExponentLimit = 1'000'000'000'000'000'000ULL;
constexpr u64 kModulus = 40'353'607ULL;  // 7^9

struct Options {
    bool run_checks = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

struct ValidationTask {
    u64 n = 0;
    u64 expected = 0;
    std::string label;
};

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0U) return false;

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) return false;

    u64 parsed = 0;
    for (const char c : tail) {
        if (c < '0' || c > '9') return false;
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) return false;
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u64 parsed = 0;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) return false;
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) return false;
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(const int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checks") {
            options.run_checks = false;
            continue;
        }
        if (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return true;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const std::size_t workload_units) {
    if (!allow_multithreading || workload_units < 2ULL) return 1U;

    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) threads = 1U;
    }

    return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload_units)));
}

u64 pow_u64(u64 base, unsigned exp) {
    u64 out = 1ULL;
    for (unsigned i = 0U; i < exp; ++i) out *= base;
    return out;
}

u64 simulate_M(const u64 bowls_count) {
    std::vector<u64> bowls(static_cast<std::size_t>(bowls_count), 1ULL);
    u64 pos = 0ULL;

    for (u64 steps = 1ULL;; ++steps) {
        const u64 beans = bowls[static_cast<std::size_t>(pos)];
        bowls[static_cast<std::size_t>(pos)] = 0ULL;

        for (u64 i = 0ULL; i < beans; ++i) {
            ++pos;
            if (pos == bowls_count) pos = 0ULL;
            ++bowls[static_cast<std::size_t>(pos)];
        }

        bool all_ones = true;
        for (const u64 value : bowls) {
            if (value != 1ULL) {
                all_ones = false;
                break;
            }
        }
        if (all_ones) return steps;
    }
}

u64 closed_form_M_for_power_of_two_plus_one(unsigned k) {
    // M(2^k + 1) = 2^(k+1) + 4^k - 3^k
    return pow_u64(2ULL, k + 1U) + pow_u64(4ULL, k) - pow_u64(3ULL, k);
}

bool run_validations(const bool allow_multithreading, const unsigned requested_threads) {
    std::vector<ValidationTask> tasks;
    tasks.reserve(12);

    tasks.push_back(ValidationTask{2ULL, 2ULL, "statement: M(2)"});
    tasks.push_back(ValidationTask{5ULL, 15ULL, "statement: M(5)"});
    tasks.push_back(ValidationTask{100ULL, 10'920ULL, "statement: M(100)"});

    for (unsigned k = 0U; k <= 8U; ++k) {
        const u64 n = (1ULL << k) + 1ULL;
        tasks.push_back(ValidationTask{
            n,
            closed_form_M_for_power_of_two_plus_one(k),
            "closed form at k=" + std::to_string(k),
        });
    }

    const unsigned thread_count =
        choose_thread_count(allow_multithreading, requested_threads, tasks.size());

    std::atomic<std::size_t> next_task{0ULL};
    std::vector<std::string> failures;
    std::mutex failure_mutex;
    std::vector<std::thread> pool;
    pool.reserve(thread_count);

    for (unsigned t = 0U; t < thread_count; ++t) {
        pool.emplace_back([&]() {
            while (true) {
                const std::size_t id = next_task.fetch_add(1ULL, std::memory_order_relaxed);
                if (id >= tasks.size()) break;

                const ValidationTask& task = tasks[id];
                const u64 got = simulate_M(task.n);
                if (got == task.expected) continue;

                std::ostringstream oss;
                oss << "Validation failed (" << task.label << "): M(" << task.n << ") = " << got
                    << ", expected " << task.expected;

                std::lock_guard<std::mutex> lock(failure_mutex);
                failures.push_back(oss.str());
            }
        });
    }

    for (std::thread& th : pool) th.join();

    if (failures.empty()) return true;
    std::sort(failures.begin(), failures.end());
    for (const std::string& line : failures) std::cerr << line << '\n';
    return false;
}

u64 mul_mod(const u64 a, const u64 b, const u64 mod) {
    return static_cast<u64>((static_cast<u128>(a) * static_cast<u128>(b)) % mod);
}

u64 pow_mod(u64 base, u64 exp, const u64 mod) {
    u64 out = 1ULL % mod;
    base %= mod;

    while (exp > 0ULL) {
        if ((exp & 1ULL) != 0ULL) out = mul_mod(out, base, mod);
        base = mul_mod(base, base, mod);
        exp >>= 1ULL;
    }

    return out;
}

u64 mod_inverse(const u64 a, const u64 mod) {
    std::int64_t t = 0;
    std::int64_t new_t = 1;
    std::int64_t r = static_cast<std::int64_t>(mod);
    std::int64_t new_r = static_cast<std::int64_t>(a % mod);

    while (new_r != 0) {
        const std::int64_t q = r / new_r;
        const std::int64_t temp_t = t - q * new_t;
        t = new_t;
        new_t = temp_t;

        const std::int64_t temp_r = r - q * new_r;
        r = new_r;
        new_r = temp_r;
    }

    if (r != 1) return 0;
    if (t < 0) t += static_cast<std::int64_t>(mod);
    return static_cast<u64>(t);
}

u64 solve() {
    const u64 inv2 = mod_inverse(2ULL, kModulus);
    const u64 inv3 = mod_inverse(3ULL, kModulus);

    const u64 p2 = pow_mod(2ULL, kExponentLimit + 1ULL, kModulus);
    const u64 p3 = pow_mod(3ULL, kExponentLimit + 1ULL, kModulus);
    const u64 p4 = pow_mod(4ULL, kExponentLimit + 1ULL, kModulus);

    const u64 sum_2 = mul_mod(2ULL, (p2 + kModulus - 1ULL) % kModulus, kModulus);
    const u64 sum_4 = mul_mod((p4 + kModulus - 1ULL) % kModulus, inv3, kModulus);
    const u64 sum_3 = mul_mod((p3 + kModulus - 1ULL) % kModulus, inv2, kModulus);

    return (sum_2 + sum_4 + kModulus - sum_3) % kModulus;
}

}  // namespace

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

    if (options.run_checks && !run_validations(options.allow_multithreading, options.requested_threads)) {
        return 1;
    }

    std::cout << solve() << '\n';
    return 0;
}

Python

def mul_mod(a, b, mod):
    return (a * b) % mod

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

def mod_inverse(a, mod):
    return pow(a, -1, mod)

def solve():
    kExponentLimit = 1000000000000000000
    kModulus = 40353607
    
    inv2 = mod_inverse(2, kModulus)
    inv3 = mod_inverse(3, kModulus)
    
    p2 = pow_mod(2, kExponentLimit + 1, kModulus)
    p3 = pow_mod(3, kExponentLimit + 1, kModulus)
    p4 = pow_mod(4, kExponentLimit + 1, kModulus)
    
    sum_2 = mul_mod(2, (p2 + kModulus - 1) % kModulus, kModulus)
    sum_4 = mul_mod((p4 + kModulus - 1) % kModulus, inv3, kModulus)
    sum_3 = mul_mod((p3 + kModulus - 1) % kModulus, inv2, kModulus)
    
    ans = (sum_2 + sum_4 + kModulus - sum_3) % kModulus
    return str(ans)

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

Java

import java.util.*;

public class Euler335 {

    static final long EXPONENT_LIMIT = 1000000000000000000L;
    static final long MODULUS = 40353607L;

    static long mulMod(long a, long b, long mod) {
        return (a * b) % mod;
    }

    static long powMod(long base, long exp, long mod) {
        long out = 1L % mod;
        base %= mod;
        while (exp > 0) {
            if ((exp & 1) != 0)
                out = mulMod(out, base, mod);
            base = mulMod(base, base, mod);
            exp >>= 1;
        }
        return out;
    }

    static long modInverse(long a, long mod) {
        return java.math.BigInteger.valueOf(a).modInverse(java.math.BigInteger.valueOf(mod)).longValue();
    }

    public static String solve() {
        long inv2 = modInverse(2, MODULUS);
        long inv3 = modInverse(3, MODULUS);

        long p2 = powMod(2, EXPONENT_LIMIT + 1, MODULUS);
        long p3 = powMod(3, EXPONENT_LIMIT + 1, MODULUS);
        long p4 = powMod(4, EXPONENT_LIMIT + 1, MODULUS);

        long sum_2 = mulMod(2, (p2 + MODULUS - 1) % MODULUS, MODULUS);
        long sum_4 = mulMod((p4 + MODULUS - 1) % MODULUS, inv3, MODULUS);
        long sum_3 = mulMod((p3 + MODULUS - 1) % MODULUS, inv2, MODULUS);

        long ans = (sum_2 + sum_4 + MODULUS - sum_3) % MODULUS;
        return String.valueOf(ans);
    }

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