Project Euler

Project Euler Problem 409: Nim Extreme

Consider nim positions where: 1. There are n non-empty piles. 2. Each pile has size less than 2^n. 3. No two piles share the same size. Define W(n) as the number of winning positions (the first pla...

Source sync Apr 19, 2026

Problem #0409

Level Level 22

Solved By 535

Languages C++, Python

Answer 253223948

Length 477 words

modular_arithmeticgame_theorybrute_force

Let $n$ be a positive integer. Consider nim positions where:

There are $n$ non-empty piles.
Each pile has size less than $2^n$.
No two piles have the same size.

Let $W(n)$ be the number of winning nim positions satisfying the above conditions (a position is winning if the first player has a winning strategy). For example, $W(1) = 1$, $W(2) = 6$, $W(3) = 168$ and $W(5) = 19764360$ and $W(100) \mod \num {1000000007} = 384777056$.

Find $W(\num {10000000} \mod \num {1000000007})$.

Project Euler Problem 409: Nim Extreme

Solution

Answer

W(10,000,000) mod 10^9 + 7 = 253223948

Mathematical Derivation

Setting up the counting. The heap sizes are distinct values from {1, 2, …, 2^n - 1}, and the piles are distinguishable (ordered). A position is losing if and only if the XOR (nim-sum) of all heap sizes is zero.

Let L(n) be the number of unordered n-element subsets of {1, …, 2^n - 1} with XOR equal to 0. Then:

W(n) = n! * (C(2^n - 1, n) - L(n))

Computing L(n) via Fourier analysis on GF(2)^n. The set {1, …, 2^n - 1} is exactly the set of nonzero elements of the vector space GF(2)^n.

Using the orthogonality of characters on GF(2)^n:

L(n) = (1 / 2^n) * [ C(2^n - 1, n) + (2^n - 1) * E_n ]

where

E_n = sum_{j=0}^{n} C(p, j) * C(q, n-j) * (-1)^{n-j}

with p = 2^{n-1} - 1 and q = 2^{n-1}. This comes from the fact that for every nonzero character a in GF(2)^n, exactly p = 2^{n-1} - 1 nonzero elements v satisfy a . v = 0 and q = 2^{n-1} satisfy a . v = 1.

Simplifying E_n via generating functions. The key insight is that E_n is the coefficient [x^n] of the product:

(1+x)^p * (1-x)^q = (1+x)^{2^{n-1}-1} * (1-x)^{2^{n-1}}
                   = (1-x) * (1-x^2)^{2^{n-1}-1}

Setting m = 2^{n-1} - 1 and h = floor(n/2):

n even: E_n = (-1)^h * C(m, h)
n odd: E_n = (-1)^{(n+1)/2} * C(m, (n-1)/2)

Final formula for W(n). Substituting and simplifying:

W(n) = n! * (2^n - 1) / 2^n * (C(2^n - 1, n) - E_n)   (mod 10^9 + 7)

where division by 2^n is performed via modular inverse.

Complexity

Time: O(n) — one loop for factorial and falling factorial, plus O(log M) for modular exponentiations.
Space: O(1) — only a constant number of accumulators.

Verification

n	W(n)	Matches Expected
1	1	Yes
2	6	Yes
3	168	Yes
5	19,764,360	Yes
100 mod 10^9+7	384,777,056	Yes

Correctness

Theorem. The method described above computes exactly the quantity requested in the problem statement.

Proof. The preceding analysis identifies the admissible objects and derives the formula, recurrence, or exhaustive search carried out by the algorithm. The computation evaluates exactly that specification, so every valid contribution is included once and no invalid contribution is counted. Therefore the returned value is the required answer. $\square$

Answer

\boxed{253223948}

C++ project_euler/problem_409/solution.cpp

#include <bits/stdc++.h>
using namespace std;

/*
 * Project Euler Problem 409: Nim Extreme
 *
 * W(n) = n! * (2^n - 1) / 2^n * (C(2^n-1, n) - E_n)  mod M
 *
 * where E_n = (-1)^h * C(2^{n-1}-1, h),  h = n/2      (n even)
 *       E_n = (-1)^{(n+1)/2} * C(2^{n-1}-1, (n-1)/2)  (n odd)
 *
 * All modular arithmetic with M = 10^9 + 7.
 */

static const long long M = 1000000007LL;

long long power(long long base, long long exp, long long mod) {
    long long result = 1;
    base %= mod;
    if (base < 0) base += mod;
    while (exp > 0) {
        if (exp & 1) result = result * base % mod;
        base = base * base % mod;
        exp >>= 1;
    }
    return result;
}

long long modinv(long long a, long long mod) {
    return power(a, mod - 2, mod);
}

long long solve(int n) {
    long long pow2n  = power(2, n, M);
    long long pow2n1 = power(2, n - 1, M);
    long long Nm1    = (pow2n - 1 + M) % M;    // 2^n - 1
    long long m_mod  = (pow2n1 - 1 + M) % M;   // 2^{n-1} - 1

    // Compute n! and C(2^n-1, n) = falling_factorial(Nm1, n) / n!
    long long fact_n = 1, falling_Nm1 = 1;
    for (int i = 0; i < n; i++) {
        fact_n = fact_n * ((i + 1) % M) % M;
        falling_Nm1 = falling_Nm1 % M * ((Nm1 - i % M + M) % M) % M;
    }
    long long inv_fact_n = modinv(fact_n, M);
    long long C_Nm1_n = falling_Nm1 % M * inv_fact_n % M;

    // Compute C(m, h) where m = 2^{n-1}-1, h = n/2
    int h = n / 2;
    long long fact_h = 1, falling_m = 1;
    for (int i = 0; i < h; i++) {
        fact_h = fact_h * ((i + 1) % M) % M;
        falling_m = falling_m % M * ((m_mod - i % M + M) % M) % M;
    }
    long long inv_fact_h = modinv(fact_h, M);
    long long C_m_h = falling_m % M * inv_fact_h % M;

    // E_n with sign
    int sign;
    if (n % 2 == 0) {
        sign = (h % 2 == 0) ? 1 : -1;
    } else {
        sign = (((n + 1) / 2) % 2 == 0) ? 1 : -1;
    }
    long long E_n_mod = (sign == 1) ? C_m_h : (M - C_m_h) % M;

    // W(n) = n! * (2^n - 1) * inv(2^n) * (C(2^n-1,n) - E_n) mod M
    long long inv_pow2n = modinv(pow2n, M);
    long long W = fact_n;
    W = W % M * (Nm1 % M) % M;
    W = W % M * (inv_pow2n % M) % M;
    W = W % M * ((C_Nm1_n - E_n_mod + M) % M) % M;
    return W;
}

int main() {
    ios_base::sync_with_stdio(false);

    // Verification
    struct { int n; long long expected; } tests[] = {
        {1, 1}, {2, 6}, {3, 168}, {5, 19764360}, {100, 384777056}
    };
    cout << "=== Verification ===" << endl;
    for (auto& t : tests) {
        long long w = solve(t.n);
        cout << "  W(" << t.n << ") mod M = " << w
             << "  expected " << t.expected
             << "  [" << (w == t.expected ? "OK" : "FAIL") << "]" << endl;
    }

    // Main computation
    cout << "\n=== Main computation ===" << endl;
    auto t0 = chrono::steady_clock::now();
    long long answer = solve(10000000);
    auto t1 = chrono::steady_clock::now();
    double elapsed = chrono::duration<double>(t1 - t0).count();

    cout << "  W(10,000,000) mod 10^9+7 = " << answer << endl;
    cout << "  Time: " << fixed << setprecision(2) << elapsed << "s" << endl;

    return 0;
}

Python project_euler/problem_409/solution.py

"""
Project Euler Problem 409: Nim Extreme
======================================
n piles with distinct sizes from {1, ..., 2^n - 1}. W(n) counts winning
positions (XOR != 0) for the first player, with piles ordered.

Formula:
  W(n) = n! * (2^n - 1) / 2^n * (C(2^n-1, n) - E_n)   mod M

where E_n = (-1)^h * C(2^{n-1}-1, h),  h = n//2   (n even)
      E_n = (-1)^{(n+1)/2} * C(2^{n-1}-1, (n-1)/2) (n odd)

All arithmetic mod M = 10^9 + 7.
"""

import time
import sys

M = 10**9 + 7

def solve(n):
    """Compute W(n) mod M in O(n) time."""
    pow2n = pow(2, n, M)
    pow2n1 = pow(2, n - 1, M)
    Nm1 = (pow2n - 1) % M          # 2^n - 1 mod M
    m_mod = (pow2n1 - 1) % M        # 2^{n-1} - 1 mod M

    # Compute n! mod M  and  C(2^n-1, n) mod M simultaneously
    fact_n = 1
    falling_Nm1 = 1
    for i in range(n):
        fact_n = fact_n * (i + 1) % M
        falling_Nm1 = falling_Nm1 * ((Nm1 - i) % M) % M
    inv_fact_n = pow(fact_n, M - 2, M)
    C_Nm1_n = falling_Nm1 * inv_fact_n % M

    # Compute C(m, h) mod M where m = 2^{n-1}-1, h = n//2
    h = n // 2
    fact_h = 1
    falling_m = 1
    for i in range(h):
        fact_h = fact_h * (i + 1) % M
        falling_m = falling_m * ((m_mod - i) % M) % M
    inv_fact_h = pow(fact_h, M - 2, M)
    C_m_h = falling_m * inv_fact_h % M

    # E_n with sign
    if n % 2 == 0:
        sign = 1 if h % 2 == 0 else -1
    else:
        sign = 1 if ((n + 1) // 2) % 2 == 0 else -1
    E_n_mod = sign * C_m_h % M

    # W(n) = n! * (2^n - 1) * inv(2^n) * (C(2^n-1,n) - E_n) mod M
    inv_pow2n = pow(pow2n, M - 2, M)
    W = fact_n * Nm1 % M * inv_pow2n % M * ((C_Nm1_n - E_n_mod) % M) % M
    return W

def brute_force(n):
    """Brute-force W(n) for small n (for verification)."""
    from itertools import combinations
    from math import factorial
    max_val = (1 << n) - 1
    winning = 0
    for combo in combinations(range(1, max_val + 1), n):
        xor_val = 0
        for v in combo:
            xor_val ^= v
        if xor_val != 0:
            winning += 1
    return winning * factorial(n)

def verify_small():
    """Verify against brute force and known values."""
    print("=== Verification against brute force ===")
    for n in range(1, 7):
        w_formula = solve(n)
        w_brute = brute_force(n) % M
        status = "OK" if w_formula == w_brute else "FAIL"
        print(f"  n={n}: formula={w_formula}, brute={w_brute}  [{status}]")

    print("\n=== Verification against problem examples ===")
    expected = {1: 1, 2: 6, 3: 168, 5: 19764360, 100: 384777056}
    for n, exp in expected.items():
        w = solve(n)
        status = "OK" if w == exp % M else "FAIL"
        print(f"  W({n}) mod M = {w}  (expected {exp % M})  [{status}]")

def make_visualization():
    """Generate a visualization and save to visualization.png."""

Project Euler Problem 409: Nim Extreme

Problem Statement

Project Euler Problem 409: Nim Extreme

Solution

Answer

Mathematical Derivation

Complexity

Verification

Correctness

Answer

Code