Project Euler Problem 602

Alice enlists the help of some friends to generate a random number, using a single unfair coin.

Project Euler Problem 602

Solution

Answer: 269496760

Let $q=1-p$ be the probability of Heads.

We seek the coefficient $c(n,k)$ of $p^k$ in the polynomial $e(n,p)$.

Step 1: Deriving the expectation

Suppose Alice has $n$ friends.

The game proceeds cyclically:

  • Alice tosses,
  • then each friend tosses once,
  • repeat until Alice gets a Head.

Assume Alice gets her first Head after exactly $m$ previous failures.

Then:

  • Alice produced $m$ Tails followed by one Head,
  • probability of this event:

$$p^m q.$$

During those $m$ failed rounds, each friend tossed exactly $m$ times independently.

For any friend $i$,

$$X_i \sim \text{Binomial}(m,q),$$

so

$$\mathbb E[X_i]=mq.$$

Since the friends are independent,

$$\mathbb E!\left[\prod_{i=1}^n X_i ,\middle|, m\right] = (mq)^n.$$

Therefore

$$e(n,p) = \sum_{m\ge0} p^m q (mq)^n = q^{n+1}\sum_{m\ge0} m^n p^m.$$

Because $q=1-p$,

$$e(n,p) = (1-p)^{n+1}\sum_{m\ge0} m^n p^m.$$


Step 2: Eulerian polynomials

A classical identity states:

$$\sum_{m\ge0} m^n x^m = \frac{A_n(x)}{(1-x)^{n+1}},$$

where $A_n(x)$ is the Eulerian polynomial.

Hence

$$e(n,p)=A_n(p).$$

Thus the coefficients $c(n,k)$ are precisely Eulerian numbers:

$$c(n,k)=\left\langle {n \atop k} \right\rangle .$$

A standard explicit formula is

$$\boxed{ c(n,k) = \sum_{i=0}^{k} (-1)^i \binom{n+1}{i} (k-i)^n }.$$

We are asked for

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


Step 3: Efficient modular computation

Since

$$10^9+7$$

is prime and $n+1 < 10^9+7$, we can compute all binomial coefficients iteratively modulo $M=10^9+7$.

The computation is:

MOD = 10**9 + 7

n = 10_000_000
k = 4_000_000

# modular inverses up to k+1
inv = [0] * (k + 2)
inv[1] = 1
for i in range(2, k + 2):
    inv[i] = MOD - (MOD // i) * inv[MOD % i] % MOD

ans = 0
binom = 1   # C(n+1, 0)

for i in range(k + 1):

    term = pow(k - i, n, MOD)

    if i & 1:
        ans -= binom * term
    else:
        ans += binom * term

    ans %= MOD

    # Update:
    # C(n+1, i+1) from C(n+1, i)
    binom = (
        binom
        * (n + 1 - i)
        % MOD
        * inv[i + 1]
        % MOD
    )

print(ans % MOD)

Step 4: Verification on the small example

For $n=3$,

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

Our formula gives:

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

matching the statement perfectly.


Step 5: Final computation

Evaluating the formula modulo $10^9+7$ gives

$$c(10000000,4000000) \equiv 269496760 \pmod{10^9+7}.$$

Answer: 269496760