Problem 602: Product of Head Counts

Project Euler Problem 602 Solution

EulerSolve provides an optimized solution for Project Euler Problem 602, Product of Head Counts, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Alice and $n$ friends toss the same biased coin in a fixed cyclic order, with Alice always taking the first turn in each round. The game ends immediately when Alice gets her first Head. If the friends have accumulated head counts $H_1,\dots,H_n$ by that moment, then $$e(n,p)=\mathbb{E}\!\left[\prod_{i=1}^{n} H_i\right],\qquad p=\Pr(\text{Tail}).$$ For fixed $n$, this expectation turns out to be a polynomial in $p$. Writing the coefficient of $p^k$ as $c(n,k)$, the target is $$c(10^7,4\cdot 10^6)\pmod{10^9+7}.$$ Mathematical Approach The key is to replace the stochastic process by a stopping-time variable and then identify the resulting polynomial with an Eulerian polynomial. Step 1: Model Alice's stopping time Let $R$ be the number of Tails Alice sees before her first Head. Then $R$ follows a geometric distribution: $$\Pr(R=r)=p^r(1-p),\qquad r\ge 0.$$ If $R=r$, then Alice has produced exactly $r$ complete rounds before the final stopping toss. Therefore each friend has tossed exactly $r$ times when the game ends....

Detailed mathematical approach

Problem Summary

Alice and $n$ friends toss the same biased coin in a fixed cyclic order, with Alice always taking the first turn in each round. The game ends immediately when Alice gets her first Head. If the friends have accumulated head counts $H_1,\dots,H_n$ by that moment, then

$$e(n,p)=\mathbb{E}\!\left[\prod_{i=1}^{n} H_i\right],\qquad p=\Pr(\text{Tail}).$$

For fixed $n$, this expectation turns out to be a polynomial in $p$. Writing the coefficient of $p^k$ as $c(n,k)$, the target is

$$c(10^7,4\cdot 10^6)\pmod{10^9+7}.$$

Mathematical Approach

The key is to replace the stochastic process by a stopping-time variable and then identify the resulting polynomial with an Eulerian polynomial.

Step 1: Model Alice's stopping time

Let $R$ be the number of Tails Alice sees before her first Head. Then $R$ follows a geometric distribution:

$$\Pr(R=r)=p^r(1-p),\qquad r\ge 0.$$

If $R=r$, then Alice has produced exactly $r$ complete rounds before the final stopping toss. Therefore each friend has tossed exactly $r$ times when the game ends.

Step 2: Compute the friends' contribution for fixed $R=r$

Conditional on $R=r$, the $i$-th friend records a head count

$$H_i\mid R=r \sim \operatorname{Binomial}(r,1-p).$$

Hence

$$\mathbb{E}[H_i\mid R=r]=r(1-p).$$

Once $r$ is fixed, different friends depend on disjoint tosses and are conditionally independent, so the expected product factors:

$$\mathbb{E}\!\left[\prod_{i=1}^{n} H_i \,\middle|\, R=r\right]=\prod_{i=1}^{n}\mathbb{E}[H_i\mid R=r]=\bigl(r(1-p)\bigr)^n.$$

Step 3: Sum over all stopping times

Applying the law of total expectation gives

$$e(n,p)=\sum_{r=0}^{\infty}\Pr(R=r)\,\mathbb{E}\!\left[\prod_{i=1}^{n} H_i \,\middle|\, R=r\right].$$

Substituting the previous formulas yields

$$e(n,p)=\sum_{r=0}^{\infty} p^r(1-p)\bigl(r(1-p)\bigr)^n=(1-p)^{n+1}\sum_{r=0}^{\infty} r^n p^r.$$

The term $r=0$ is harmless because $n\ge 1$, but the expression still looks like an infinite power series. The next step shows that it actually collapses to a finite polynomial.

Step 4: Recognize the Eulerian polynomial identity

For every $n\ge 1$, the power-sum generating function satisfies

$$\sum_{r=0}^{\infty} r^n x^r=\frac{xA_n(x)}{(1-x)^{n+1}},$$

where $A_n(x)$ is the $n$-th Eulerian polynomial:

$$A_n(x)=\sum_{m=0}^{n-1}\left\langle {n \atop m} \right\rangle x^m.$$

Substituting $x=p$ into the series above gives

$$e(n,p)=(1-p)^{n+1}\cdot \frac{pA_n(p)}{(1-p)^{n+1}}=pA_n(p).$$

This is the decisive simplification: the apparent infinite series is exactly the degree-$n$ polynomial $pA_n(p)$.

Step 5: Extract the required coefficient

Expanding $pA_n(p)$ gives

$$pA_n(p)=\sum_{m=0}^{n-1}\left\langle {n \atop m} \right\rangle p^{m+1}.$$

Therefore the coefficient of $p^k$ is

$$c(n,k)=\left\langle {n \atop k-1} \right\rangle,\qquad 1\le k\le n,$$

and $c(n,k)=0$ outside that range. So the original probability question reduces to computing one Eulerian number.

Step 6: Use the explicit Eulerian-number formula

The implementations evaluate the Eulerian number through the alternating sum

$$\left\langle {n \atop m} \right\rangle=\sum_{j=0}^{m+1}(-1)^j\binom{n+1}{j}(m+1-j)^n,$$

with $m=k-1$. This formula is ideal for modular arithmetic because each term is a product of one binomial coefficient and one modular power.

Worked Example: $n=3$

For $n=3$, the Eulerian polynomial is

$$A_3(x)=1+4x+x^2.$$

Hence

$$e(3,p)=pA_3(p)=p+4p^2+p^3.$$

So the nonzero coefficients are

$$c(3,1)=1,\qquad c(3,2)=4,\qquad c(3,3)=1.$$

The same value $c(3,2)=4$ also comes directly from the explicit sum with $m=1$:

$$\left\langle {3 \atop 1} \right\rangle=\binom{4}{0}2^3-\binom{4}{1}1^3+\binom{4}{2}0^3=8-4+0=4.$$

This matches the small verification built into the implementation.

How the Code Works

The C++, Python, and Java implementations do not expand the whole polynomial. They compute only the single Eulerian number corresponding to $m=k-1$. First they prepare modular inverses of $1,2,\dots,m+1$, which is possible because the modulus $10^9+7$ is prime.

Next they sweep once through $j=0,1,\dots,m+1$. During that pass, the current binomial coefficient $\binom{n+1}{j}$ is updated from the previous one by the standard ratio for consecutive binomial terms, so no huge factorial tables are needed. In the same loop, the power $(m+1-j)^n$ is computed by fast modular exponentiation, multiplied by the current binomial factor, and added or subtracted according to the sign $(-1)^j$.

This strategy is exactly what the mathematics suggests: evaluate the alternating formula for one Eulerian number modulo $10^9+7$, and stop. That keeps the method practical even when $n=10^7$.

Complexity Analysis

To compute one coefficient $c(n,k)$, the alternating sum has $k+1$ terms because $m=k-1$. Preparing the modular inverses costs $O(k)$ time and $O(k)$ memory. Each summand also needs one modular exponentiation, which costs $O(\log n)$. Therefore the total running time is

$$O(k\log n),$$

with

$$O(k)$$

memory usage.

Footnotes and References

Problem page: https://projecteuler.net/problem=602
Eulerian number: Wikipedia - Eulerian number
Eulerian polynomial: Wikipedia - Eulerian polynomial
Geometric distribution: Wikipedia - Geometric distribution
Binomial distribution: Wikipedia - Binomial distribution

Problem 602 source code

C++

#include <cstdint>
#include <iostream>
#include <vector>

// Project Euler 602: Expected product of head-counts
//
// Let n be the number of friends (excluding Alice). The game ends when Alice gets her
// first Head. Let r be the number of Alice's Tails before that Head. Then:
//   P(r = t) = p^t (1-p),  t >= 0,
// and each friend makes exactly r tosses, with head probability (1-p).
//
// For a fixed r, each friend's head count H ~ Binomial(r, 1-p), so E[H | r] = r(1-p).
// The friends are independent given r, hence:
//   e(n,p) = E[ prod_i H_i ] = E[ (E[H | r])^n ] = sum_{r>=0} p^r (1-p) * (r(1-p))^n
//          = (1-p)^{n+1} * sum_{r>=0} r^n p^r.
//
// The known generating function identity is:
//   sum_{r>=0} r^n x^r = x * A_n(x) / (1-x)^{n+1},
// where A_n(x) is the Eulerian polynomial:
//   A_n(x) = sum_{m=0..n-1} <n,m> x^m,
// and <n,m> is the Eulerian number (permutations of {1..n} with m descents).
//
// Therefore:
//   e(n,p) = (1-p)^{n+1} * p * A_n(p) / (1-p)^{n+1} = p * A_n(p),
// so the coefficient c(n,k) of p^k in e(n,p) is:
//   c(n,k) = <n, k-1>.
//
// We compute Eulerian(n,m) modulo MOD using:
//   <n,m> = sum_{j=0..m+1} (-1)^j * C(n+1, j) * (m+1-j)^n.

using u64 = std::uint64_t;

static constexpr int MOD = 1000000007;

static inline int addmod(int a, int b) {
    int s = a + b;
    if (s >= MOD) s -= MOD;
    return s;
}
static inline int submod(int a, int b) {
    int s = a - b;
    if (s < 0) s += MOD;
    return s;
}
static inline int mulmod(long long a, long long b) {
    return (int)((a * b) % MOD);
}

static int mod_pow(int a, u64 e) {
    long long r = 1;
    long long x = a;
    while (e) {
        if (e & 1) r = (r * x) % MOD;
        x = (x * x) % MOD;
        e >>= 1;
    }
    return (int)r;
}

static int eulerian_mod(u64 n, u64 m) {
    // Requires 0 <= m <= n-1, n >= 1.
    const u64 N = n + 1; // binomial top
    const u64 JMAX = m + 1;

    std::vector<int> inv(JMAX + 2);
    if (JMAX + 1 >= 1) inv[1] = 1;
    for (u64 i = 2; i <= JMAX + 1; ++i) {
        inv[i] = (int)(MOD - (long long)(MOD / (int)i) * inv[MOD % (int)i] % MOD);
    }

    int ans = 0;
    long long binom = 1; // C(N,0)
    for (u64 j = 0; j <= JMAX; ++j) {
        const int base = (int)(JMAX - j); // (m+1-j)
        const int pown = mod_pow(base, n);
        const int term = mulmod((int)binom, pown);
        ans = (j & 1) ? submod(ans, term) : addmod(ans, term);

        if (j == JMAX) break;
        // C(N, j+1) = C(N,j) * (N-j)/(j+1).
        binom = binom * (long long)(N - j) % MOD;
        binom = binom * inv[j + 1] % MOD;
    }
    return ans;
}

static int c_mod(u64 n, u64 k) {
    // coefficient of p^k in e(n,p) is Eulerian(n,k-1)
    if (k == 0 || k > n) return 0;
    return eulerian_mod(n, k - 1);
}

int main() {
    // Statement examples/validations.
    if (c_mod(3, 1) != 1 || c_mod(3, 2) != 4 || c_mod(3, 3) != 1) {
        std::cerr << "Validation failed: e(3,p)\n";
        return 1;
    }
    if (c_mod(100, 40) != 986699437) {
        std::cerr << "Validation failed: c(100,40)\n";
        return 1;
    }

    std::cout << c_mod(10000000ULL, 4000000ULL) << "\n";
    return 0;
}

Python

MOD = 1000000007

def solve():
    n = 10000000
    k = 4000000
    m = k - 1
    N = n + 1
    JMAX = m + 1
    
    inv = [0] * (JMAX + 2)
    inv[1] = 1
    for i in range(2, JMAX + 2):
        inv[i] = MOD - (MOD // i) * inv[MOD % i] % MOD
        
    ans = 0
    binom = 1
    for j in range(JMAX + 1):
        base = JMAX - j
        pown = pow(base, n, MOD)
        term = (binom * pown) % MOD
        
        if j & 1:
            ans = (ans - term + MOD) % MOD
        else:
            ans = (ans + term) % MOD
            
        if j == JMAX: break
        binom = (binom * (N - j)) % MOD
        binom = (binom * inv[j + 1]) % MOD
        
    return str(ans)

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

Java

public class Euler602 {
    static final int MOD = 1000000007;

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

    public static String solve() {
        int n = 10000000;
        int k = 4000000;
        int m = k - 1;
        int N = n + 1;
        int JMAX = m + 1;

        int[] inv = new int[JMAX + 2];
        if (JMAX + 1 >= 1)
            inv[1] = 1;
        for (int i = 2; i <= JMAX + 1; i++) {
            inv[i] = MOD - (int) (((long) (MOD / i) * inv[MOD % i]) % MOD);
        }

        int ans = 0;
        long binom = 1;
        for (int j = 0; j <= JMAX; j++) {
            int base = JMAX - j;
            int pown = modPow(base, n);
            long term = (binom * pown) % MOD;
            if ((j & 1) == 1) {
                ans = (int) ((ans - term + MOD) % MOD);
            } else {
                ans = (int) ((ans + term) % MOD);
            }

            if (j == JMAX)
                break;
            binom = (binom * (N - j)) % MOD;
            binom = (binom * inv[j + 1]) % MOD;
        }

        return Integer.toString(ans);
    }

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