Problem 624: Two Heads Are Better Than One

Project Euler Problem 624 Solution

EulerSolve provides an optimized solution for Project Euler Problem 624, Two Heads Are Better Than One, with C++, Python, Java, and a step-by-step mathematical explanation.

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....

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}.$$

$$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,$$

$$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

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

Problem 624 source code

C++

#include <cassert>
#include <cstdint>
#include <iostream>

using u64 = unsigned long long;
using u128 = __uint128_t;

struct Mat2 {
    u64 a00, a01, a10, a11;
};

static u64 mod_mul(u64 a, u64 b, u64 mod) { return (u128)a * b % mod; }

static u64 mod_pow(u64 a, u64 e, u64 mod) {
    u64 r = 1 % mod;
    while (e) {
        if (e & 1) r = mod_mul(r, a, mod);
        a = mod_mul(a, a, mod);
        e >>= 1;
    }
    return r;
}

static Mat2 mat_mul(const Mat2 &x, const Mat2 &y, u64 mod) {
    Mat2 r;
    r.a00 = (mod_mul(x.a00, y.a00, mod) + mod_mul(x.a01, y.a10, mod)) % mod;
    r.a01 = (mod_mul(x.a00, y.a01, mod) + mod_mul(x.a01, y.a11, mod)) % mod;
    r.a10 = (mod_mul(x.a10, y.a00, mod) + mod_mul(x.a11, y.a10, mod)) % mod;
    r.a11 = (mod_mul(x.a10, y.a01, mod) + mod_mul(x.a11, y.a11, mod)) % mod;
    return r;
}

static Mat2 mat_pow(Mat2 base, u64 e, u64 mod) {
    Mat2 r{1, 0, 0, 1};
    while (e) {
        if (e & 1) r = mat_mul(r, base, mod);
        base = mat_mul(base, base, mod);
        e >>= 1;
    }
    return r;
}

static u64 P_mod(u64 n, u64 mod) {
    assert(mod > 2);
    const u64 inv2 = (mod + 1) / 2;
    const Mat2 A{inv2, inv2, inv2, 0};

    const Mat2 An = mat_pow(A, n, mod);
    const Mat2 An1 = mat_pow(A, n - 1, mod);

    const u64 w0 = An1.a00;
    const u64 w1 = An1.a10;

    const u64 m00 = (1 + mod - An.a00) % mod;
    const u64 m01 = (0 + mod - An.a01) % mod;
    const u64 m10 = (0 + mod - An.a10) % mod;
    const u64 m11 = (1 + mod - An.a11) % mod;

    const u64 det = (mod_mul(m00, m11, mod) + mod - mod_mul(m01, m10, mod)) % mod;
    const u64 inv_det = mod_pow(det, mod - 2, mod);

    const u64 s1 = (mod_mul((mod - m10) % mod, w0, mod) + mod_mul(m00, w1, mod)) % mod;
    return mod_mul(inv2, mod_mul(s1, inv_det, mod), mod);
}

static u64 Q_P(u64 n, u64 mod) {
    u64 r = P_mod(n, mod);
    return (r == 0) ? mod : r;
}

int main() {
    assert(Q_P(2, 109) == 66);
    assert(Q_P(3, 109) == 46);
    std::cout << Q_P(1'000'000'000'000'000'000ULL, 1'000'000'009ULL) << "\n";
    return 0;
}

Python

def solve():
    MOD = 1000000009
    
    def mod_mul(a, b):
        return (a * b) % MOD
        
    def mat_mul(x, y):
        r00 = (mod_mul(x[0], y[0]) + mod_mul(x[1], y[2])) % MOD
        r01 = (mod_mul(x[0], y[1]) + mod_mul(x[1], y[3])) % MOD
        r10 = (mod_mul(x[2], y[0]) + mod_mul(x[3], y[2])) % MOD
        r11 = (mod_mul(x[2], y[1]) + mod_mul(x[3], y[3])) % MOD
        return (r00, r01, r10, r11)
        
    def mat_pow(base, e):
        r = (1, 0, 0, 1)
        cur = base
        while e > 0:
            if e & 1:
                r = mat_mul(r, cur)
            cur = mat_mul(cur, cur)
            e >>= 1
        return r
        
    def mod_pow(a, e):
        return pow(a, e, MOD)
        
    def P_mod(n):
        inv2 = (MOD + 1) // 2
        A = (inv2, inv2, inv2, 0)
        
        An = mat_pow(A, n)
        An1 = mat_pow(A, n - 1)
        
        w0 = An1[0]
        w1 = An1[2]
        
        m00 = (1 + MOD - An[0]) % MOD
        m01 = (0 + MOD - An[1]) % MOD
        m10 = (0 + MOD - An[2]) % MOD
        m11 = (1 + MOD - An[3]) % MOD
        
        det = (mod_mul(m00, m11) + MOD - mod_mul(m01, m10)) % MOD
        inv_det = mod_pow(det, MOD - 2)
        
        s1 = (mod_mul((MOD - m10) % MOD, w0) + mod_mul(m00, w1)) % MOD
        return mod_mul(inv2, mod_mul(s1, inv_det))
        
    def Q_P(n):
        r = P_mod(n)
        return MOD if r == 0 else r
        
    return str(Q_P(1000000000000000000))

if __name__ == '__main__':
    print(solve())

Java

public class Euler624 {
    static final long MOD = 1000000009L;

    static long modPow(long a, long e) {
        long r = 1 % MOD;
        a %= MOD;
        while (e > 0) {
            if ((e & 1) == 1)
                r = (r * a) % MOD;
            a = (a * a) % MOD;
            e >>= 1;
        }
        return r;
    }

    static class Mat2 {
        long a00, a01, a10, a11;

        Mat2(long a00, long a01, long a10, long a11) {
            this.a00 = a00;
            this.a01 = a01;
            this.a10 = a10;
            this.a11 = a11;
        }
    }

    static Mat2 matMul(Mat2 x, Mat2 y) {
        long r00 = (x.a00 * y.a00 % MOD + x.a01 * y.a10 % MOD) % MOD;
        long r01 = (x.a00 * y.a01 % MOD + x.a01 * y.a11 % MOD) % MOD;
        long r10 = (x.a10 * y.a00 % MOD + x.a11 * y.a10 % MOD) % MOD;
        long r11 = (x.a10 * y.a01 % MOD + x.a11 * y.a11 % MOD) % MOD;
        return new Mat2(r00, r01, r10, r11);
    }

    static Mat2 matPow(Mat2 base, long e) {
        Mat2 r = new Mat2(1, 0, 0, 1);
        Mat2 cur = base;
        while (e > 0) {
            if ((e & 1) == 1)
                r = matMul(r, cur);
            cur = matMul(cur, cur);
            e >>= 1;
        }
        return r;
    }

    static long pMod(long n) {
        long inv2 = (MOD + 1) / 2;
        Mat2 A = new Mat2(inv2, inv2, inv2, 0);

        Mat2 An = matPow(A, n);
        Mat2 An1 = matPow(A, n - 1);

        long w0 = An1.a00;
        long w1 = An1.a10;

        long m00 = (1 + MOD - An.a00) % MOD;
        long m01 = (0 + MOD - An.a01) % MOD;
        long m10 = (0 + MOD - An.a10) % MOD;
        long m11 = (1 + MOD - An.a11) % MOD;

        long det = (m00 * m11 % MOD + MOD - m01 * m10 % MOD) % MOD;
        long invDet = modPow(det, MOD - 2);

        long s1 = ((MOD - m10) % MOD * w0 % MOD + m00 * w1 % MOD) % MOD;
        return inv2 * (s1 * invDet % MOD) % MOD;
    }

    static long qP(long n) {
        long r = pMod(n);
        return (r == 0) ? MOD : r;
    }

    public static String solve() {
        return Long.toString(qP(1000000000000000000L));
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}