Detailed mathematical approach

Problem Summary

A fair coin is tossed until two consecutive heads appear for the first time. Let \(M\) be the index of the second toss in that first \(HH\) block, so \(M\ge 2\). For a positive integer \(n\), define \(P(n)=\Pr(n\mid M)\), the probability that \(n\) divides the stopping time. If \(P(n)=a/b\) in lowest terms and \(p\) is prime, then \(Q(P(n),p)\) is the least positive residue \(q\) such that \(a\equiv bq\pmod p\).

The target is \(Q(P(10^{18}),10^9+9)\). The number \(10^{18}\) is far too large for direct summation over all multiples of \(n\), so the solution turns the stopping-time question into a 2x2 linear algebra problem and evaluates it directly in modular arithmetic.

Mathematical Approach

Before the first appearance of \(HH\), only one piece of history matters: whether the last toss was a head. That observation collapses the probabilistic process to two transient states.

Step 1: Build the transient-state model

Use two non-absorbing states:

\(S\): the last toss was not a head, or no toss has happened yet.

\(H\): the last toss was a head, but the process has not stopped yet.

From \(S\), a tail keeps the process in \(S\), while a head moves it to \(H\). From \(H\), a tail returns to \(S\), while a head creates the first \(HH\) and ends the experiment. Therefore the transient transition matrix is

$T=\begin{bmatrix}\frac12 & \frac12\\[4pt]\frac12 & 0\end{bmatrix}.$

If

$v_t=\begin{bmatrix}s_t\\ h_t\end{bmatrix}$

is the transient-state vector after \(t\) tosses, then

$v_{t+1}=T\,v_t,\qquad v_0=\begin{bmatrix}1\\0\end{bmatrix}.$

Step 2: Express the exact stopping distribution

The process ends at toss \(t+1\) exactly when, after \(t\) tosses, it is in state \(H\) and the next toss is another head. Hence

$\Pr(M=t+1)=\frac12\,h_t=\frac12\begin{bmatrix}0&1\end{bmatrix}v_t.$

Substituting \(v_t=T^t v_0\) gives

$\Pr(M=t+1)=\frac12\begin{bmatrix}0&1\end{bmatrix}T^t\begin{bmatrix}1\\0\end{bmatrix}.$

The matrix \(T\) is one half of the standard Fibonacci companion matrix, so its powers satisfy

$T^t=\frac{1}{2^t}\begin{bmatrix}F_{t+1}&F_t\\ F_t&F_{t-1}\end{bmatrix}.$

Therefore the exact stopping distribution is

$\Pr(M=m)=\frac{F_{m-1}}{2^m}\qquad (m\ge 2).$

The first values are \(1/4,1/8,1/8,3/32,\dots\), which agrees with a direct case split.

Step 3: Restrict the sum to multiples of \(n\)

By definition,

$P(n)=\sum_{k\ge 1}\Pr(M=kn)=\frac12\sum_{k\ge 1}\begin{bmatrix}0&1\end{bmatrix}T^{kn-1}\begin{bmatrix}1\\0\end{bmatrix}.$

Using the Fibonacci form, the same quantity can be written as

$P(n)=\sum_{k\ge 1}\frac{F_{kn-1}}{2^{kn}}.$

This identity is useful conceptually, but for \(n=10^{18}\) it is still not computationally practical. The crucial step is to rewrite the infinite sum as a matrix geometric series.

Step 4: Turn the infinite series into one inverse

Set

$B=T^n,\qquad w=T^{n-1}\begin{bmatrix}1\\0\end{bmatrix}.$

Then

$T^{kn-1}=T^{(k-1)n}T^{n-1}=B^{k-1}w,$

so

$P(n)=\frac12\begin{bmatrix}0&1\end{bmatrix}\left(\sum_{j\ge 0}B^j\right)w.$

Over the real numbers, every eigenvalue of \(T\) has absolute value less than \(1\), so the series converges and

$\sum_{j\ge 0}B^j=(I-B)^{-1}.$

This yields the exact closed form used by the implementations:

$\boxed{P(n)=\frac12\begin{bmatrix}0&1\end{bmatrix}(I-T^n)^{-1}T^{n-1}\begin{bmatrix}1\\0\end{bmatrix}.}$

Step 5: Pass to arithmetic modulo a prime

If \(P(n)=a/b\) in lowest terms, then \(Q(P(n),p)\) is simply the fraction \(a/b\) interpreted in the field \(\mathbb F_p\):

$Q(P(n),p)\equiv a\,b^{-1}\pmod p.$

So the computation never has to reconstruct the rational number explicitly. Every division by \(2\) or by a 2x2 determinant can be replaced by a modular inverse. That is why the same matrix formula can be evaluated directly modulo \(10^9+9\).

Worked Example: \(n=2\)

For \(n=2\),

$T^2=\begin{bmatrix}\frac12&\frac14\\[4pt]\frac14&\frac14\end{bmatrix},\qquad T\begin{bmatrix}1\\0\end{bmatrix}=\begin{bmatrix}\frac12\\[4pt]\frac12\end{bmatrix}.$

Hence

$I-T^2=\begin{bmatrix}\frac12&-\frac14\\[4pt]-\frac14&\frac34\end{bmatrix}.$

Inverting this matrix and multiplying by \(T(1,0)^T\) gives

$ (I-T^2)^{-1}T\begin{bmatrix}1\\0\end{bmatrix}=\begin{bmatrix}\frac85\\[4pt]\frac65\end{bmatrix}. $

Taking the second component and multiplying by \(1/2\) yields

$P(2)=\frac12\cdot\frac65=\frac35.$

Modulo \(109\), the inverse of \(5\) is \(22\), so

$Q(P(2),109)\equiv 3\cdot 22\equiv 66\pmod{109}.$

A second small checkpoint is \(P(3)=9/31\), which reduces to \(46\) modulo \(109\).

How the Code Works

The C++, Python, and Java implementations all encode the same 2x2 matrix over the prime field \(\mathbb F_{10^9+9}\). The value \(1/2\) is represented by the modular inverse of \(2\), so no floating-point arithmetic is used anywhere.

They compute \(T^n\) and \(T^{n-1}\) with binary exponentiation, form the matrix \(I-T^n\), invert that 2x2 matrix with the determinant formula, and multiply the result by the starting-state column vector. The component corresponding to state \(H\) is then multiplied by \(1/2\) to obtain the stopping probability on the next toss. Because the whole computation already lives in a prime field, the final residue is exactly \(Q(P(n),10^9+9)\).

Complexity Analysis

Every matrix is 2x2, so each multiplication and inversion costs constant time. The dominant work is binary exponentiation for \(T^n\) and \(T^{n-1}\), which takes \(O(\log n)\) matrix multiplications. The memory usage is \(O(1)\).

Footnotes and References

1. Problem page: https://projecteuler.net/problem=624
2. Absorbing Markov chains: Wikipedia — Absorbing Markov chain
3. Matrix geometric series: Wikipedia — Matrix geometric series
4. Fibonacci numbers: Wikipedia — Fibonacci number
5. Modular multiplicative inverse: Wikipedia — Modular multiplicative inverse