Detailed mathematical approach

Problem Summary

There are \(k\) bacterial types. Type \(i\) does not have just one reproduction rule; it has \(m\) equally likely rules extracted from the deterministic pseudo-random sequence

$r_0=306,\qquad r_{t+1}\equiv r_t^2 \pmod{10007}.$

The block \(r_{im},r_{im+1},\dots,r_{im+m-1}\) belongs to type \(i\). Starting from a single bacterium of type \(0\), we repeatedly apply these rules independently to all descendants and ask for the probability that the entire population eventually disappears.

Mathematical Approach

This is a multitype branching process. The key simplification is that each type’s \(m\) rules can be compressed into five residue counts, after which the extinction problem becomes a fixed-point system on \([0,1]^k\).

Step 1: Compress the pseudo-random rules into five counts

For each type \(i\) and each residue \(q\in\{0,1,2,3,4\}\), define

$c_{i,q}=\#\left\{j\in\{0,\dots,m-1\}: r_{im+j}\equiv q \pmod{5}\right\}.$

These counts satisfy

$\sum_{q=0}^4 c_{i,q}=m.$

Once the counts are known, the individual rule values are irrelevant. Every bacterium of type \(i\) simply chooses uniformly from a multiset containing \(c_{i,0}\) rules of kind \(0\), \(c_{i,1}\) rules of kind \(1\), and so on.

Step 2: Translate each residue class into an offspring pattern

The five residue classes encode five possible outcomes for one bacterium of type \(i\).

Residue \(0\): no offspring are produced.

Residue \(1\): two offspring of type \(i\) are produced.

Residue \(2\): one offspring of type \(2i \bmod k\) is produced.

Residue \(3\): three offspring of type \((i^2+1) \bmod k\) are produced.

Residue \(4\): one offspring of type \(i\) and one offspring of type \((i+1) \bmod k\) are produced.

If \(s_j\) denotes the extinction contribution of a single child of type \(j\), these five cases contribute the monomials

$1,\qquad s_i^2,\qquad s_{2i \bmod k},\qquad s_{(i^2+1) \bmod k}^3,\qquad s_i s_{(i+1) \bmod k}.$

The squares, cubes, and products come from independence: all descendants created by the chosen rule must die out.

Step 3: Build the generating polynomial and the fixed-point system

Let \(\mathbf{s}=(s_0,\dots,s_{k-1})\). Averaging over the \(m\) equally likely rules of type \(i\) gives the offspring generating polynomial

$f_i(\mathbf{s})=\frac{c_{i,0}+c_{i,1}s_i^2+c_{i,2}s_{2i \bmod k}+c_{i,3}s_{(i^2+1) \bmod k}^3+c_{i,4}s_i s_{(i+1) \bmod k}}{m}.$

Now let \(p_i\) be the probability that the process eventually becomes extinct when it starts from one bacterium of type \(i\). Standard branching-process reasoning gives

$p_i=f_i(\mathbf{p})=\frac{c_{i,0}+c_{i,1}p_i^2+c_{i,2}p_{2i \bmod k}+c_{i,3}p_{(i^2+1) \bmod k}^3+c_{i,4}p_i p_{(i+1) \bmod k}}{m}.$

So the desired answer is the first coordinate of the extinction vector \(\mathbf{p}\).

Step 4: Solve the system by monotone iteration

Define the map

$T(\mathbf{x})=\left(f_0(\mathbf{x}),f_1(\mathbf{x}),\dots,f_{k-1}(\mathbf{x})\right).$

Because all coefficients are non-negative and each row is averaged by \(m\), the map sends \([0,1]^k\) into itself and is monotone in every coordinate. Therefore, if we start from

$\mathbf{p}^{(0)}=(0,\dots,0),\qquad \mathbf{p}^{(t+1)}=T\left(\mathbf{p}^{(t)}\right),$

the iterates form an increasing sequence bounded above by \((1,\dots,1)\). They converge to the minimal fixed point of the system, which is exactly the extinction vector for the branching process.

Worked Example: \(k=2\), \(m=2\)

The first four pseudo-random values are

$306,\qquad 3573,\qquad 7404,\qquad 870,$

so the residues modulo \(5\) are

$1,\qquad 3,\qquad 4,\qquad 0.$

Thus type \(0\) receives one residue-\(1\) rule and one residue-\(3\) rule, while type \(1\) receives one residue-\(4\) rule and one residue-\(0\) rule. The extinction equations become

$p_0=\frac{p_0^2+p_1^3}{2},\qquad p_1=\frac{1+p_0p_1}{2}.$

Starting from \((0,0)\), the first synchronous iterates are

$\begin{aligned} \mathbf{p}^{(1)}&=\left(0,\frac12\right),\\ \mathbf{p}^{(2)}&=\left(\frac1{16},\frac12\right),\\ \mathbf{p}^{(3)}&=\left(\frac{33}{512},\frac{33}{64}\right). \end{aligned}$

Continuing this process gives

$p_0\approx 0.07243802,$

which matches the small verification case used by the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same structure. First they generate the full pseudo-random stream of length \(km\), then they split it into \(k\) consecutive blocks of size \(m\) and tally the five residue counts for each type.

Next they keep two probability vectors. During one sweep over all types, every new coordinate is computed only from the previous sweep’s vector, so the update is synchronous rather than in-place. After the sweep, the vectors are swapped.

Each sweep also records the largest absolute change between old and new coordinates. The C++ implementation uses extended precision and stops once this change is below \(10^{-20}\); the Python and Java implementations use double precision and stop at \(10^{-15}\). All three impose a hard cap of \(10^6\) iterations.

After convergence, the implementation prints the extinction probability of type \(0\) rounded to eight decimal places.

Complexity Analysis

Generating the pseudo-random stream and building the five-count histogram for every type costs \(O(km)\) time. If the fixed-point solver needs \(I\) sweeps, the nonlinear iteration costs \(O(Ik)\) time, so the total running time is

$O(km+Ik).$

As written, the implementations store the entire stream of length \(km\), the \(5k\) histogram counts, and two probability vectors of length \(k\). Therefore the memory usage is \(O(km)\).

Footnotes and References

1. Problem page: https://projecteuler.net/problem=666
2. Branching process: Wikipedia — Branching process
3. Multi-type branching process: Wikipedia — Multi-type branching process
4. Probability-generating function: Wikipedia — Probability-generating function
5. Fixed-point iteration: Wikipedia — Fixed-point iteration