Project Euler Problem 329
Susan has a prime frog.
Solution
Answer: 199740353/29386561536000
Let the observed croak sequence be
$$S=\texttt{PPPPNNPPPNPPNPN}$$
of length $15$.
We must compute the exact probability that the frog produces this sequence.
Mathematical analysis
The frog moves on the integers $1,\dots,500$.
- At square $1$, it must jump to $2$.
- At square $500$, it must jump to $499$.
- Otherwise it jumps left/right with probability $1/2$.
Before each jump, it croaks:
-
On a prime square:
-
$P$ with probability $2/3$,
-
$N$ with probability $1/3$.
-
On a non-prime square:
-
$P$ with probability $1/3$,
-
$N$ with probability $2/3$.
The initial position is uniformly distributed over all $500$ squares.
Step 1: Emission probability
Define
$$e(x,c)= \begin{cases} \frac23,& \text{if }x\text{ is prime and }c=P,\[4pt] \frac13,& \text{if }x\text{ is prime and }c=N,\[4pt] \frac13,& \text{if }x\text{ is nonprime and }c=P,\[4pt] \frac23,& \text{if }x\text{ is nonprime and }c=N. \end{cases}$$
If the frog is currently at square $x$, and the next heard croak is $c$, then the probability contribution is $e(x,c)$.
Step 2: Dynamic programming
Let
$$dp_k(x)$$
be the probability that after processing the first $k$ croaks, the frog is at square $x$.
Initially:
$$dp_0(x)=\frac1{500}.$$
For each character $c_k$ in the sequence:
- multiply by the croak probability $e(x,c_k)$,
- then move according to the jump rules.
Transition:
- from $1$ to $2$ with probability $1$,
- from $500$ to $499$ with probability $1$,
- otherwise split equally left/right.
After all $15$ croaks, summing all states gives the desired probability.
Because every transition and emission probability is rational, the final answer is an exact fraction.
Python implementation
from fractions import Fraction
# ---------------------------------------------------------
# Problem data
# ---------------------------------------------------------
N = 500
sequence = "PPPPNNPPPNPPNPN"
# ---------------------------------------------------------
# Sieve of Eratosthenes
# ---------------------------------------------------------
is_prime = [False] * (N + 1)
prime = [True] * (N + 1)
prime[0] = prime[1] = False
for p in range(2, int(N**0.5) + 1):
if prime[p]:
for multiple in range(p * p, N + 1, p):
prime[multiple] = False
for i in range(2, N + 1):
is_prime[i] = prime[i]
# ---------------------------------------------------------
# Initial distribution
# ---------------------------------------------------------
dp = {x: Fraction(1, N) for x in range(1, N + 1)}
# ---------------------------------------------------------
# Process croaks one by one
# ---------------------------------------------------------
for step, ch in enumerate(sequence):
next_dp = {}
for pos, prob in dp.items():
# Probability of hearing this croak at this position
if is_prime[pos]:
emit = Fraction(2, 3) if ch == 'P' else Fraction(1, 3)
else:
emit = Fraction(1, 3) if ch == 'P' else Fraction(2, 3)
prob *= emit
# Last croak: no need to move afterward
if step == len(sequence) - 1:
next_dp[pos] = next_dp.get(pos, Fraction(0)) + prob
continue
# Movement
if pos == 1:
next_dp[2] = next_dp.get(2, Fraction(0)) + prob
elif pos == N:
next_dp[N - 1] = next_dp.get(N - 1, Fraction(0)) + prob
else:
next_dp[pos - 1] = (
next_dp.get(pos - 1, Fraction(0))
+ prob / 2
)
next_dp[pos + 1] = (
next_dp.get(pos + 1, Fraction(0))
+ prob / 2
)
dp = next_dp
# ---------------------------------------------------------
# Final probability
# ---------------------------------------------------------
answer = sum(dp.values())
print(answer)
print(answer.numerator)
print(answer.denominator)
Code walkthrough
Prime sieve
We first compute which integers from $1$ to $500$ are prime.
Initial state
dp = {x: Fraction(1, N) for x in range(1, N + 1)}
The frog starts uniformly at random.
Processing each croak
For each current position:
- multiply by the probability of hearing the observed croak there;
- distribute probability mass to neighboring squares.
At interior squares:
prob / 2
goes left and right.
At endpoints:
1 -> 2
500 -> 499
with probability $1$.
Exact arithmetic
Using Fraction guarantees no floating-point error.
Final verification
The denominator should contain powers of:
- $500$ from the initial distribution,
- $3^{15}$ from the croak probabilities,
- $2^{14}$ from the jumps.
Indeed,
$$29386561536000 = 500 \cdot 3^{15} \cdot 2^{14},$$
which is exactly the expected structure.
The reduced fraction obtained is
$$\frac{199740353}{29386561536000}.$$
Thus the computation is internally consistent.
Answer: 199740353/29386561536000