Problem 339: Peredur Fab Efrawg

Detailed mathematical approach

Problem Summary

A river separates two flocks, each starting with \(n\) sheep: one white, one black. Every sheep is equally likely to bleat next, and a bleat makes one sheep of the other flock cross the river and change colour, so a white bleat moves the position \((w,b)\) of white and black sheep to \((w+1,b-1)\) and a black bleat to \((w-1,b+1)\), and the flock never changes size. After each bleat and crossing, Peredur may remove any number of white sheep, to maximise the expected final number of black sheep. Once one colour is gone nothing can change, and the final count is the remaining flock if the white sheep are gone, \(0\) otherwise. Let \(E(n)\) be that count in expectation under an optimal strategy. The statement gives \(E(5)=6.871346\) and asks for \(E(10\,000)\), both to six decimals. ^[1] Its opening passage comes from the Welsh Mabinogion. ^[2]

The whole solution in one sentence: whenever both colours remain it is optimal to remove white sheep until they are strictly outnumbered by the black ones, and the expectation then follows from a sweep over the flock size \(s=1,\dots,2n\) in which each \(s\) costs one small system of linear equations linking neighbouring white counts.

Mathematical Approach

1. The moves, the two ends and the drift

Let \(s=w+b\) be the flock size. In \((w,b)\) the next bleat is white with probability \(w/s\) and black with probability \(b/s\); neither bleat changes \(s\), since sheep only change colour. Removing \(k\) white sheep is the only move that changes \(s\), to \(s-k\).

The smallest case is settled by hand: from \((1,1)\) the first bleat leads to \((2,0)\) or \((0,2)\), each with probability \(\tfrac12\), and both are final, so \(E(1)=\tfrac12(0+2)=1\); removals cannot help, since Peredur may act only after a bleat.

One line already explains the strategy: the expected change of the number of black sheep in one bleat is

$$\frac{w}{s}(-1)+\frac{b}{s}(+1)=\frac{b-w}{s},$$

so the black flock grows in expectation exactly while it is in the majority: the drift is zero when the colours are level and negative when white sheep outnumber black ones. White sheep are not worthless — each may still cross and become black — but they stop paying once they catch up.

2. The Bellman equation

Let \(V(w,b)\) be the largest expected final number of black sheep from \((w,b)\) when a removal is allowed right now. Keeping \(k\ge1\) white sheep has the one-step value

$$C(k,b)=\frac{k}{k+b}V(k+1,b-1)+\frac{b}{k+b}V(k-1,b+1),$$

while \(k=0\) removes every white sheep and ends the game at once with \(b\) black sheep. Hence

$$V(w,b)=\max\Bigl(b,\ \max_{1\le k\le w}C(k,b)\Bigr),$$

with \(V(w,0)=0\) and \(V(0,b)=b\). The inner maximum over \(k\le w\) contains the one over \(k\le w-1\), so this is the same as

$$V(w,b)=\max\bigl(V(w-1,b),\,C(w,b)\bigr),$$

which reads: keep the white flock for one more bleat, or throw one sheep away and ask again. Everything is finite because the game ends almost surely: at most \(2n\) removals are possible, and between two of them the white count performs a random walk on \(\{0,1,\dots,s\}\), which reaches an absorbing end with probability \(1\).

3. Layers of constant flock size

Since a bleat freezes \(s\), group the positions by flock size: put \(U_s(w)=V(w,s-w)\), so \(U_s(0)=s\) and \(U_s(s)=0\). With

$$C_s(w)=\frac{w}{s}U_s(w+1)+\frac{s-w}{s}U_s(w-1)$$

the equation of Section 2 becomes, for \(1\le w\le s-1\),

$$U_s(w)=\max\bigl(U_{s-1}(w-1),\,C_s(w)\bigr).$$

This is an obstacle problem: inside the layer \(s\) we face an optimal stopping problem for the walk, where "stopping" means removing one white sheep and dropping into the layer \(s-1\) for the reward \(U_{s-1}(w-1)\). Its solution is the Snell envelope: the smallest process that dominates the reward and whose expected next value never exceeds its current one, which is what the equation says. ^[6] A layer needs only the one below it, so they can be built in the order \(s=1,2,3,\dots\).

4. If nothing is removed: a gambler's ruin with binomial weights

Suppose no sheep is ever removed, and write \(P_s(w)\) for the expected final count then. Conditioning on the next bleat gives \(P_s(0)=s\), \(P_s(s)=0\) and

$$P_s(w)=\frac{w}{s}P_s(w+1)+\frac{s-w}{s}P_s(w-1).$$

Multiplying by \(s\) and rearranging gives a statement about the differences \(D(w)=P_s(w)-P_s(w-1)\):

$$w\,D(w+1)=(s-w)\,D(w),$$

so \(D(w)=D(1)\binom{s-1}{w-1}\) by induction — the ratios \((s-1)/1,(s-2)/2,\dots\) build up a binomial coefficient. Summing the differences gives the general solution

$$P_s(w)=s+D(1)\sum_{j=0}^{w-1}\binom{s-1}{j},$$

and \(P_s(s)=0\) forces \(D(1)=-s/2^{\,s-1}\). The value of an untouched flock is therefore

$$P_s(w)=\frac{s}{2^{\,s-1}}\sum_{i=w}^{s-1}\binom{s-1}{i},$$

which is \(s\) times the probability that a binomial variable with \(s-1\) trials and probability \(\tfrac12\) is at least \(w\): the classical gambler's-ruin computation with state-dependent step probabilities. ^[5] For instance \(P_5(2)=5\cdot\tfrac{11}{16}=3.4375\), against \(\tfrac{125}{33}\approx3.787879\) for the optimal play of that position \((2,3)\). Most striking is the balanced start: half of the coefficients \(\binom{2n-1}{i}\) lie above the middle, so \(P_{2n}(n)=n\). Doing nothing yields exactly half the flock, for every \(n\); this formula and that special case also appear in the literature. ^[4] Everything above \(n\) must be earned.

5. The key claim: keep the white sheep in the minority

Claim. In a non-terminal position \((w,b)\), so with \(b\ge1\), it is optimal to reduce the white flock to \(\min(w,b-1)\): do nothing while \(w\le b-1\), otherwise remove white sheep until one fewer white than black remains.

This is the optimal policy of the Mabinogion sheep problem, which David Williams introduced and proved optimal by martingale methods in Probability with Martingales (1991); for a general flock that theorem is what we rely on. ^{[3] [4]} The drift of Section 1 is only the intuition: \(w=b-1\) is the largest white flock compatible with a black lead, and removing more wastes sheep that might still cross over.

Any finite set of layers can also be certified on its own, because the Bellman equation has only one solution.

Lemma. Fix \(s\), suppose layer \(s-1\) is known, and let \(C_h\) be built from a function \(h\) as \(C_s\) is from \(U_s\),

$$C_h(w)=\frac{w}{s}h(w+1)+\frac{s-w}{s}h(w-1).$$

If \(h(0)=s\), \(h(s)=0\) and \(h(w)=\max(U_{s-1}(w-1),C_h(w))\) for \(1\le w\le s-1\), then \(h=U_s\).

Proof. Let \(m\) be the largest value of \(|h-U_s|\), attained at some \(w\). If \(1\le w\le s-1\) then, since \(x\mapsto\max(c,x)\) never increases distances,

$$\begin{aligned}m&\le\frac{w}{s}\bigl|h-U_s\bigr|(w+1)\\&\quad+\frac{s-w}{s}\bigl|h-U_s\bigr|(w-1).\end{aligned}$$

Both weights are positive and both differences are at most \(m\), so \(m\) is attained at \(w-1\) too; repeating walks down to \(w=0\), where the difference is \(0\), so \(m=0\). \(\square\)

Checking that equation in exact rational arithmetic succeeds at every position of every flock with \(s\le160\), which proves the claim for those flocks: waiting is strictly better throughout \(w\le b-1\), removing throughout \(w\ge b\), the only tie being \((1,1)\). The last section reports the same test in floating point at full size — evidence, not an extension of the proof.

6. A consequence: the black flock never dies out

For \(n\ge2\) the first crossing leaves a black sheep, and after that removal and every later one the position has \(w\le b-1\), which for \(w\ge1\) forces \(b\ge2\): a white bleat still leaves a black sheep, a black bleat adds one. So the black sheep never run out, the game ends at \(w=0\), and every remaining sheep is black. The final count is the final flock size:

$$E(n)=2n-R(n),$$

where \(R(n)\) is the expected number of sheep removed under this policy. Among policies certain to end all black, maximising black sheep is the same as minimising removals; a policy that can lose every black sheep does not satisfy the identity, and removing nothing gives only \(n\). For \(n=2\) one gets \(E(2)=\tfrac{55}{24}\) and \(R(2)=\tfrac{41}{24}\approx1.708333\).

7. The recurrence the program solves

Fix a layer \(s\) and let \(t=\lfloor(s-1)/2\rfloor\), the largest \(w\) with \(w\le s-w-1\): the largest white flock still in the minority. The claim splits the layer in two.

• Removal region \(t+1\le w\le s-1\): one sheep goes, so \(U_s(w)=U_{s-1}(w-1)\); repeating this cuts the white flock down to \(b-1\).
• Waiting region \(1\le w\le t\): nothing goes, so \(U_s(w)=C_s(w)\), that is

$$-\frac{s-w}{s}U_s(w-1)+U_s(w)-\frac{w}{s}U_s(w+1)=0.$$

The waiting region is a linear system in \(U_s(1),\dots,U_s(t)\), with the known values \(U_s(0)=s\) and \(U_s(t+1)=U_{s-1}(t)\) at its two ends. Its matrix is tridiagonal with diagonal \(1\) and off-diagonal entries \(-(s-w)/s\) and \(-w/s\), whose absolute values sum to \(1\) in an inner row and to less in the first and last row, where a neighbour has moved to the right-hand side. Such a matrix — irreducible, weakly diagonally dominant, strictly dominant in at least one row — is non-singular, ^[8] so each layer has one solution, and row dominance also makes the Thomas algorithm, a forward elimination plus a back substitution in \(O(t)\) operations, stable without pivoting. ^[7]

The difference equation of Section 4 also solves the waiting region in closed form, \(U_s(w)=s+D\sum_{j\lt w}\binom{s-1}{j}\) with \(D\) fixed by the right boundary, which proves unique solvability a second time. Forming those coefficients is hopeless in machine arithmetic — \(\binom{19999}{9999}\) has \(6019\) digits — and a scaled evaluation would avoid that, but the code solves the small system instead.

8. Where the clock starts

Peredur may act only after a bleat, so he cannot touch the balanced start \((n,n)\). The first bleat is white or black with probability \(\tfrac12\) each, and both outcomes lie in the layer \(s=2n\):

$$E(n)=\tfrac12\bigl(U_{2n}(n+1)+U_{2n}(n-1)\bigr).$$

The restriction is real: with a removal allowed first the value would be the larger \(V(n,n)=U_{2n}(n)\), for instance \(\tfrac73\) against \(\tfrac{55}{24}\) at \(n=2\). The gap shrinks like \(1/(4n)\) — \(2.478\cdot10^{-5}\) at \(n=10\,000\) — still far more than the six decimals asked for, so the two states must be averaged.

9. Worked example: \(E(3)\) by hand

The layers up to \(s=5\) give \(U_5(2)=V(2,3)=\tfrac{125}{33}\), \(U_5(3)=V(3,2)=\tfrac73\) and \(U_5(4)=V(4,1)=1\). In the layer \(s=6\), \(t=2\), so \(w=3,4,5\) copy those three values. With \(U_6(0)=6\) and the right boundary \(U_6(3)=\tfrac{125}{33}\), the waiting region \(w=1,2\) has the equations

$$U_6(1)-\tfrac16U_6(2)=5,$$

$$-\tfrac23U_6(1)+U_6(2)=\tfrac{125}{99}.$$

Substituting the first into the second gives \(\tfrac89U_6(2)=\tfrac{455}{99}\), hence \(U_6(2)=\tfrac{455}{88}\approx5.170455\) and \(U_6(1)=\tfrac{3095}{528}\approx5.861742\). With \(U_6(4)=\tfrac73\) the average of Section 8 is

$$\begin{aligned}E(3)&=\tfrac12\Bigl(\tfrac{455}{88}+\tfrac73\Bigr)\\&=\tfrac{1981}{528}\approx3.751894.\end{aligned}$$

The statement's checkpoint comes out of the layer \(s=10\) the same way: \((6,4)\) is a removal position and three removals lead to \((3,4)\), so \(U_{10}(6)=\tfrac{335}{63}\approx5.317460\), while the waiting position \((4,6)\) has \(U_{10}(4)=\tfrac{3691465}{438144}\approx8.425232\). Their average is

$$E(5)=\frac{18063835}{2628864}=6.871346\ldots,$$

which rounds to the \(6.871346\) of the statement.

10. A faster route along the diagonal

Along the diagonal the values obey a one-term recursion. With \(p_k=4^{-k}\binom{2k}{k}\) and \(v_k=V(k,k)\), the published analysis gives ^[4]

$$v_{k+1}=\frac{1-p_k}{1+p_k}v_k+\frac{2(2k+1)p_k}{1+p_k},$$

with \(v_1=1\) and \(p_1=\tfrac12\), \(p_{k+1}=p_k(2k+1)/(2k+2)\); we checked it against the exact rational values for all \(k\le14\). It reaches \(E(n)\) in \(O(n)\) steps: the removal rule makes \(U_{2n}(n+1)=v_{n-1}\), since \((n+1,n-1)\) and \((n-1,n-1)\) are both cut back to \((n-2,n-1)\), and \(U_{2n}(n)=v_n\) is the right boundary of that layer, so the closed form of Section 7 gives \(U_{2n}(n-1)=v_n+2p_n(2n-v_n)\). Averaging,

$$E(n)=\frac{v_{n-1}+v_n}{2}+p_n(2n-v_n)\qquad(n\ge2).$$

It also says what the answer means: the loss \(2n-E(n)\) of Section 6 stays numerically close to \(\sqrt{\pi n}\) — \(24.266415\) against \(25.066283\) at \(n=200\), \(395.546346\) against \(396.332730\) at \(n=50\,000\) — so nearly all \(2n\) sheep end up black, against \(n\) with no removals at all.

How the Code Works

solve_fast implements Sections 7 and 8 directly. ^[9] Two buffers prev and cur hold the layers \(s-1\) and \(s\). For each \(s\) up to \(2n\) it writes the ends cur[0] = s, cur[s] = 0, computes t = (s - 1) / 2 in integers and fills the removal region by cur[w] = prev[w-1]. If \(t\ge1\) it eliminates the waiting system in the scratch arrays b, c, d: diagonal, superdiagonal, right-hand side. Both boundaries enter through d: d[1] = s - 1 is \(\frac{s-1}{s}U_s(0)\), and the term carrying right_boundary = prev[t] goes to d[1] if \(t=1\) and to d[t] otherwise. The forward sweep uses mult = a / b[w-1] with a = -(s - w) / s, the back substitution runs down from cur[t], the buffers swap, and the function finally returns 0.5 * (prev[n-1] + prev[n+1]).

solve_exact_table is an independent solver that assumes nothing about the policy, in long double despite its name. For every layer it stores the obstacle obstacle[w] = value[s-1][w-1], starts there, and sweeps \(w=1,\dots,s-1\) replacing each entry by the larger of the obstacle and the continuation value built from the already updated left neighbour. That projected Gauss–Seidel iteration climbs monotonically to the solution of Section 3. It stops when a sweep moves nothing by more than \(10^{-18}\) — a stopping rule, not an error bound — and throws after a million sweeps; a flock of \(80\) needs \(309\) sweeps in its slowest layer on the build measured here, so only small flocks go through it. solve_exact_small reads \(E(n)\) from such a table with expectation_from_table.

verify_cap_policy_on_small_domain tests the claim of Section 5 by brute force: for every position with at most \(80\) sheep it reads \(b\) and every candidate \(C(k,b)\) from that table, keeps a candidate only if it beats the best by more than \(10^{-12}\) so that ties fall to the smaller \(k\), and fails unless the winner is \(\min(w,b-1)\). run_checkpoints runs that test, compares solve_fast with solve_exact_small for \(n=1,5,10,20,30\) within \(10^{-11}\), and checks \(E(5)\) against the statement within \(5\cdot10^{-7}\). main accepts --n= and --skip-checkpoints and prints six decimals. Its parsing rejects overflowing arguments, though the later 2 * n is an unchecked int, so the accepted range is wider than the usable one.

The Python and Java files hold the fast solver only, with \(n\) fixed at \(10\,000\) and no self-checks; Java formats with String.format(Locale.US, "%.6f", ans), Python with an f-string. The C++ computes in long double, whose width depends on the platform, while Java and Python work in double; the runs below agree to all six printed decimals. The C++ reuses its buffers, whereas Java and Python allocate a fresh layer array each time and advance with prev = cur, costing time but no accuracy.

Complexity and Verification

Layer \(s\) is filled in \(O(s)\) operations, so the sweep costs \(\sum_{s\le2n}O(s)=O(n^2)\); its interior states number \(\sum_{s\le2n}(s-1)=199\,990\,000\) for \(n=10\,000\). Only two layers and three scratch arrays are alive, some \(70\,000\) entries, so the working set is \(O(n)\) and under two megabytes. On one laptop, with clang at -O2, the C++ needs \(0.88\) seconds at the target size and its checkpoints less than \(0.05\) more, Java \(0.96\) and Python \(19\) seconds.

The C++ runs the checkpoints above and Python and Java none; all three print the same six decimals, the policy test passes up to \(80\) sheep, and the two solvers agree to \(6.4\cdot10^{-14}\) or better at \(n=1,5,10,20,30\). Further checks made for this article:

• A solver in exact rational arithmetic, built from the closed form of Section 7, matches the published program for \(n\le32\) to within \(2.8\cdot10^{-14}\), and gives \(E(2)=\tfrac{55}{24}\), \(E(3)=\tfrac{1981}{528}\), \(E(4)=\tfrac{468905}{88704}\), \(E(5)=\tfrac{18063835}{2628864}\).
• In the same arithmetic the lemma of Section 5 holds at every position of every flock with at most \(160\) sheep, which proves the policy optimal there.
• That test in floating point at \(n=10\,000\) passes at all \(199\,990\,000\) interior positions. Apart from the tie at \((1,1)\), the smallest margins are \(5.0\cdot10^{-5}\) for waiting over removing, at \((9999,10000)\), and \(2.5\cdot10^{-5}\) for removing over waiting, at the balanced \((10000,10000)\) — six orders of magnitude above the rounding error, though not a rigorous bound on it. The elimination is sound too: its smallest pivot is \(0.5056\) and every solved row satisfies its equation to \(6.3\cdot10^{-16}\) relative to the flock size.
• Rerunning the recurrence in double-double arithmetic, about \(31\) digits, moves the result by \(8.7\cdot10^{-12}\) at \(n=2000\) and \(3.3\cdot10^{-11}\) at \(n=10\,000\), far below half a unit in the sixth decimal; that supports the printed digits without proving them, since a value can in principle sit that close to a rounding boundary.
• The diagonal recursion of Section 10, a different route costing \(O(n)\), gives the same rational values for \(n\le14\) and agrees with the layer sweep to \(3.6\cdot10^{-11}\) at \(n=10\,000\).
• In \(200\,000\) simulated games at \(n=5\) the mean harvest was \(6.873\) black sheep, and harvest plus loss came to \(10\) in every game, as Section 6 requires.

Footnotes and References

The references below supply the background facts used above. The layer recurrence, the closed form of the waiting region, the verification lemma and the link from the diagonal values to \(E(n)\) are derived in this article.

1. Project Euler 339 — Peredur Fab Efrawg. The statement starts both flocks with n sheep, lets Peredur remove white sheep after a bleat and a crossing, defines E(n) as the expected final number of black sheep under optimal play, gives E(5) = 6.871346 and asks for E(10000) to six decimals.
2. The Mabinogion (Guest 1877) — Peredur the Son of Evrawc. The translation the statement quotes: in a wooded valley with a river, white sheep stand on one side and black sheep on the other, and whenever one bleats a sheep of the other flock crosses over and changes colour.
3. Wikipedia — Mabinogion sheep problem. The same dynamics as a control problem of David Williams (Probability with Martingales, 1991); it states the rule of Section 5, to remove just enough white sheep so that there are more black sheep than white sheep.
4. D. Stenlund, On the Mabinogion urn model, Advances in Applied Probability 50 (2018), 327–346. It restates Williams' policy exactly as in Section 5 and notes that he proved it optimal with martingale theory (Probability with Martingales, pp. 159–163). Its Proposition 1 is the uncontrolled value P of Section 4, with the balanced case P = k for k white and k black sheep (the paper calls that V, which here means the optimal value); its Proposition 3 is the diagonal recursion of Section 10.
5. Wikipedia — Gambler's ruin. It solves the classical walk with a fixed win probability p, where the ruin probabilities are built from powers of q/p; Section 4 repeats that calculation for the state-dependent probabilities w/s, whose ratios multiply to binomial coefficients.
6. Wikipedia — Snell envelope. It defines the envelope as the smallest supermartingale dominating a payoff process, computed backwards as the maximum of the current payoff and the conditional expectation of the next value — the equation of Section 3.
7. Wikipedia — Tridiagonal matrix algorithm. It gives the Thomas algorithm as a forward sweep of modified coefficients plus a back substitution in O(n) operations, and states that it is stable without pivoting for matrices diagonally dominant by rows or columns, the case of Section 7.
8. Wikipedia — Diagonally dominant matrix. It calls a matrix irreducibly diagonally dominant when it is irreducible, weakly dominant and strictly dominant in at least one row, and states that such a matrix is non-singular, which is the waiting region’s case.
9. C++ — Euler339.cpp. The linked immutable revision contains solve_fast, solve_exact_table, verify_cap_policy_on_small_domain and run_checkpoints, and the comment stating the removal rule as white = black − 1; the independent checks above are not part of it.