Problem 409: Nim Extreme

Project Euler Problem 409 Solution

EulerSolve provides an optimized solution for Project Euler Problem 409, Nim Extreme, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $W(n)$ denote the number of winning positions in the constrained $n$-heap Nim variant where each heap size is a distinct positive integer smaller than $2^n$. In ordinary Nim, a position is losing exactly when the bitwise xor of all heap sizes is zero, so the task is to count ordered $n$-tuples $(a_1,\dots,a_n)$ with $1 \le a_i \lt 2^n$, all $a_i$ different, and $a_1 \oplus \cdots \oplus a_n \ne 0$, then evaluate $W(10^7)\bmod 10^9+7$. Mathematical Approach Step 1: Rewrite the game in $\mathbb{F}_2^n$ Write each heap size in binary and identify it with a vector in $G=\mathbb{F}_2^n$. The admissible heap sizes are the nonzero vectors $G^\ast=G\setminus\{0\}$, so $|G^\ast|=2^n-1$. Because xor is addition in $\mathbb{F}_2^n$, a position is losing precisely when $$a_1+\cdots+a_n=0 \quad \text{in } \mathbb{F}_2^n.$$ This turns the Nim condition into a clean linear-algebraic constraint. Step 2: Count all admissible ordered tuples Set $q=2^n$. Ignoring the xor condition, we choose $n$ distinct nonzero vectors in order. Therefore the total number of admissible tuples is the falling factorial $$N_{\mathrm{all}}=(q-1)_n=\prod_{i=0}^{n-1}(q-1-i).$$ This is the count of all candidate positions before separating losing from winning ones....

Detailed mathematical approach

Problem Summary

Nim is played with several piles of stones. A move removes any positive number of stones from exactly one pile, and the player who cannot move, every pile being empty, loses. For a positive integer $n$, the problem looks only at positions of one special shape: exactly $n$ non-empty piles, each holding fewer than $2^n$ stones, no two of the same size. Let $W(n)$ be the number of such positions that are winning, meaning that the player who moves first can force a win. The statement gives $W(1)=1$, $W(2)=6$, $W(3)=168$ and $W(5)=19764360$, says that $W(100)$ reduced modulo $1\,000\,000\,007$ is $384777056$, and asks for ^[1]

$$W(10\,000\,000)\bmod 1\,000\,000\,007.$$

The piles are counted as distinguishable. For $n=2$ the allowed sizes are $1,2,3$, and the three sets of two different sizes, $\{1,2\}$, $\{1,3\}$ and $\{2,3\}$, are all winning, which would give $3$; the stated $W(2)=6$ is twice that. A position is therefore an ordered list $(a_1,\dots,a_n)$, and two lists that differ only in the order of the piles count twice.

The whole solution in one sentence: a nim position is lost by the player to move exactly when the bitwise xor of all pile sizes is zero, so $W(n)$ counts the ordered $n$-tuples of distinct non-zero $n$-bit numbers with non-zero xor, and rewriting the condition "xor is zero" as an average of signs turns that count into one binomial coefficient, which a program evaluates with $O(n)$ modular multiplications.

Mathematical Approach

1. Nim is decided by one xor

Write $a\oplus b$ for the bitwise exclusive or: compare the two numbers bit by bit and write $1$ where the bits differ, $0$ where they agree. For instance $3=011_2$ and $5=101_2$ give $3\oplus5=110_2=6$. Equal bits cancel, so $a\oplus a=0$. Call $s=a_1\oplus\cdots\oplus a_n$ the nim-sum of a position; Bouton proved in 1901 that the player to move loses exactly when $s=0$. ^[2] The restrictions on the sizes describe only the positions being counted: play follows the ordinary rules, so a move may empty a pile or make two piles equal. Two short lemmas give the whole proof.

From $s=0$ every move leads to $s\ne0$. A move replaces one size $a$ by some $a'\lt a$. The new nim-sum is $s\oplus a\oplus a'$, because $a$ leaves the xor and $a'$ enters it. With $s=0$ this is $a\oplus a'$, which is non-zero since $a\ne a'$.
From $s\ne0$ some move leads to $s=0$. Let $b$ be the highest set bit of $s$. At least one pile $a$ has bit $b$ set, otherwise that bit of $s$ would be $0$. Put $a'=a\oplus s$. Above bit $b$ nothing changes, and at bit $b$ a $1$ becomes a $0$, so $a'\lt a$ and removing $a-a'$ stones is a legal move. The new nim-sum is $s\oplus a\oplus a'=s\oplus s=0$.

Every move removes at least one stone, so the game is finite, and the empty position has $s=0$ and loses. Induction on the number of stones finishes the proof: from $s=0$ every move hands the opponent a position with $s\ne0$, a win for them, while from $s\ne0$ one move reaches $s=0$. For example $(3,4,5)$ has $s=2$, whose highest set bit is bit $1$; the pile $3$ has that bit, and taking two stones there gives $(1,4,5)$ with nim-sum $0$.

2. Pile sizes as vectors over the two-element field

Let $\mathbb F_2=\{0,1\}$ with addition modulo $2$, and let $V=\mathbb F_2^{\,n}$ be the space of bit strings of length $n$. Addition in $V$ is componentwise addition modulo $2$, which is exactly the bitwise xor. ^[3] Writing a number as its string of $n$ bits identifies $0,1,\dots,2^n-1$ with the $2^n$ elements of $V$, and the admissible pile sizes $1\le a\lt 2^n$ with the non-zero vectors

$$V^{*}=V\setminus\{0\},\quad |V^{*}|=2^n-1.$$

By Section 1, a position is losing for the player to move exactly when $a_1+\cdots+a_n=0$ in $V$. So $W(n)$ is the number of ordered $n$-tuples of distinct elements of $V^{*}$ whose sum is not $0$. Nothing about the game is left: this is a pure counting problem in a vector space of dimension $n$.

3. What has to be counted

Let $T$ count the ordered $n$-tuples of distinct elements of $V^{*}$, with no condition on the sum. The first pile can be any of the $2^n-1$ non-zero vectors, the second any of the remaining $2^n-2$, and so on, so $T$ is a falling factorial of length $n$ ^[4]

$$T=(2^n-1)(2^n-2)\cdots(2^n-n).$$

Let $N_0$ be the number of those tuples whose sum is $0$. Then $W(n)=T-N_0$. For $n=1$ a single non-zero vector is never $0$, and for $n=2$ the equation $a+b=0$ forces $a=b$, which is excluded; so $N_0=0$ in both cases, and $T=1$ and $T=6$ already give the stated $W(1)=1$ and $W(2)=6$. The work is entirely in $N_0$.

4. Replacing "sum is zero" by an average of signs

For $u,x\in V$ let $u\cdot x=u_1x_1+\cdots+u_nx_n$ in $\mathbb F_2$: it is $1$ exactly when the two vectors share an odd number of $1$ bits, so $(-1)^{u\cdot x}$ is $+1$ or $-1$. Try $n=2$: for $x=01$, letting $u$ run through $00,01,10,11$ gives signs $+1,-1,+1,-1$ adding up to $0$, while for $x=00$ all four are $+1$ and add up to $4$. In general the key identity is

$$\sum_{u\in V}(-1)^{u\cdot x}= \begin{cases} 2^n, & x=0,\\ 0, & x\ne 0. \end{cases}$$

If $x=0$ all $2^n$ terms equal $1$. If $x\ne0$, pick $i$ with $x_i=1$ and let $e_i$ be the vector with a single $1$ in position $i$. The map $u\mapsto u+e_i$ splits $V$ into $2^{n-1}$ pairs, and $(u+e_i)\cdot x=u\cdot x+1$, so the two terms of each pair cancel. The functions $x\mapsto(-1)^{u\cdot x}$ are the characters of $V$, and this is their orthogonality relation against the trivial character. ^[5]

Applied to $x=a_1+\cdots+a_n$, the identity gives $2^n$ exactly for the tuples we want and $0$ for all others. Dividing by $2^n$ and exchanging the summations gives

$$\begin{aligned} N_0&=\frac{1}{2^n}\sum_{u\in V}\Sigma(u),\\ \Sigma(u)&=\sum_{(a_1,\dots,a_n)}\ \prod_{i=1}^{n}(-1)^{u\cdot a_i}, \end{aligned}$$

where the inner sum runs over all $T$ ordered tuples of distinct non-zero vectors; the product form uses $u\cdot(a+b)=u\cdot a+u\cdot b$. For $u=0$ every factor is $1$, so $\Sigma(0)=T$. The gain is that $\Sigma(u)$ no longer mentions the sum of the tuple: each tuple contributes one product of signs, one per chosen size, and Section 6 is where distinctness is enforced.

5. Every non-zero $u$ splits the sizes in half

Fix $u\ne0$ and pick $i$ with $u_i=1$. The map $x\mapsto x+e_i$ is a bijection from $\{x\in V:u\cdot x=0\}$ onto $\{x\in V:u\cdot x=1\}$, so each of the two sets has $r=2^{n-1}$ elements. The zero vector lies in the first set, so among the $2^n-1$ admissible pile sizes exactly $r-1$ have sign $+1$ and exactly $r$ have sign $-1$. For $n=3$ and $u=001$ the sign is $+1$ for the sizes $2,4,6$, which is $r-1=3$ of them, and $-1$ for $1,3,5,7$, which is $r=4$.

6. From ordered tuples to a polynomial coefficient

In $\Sigma(u)$ the product $\prod_i(-1)^{u\cdot a_i}$ does not depend on the order of the tuple, and every $n$-element subset of $V^{*}$ can be ordered in $n!$ ways, so $\Sigma(u)=n!\,e_n$. Here $e_n$ is the elementary symmetric polynomial of degree $n$ in the $2^n-1$ signs: the sum, over all $n$-element subsets of $V^{*}$, of the product of their signs. Such polynomials are read off from a product of linear factors, ^[6]

$$\prod_{a\in V^{*}}\bigl(1+(-1)^{u\cdot a}x\bigr)=\sum_{k\ge0}e_kx^k,$$

because taking $(-1)^{u\cdot a}x$ from $k$ factors and $1$ from the rest contributes exactly the product of one $k$-element subset. By Section 5 the product has $r-1$ factors $(1+x)$ and $r$ factors $(1-x)$, so with $[x^n]$ for "coefficient of $x^n$",

$$\Sigma(u)=n!\,[x^n](1-x)^{r}(1+x)^{r-1}.$$

This depends only on the two counts, not on which $u\ne0$ was chosen, so all $2^n-1$ non-zero $u$ give the same value, which we call $C_n$.

7. The polynomial collapses

Pairing one $(1-x)$ with each $(1+x)$ uses the difference of two squares:

$$\begin{aligned} &(1-x)^{r}(1+x)^{r-1}\\ &\quad=(1-x)(1-x^2)^{r-1}. \end{aligned}$$

The binomial theorem gives

$$(1-x^2)^{r-1}=\sum_{j=0}^{r-1}(-1)^j\binom{r-1}{j}x^{2j},$$

a polynomial with even powers only. Multiplying by $1-x$ keeps those coefficients and adds their negatives one degree higher, so $x^{2m}$ keeps the coefficient of the even series while $x^{2m+1}$ takes its negative. With $m=\lfloor n/2\rfloor$, the coefficient of $x^n$ in $(1-x)(1-x^2)^{r-1}$ is therefore

$$\begin{cases} (-1)^{m}\binom{r-1}{m}, & n=2m,\\ (-1)^{m+1}\binom{r-1}{m}, & n=2m+1, \end{cases}$$

and both cases are the single expression $(-1)^{\lceil n/2\rceil}\binom{r-1}{m}$. Hence

$$C_n=(-1)^{\lceil n/2\rceil}\,n!\,\binom{r-1}{m}.$$

A more pedestrian route reaches the same number: choosing $j$ of the $n$ sizes among the $r$ negative ones and $n-j$ among the $r-1$ positive ones gives

$$e_n=\sum_{j=0}^{n}(-1)^j\binom{r}{j}\binom{r-1}{n-j},$$

with the convention that a binomial coefficient is $0$ when more elements are asked for than exist. The collapse replaces those $n+1$ terms by one.

8. The closed form for $W(n)$

Of the $2^n$ vectors $u$, one is $0$ and contributes $T$, and the other $2^n-1$ each contribute $C_n$. Section 4 therefore gives

$$\begin{aligned} N_0&=\frac{T+(2^n-1)\,C_n}{2^n},\\ W(n)&=T-N_0. \end{aligned}$$

Eliminating $N_0$ simplifies this to a single product:

$$W(n)=\frac{2^n-1}{2^n}\bigl(T-C_n\bigr).$$

Since $W(n)$ is an integer and the consecutive numbers $2^n-1$ and $2^n$ are coprime, the second form forces $2^n$ to divide $T-C_n$; direct checks confirm this for all $n\le20$.

9. A worked example and the first values

Take $n=3$. The admissible sizes are $1,\dots,7$, and $T=7\cdot6\cdot5=210$. Here $r=4$, so for every $u\ne0$ the generating product is $(1-x)^4(1+x)^3=(1-x)(1-x^2)^3$. Since $(1-x^2)^3=1-3x^2+3x^4-x^6$, multiplying by $1-x$ gives the coefficient $3$ at $x^3$, hence $C_3=3!\cdot3=18$, as the closed form $(-1)^2\cdot6\cdot\binom{3}{1}$ also gives. Therefore

$$\begin{aligned} N_0&=\frac{210+7\cdot18}{8}=42,\\ W(3)&=210-42=168, \end{aligned}$$

and the simplified form agrees: $\tfrac78(210-18)=168$. The count $N_0=42$ can be checked by hand: the triples of distinct sizes with xor $0$ are $\{1,2,3\}$, $\{1,4,5\}$, $\{1,6,7\}$, $\{2,4,6\}$, $\{2,5,7\}$, $\{3,4,7\}$ and $\{3,5,6\}$ — the seven lines of the Fano plane — and $7\cdot3!=42$.

The same computation at $n=4$ gives $C_4=4!\binom72=504$, $T=32760$, $N_0=2520$ and $W(4)=30240$, where $2520$ is $105$ losing sets times $4!$ and those $105$ sets were counted directly; at $n=5$ it gives $C_5=-5!\binom{15}{2}=-12600$, $T=20389320$, $N_0=624960$ and $W(5)=19764360$, the value quoted in the statement.

10. Evaluating the formula modulo a prime

Building the exact integers is out of the question at $n=10^7$: each of the $n$ factors of $T$ has about three million decimal digits, so $T$ itself has some $3\cdot10^{13}$. Only remainders are wanted, so all arithmetic runs modulo the prime $p=10^9+7$. Three ingredients are needed.

Reduce the factors before multiplying. Taking remainders commutes with addition, subtraction and multiplication, so reducing after every step leaves the final remainder unchanged. Each factor $2^n-i$ of $T$ is reduced first: compute $2^n\bmod p$ by repeated squaring, subtract $i$, and multiply the $n$ reduced factors modulo $p$.
Division becomes multiplication by an inverse. For $a$ not divisible by $p$, Fermat's little theorem gives $a^{p-1}\equiv1\pmod p$, so $a^{p-2}$ is the inverse of $a$. ^[7] This is how $m!$ and $2^n$ are divided out; $2^n\bmod p$ is never $0$ because $p$ is odd.
A binomial coefficient with an enormous top. The top $M=r-1=2^{n-1}-1$ never has to be built: its remainder is enough, and only $m=\lfloor n/2\rfloor$ factors are needed:

$$\binom{M}{m}=\frac{M(M-1)\cdots(M-m+1)}{m!}.$$

Reducing $M$ modulo $p$ and multiplying $m$ consecutive reduced factors gives the numerator; the inverse of $m!$ finishes it. The identity $\binom{M}{m}\,m!=M(M-1)\cdots(M-m+1)$ holds over the integers, and $m!$ is invertible modulo $p$ because $m=5\cdot10^6\lt p$. It is the only factorial we ever divide by: $n!$ is multiplied in, never inverted.

Lucas's theorem, the standard tool for binomial coefficients modulo a prime, agrees with this shortcut rather than replacing it. ^[8] It multiplies binomial coefficients built from the base-$p$ digits of the two arguments; since $m\lt p$, every digit of $m$ above the last is $0$ and those factors are $1$, leaving

$$\binom{M}{m}\equiv\binom{M\bmod p}{m}\pmod p.$$

The falling product computes exactly that, including the case $M\bmod p\lt m$, where it picks up a zero factor and both sides vanish. We checked all three expressions against each other for the primes $7,11,13,101$.

How the Code Works

All three programs implement Sections 8 and 10 directly, with no array and no search. ^[9] In C++ Options holds $n$ and whether to run the built-in checks, while parse_arguments and parse_int_after_prefix handle --n= and --skip-checkpoints, reject an unknown option and require $n\ge1$. Their std::stoi call does not check that the whole tail was consumed, so a malformed --n=3junk is silently accepted as $3$. mod_pow is square-and-multiply exponentiation modulo MOD, used for $2^n$ and for the Fermat inverses.

The mathematics all sits in W_mod(n), in the order of the derivation. It first reduces pow2n $=2^n$ and total_base $=2^n-1$. One loop of $n$ steps multiplies T by total_base - i for $i=0,\dots,n-1$, giving $T$ of Section 3. A second loop to $n$ builds fact_n $=n!$ and, on the way, fact_m $=m!$ for m = n / 2, inverted as inv_fact_m. A third loop of $m$ steps multiplies falling_q by q_minus_1 - i, with q_minus_1 $=2^{n-1}-1$, so that comb $=$ falling_q * inv_fact_m is $\binom{r-1}{m}$. The parity test of Section 7 supplies the sign — negate when $n$ is even and $m$ odd, or $n$ odd and $m$ even, which is $(-1)^{\lceil n/2\rceil}$ — and fact_n turns comb into S $=C_n$. Then losing is $(T+(2^n-1)S)$ times inv_pow2n, that is $N_0$, and winning $=$ T - losing is returned. The additions of MOD before each % keep intermediates non-negative, since total_base - i and q_minus_1 - i turn negative once $i$ passes the reduced residue.

The counting residues in W_mod are kept in int, and each multiplication starts with a 1LL * cast so that the product is formed in long long: a product of two residues approaches $10^{18}$ and would overflow $32$ bits, and the expression for losing keeps the widened intermediate below $2^{63}$. Inside mod_pow no cast is needed, as result and cur already use the long long alias i64. run_checkpoints recomputes $W(1)$, $W(2)$, $W(3)$, $W(5)$ and $W(100)\bmod p$ and compares them with the statement; on a mismatch it writes a message and returns false, whereupon main exits with status $2$. Otherwise main prints $W(n)\bmod p$.

The Java port is a line-by-line translation: modPow, the same three loops, the same sign test and the same 1L * widening, with counting residues in int and products in long, while modPow keeps its own residues in long. Python is shorter because the language does the bookkeeping: pow(base, exp, MOD) is built in, integers are unbounded, and the remainder operator already returns a non-negative value, so total_base - i needs no correction. Both ports compute in solve, fix $n=10^7$ there and return a decimal string, printed from main in Java and under the __main__ guard in Python; neither has checkpoints, so the five reference values are verified only by the C++ program.

Complexity and Verification

For one value of $n$ the three loops perform $n$, $n+m$ and $m$ modular multiplications, where $m=\lfloor n/2\rfloor$; the total $2n+2m$ is $3n$ for even $n$ and $3n-1$ for odd $n$, so $3\cdot10^7$ at the target. Added to that are a few final products, four modular exponentiations of at most $30$ squarings and, unless skipped, the five checkpoints. The running time is $O(n)$ and the auxiliary storage $O(1)$: a fixed number of integer variables, no array at all. Measured on one machine, with the usual dependence on hardware, compiler and runtime, the C++ program needs $0.09$ s including the checkpoints, the Java port $0.13$ s and Python $0.9$ s, mostly the interpreted loop. Reading $\binom{M\bmod p}{m}$ off a factorial table would need $10^9$ entries; the falling product uses only $m$ of them, which keeps the cost linear in $n$.

Beyond the five checkpoints of the C++ program, the following was verified for this article, outside the published sources:

An exhaustive game solver that recurses over all legal moves, knowing nothing about xor, agrees with Section 1 on all $216$ ordered triples of sizes in $\{0,1,2,3,4,5\}$, empty piles included.
A dynamic programme over subsets, counting for every subset size and partial xor how many subsets of $V^{*}$ realise it, reproduces $T$, $N_0$ and $W(n)$ exactly for all $n\le10$ and matches Section 8 every time, $W(1)$, $W(2)$, $W(3)$ and $W(5)$ among them.
Direct enumeration of all ordered tuples at $n=3,4,5$ confirms $\Sigma(0)=T$ and $\Sigma(u)=C_n$ for every $u\ne0$; the split of Section 5 holds for every non-zero $u$ at $n=3,4,8$.
Exact polynomial multiplication confirms the identity of Section 7, the coefficient and the alternating sum for all $n\le12$, and $2^n$ divides $T-C_n$ for all $n\le20$.
An implementation written only from the formula reproduces $W(100)\bmod p=384777056$ and agrees with the published C++ program at $n=1,\dots,7,10,100,1000$ and $12345$.
The three published programs, run exactly as they stand, produce the same value at $n=10^7$; their outputs were compared through cryptographic hashes, so the value never appeared.
With $p=7$ and $p=11$ in place of $10^9+7$, the published arithmetic reproduces $W(n)\bmod p$ for all $n\le14$, which covers $n\ge p$ as well.

That last check also settles what happens beyond this problem's range: a winning position has $n!$ orderings of its distinct sizes, so $n!$ divides $W(n)$ and $W(n)\equiv0\pmod p$ once $n\ge p$ — exactly what the program returns there, both $T$ and its factor $n!$ reducing to $0$.

Footnotes and References

The references below supply the background facts used above. The xor criterion is proved in Section 1, and the character identity, the polynomial collapse and the closed form for $W(n)$ in Sections 4 to 8.

Project Euler 409 — Nim Extreme. The official statement considers nim positions with n non-empty piles, each of size less than 2ⁿ and no two of the same size, defines W(n) as the number of them that are winning for the first player, lists W(1) = 1, W(2) = 6, W(3) = 168, W(5) = 19764360 and W(100) mod 1000000007 = 384777056, and asks for W(10000000) modulo 1000000007.
Wikipedia — Nim. The article states that the nim-sum, the binary digital sum of the heap sizes, is the key to the theory and that the first player has a winning strategy if and only if it is not zero; it credits Charles L. Bouton with the complete theory in 1901. MathWorld — Nim cites Bouton's paper, "Nim, A Game with a Complete Mathematical Theory", Annals of Mathematics 3, pages 35–39, 1901–1902.
Wikipedia — GF(2). The article states that addition in the two-element field is the logical XOR operation and that fixed-length bit strings carry a vector-space structure over GF(2) whose addition is the bitwise exclusive or, the identification used in Section 2.
Wikipedia — Falling and rising factorials. The article defines the falling factorial as the product x(x − 1)⋯(x − n + 1) of n factors, written (x)ₙ, with the value 1 for n = 0: the shape of T in Section 3 and of the binomial numerator in Section 10.
Wikipedia — Character group. The article states that a finite abelian group of order n has exactly n distinct characters and proves the orthogonality relation AA* = A*A = nI, whose column sums give the identity of Section 4: the group order at the identity element, 0 elsewhere. Section 4 proves that case directly.
Wikipedia — Elementary symmetric polynomial. The article defines eₖ as the sum of all products of k distinct variables and expands ∏(T − Xⱼ) = Tⁿ − e₁Tⁿ⁻¹ + ⋯ + (−1)ⁿeₙ; substituting T = −1/x and multiplying by (−x)ⁿ gives the form ∏(1 + Xⱼx) = Σ eₖxᵏ of Section 6.
Wikipedia — Fermat's little theorem. The article states that aᵖ − a is a multiple of p for a prime p and any integer a, and that for a not divisible by p this is equivalent to aᵖ⁻¹ ≡ 1 (mod p); multiplying by the inverse of a shows that aᵖ⁻² is that inverse, which is how Section 10 divides by m! and by 2ⁿ.
Wikipedia — Lucas's theorem. The theorem writes a binomial coefficient modulo a prime p as a product of binomial coefficients built from the base-p digits of its two arguments. Since the lower argument m is smaller than p, only the lowest digit of the upper argument survives, which is the agreement shown in Section 10.
C++ — Euler409.cpp. The linked immutable revision contains mod_pow, W_mod, the argument parsing and run_checkpoints with the five reference values; it is the implementation record for this article. The extra checks listed above are not packaged in it.

Problem 409 source code

C++

#include <cstdint>
#include <iostream>
#include <string>

namespace {

using i64 = long long;
using u64 = std::uint64_t;

constexpr int MOD = 1000000007;

struct Options {
    int n = 10000000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        value = std::stoi(tail);
    } catch (...) {
        return false;
    }
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1;
}

int mod_pow(i64 base, i64 exp) {
    i64 result = 1;
    i64 cur = base % MOD;
    i64 e = exp;
    while (e > 0) {
        if (e & 1LL) {
            result = result * cur % MOD;
        }
        cur = cur * cur % MOD;
        e >>= 1LL;
    }
    return static_cast<int>(result);
}

int W_mod(const int n) {
    const int pow2n = mod_pow(2, n);
    const int total_base = (pow2n - 1 + MOD) % MOD;  // 2^n - 1

    // Total ordered distinct non-zero tuples: T = (2^n-1)_n.
    int T = 1;
    for (int i = 0; i < n; ++i) {
        T = static_cast<int>(1LL * T * ((total_base - i + MOD) % MOD) % MOD);
    }

    const int m = n / 2;

    int fact_n = 1;
    int fact_m = 1;
    for (int i = 1; i <= n; ++i) {
        fact_n = static_cast<int>(1LL * fact_n * i % MOD);
        if (i <= m) {
            fact_m = static_cast<int>(1LL * fact_m * i % MOD);
        }
    }
    const int inv_fact_m = mod_pow(fact_m, MOD - 2);

    const int q_minus_1 = (mod_pow(2, n - 1) - 1 + MOD) % MOD;  // 2^(n-1)-1
    int falling_q = 1;
    for (int i = 0; i < m; ++i) {
        falling_q = static_cast<int>(1LL * falling_q * ((q_minus_1 - i + MOD) % MOD) % MOD);
    }
    int comb = static_cast<int>(1LL * falling_q * inv_fact_m % MOD);  // C(2^(n-1)-1, floor(n/2))

    // Coefficient [x^n] of (1-x)^(2^(n-1)) (1+x)^(2^(n-1)-1):
    // if n even:  (+/-) C(2^(n-1)-1, n/2) with sign (-1)^(n/2)
    // if n odd :  (+/-) C(2^(n-1)-1, (n-1)/2) with sign -(-1)^((n-1)/2)
    if ((n & 1) == 0) {
        if ((m & 1) == 1) {
            comb = (MOD - comb) % MOD;
        }
    } else {
        if (((m) & 1) == 0) {
            comb = (MOD - comb) % MOD;
        }
    }

    const int S = static_cast<int>(1LL * fact_n * comb % MOD);

    const int inv_pow2n = mod_pow(pow2n, MOD - 2);
    const int losing = static_cast<int>(
        1LL * ((T + 1LL * (pow2n - 1 + MOD) % MOD * S) % MOD) * inv_pow2n % MOD);
    const int winning = (T - losing + MOD) % MOD;
    return winning;
}

bool run_checkpoints() {
    if (W_mod(1) != 1) {
        std::cerr << "Checkpoint failed: W(1)\n";
        return false;
    }
    if (W_mod(2) != 6) {
        std::cerr << "Checkpoint failed: W(2)\n";
        return false;
    }
    if (W_mod(3) != 168) {
        std::cerr << "Checkpoint failed: W(3)\n";
        return false;
    }
    if (W_mod(5) != 19764360) {
        std::cerr << "Checkpoint failed: W(5)\n";
        return false;
    }
    if (W_mod(100) != 384777056) {
        std::cerr << "Checkpoint failed: W(100)\n";
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << W_mod(options.n) << '\n';
    return 0;
}

Python

MOD = 1000000007

def mod_pow(base, exp):
    return pow(base, exp, MOD)

def solve():
    n = 10000000
    pow2n = mod_pow(2, n)
    total_base = (pow2n - 1) % MOD
    
    T = 1
    for i in range(n):
        T = (T * (total_base - i)) % MOD
        
    m = n // 2
    
    fact_n = 1
    fact_m = 1
    for i in range(1, n + 1):
        fact_n = (fact_n * i) % MOD
        if i <= m:
            fact_m = (fact_m * i) % MOD
            
    inv_fact_m = mod_pow(fact_m, MOD - 2)
    
    q_minus_1 = (mod_pow(2, n - 1) - 1) % MOD
    falling_q = 1
    for i in range(m):
        falling_q = (falling_q * (q_minus_1 - i)) % MOD
        
    comb = (falling_q * inv_fact_m) % MOD
    
    if (n & 1) == 0:
        if (m & 1) == 1:
            comb = (MOD - comb) % MOD
    else:
        if (m & 1) == 0:
            comb = (MOD - comb) % MOD
            
    S = (fact_n * comb) % MOD
    
    inv_pow2n = mod_pow(pow2n, MOD - 2)
    losing = ((T + (total_base * S) % MOD) * inv_pow2n) % MOD
    winning = (T - losing + MOD) % MOD
    
    return str(winning)

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

Java

public class Euler409 {
    private static final int MOD = 1000000007;

    private static int modPow(long base, long exp) {
        long result = 1;
        long cur = base % MOD;
        long e = exp;
        while (e > 0) {
            if ((e & 1L) != 0) {
                result = (result * cur) % MOD;
            }
            cur = (cur * cur) % MOD;
            e >>= 1L;
        }
        return (int) result;
    }

    public static String solve() {
        int n = 10000000;
        int pow2n = modPow(2, n);
        int totalBase = (pow2n - 1 + MOD) % MOD;

        int T = 1;
        for (int i = 0; i < n; ++i) {
            T = (int) (1L * T * ((totalBase - i + MOD) % MOD) % MOD);
        }

        int m = n / 2;

        int factN = 1;
        int factM = 1;
        for (int i = 1; i <= n; ++i) {
            factN = (int) (1L * factN * i % MOD);
            if (i <= m) {
                factM = (int) (1L * factM * i % MOD);
            }
        }
        int invFactM = modPow(factM, MOD - 2);

        int qMinus1 = (modPow(2, n - 1) - 1 + MOD) % MOD;
        int fallingQ = 1;
        for (int i = 0; i < m; ++i) {
            fallingQ = (int) (1L * fallingQ * ((qMinus1 - i + MOD) % MOD) % MOD);
        }
        int comb = (int) (1L * fallingQ * invFactM % MOD);

        if ((n & 1) == 0) {
            if ((m & 1) == 1) {
                comb = (MOD - comb) % MOD;
            }
        } else {
            if ((m & 1) == 0) {
                comb = (MOD - comb) % MOD;
            }
        }

        int S = (int) (1L * factN * comb % MOD);

        int invPow2n = modPow(pow2n, MOD - 2);
        int losing = (int) (1L * ((T + 1L * totalBase * S) % MOD) * invPow2n % MOD);
        int winning = (T - losing + MOD) % MOD;

        return String.valueOf(winning);
    }

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