Problem 329: Prime Frog

Detailed mathematical approach

Problem Summary

A frog hops over \(n=500\) squares numbered \(1\) to \(500\), each move going one square left or right with probability \(\frac12\) and never off the board: from square \(1\) the move is forced to \(2\), from \(500\) to \(499\). Just before every move he croaks once: on a prime square P with probability \(\frac23\) and N with \(\frac13\), on a non-prime square the other way round. His starting square is uniform over all \(500\). Listening to his first \(15\) croaks, we want the exact probability of PPPPNNPPPNPPNPN — ten Ps and five Ns — as a fraction \(p/q\) in lowest terms. ^[1] The statement gives no numerical checkpoints, so Sections 7 and 9 build some.

The whole solution in one sentence: keep one exact probability per square, multiply by the probability of each required croak, spread the weights to the neighbouring squares between successive croaks, and sum the entries after the last croak.

Mathematical Approach

1. The frog's walk, written out

First, why this is easy. On a five-square board the croaks PN reach square \(2\) by the route \((1,2)\) or \((3,2)\), each of probability \(\frac1{45}\); their sum \(\frac2{45}\) may be kept as one number there, since both routes face the same future from that square. Merging histories square by square is the whole idea.

Write \(X_t\) for the square the frog occupies at his \(t\)-th croak, and \(c_t\) for that croak. Only the fourteen moves between croaks matter, because the croak comes first: croak \(1\) happens on the random starting square, and the move after croak \(15\) changes nothing Susan hears. The move probabilities form the matrix

$$T_{i,j}=\begin{cases}1,&i\in\{1,n\},\ |i-j|=1,\\ \frac12,&2\le i\le n-1,\ |i-j|=1,\\ 0,&|i-j|\ne1.\end{cases}$$

The move rule looks only at the frog's current square, so the walk is a Markov chain: the next square does not depend on the earlier ones. Every entry of \(T\) is non-negative and every row sums to \(1\), so \(T\) is stochastic. ^[4]

2. What exactly must be added up

Call a route a list \(x=(x_1,\dots,x_L)\) of squares with \(|x_{t+1}-x_t|=1\), where \(L=15\), and let \(q_t(i)\in\{\frac13,\frac23\}\) be the probability that a frog on square \(i\) croaks \(c_t\). The frog takes the route \(x\) and croaks our sequence with probability

$$\begin{aligned}w(x)&=\frac1n\prod_{t=1}^{L}q_t(x_t)\\&\quad\times\prod_{t=1}^{L-1}T_{x_t,x_{t+1}},\end{aligned}$$

and we want \(\sum_x w(x)\) over every route. Each move has two choices, only one off an end, so the routes number at most \(n2^{L-1}\); here exactly \(8\,108\,856\). Enumerating them would work here, but Section 3 instead keeps one accumulated number per square and croak.

3. A hidden Markov model, and the forward recursion

Susan never sees the frog, and each croak depends only on the square he is on. A Markov chain whose states are invisible and which emits one observation per step, depending only on the current state, is a hidden Markov model; here the square is the hidden state, the croak the observation, and the problem asks the standard likelihood question for such a model. ^[2] Its classical answer is this recursion. For \(1\le t\le L\) let \(\alpha_t(i)\) be the probability that the first \(t\) croaks are \(c_1,\dots,c_t\) and that \(X_t=i\). Then

$$\alpha_1(i)=\frac{q_1(i)}{n},$$

$$\alpha_{t+1}(j)=q_{t+1}(j)\sum_{i}\alpha_t(i)\,T_{i,j}.$$

The first holds because the frog starts on square \(i\) with probability \(\frac1n\) and then croaks \(c_1\) with probability \(q_1(i)\). For the second, split the partial routes counted by \(\alpha_{t+1}(j)\) by the square \(i=x_t\) held at the previous croak: each is one counted by \(\alpha_t(i)\), then the move \(i\to j\) of probability \(T_{i,j}\), then the croak \(c_{t+1}\) of probability \(q_{t+1}(j)\); since \(x_t\) is fixed by the route, nothing is counted twice. Full routes are counted once under their last square, so

$$\Pr(c_1c_2\ldots c_L)=\sum_{i=1}^{n}\alpha_L(i).$$

This is the forward algorithm. ^[3] A step costs \(O(n^2)\) for a general model, the inner sum running over all states; here \(T_{i,j}\ne0\) only when \(|i-j|=1\), so a step costs \(O(n)\) and the pass \(O(Ln)\).

4. Which squares are prime

A sieve of Eratosthenes settles all \(500\) squares at once: mark \(0\) and \(1\) as not prime — neither is composite — then for every still-marked \(p\) with \(p^2\le500\) strike out \(p^2,\,p^2+p,\,p^2+2p,\dots\). Starting at \(p^2\) loses nothing, since every multiple \(kp\) with \(2\le k\lt p\) has a prime factor below \(p\) and went earlier, and stopping once \(p^2\gt500\) is safe, since every composite \(m\le500\) has a prime factor at most \(\sqrt m\). ^[5] There are \(95\) primes in \(1,\dots,500\) — the \(95\)-th is \(499\) and the \(96\)-th is \(503\) ^[6] — so \(405\) squares are not prime. The sieve runs once per probability.

5. From fractions to whole numbers

Apart from the opening \(\frac1n\), every factor in \(w(x)\) is one of \(\frac23,\frac13,\frac12,1\), so clearing denominators turns the work into integer arithmetic, with nothing rounded. Let \(\chi(i)=1\) when \(i\) is prime and \(0\) otherwise, and \(\sigma_t=1\) when \(c_t\) is P and \(0\) when it is N. The croak \(c_t\) then matches square \(i\) precisely when \(\chi(i)=\sigma_t\), and \(q_t(i)=e_t(i)/3\) for

$$e_t(i)=\begin{cases}2,&\chi(i)=\sigma_t,\\1,&\chi(i)\ne\sigma_t.\end{cases}$$

Likewise put \(U_{i,j}=2T_{i,j}\), so \(U_{i,j}=1\) for an ordinary move and \(U_{i,j}=2\) for a forced move off an end. With the scaled quantities \(v_t(i)=n\,3^{t}\,2^{t-1}\,\alpha_t(i)\), which Section 6 shows to be non-negative integers,

$$\Pr(c_1\ldots c_L)=\frac{\sum_i v_L(i)}{n\,3^{L}\,2^{L-1}}.$$

For \(n=500\) and \(L=15\) that denominator is \(500\cdot3^{15}\cdot2^{14}\), and since \(500=2^2\cdot5^3\) it equals \(2^{16}3^{15}5^{3}\). Whatever cancels, the printed denominator is built from \(2\), \(3\) and \(5\) alone — a cheap check on the output.

6. The integer recursion the code runs

Substituting \(q_t=e_t/3\), \(T=U/2\) and \(\alpha_t=v_t/(n3^{t}2^{t-1})\) into Section 3 cancels every power of \(2\), \(3\) and \(n\) and leaves

$$\begin{aligned}v_1(i)&=e_1(i),\\ v_{t+1}(j)&=e_{t+1}(j)\sum_{i}v_t(i)\,U_{i,j}.\end{aligned}$$

Both right-hand sides are integers built from integers, so by induction every \(v_t(i)\) is a non-negative integer. This is what the programs do: they hold \(v_t\) in whole numbers and keep the denominator separately, multiplying it by \(3\) at each croak and by \(2\) at each move, then divide both by their greatest common divisor, found by the Euclidean algorithm. ^[7] Expanding the recursion back into a sum over routes describes the numerator neatly,

$$\sum_i v_L(i)=\sum_x 2^{\,m(x)+f(x)},$$

where \(m(x)\) counts the route's matching croaks and \(f(x)\) its forced moves: every route contributes a power of two.

7. A complete worked example

Shrink the board to \(n=5\) and listen to \(L=2\) croaks, PN. The primes are \(2,3,5\), the non-primes \(1,4\), the denominator \(5\cdot3^2\cdot2=90\). Croak \(1\) is P, so \(e_1\) doubles the prime squares: \(v_1=(1,2,2,1,2)\), total \(8\). Now push along \(U\): the ends send twice their value to their one neighbour, \(2\) from square \(1\) into \(2\) and \(4\) from square \(5\) into \(4\), while each interior square sends its value into both neighbours. The arrivals are \((2,4,3,6,1)\), total \(16\), exactly twice the previous total — and the \(4\) on square \(2\) is the \(\frac2{45}\) of Section 1. Croak \(2\) is N, doubling the non-prime squares \(1\) and \(4\):

$$v_2=(4,4,3,12,1),\qquad \textstyle\sum_i v_2(i)=24.$$

So hearing PN here has probability \(24/90=4/15\). The eight routes give the same answer directly, with weights \(2^{m+f}\) of \(2,4,2,2,4,1,1,8\) for \((1,2),(2,1),(2,3),(3,2)\) and \((3,4),(4,3),(4,5),(5,4)\): the last has a matching P on the prime \(5\), the forced move \(5\to4\) and a matching N on the non-prime \(4\), weight \(2^3\). They sum to \(24\), as they must.

8. Why the croaks are not independent

A tempting shortcut is to multiply one-croak probabilities: \(\frac{119}{300}\) ten times and \(\frac{181}{300}\) five times (Section 9). That assumes independence, which fails here. A marginal argument cannot settle it: on the five-square board a croak is P with probability \(\frac8{15}\) and N with \(\frac7{15}\), and the second croak has those very marginals too. But Section 7 gives the joint value \(\Pr(\mathtt{PN})=\frac4{15}\), so dividing by \(\frac8{15}\) yields

$$\Pr(c_2=\mathtt N\mid c_1=\mathtt P)=\frac12,$$

which is not \(\frac7{15}\), and the joint \(\frac4{15}\) differs from the product \(\frac{56}{225}\) of the marginals. Hearing P reweights the squares the frog may be on, and those weights steer the next croak; \(v_t\) carries exactly that information.

9. Two checks that cost nothing

First, the case \(L=1\): with no move, the answer is a plain average over the starting squares,

$$\begin{aligned}\Pr(\mathtt P)&=\frac{95\cdot\frac23+405\cdot\frac13}{500}\\&=\frac{595}{1500}=\frac{119}{300},\end{aligned}$$

whence \(\Pr(\mathtt N)=\frac{181}{300}\), and the two add to \(1\). Second, a normalisation identity. Every row of \(U\) sums to \(2\), so a move doubles the total weight exactly — Section 7 shows \(8\to16\) — and the factors \(e_t(i)\) for the two possible croaks at step \(t\) are \(2\) and \(1\) in some order, summing to \(3\). Adding the numerators over all \(2^{L}\) croak sequences therefore replaces every emission factor by \(3\) and gives \(n\,3^{L}\,2^{L-1}\), the common denominator: the \(2^{15}=32768\) probabilities add to \(1\).

How the Code Works

In the C++ program ^[9] parse_arguments fills an Options struct from --squares=, --sequence= and --skip-checkpoints, defaulting to the statement's board and sequence; it accumulates digits in an int with no overflow check, so an oversized --squares= value is not safely rejected. sieve_prime_flags is Section 4.

probability_fraction is Sections 5 and 6 in a single loop. It fills state with ones and sets denom to the number of squares, standing together for the uniform start before any croak; state equals \(v_t\) just after the emission step for croak \(t\). For each croak it doubles state[pos] wherever the boolean good, which is \(\chi(i)=\sigma_t\), holds, and multiplies denom by \(3\). For every croak but the last it builds a fresh vector next: squares \(1\) and squares push twice their value into their one neighbour — the entries \(U_{i,j}=2\) — interior squares push theirs into both neighbours, and denom gains a factor \(2\). At the end it sums the entries and divides both parts by gcd_cpp, the Euclidean loop.

run_checkpoints checks \(\Pr(\mathtt P)=\frac{119}{300}\), then recomputes that value with \(\Pr(\mathtt N)\) and confirms their sum by cross-multiplication, so that no division is needed. Both use \(500\) squares whatever --squares says, being fixed reference values rather than a test of the options. These are the one-croak identities of Section 9 only: its normalisation identity goes untested and no move step is exercised. Each call sieves afresh, so a default C++ run sieves four times where Python and Java sieve once.

Python and Java run the same recursion and loop structure without the options or checkpoints; only the arithmetic differs. Python's integers are unlimited by nature, so state is a list of int and math.gcd finishes; Java uses immutable BigInteger values, assigning each multiply and add back into the array, and reduces with numer.gcd(denom); C++ borrows cpp_int from Boost.Multiprecision. ^[8] One asymmetry is harmless: C++ initialises state[0] to \(1\) where Python clears it and Java leaves it null, but every loop starts at index \(1\).

Arbitrary precision is in fact not needed. Dividing the total weight after \(t\) croaks by \(n\,3^{t}2^{t-1}\) turns it into a probability, so that product bounds it; for \(n=500\) the bound stays below \(2^{63}\) up to \(t=21\), and the largest entry of any \(v_t\) has eight digits. Signed 64-bit arithmetic would serve; the big-integer types keep the C++ program correct for an arbitrary --sequence.

Complexity and Verification

Each probability computed costs one sieve, \(O(n\log\log n)\) steps. ^[5] The forward pass visits \(L\cdot n=7500\) square-and-croak pairs, but only matching ones double an entry: the ten Ps match \(95\) squares each and the five Ns match \(405\), giving \(2975\) doublings. The move loops visit \(7000\) sources and make \(14\cdot(2\cdot498+2)=13972\) additions. That is \(O(Ln)\) big-integer operations and two arrays of \(n+1\) entries — counts of operations and stored numbers, not of seconds and bytes, since the entries grow to \(O(L+\log n)\) bits. The work sits more than two orders of magnitude below the route count of Section 2, and all three programs finish in about \(0.01\) s for C++ and Python, \(0.04\) s for Java.

Beyond the C++ checkpoints, the following checks all agreed exactly. The integer recursion, a forward pass using exact rational probabilities and a direct enumeration of every route agree on all \(434\) cases given by every P/N sequence of length \(1\) to \(5\) on every board of \(2\) to \(8\) squares, and on larger ones: the C++ program prints \(\frac4{15}\) for \(n=5\) with PN, \(\frac{1595}{27216}\) for \(n=7\) with PNPNP and \(\frac{23821}{37791360}\) for \(n=20\) with PPPPNNPPPN. On the full problem the published fraction comes back from the exact rational implementation and from a backward recursion that sums the routes by their first square, not their last. A floating-point version using an explicit \(500\times500\) transition matrix, its rows confirmed to sum to \(1\), agrees to a relative error of \(9\cdot10^{-16}\). The route-weight formula of Section 6 reproduces the numerator in five cases; replacing every emission factor by \(3\) returns \(500\cdot3^{15}\cdot2^{14}\), and on a \(12\)-square board the \(2^{7}\) probabilities of all length-\(7\) sequences sum to \(1\), as Section 9 requires. Finally, \(200\) million simulated frogs produced \(1380\) that croaked the requested sequence, \(0.6\) standard deviations from the exact probability.

Footnotes and References

The references below supply the background facts; the recursion, the rescaling, the route weights and the normalisation identity are derived here.

1. Project Euler 329 — Prime Frog. The official statement moves a frog over squares 1 to 500, one square left or right with equal probability, forced inwards from either end, croaking P with probability 2/3 on a prime square and 1/3 otherwise and the reverse for N; from a random start it asks the probability that his first 15 croaks are PPPPNNPPPNPPNPN, as a reduced fraction.
2. Wikipedia — Hidden Markov model. Defines the model as a Markov process whose states cannot be observed directly, each observation depending only on the state at that step, and gives an observation sequence's likelihood as the sum of P(Y | X) P(X) over all hidden X.
3. Wikipedia — Forward algorithm. Gives the forward variable α(xₜ) = p(xₜ, y₁:ₜ) and the recursion α(xₜ) = p(yₜ | xₜ) Σ p(xₜ | xₜ₋₁) α(xₜ₋₁) of Section 3, at a cost of Θ(n m²) for m states.
4. Wikipedia — Markov chain. States the Markov property, that the next state depends only on the present one, and describes a finite chain's transition matrix as right stochastic, its rows summing to one — the form of \(T\) in Section 1.
5. Wikipedia — Sieve of Eratosthenes. Gives the variant of Section 4, crossing out i², i²+i, i²+2i, … for each i up to √n still marked prime, at a running time of O(n log log n).
6. Wikipedia — List of prime numbers. Tabulates the first 1000 primes in rows labelled by rank; ranks 81–100 show 499 as the 95th prime and 503 as the 96th, so 95 of the 500 squares carry a prime.
7. Wikipedia — Euclidean algorithm. Describes the remainder form used for the final reduction, replacing the larger number by its remainder on division by the smaller until zero, and notes that it reduces fractions to lowest terms.
8. Python — numeric types gives integers unlimited precision; Java API — BigInteger documents immutable arbitrary-precision integers and a gcd of its operands' absolute values; Boost.Multiprecision — cpp_int documents the signed arbitrary-precision type from boost/multiprecision/cpp_int.hpp.
9. C++ — Euler329.cpp. The linked immutable revision holds sieve_prime_flags, gcd_cpp, probability_fraction and run_checkpoints; it is this article's implementation record.