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

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)$.

The Key Closed Form for $n=2^k+1$

The nontrivial combinatorial ingredient used by the solver is the identity

$$\boxed{M(2^k+1)=2^{k+1}+4^k-3^k\qquad(k\ge 0).}$$

The code treats this as the central lemma and verifies it against explicit simulation for small $k$. The first values are:

$$M(2)=2,\qquad M(3)=5,\qquad M(5)=15,\qquad M(9)=53,$$

which match

$$2^{1}+4^0-3^0=2,\qquad 2^{2}+4^1-3^1=5,\qquad 2^{3}+4^2-3^2=15,\qquad 2^{4}+4^3-3^3=53.$$

The implementation also checks the statement values $M(5)=15$ and $M(100)=10920$ by direct simulation before computing the final answer.

Reducing the Required Sum

Set

$$E=10^{18},\qquad S(E)=\sum_{k=0}^{E} M(2^k+1).$$

Substituting the closed form and splitting the sum termwise gives

$$S(E)=\sum_{k=0}^{E}\left(2^{k+1}+4^k-3^k\right) =\sum_{k=0}^{E}2^{k+1}+\sum_{k=0}^{E}4^k-\sum_{k=0}^{E}3^k.$$

So the original process problem is reduced to three standard geometric sums.

Evaluating the Geometric Series

For any ratio $r\neq 1$,

$$\sum_{k=0}^{E} r^k=\frac{r^{E+1}-1}{r-1}.$$

Applying this identity to the three terms:

$$\sum_{k=0}^{E}2^{k+1}=2\sum_{k=0}^{E}2^k=2(2^{E+1}-1),$$

$$\sum_{k=0}^{E}4^k=\frac{4^{E+1}-1}{3},\qquad \sum_{k=0}^{E}3^k=\frac{3^{E+1}-1}{2}.$$

Therefore

$$\boxed{S(E)=2(2^{E+1}-1)+\frac{4^{E+1}-1}{3}-\frac{3^{E+1}-1}{2}.}$$

This is the exact closed form the code evaluates, with $E=10^{18}$.

Modular Evaluation

The required modulus is

$$m=7^9=40{,}353{,}607.$$

Because $m$ is odd and not divisible by $3$, both $2$ and $3$ are invertible modulo $m$. Thus the divisions in the geometric-series formula can be replaced by multiplication with modular inverses:

$$S(E)\equiv 2(2^{E+1}-1)+(4^{E+1}-1)\cdot 3^{-1}-(3^{E+1}-1)\cdot 2^{-1}\pmod m.$$

The only large quantities left are $2^{E+1}$, $3^{E+1}$, and $4^{E+1}$ modulo $m$. These are computed by binary exponentiation, so the huge exponent $10^{18}+1$ is handled in $O(\log E)$ multiplications rather than $O(E)$.

How the Code Works

The C++ file has two clearly separated parts. First, simulate_M and run_validations perform checkpoint verification on small inputs, including the statement values and the formula $M(2^k+1)=2^{k+1}+4^k-3^k$ for $k=0,\dots,8$. Second, solve() computes the modular inverses of $2$ and $3$, evaluates

$$2^{E+1}\bmod m,\qquad 3^{E+1}\bmod m,\qquad 4^{E+1}\bmod m,$$

forms the three partial sums sum_2, sum_4, and sum_3, and combines them as (sum_2 + sum_4 - sum_3) mod m. The Python version mirrors the same mathematics using pow(base, exp, mod) and Python's built-in modular inverse. The Java version uses the same repeated squaring logic and obtains inverses through BigInteger.modInverse.

Complexity Analysis

The final computation needs only three modular exponentiations and a constant number of modular arithmetic operations, so it runs in $O(\log E)$ time and $O(1)$ memory. Directly simulating the bean process up to $10^{18}$ terms would be hopelessly too slow; simulation appears only in the fixed-size validation phase.

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());
    }
}