Project Euler Problem 329

Susan has a prime frog.

Project Euler Problem 329

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:

  1. multiply by the croak probability $e(x,c_k)$,
  2. 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:

  1. multiply by the probability of hearing the observed croak there;
  2. 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