Project Euler Problem 666

Solution

Answer: 0.48023168

This is a classic multi-type Galton–Watson branching process.

Let

$$x_i=\Pr(\text{eventual extinction}\mid \text{start with one bacterium of type }\alpha_i).$$

The required answer is $P_{500,10}=x_0$.

For a branching process, the extinction probabilities are the minimal nonnegative fixed point of the offspring generating functions.

For each type $i$, define

$$f_i(x_0,\dots,x_{k-1})$$

as the generating function corresponding to one minute of evolution of a single bacterium of type $\alpha_i$.

Because each choice of $j\in[0,m-1]$ is equally likely, we average over the $m$ possible rules.

For a given $q=r_{im+j}\bmod 5$:

$q=0$: death

contribution $=1$

$q=1$: clone itself

offspring: two $\alpha_i$

contribution $=x_i^2$

$q=2$: mutate to $\alpha_{2i\bmod k}$

contribution $=x_{2i\bmod k}$

$q=3$: split into three $\alpha_{i^2+1\bmod k}$

contribution $=x_{i^2+1\bmod k}^3$

$q=4$: keep itself and spawn $\alpha_{i+1\bmod k}$

contribution $=x_i x_{i+1\bmod k}$

Hence

$$x_i=f_i(x)$$

with

$$f_i(x)=\frac1m\sum_{j=0}^{m-1} g_{i,j}(x),$$

where $g_{i,j}$ is the corresponding contribution above.

Starting from the zero vector and iterating

$$x^{(n+1)} = f(x^{(n)})$$

converges monotonically to the minimal fixed point, which is exactly the extinction probability vector.

Python implementation

from math import *

def extinction_probability(k, m):
    MOD = 10007

    # Generate q-values
    r = 306
    qs = []
    for _ in range(k * m):
        qs.append(r % 5)
        r = (r * r) % MOD

    # qs for each type i
    rules = [qs[i*m:(i+1)*m] for i in range(k)]

    # Fixed-point iteration
    x = [0.0] * k

    for _ in range(1000):
        nx = [0.0] * k

        for i in range(k):
            s = 0.0

            for q in rules[i]:

                if q == 0:
                    # death
                    s += 1.0

                elif q == 1:
                    # two copies of type i
                    s += x[i] ** 2

                elif q == 2:
                    # mutate
                    s += x[(2 * i) % k]

                elif q == 3:
                    # split into three
                    t = (i * i + 1) % k
                    s += x[t] ** 3

                else:  # q == 4
                    # keep itself and spawn one new bacterium
                    s += x[i] * x[(i + 1) % k]

            nx[i] = s / m

        # convergence check
        diff = max(abs(nx[i] - x[i]) for i in range(k))
        x = nx

        if diff < 1e-16:
            break

    return x[0]

# Verification on the examples
print(round(extinction_probability(2, 2), 8))
print(round(extinction_probability(4, 3), 8))
print(round(extinction_probability(10, 5), 8))

# Required value
ans = extinction_probability(500, 10)
print(round(ans, 8))

Code walkthrough

1. Generate the pseudo-random sequence

r = 306
r = (r * r) % 10007

This produces the sequence $r_n$, and each rule uses

$$q = r_n \bmod 5.$$

2. Build the rule table

rules = [qs[i*m:(i+1)*m] for i in range(k)]

For each type $i$, we store its $m$ possible transformations.

3. Fixed-point iteration

We begin with

x = [0.0] * k

which corresponds to “certain survival” nowhere.

Repeatedly applying the generating functions converges to the smallest fixed point, which is the extinction probability vector.

4. Contributions for each rule

Example:

s += x[i] ** 2

corresponds to producing two independent descendants of type $i$, so both lineages must die out.

Similarly,

s += x[t] ** 3

means three independent descendants.

Verification against the examples

The program reproduces:

$$P_{2,2}=0.07243802$$

$$P_{4,3}=0.18554021$$

$$P_{10,5}=0.53466253$$

matching the statement exactly (to 8 decimal places).

For $(k,m)=(500,10)$, the iteration converges to

$$P_{500,10}=0.48023167580091586\ldots$$

Rounded to 8 decimal places:

$$0.48023168$$

Answer: 0.48023168