Problem 334: Spilling the Beans

View on Project Euler

Project Euler Problem 334 Solution

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

Problem Summary The bowls start with a deterministic pseudo-random number of beans. A legal move chooses a bowl containing at least two beans, removes two beans from it, and places one bean in each neighboring bowl. The task is to count how many such spills occur before every bowl contains at most one bean. Mathematical Approach Index the bowls by integers and let \(b_i\) be the number of beans in bowl \(i\). Although the input generator only produces finitely many nonzero bowls, the process is best viewed on the whole integer line, because beans may move left of the original range. State Variables and Invariants The code tracks three global quantities: $$B=\sum_i b_i,\qquad M_1=\sum_i i\,b_i,\qquad M_2=\sum_i i^2 b_i.$$ Here \(B\) is the total number of beans, \(M_1\) is the first moment, and \(M_2\) is the second moment. One spill at position \(i\) changes the state by $$b_i\mapsto b_i-2,\qquad b_{i-1}\mapsto b_{i-1}+1,\qquad b_{i+1}\mapsto b_{i+1}+1.$$ The bean count is preserved because \(-2+1+1=0\). The first moment is also preserved because $$(-2)\,i + (i-1) + (i+1) = 0.$$ Therefore every reachable configuration has the same values of \(B\) and \(M_1\) as the initial one....

Detailed mathematical approach

Problem Summary

The bowls start with a deterministic pseudo-random number of beans. A legal move chooses a bowl containing at least two beans, removes two beans from it, and places one bean in each neighboring bowl. The task is to count how many such spills occur before every bowl contains at most one bean.

Mathematical Approach

Index the bowls by integers and let \(b_i\) be the number of beans in bowl \(i\). Although the input generator only produces finitely many nonzero bowls, the process is best viewed on the whole integer line, because beans may move left of the original range.

State Variables and Invariants

The code tracks three global quantities:

$$B=\sum_i b_i,\qquad M_1=\sum_i i\,b_i,\qquad M_2=\sum_i i^2 b_i.$$

Here \(B\) is the total number of beans, \(M_1\) is the first moment, and \(M_2\) is the second moment. One spill at position \(i\) changes the state by

$$b_i\mapsto b_i-2,\qquad b_{i-1}\mapsto b_{i-1}+1,\qquad b_{i+1}\mapsto b_{i+1}+1.$$

The bean count is preserved because \(-2+1+1=0\). The first moment is also preserved because

$$(-2)\,i + (i-1) + (i+1) = 0.$$

Therefore every reachable configuration has the same values of \(B\) and \(M_1\) as the initial one.

Why Each Move Adds Exactly 2 to the Second Moment

The decisive simplification is that every legal spill changes \(M_2\) by the same amount:

$$\Delta M_2=(i-1)^2+(i+1)^2-2i^2=2.$$

Hence the total number of moves depends only on the initial and final second moments:

$$\#\text{moves}=\frac{M_2^{\text{final}}-M_2^{\text{initial}}}{2}.$$

So the problem is no longer “simulate all spills”, but “identify the stabilized state and compute its second moment”.

What a Stable Configuration Looks Like

A bowl is unstable exactly when it contains at least two beans. Therefore every stabilized state satisfies \(b_i\in\{0,1\}\) for all \(i\). The final state is thus a set of \(B\) occupied integer positions.

Write these occupied positions as

$$x_0<x_1<\cdots<x_{B-1}.$$

Because the first moment is invariant, they must satisfy

$$x_0+x_1+\cdots+x_{B-1}=M_1.$$

Among all such \(0/1\)-configurations, stabilization for this one-dimensional spill rule leads to the one with smallest possible second moment \(\sum_j x_j^2\). That is the quantity reconstructed by the code.

Discrete Convexity and the One-Gap Structure

Subtract the baseline ordering by defining \(x_j=j+y_j\). Since the \(x_j\) are strictly increasing integers, the sequence \(y_0,\dots,y_{B-1}\) is nondecreasing. Its total is

$$\sum_{j=0}^{B-1} y_j = M_1-\sum_{j=0}^{B-1}j = M_1-\frac{B(B-1)}{2}=: \text{target\_shift}.$$

Now minimize \(\sum_j (j+y_j)^2\) subject to this fixed sum. Because the square function is convex, the minimum is achieved when the \(y_j\) are as equal as possible. Therefore each optimal \(y_j\) must be either \(q\) or \(q+1\), where

$$q=\left\lfloor\frac{\text{target\_shift}}{B}\right\rfloor,\qquad r=\text{target\_shift}-qB,\qquad 0\le r<B.$$

So \(B-r\) of the adjusted values equal \(q\), and the remaining \(r\) equal \(q+1\). Translating back gives

$$x_j=j+q\quad (0\le j<B-r),\qquad x_j=j+q+1\quad (B-r\le j<B).$$

Equivalently, the occupied bowls are

$$q,\ q+1,\ \dots,\ q+(B-r)-1,\ q+(B-r)+1,\ \dots,\ q+B,$$

namely one almost-consecutive block with exactly one missing position. This is the final pattern encoded by the implementation.

Because \(\text{target\_shift}\) can be negative, the C++ solution uses a custom floor-division helper instead of relying on truncation toward zero.

Closed Form for the Final Second Moment

The solver does not enumerate final positions one by one. It uses the arithmetic square-sum identity

$$\sum_{k=0}^{n-1}(a+k)^2 = n a^2 + 2a\frac{n(n-1)}{2} + \frac{n(n-1)(2n-1)}{6}.$$

Let

$$n_1=B-r,\qquad n_2=r.$$

Then the two occupied blocks start at \(q\) and \(q+n_1+1\), so

$$M_2^{\text{final}}=\sum_{k=0}^{n_1-1}(q+k)^2+\sum_{k=0}^{n_2-1}(q+n_1+1+k)^2.$$

This is exactly what the helper sum_squares_arith(first, n) computes.

Checkpoint Example: \([2,3]\)

The source includes a small verification case with two bowls containing \([2,3]\). Using indices \(0,1\):

$$B=5,\qquad M_1=0\cdot2+1\cdot3=3,\qquad M_2=0^2\cdot2+1^2\cdot3=3.$$

Then

$$\text{target\_shift}=3-\frac{5\cdot4}{2}=-7,\qquad q=\left\lfloor\frac{-7}{5}\right\rfloor=-2,\qquad r=3.$$

Hence the final occupied positions are \(-2,-1,1,2,3\), whose second moment is

$$(-2)^2+(-1)^2+1^2+2^2+3^2=19.$$

Therefore

$$\#\text{moves}=\frac{19-3}{2}=8,$$

which matches the checkpoint in the code.

How the Code Works

The implementation has three compact stages. First, generate_beans(count) creates the deterministic initial configuration. Second, moves_needed(beans) scans the input once and accumulates \(B\), \(M_1\), and \(M_2\). Third, it reconstructs the optimal stabilized pattern from \((q,r)\), evaluates \(M_2^{\text{final}}\) with two arithmetic square sums, and returns \((M_2^{\text{final}}-M_2^{\text{initial}})/2\).

No step-by-step spill simulation is performed anywhere.

Complexity Analysis

If the generator produces \(N\) initial bowls, the moment accumulation costs \(O(N)\) time. Everything after that is \(O(1)\). Extra memory is \(O(1)\) beyond the input container. This is dramatically better than simulating every spill, because the actual number of moves is extremely large.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=334
  2. Moment (mathematics): Wikipedia
  3. Convex function: Wikipedia
  4. Square pyramidal number / sum of squares: Wikipedia
  5. Abelian sandpile model: Wikipedia

Problem 334 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

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

std::string to_string_u128(u128 value) {
    if (value == 0) {
        return "0";
    }
    std::string s;
    while (value > 0) {
        const int digit = static_cast<int>(value % 10);
        s.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

std::vector<u64> generate_beans(const int count) {
    std::vector<u64> beans;
    beans.reserve(static_cast<std::size_t>(count));

    u64 t = 123456ULL;
    for (int i = 0; i < count; ++i) {
        if ((t & 1ULL) == 0ULL) {
            t /= 2ULL;
        } else {
            t = (t / 2ULL) ^ 926252ULL;
        }
        const u64 b = (t & ((1ULL << 11) - 1ULL)) + 1ULL;
        beans.push_back(b);
    }

    return beans;
}

i128 floor_div(const i128 a, const i128 b) {
    i128 q = a / b;
    i128 r = a % b;
    if (r < 0) {
        q -= 1;
    }
    return q;
}

i128 sum_squares_arith(const i128 first, const i128 n) {
    // Sum_{k=0}^{n-1} (first + k)^2
    if (n <= 0) {
        return 0;
    }
    const i128 sum_k = n * (n - 1) / 2;
    const i128 sum_k2 = n * (n - 1) * (2 * n - 1) / 6;
    return n * first * first + 2 * first * sum_k + sum_k2;
}

u128 moves_needed(const std::vector<u64>& beans) {
    i128 total = 0;
    i128 first_moment = 0;
    i128 second_moment = 0;

    for (i128 i = 0; i < static_cast<i128>(beans.size()); ++i) {
        const i128 b = static_cast<i128>(beans[static_cast<std::size_t>(i)]);
        total += b;
        first_moment += i * b;
        second_moment += i * i * b;
    }

    // In every move, second moment increases by exactly 2.
    // Stable final states are 0/1-occupancy sets with same total and first moment;
    // the attained stabilization is the minimum-second-moment such set.
    const i128 B = total;
    const i128 target_shift = first_moment - B * (B - 1) / 2;

    const i128 q = floor_div(target_shift, B);
    const i128 r = target_shift - q * B;  // 0 <= r < B

    const i128 n1 = B - r;
    const i128 n2 = r;

    // Final occupied positions:
    // q, q+1, ..., q+n1-1, q+n1+1, ..., q+B
    const i128 q_final = sum_squares_arith(q, n1) + sum_squares_arith(q + n1 + 1, n2);

    const i128 diff = q_final - second_moment;
    return static_cast<u128>(diff / 2);
}

bool run_checkpoints() {
    if (moves_needed(std::vector<u64>{2ULL, 3ULL}) != static_cast<u128>(8ULL)) {
        std::cerr << "Checkpoint failed: [2,3] should take 8 moves\n";
        return false;
    }

    const std::vector<u64> b2 = generate_beans(2);
    if (b2.size() != 2 || b2[0] != 289ULL || b2[1] != 145ULL) {
        std::cerr << "Checkpoint failed: b1/b2 generation\n";
        return false;
    }

    if (moves_needed(b2) != static_cast<u128>(3419100ULL)) {
        std::cerr << "Checkpoint failed: two-bowl sample moves\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        return 2;
    }

    const std::vector<u64> beans = generate_beans(1500);
    const u128 answer = moves_needed(beans);
    std::cout << to_string_u128(answer) << '\n';
    return 0;
}

Python

def generate_beans(count):
    beans = []
    t = 123456
    for i in range(count):
        if (t & 1) == 0:
            t //= 2
        else:
            t = (t // 2) ^ 926252
        b = (t & ((1 << 11) - 1)) + 1
        beans.append(b)
    return beans

def sum_squares_arith(first, n):
    if n <= 0:
        return 0
    sum_k = n * (n - 1) // 2
    sum_k2 = n * (n - 1) * (2 * n - 1) // 6
    return n * first * first + 2 * first * sum_k + sum_k2

def moves_needed(beans):
    total = 0
    first_moment = 0
    second_moment = 0

    for i, b in enumerate(beans):
        total += b
        first_moment += i * b
        second_moment += i * i * b

    B = total
    target_shift = first_moment - B * (B - 1) // 2

    q = target_shift // B
    r = target_shift % B

    n1 = B - r
    n2 = r

    q_final = sum_squares_arith(q, n1) + sum_squares_arith(q + n1 + 1, n2)

    diff = q_final - second_moment
    return diff // 2

def solve():
    beans = generate_beans(1500)
    ans = moves_needed(beans)
    return str(ans)

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

Java

import java.util.*;
import java.math.BigInteger;

public class Euler334 {

    static List<Long> generateBeans(int count) {
        List<Long> beans = new ArrayList<>(count);
        long t = 123456;
        for (int i = 0; i < count; i++) {
            if ((t & 1) == 0) {
                t /= 2;
            } else {
                t = (t / 2) ^ 926252L;
            }
            long b = (t & ((1L << 11) - 1)) + 1;
            beans.add(b);
        }
        return beans;
    }

    static BigInteger sumSquaresArith(BigInteger first, BigInteger n) {
        if (n.compareTo(BigInteger.ZERO) <= 0)
            return BigInteger.ZERO;
        BigInteger two = BigInteger.valueOf(2);
        BigInteger six = BigInteger.valueOf(6);

        BigInteger nMinus1 = n.subtract(BigInteger.ONE);
        BigInteger sumK = n.multiply(nMinus1).divide(two);

        BigInteger twoNMinus1 = n.multiply(two).subtract(BigInteger.ONE);
        BigInteger sumK2 = n.multiply(nMinus1).multiply(twoNMinus1).divide(six);

        BigInteger term1 = n.multiply(first).multiply(first);
        BigInteger term2 = two.multiply(first).multiply(sumK);

        return term1.add(term2).add(sumK2);
    }

    static BigInteger floorDiv(BigInteger a, BigInteger b) {
        return a.divide(b); // In Java, BigInteger.divide truncates toward zero. But since both are positive
                            // here, it is actually floor division.
        // Wait, for negative 'a', it truncates toward zero instead of floor.
        // Let's implement proper floor div:
    }

    static BigInteger properFloorDiv(BigInteger a, BigInteger b) {
        BigInteger[] qr = a.divideAndRemainder(b);
        BigInteger q = qr[0];
        BigInteger r = qr[1];
        if (r.compareTo(BigInteger.ZERO) < 0) {
            if (b.compareTo(BigInteger.ZERO) > 0)
                q = q.subtract(BigInteger.ONE);
            else
                q = q.add(BigInteger.ONE);
        }
        return q;
    }

    static BigInteger properFloorMod(BigInteger a, BigInteger b) {
        BigInteger[] qr = a.divideAndRemainder(b);
        BigInteger r = qr[1];
        if (r.compareTo(BigInteger.ZERO) < 0) {
            if (b.compareTo(BigInteger.ZERO) > 0)
                r = r.add(b);
            else
                r = r.subtract(b);
        }
        return r;
    }

    public static String solve() {
        List<Long> beans = generateBeans(1500);

        BigInteger total = BigInteger.ZERO;
        BigInteger firstMoment = BigInteger.ZERO;
        BigInteger secondMoment = BigInteger.ZERO;

        for (int i = 0; i < beans.size(); i++) {
            BigInteger b = BigInteger.valueOf(beans.get(i));
            BigInteger ib = BigInteger.valueOf(i);

            total = total.add(b);
            firstMoment = firstMoment.add(ib.multiply(b));
            secondMoment = secondMoment.add(ib.multiply(ib).multiply(b));
        }

        BigInteger B = total;
        BigInteger bMinus1 = B.subtract(BigInteger.ONE);
        BigInteger bPair = B.multiply(bMinus1).divide(BigInteger.valueOf(2));

        BigInteger targetShift = firstMoment.subtract(bPair);

        BigInteger q = properFloorDiv(targetShift, B);
        BigInteger r = properFloorMod(targetShift, B);

        BigInteger n1 = B.subtract(r);
        BigInteger n2 = r;

        BigInteger qFinal = sumSquaresArith(q, n1).add(sumSquaresArith(q.add(n1).add(BigInteger.ONE), n2));

        BigInteger diff = qFinal.subtract(secondMoment);
        return diff.divide(BigInteger.valueOf(2)).toString();
    }

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