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

1. Problem page: https://projecteuler.net/problem=602
2. Eulerian number: Wikipedia - Eulerian number
3. Eulerian polynomial: Wikipedia - Eulerian polynomial
4. Geometric distribution: Wikipedia - Geometric distribution
5. Binomial distribution: Wikipedia - Binomial distribution