Problem 344: Silver Dollar Game

Project Euler Problem 344 Solution

EulerSolve provides an optimized solution for Project Euler Problem 344, Silver Dollar Game, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary There are $m=c+1$ coins on a strip of $n$ squares, and exactly one of them is the silver dollar. A legal move is either a regular left move of one coin or the special move of pocketing the leftmost coin. The first player wins iff they eventually pocket the silver dollar. So we must count all configurations from which the first player can force that outcome. If we first choose the occupied squares and then choose which of the $m$ coins is silver, the total number of raw configurations is $$T=m\binom{n}{m}.$$ Mathematical Approach Gap Model Write the coin positions as $$0\le x_0 \lt x_1 \lt \cdots \lt x_{m-1}\le n-1.$$ Instead of working with positions directly, the code switches to gaps: $$g_0=x_0,\qquad g_i=x_i-x_{i-1}-1\ (1\le i\le m-1),\qquad g_m=n-1-x_{m-1}.$$ These are the empty squares before the first coin, between consecutive coins, and after the last coin. Every legal placement corresponds to exactly one nonnegative gap vector $(g_0,\dots,g_m)$, and the total number of empty squares is $$g_0+g_1+\cdots+g_m=n-m=:N.$$ This is why the implementation begins with the parameter reduction $$m=c+1,\qquad N=n-m.$$ Ordinary Silver Dollar Positions as Nim Heaps The crucial combinatorial fact is that the ordinary silver dollar game can be described by alternating gaps ....

Detailed mathematical approach

Problem Summary

A strip is divided into $n$ squares in a row. On it lie $c$ worthless coins and one coin of value, the silver dollar, at most one per square, so $m=c+1$ coins in all. Two players move alternately. A regular move slides one coin one or more squares to the left, without leaving the strip and without landing on or jumping over another coin. Instead, a player may make the special move of pocketing the leftmost coin, and a player with no regular move left is forced to pocket it. Whoever pockets the silver dollar wins. A configuration is a placement of the $m$ coins together with the choice of which one is the dollar; it is winning when the player to move can force a win, and $W(n,c)$ counts those. We are given $W(10,2)=324$ and $W(100,10)=1514704946113500$, and asked for $W(10^6,100)$ modulo the semiprime ^[1]

$$M=1000036000099=1000003\cdot1000033.$$

The game is a variant of one invented by N. G. de Bruijn. ^[2]

The whole solution in one sentence: the game is Nim in disguise, its heaps being the gaps inside a pairing of the coins, so a configuration is losing exactly when the dollar is not leftmost and those gaps, after one small correction, have bitwise exclusive-or zero; we count such gap vectors one binary column at a time and subtract them from all configurations.

Mathematical Approach

1. Coins, gaps, and how many configurations there are

Number the squares $0,1,\dots,n-1$ from the left and let the coins stand on

$$0\le x_0\lt x_1\lt\cdots\lt x_{m-1}\le n-1.$$

Positions are awkward, since a move changes one of them; the gaps are better. Put

$$g_0=x_0,\qquad g_i=x_i-x_{i-1}-1,\qquad g_m=n-1-x_{m-1},$$

for $1\le i\le m-1$: the empty squares left of the first coin, between consecutive coins, and right of the last. Gaps and placements determine each other, so placements are exactly the vectors $(g_0,\dots,g_m)$ of $m+1$ nonnegative integers with

$$g_0+g_1+\cdots+g_m=n-m=:N.$$

By the stars-and-bars count there are $\binom{N+m}{m}=\binom{n}{m}$ of them, ^[3] one for each choice of occupied squares. Each placement carries $m$ choices of dollar, so there are $T=m\binom nm$ configurations. We count the losing ones, $L$, and report $W=T-L$, since it is those that a single equation describes.

2. What a move does to the gaps

Sliding coin $i$ left by $d\ge1$ squares is legal exactly when $d\le g_i$, and it replaces $g_i$ by $g_i-d$ and $g_{i+1}$ by $g_{i+1}+d$: a slide carries $d$ empty squares into the gap on the right. Pocketing the leftmost coin leaves the others where they are, so the new first gap is $x_1=g_0+g_1+1$ and the later gaps keep their values with indices shifted down by one. The total grows by one, as it must.

3. Two coins, worked out by hand

Take $n=5$ and $c=1$, so $m=2$ and $T=2\binom52=20$. If the dollar is the left coin, the player to move pockets it at once, so those $10$ configurations are winning. Let the dollar be the right coin and $g_1$ the gap between the two coins. If $g_1\ge1$, slide the dollar left until it touches the other coin, handing the opponent a position with $g_1=0$.

And $g_1=0$ is lost for the player to move. Pocketing the left coin makes the dollar leftmost and the opponent takes it. The only other move slides the left coin left by some $d\le g_0$, and we answer by sliding the dollar left by the same $d$, restoring $g_1=0$ with $g_0$ smaller by $d$. As $g_0$ strictly decreases, the coins eventually sit on squares $0$ and $1$, no slide is left, the opponent must pocket the left coin, and we pocket the dollar.

So exactly the four configurations with touching coins and the dollar on the right are losing, $W(5,1)=20-4=16$, and for general $n$ the same argument gives $L=n-1$ and $W(n,1)=(n-1)^2$. Two features survive: one gap decides everything, and widening it is worthless, since the partner coin undoes that.

4. Pairing the coins: the heaps

Pair consecutive coins from the right: the last with the one before it, then the next two, and so on. For even $m$ every coin gets a partner; for odd $m$ the first coin is left over and the left wall partners it, as if an immovable coin stood on square $-1$. Call the gap inside a pair a heap and the other gaps free. The heaps are then

$$g_1,g_3,\dots,g_{m-1}$$

for even $m$ and

$$g_0,g_2,\dots,g_{m-1}$$

for odd $m$: $m/2$ heaps in the first case, $(m+1)/2$ in the second.

If the dollar is leftmost the player to move pockets it and wins, so let its rank from the left be $j\in\{1,\dots,m-1\}$. The nim-sum $a\oplus b$ is the bitwise exclusive-or: an output bit is $1$ exactly when an odd number of inputs carry a $1$ there. So $2\oplus3=1$, binary $10$ against $11$ giving $01$, and equal numbers cancel. ^[4] Let the heap sum $\Phi$ be the nim-sum of the heaps, with one correction for odd $m$:

$$\Phi=g_1\oplus g_3\oplus\cdots\oplus g_{m-1},$$

$$\Phi=(g_0+\varepsilon)\oplus g_2\oplus\cdots\oplus g_{m-1},$$

again for even and for odd $m$, where $\varepsilon=0$ if $j=1$ and $\varepsilon=1$ if $j\ge2$.

Theorem. A configuration with the dollar not leftmost is losing for the player to move precisely when $\Phi=0$.

For $m=2$ this is Section 3: the only heap is $g_1$. The correction says that with an odd number of coins the wall sits on square $-1$ when the dollar is the second coin and on square $-2$ when it lies further right. Section 6 shows why.

5. Every slide changes exactly one heap

The heap indices run up to $m-1$ in steps of two, so they all have one parity and the free indices the other, and two consecutive gaps are never both heaps. A slide on coin $i$ touches $g_i$ and $g_{i+1}$, and exactly one is a heap: if $i$ has the heap parity, $g_i$ is a heap and shrinks; otherwise $i\ne m-1$, since $m-1$ has the heap parity, so $i+1\le m-1$ is a heap index and $g_{i+1}$ is a heap and grows.

So a slide changes exactly one heap, to a different value, hence it changes $\Phi$. Shrinking is the useful direction: for any integer $0\le v\lt g_i$ the slide with $d=g_i-v$ resets the heap $g_i$ to $v$. The one exception is the first heap of an odd number of coins when $\varepsilon=1$: its value is $g_0+1$, and sliding the first coin brings it down only to $1$, never to $0$.

6. What pocketing does to the heaps

Let $j\ge1$, so the leftmost coin is not the dollar. Pocketing it changes the number of coins, hence the pairing and the heaps. From the new gaps of Section 2:

Either parity, $j=1$: pocketing makes the dollar leftmost and the opponent takes it at once, so this move simply loses, whatever the heaps say. For odd $m$ that is why the first heap then keeps its plain size $g_0$.
Odd $m$, $j\ge2$: the new count $m-1$ is even and the new heaps are $g_2,g_4,\dots,g_{m-1}$, the old ones with the first deleted. So pocketing is exactly the move "empty the first heap", the one move Section 5 denies to a slide. That is why the first heap counts as $g_0+1$: it holds the $g_0$ empty squares plus one token that only pocketing can take, and pocketing takes them all, making it an ordinary nim heap of size $g_0+1$. The new nim-sum is $\Phi\oplus(g_0+1)$, which for $\Phi=0$ is $g_0+1\ge1$.
Even $m$, $j\ge2$: the new count is odd, the new first heap is $g_0+g_1+1+\varepsilon'$ with $\varepsilon'=0$ for $j=2$ and $\varepsilon'=1$ for $j\ge3$, and the others are $g_3,g_5,\dots,g_{m-1}$. Only the first heap changes, growing from $g_1$ to at least $g_1+1$: if $\Phi=0$ the new nim-sum is $(g_0+g_1+1+\varepsilon')\oplus g_1\ne0$, since $g_0+g_1+1+\varepsilon'\gt g_1$.

7. Proof of the theorem

Pocketing reduces the number of coins, and a slide keeps it while reducing $x_0+\cdots+x_{m-1}$, so we may induct on (number of coins, sum of positions) ordered lexicographically. A leftmost dollar is pocketed at once for the win, so let $j\ge1$.

If $\Phi=0$, every move hands the opponent a winning position. A slide changes exactly one heap, so it makes the heap sum nonzero. Pocketing loses at once when $j=1$, and for $j\ge2$, of either parity, Section 6 computes the new heap sum and finds it nonzero. By induction the position is losing.

If $\Phi\ne0$, some move reaches $\Phi=0$. Let $b$ be the highest bit set in $\Phi$; some heap $h$ has bit $b$ set, and then $h\oplus\Phi\lt h$, the standard nim-sum argument. ^[4] If that heap is a plain gap $g_i$, slide coin $i$ left by $g_i-(g_i\oplus\Phi)\ge1$ squares: the other heaps are untouched and the new heap sum is $0$. If instead it is the first heap of an odd number of coins with $\varepsilon=1$, put $h=g_0+1$ and $t=h\oplus\Phi\lt h$. For $t\ge1$ the slide of the first coin by $h-t$ is legal, since $1\le h-t\le g_0$, and it sets that heap to $t$. Only $t=0$ is beyond a slide, and it occurs exactly when $\Phi=g_0+1$, that is when $g_2\oplus\cdots\oplus g_{m-1}=0$; then pocketing deletes the first heap, leaving heap sum $0$ with $j-1\ge1$. Either way the opponent inherits a losing position, so the position is winning.

An example with four coins: $n=9$, coins on $1,2,5,7$, the dollar third. The gaps $(1,0,2,1,1)$ give heaps $g_1=0$, $g_3=1$ and $\Phi=1$, so the mover wins by sliding the last coin to square $6$, leaving $(1,0,2,0,2)$.

8. From the theorem to three counting problems

The criterion uses $j$ only through $\varepsilon$, so losing configurations can be counted rank by rank. Write $Z_{k,r}(S)$ for the number of ways to write $S$ as an ordered sum of $k+r$ nonnegative integers, $k$ of them marked and $r$ free, with the marked ones of nim-sum $0$. For even $m$ the criterion is the same for all $m-1$ admissible ranks, with $k=m/2$ marked gaps and $r=m/2+1$ free ones, so $L=(m-1)Z_{k,r}(N)$.

For odd $m$, where $k=r=(m+1)/2$, there are two families. Rank $j=1$ contributes $Z_{k,r}(N)$. Each of the $m-2$ ranks $j\ge2$ contributes the gap vectors with $(g_0+1)\oplus g_2\oplus\cdots\oplus g_{m-1}=0$; substituting $g_0'=g_0+1\ge1$ makes that a nim-sum condition on an ordinary vector of total $N+1$ whose first marked entry is positive. Vectors of total $N+1$ with first marked entry $0$ number $Z_{k-1,r}(N+1)$, since a marked $0$ adds nothing to the nim-sum and may be dropped. Hence

$$L=Z_{k,r}(N)+(m-2)\bigl(Z_{k,r}(N+1)-Z_{k-1,r}(N+1)\bigr).$$

For $m=1$ the single coin is the dollar, so $L=0$.

9. Counting nim-sum-zero vectors column by column

As the exclusive-or acts on each bit separately, a nim-sum is zero if and only if every binary column of the marked entries holds an even number of ones. Fix a column: if $x$ of the $k$ marked and $y$ of the $r$ free entries have a one in it, the column is admissible exactly when $x$ is even, and $\binom{k}{x}\binom{r}{y}$ choices of entries give that. Grouping by the total number of ones $s=x+y$,

$$w(s)=\sum_{2\mid x}\binom{k}{x}\binom{r}{s-x}$$

for $0\le s\le k+r$, the sum running over the even $x$ from $0$ to $k$ and $\binom{r}{s-x}$ read as $0$ when $s-x$ falls outside $0,\dots,r$. For $k=r=2$ this is $w=(1,2,2,2,1)$ for $s=0,\dots,4$; for instance $w(2)=\binom20\binom22+\binom22\binom20=2$.

10. Fixing the total: dynamic programming on carries

The columns are independent except for one global condition: the entries must add up to $S$. If column $b$ holds $s_b$ ones the total is $\sum_b s_b2^b$, a sum of digits that may exceed $1$, so it must be added with carries. Write $\beta_b$ for bit $b$ of $S$ and $c_b$ for the carry entering column $b$, with $c_0=0$. Working upwards, column $b$ is consistent with $S$ when

$$\begin{aligned}s_b+c_b&\equiv\beta_b\pmod 2,\\ c_{b+1}&=\frac{s_b+c_b-\beta_b}{2}.\end{aligned}$$

Summing $s_b=\beta_b-c_b+2c_{b+1}$ against $2^b$ and telescoping gives, for every $b$,

$$\sum_{t\lt b}s_t2^t=(S\bmod 2^b)+c_b2^b.$$

Let $B$ be the number of binary digits of $S$. Every entry is at most $S$, so all ones sit below column $B$, and the identity at $b=B$ says the entries total $S$ exactly when $c_B=0$. Further columns are harmless: a one in column $b\ge B$ would make some entry at least $2^B\gt S$, so no such choice ends with carry $0$. Now let $\mathrm{dp}_b(c)$ count the assignments to columns $0,\dots,b-1$ that match the lowest $b$ bits of $S$ and leave carry $c$, so $\mathrm{dp}_0(0)=1$ and $\mathrm{dp}_0(c)=0$ for $c\gt0$. Assignments sharing a carry may be merged, since later columns depend on nothing else. The recursion

$$\mathrm{dp}_{b+1}(c_{b+1})\mathrel{+}=\mathrm{dp}_b(c_b)\,w(s_b)$$

over all admissible $s_b$ therefore gives $Z_{k,r}(S)=\mathrm{dp}_B(0)$. The carry stays bounded: if $c_b\le k+r$ then $c_{b+1}\le(s_b+c_b)/2\le k+r$, so $k+r+1$ states suffice.

A small run in full, $S=2$ in binary $10$, with $k=r=2$. Column $0$ has $\beta_0=0$ and carry $0$, so $s_0\in\{0,2,4\}$ with weights $1,2,1$ and carries $0,1,2$: the state vector becomes $(1,2,1,0,0)$. Column $1$ has $\beta_1=1$, and only $s_1=1$ from carry $0$ and $s_1=0$ from carry $1$ end at carry $0$, giving $1\cdot2+2\cdot1=4$. Indeed the vectors are those with $g_0=g_2=a$ and $g_1+g_3=2-2a$: $3+1=4$ of them.

11. The checkpoint $W(10,2)=324$ in full

Here $m=3$, $N=7$, $k=r=2$ and $T=3\binom{10}{3}=360$. As $m$ is odd the heaps are $g_0$ (or $g_0+1$) and $g_2$, and a nim-sum of two numbers vanishes exactly when they are equal, so both counts can be worked out by hand.

Rank $j=1$: $g_0=g_2=a$ and $g_1+g_3=7-2a$, with $8-2a$ solutions. Over $a=0,1,2,3$ that is $8+6+4+2=20=Z_{2,2}(7)$.
Rank $j=2$: $g_2=g_0+1$, so with $g_0=a$ we need $g_1+g_3=6-2a$, with $7-2a$ solutions, and $7+5+3+1=16$.

The second count is the difference of Section 8: $Z_{2,2}(8)=9+7+5+3+1=25$ and $Z_{1,2}(8)=9$, a single marked entry having nim-sum zero only when it is zero, leaving $g_1+g_3=8$. So $L=20+(3-2)\cdot16=36$ and $W=360-36=324$, the published value. For general $n\ge3$, rank $j=1$ contributes the positive terms of $N+1,N-1,\dots$ and rank $j=2$ those of $N,N-2,\dots$; together these are the integers $1$ to $N+1$, each once, so $L$ is an arithmetic series:

$$L=\frac{(N+1)(N+2)}{2},\qquad W(n,2)=\frac{(n-1)^2(n-2)}{2},$$

a sequence recorded as $t(t+1)^2/2$ with $t=n-2$. ^[5] At $n=10$ it reads $9^2\cdot8/2=324$.

12. Working modulo the semiprime

Only $W$ modulo $M=1000003\cdot1000033$ is wanted, and both factors are prime and exceed $10^6$, so no factorial $i!$ with $i\le10^6$ is divisible by either. Hence $\binom{n}{m}$ comes from tables of factorials and inverse factorials modulo each prime, the inverse of $a$ modulo a prime $p$ being $a^{p-2}$ by Fermat's little theorem. ^[6] As the moduli are coprime, the two least residues $a_1,a_2$ determine one residue modulo $M$: take $k$ as the least residue of $(a_2-a_1)\cdot1000003^{-1}$ modulo $1000033$, and then $a_1+1000003k$ lies below $M$ and leaves $a_1$ and $a_2$ on division by the two primes. ^[7] Everything else is addition and multiplication, which commute with reduction, with one caution: two residues below $M$ multiply to roughly $10^{24}$, so products need a type wider than 64 bits.

How the Code Works

Each program builds the Pascal triangle, then the two factorial tables, then $T$, and only then counts the losing configurations; the C++ runs its exact checkpoints in between. ^[8] The names below are the C++ and Python ones: Java builds SmallComb, FactTable and WaysEvenMod as constructors and uses binomModPrime, crt, countWithWaysMod, computeWMod. build_small_comb(52) serves Section 9, which at $c=100$ needs $k=r=51$; that bound is tied to this input. build_fact_table(n, p) fills fact upwards, inverts fact[n] by one modular exponentiation, the Fermat inverse of Section 12, then fills invfact downwards; binom_mod_prime multiplies three of their entries, and binom_mod_composite turns the two residues into one modulo $M$ through crt.

build_ways_even_mod(k, r, comb) is the array $w$ of Section 9: it loops over even x and all y, adding $\binom kx\binom ry$ to ways[x+y]; the factors reach $\binom{51}{25}=247959266474052$, so their product needs $96$ bits. count_with_ways_mod(S, ways) is the carry recursion of Section 10, with dp indexed by the carry up to max_s${}=k+r$: per column it adds dp[carry] * ways[s] into next[(s+carry)/2] for every admissible pair. It runs bit_length(S) + 7 columns and returns dp[0], the extra columns being harmless by Section 10.

compute_W_mod then assembles the answer: $m=c+1$, $N=n-m$, total${}=m\binom nm$, an early return for $m=1$, and either the even branch with one call and weight $m-1$ or the odd branch with the three calls of Section 8. The subtraction $T-L$ runs as (total + MOD - losing) % MOD, so nothing turns negative.

All three run the same recurrence and differ only in representation. C++ keeps residues in uint64_t and casts every modular product to unsigned __int128. Java keeps long[] arrays and calls BigInteger inside mulMod, since a product of two residues near $M$ overflows long; its factorial loops, binomModPrime and crt multiply directly in long, with operands below $10^{12}$. Python holds arbitrary-precision integers in lists. The C++ also carries an exact copy of the pipeline, compute_W_exact and its helpers, printed by to_string_u128; its binom_exact alternates multiplication and division, res = res * (n-k+i) / i, exact because the partial product is always $\binom{n-k+i}{i}$. main runs it on the two checkpoints first; the Python and Java ports keep only the modular path.

Complexity and Verification

One call to count_with_ways_mod runs through the binary digits of $S$ plus seven further columns, each trying every pair (carry, column total), so it costs $O(m^2\log n)$ time in $O(m)$ memory; with the factorial preprocessing the whole modular computation is $O(n+m^2\log n)$ operations in $O(n+m)$ storage. For $n=10^6$ and $c=100$ we have $m=101$, $k=r=51$, max_s${}=102$ and $N=999899$ with $20$ binary digits, hence $27$ columns: each of the three calls examines at most $27\cdot103^2=286443$ pairs. Only $w(102)$ vanishes, since $s=102$ would force $x=y=51$ and $51$ is odd.

The dominant cost is elsewhere: the four factorial arrays hold $4(n+1)$ entries, about $32$ MB in C++ and Java, and filling them takes $4n$ modular multiplications and two exponentiations, all for one binomial coefficient that $101$ factors would also give. Python's lists hold references to separate integer objects, and the same four tables measured about five times as much. Each program still finishes in well under a second, and all three return the same residue for $W(10^6,100)$. The C++ also prints the two statement checkpoints from the exact path, $W(10,2)=324$ and $W(100,10)=1514704946113500$, both as published.

The method was also checked independently. A brute force plays the real game by backward induction over positions, trying every slide and the pocket move, with no reference to gaps or nim-sums. Over all $69$ pairs $(n,c)$ with $n\le14$ and $c\le5$, covering $80757$ configurations, every value agrees with the criterion of Section 4 and every total with the formula of Section 8; for $W(5,1)=16$ and $W(5,2)=24$ the losing configurations were listed and compared. The closed forms $W(n,1)=(n-1)^2$ and $W(n,2)=(n-1)^2(n-2)/2$ match the dynamic programming for all $n\le200$, the second agreeing with the recorded sequence, ^[5] and the carry recursion reproduced the hand counts $Z_{2,2}(7)=20$, $Z_{2,2}(8)=25$, $Z_{1,2}(8)=9$ and $Z_{2,2}(2)=4$. Beyond the two checkpoints the solution files package no test.

Footnotes and References

The references below supply the background facts used above; the heap pairing, the theorem and its proof, the reduction to $Z_{k,r}$ and the carry recursion are derived here.

Project Euler 344 — Silver Dollar Game. The official statement describes one variant of N. G. de Bruijn's silver dollar game, defines the regular and special moves, names the winner as whoever pockets the silver dollar, defines a winning configuration and W(n, c), gives W(10, 2) = 324 and W(100, 10) = 1514704946113500, and asks for W(1000000, 100) modulo the semiprime 1000036000099.
Cut the Knot — the Silver Dollar game. The page attributes the Silver Dollar game to N. G. de Bruijn and describes variants with a terminal square where counters may accumulate. It is cited only for that attribution, reference [1] fixing the exact rules, and gives no winning strategy.
Wikipedia — Stars and bars. Its second theorem gives C(n + k − 1, k − 1) for the number of k-tuples of non-negative integers with sum n, used in Section 1 with k = m + 1 and n = N.
Wikipedia — Nim. Defines the nim-sum as the bitwise exclusive-or of the heap sizes and states Bouton's theorem, that under normal play the first player wins exactly when the nim-sum is not zero; its second lemma, that a nonzero nim-sum can be moved to zero by reducing a heap carrying its highest bit, is reused in Section 7.
OEIS A006002. The sequence a(t) = t(t + 1)²/2, with terms 0, 2, 9, 24, 50, 90, 147, 224, 324, 450, is the closed form of Section 11: W(n, 2) = a(n − 2), which at n = 10 gives the checkpoint 324.
Wikipedia — Fermat's little theorem. For a prime p and an integer a coprime to p it states a^(p−1) ≡ 1 (mod p), so a^(p−2) inverts a modulo p, which seeds the inverse factorials of Section 12.
Wikipedia — Chinese remainder theorem. For pairwise coprime moduli with product N and prescribed residues it states that exactly one integer x satisfies 0 ≤ x < N, and gives the two-modulus construction from Bézout's identity.
C++ — Euler344.cpp. The linked immutable revision holds build_fact_table, crt, build_ways_even_mod, count_with_ways_mod and compute_W_mod, their exact counterparts, and the block in main that checks W(10, 2) and W(100, 10). It is the implementation record for this article.

Problem 344 source code

C++

#include <algorithm>
#include <cstdint>
#include <future>
#include <iostream>
#include <string>
#include <vector>

using u128 = unsigned __int128;

namespace {

constexpr uint64_t kMod = 1000036000099ULL;
constexpr uint64_t kP1 = 1000003ULL;
constexpr uint64_t kP2 = 1000033ULL;

struct SmallComb {
    std::vector<std::vector<unsigned long long>> c;
};

SmallComb build_small_comb(int max_n) {
    SmallComb comb;
    comb.c.assign(max_n + 1, std::vector<unsigned long long>(max_n + 1, 0));
    for (int n = 0; n <= max_n; ++n) {
        comb.c[n][0] = comb.c[n][n] = 1;
        for (int k = 1; k < n; ++k) {
            comb.c[n][k] = comb.c[n - 1][k - 1] + comb.c[n - 1][k];
        }
    }
    return comb;
}

uint64_t mod_pow(uint64_t base, uint64_t exp, uint64_t mod) {
    uint64_t result = 1;
    while (exp > 0) {
        if (exp & 1) result = static_cast<uint64_t>((u128)result * base % mod);
        base = static_cast<uint64_t>((u128)base * base % mod);
        exp >>= 1;
    }
    return result;
}

uint64_t mod_inv_prime(uint64_t a, uint64_t p) {
    return mod_pow(a, p - 2, p);
}

struct FactTable {
    std::vector<uint64_t> fact;
    std::vector<uint64_t> invfact;
};

FactTable build_fact_table(int n, uint64_t p) {
    FactTable ft;
    ft.fact.resize(n + 1);
    ft.invfact.resize(n + 1);
    ft.fact[0] = 1;
    for (int i = 1; i <= n; ++i) {
        ft.fact[i] = static_cast<uint64_t>((u128)ft.fact[i - 1] * i % p);
    }
    ft.invfact[n] = mod_pow(ft.fact[n], p - 2, p);
    for (int i = n; i >= 1; --i) {
        ft.invfact[i - 1] = static_cast<uint64_t>((u128)ft.invfact[i] * i % p);
    }
    return ft;
}

uint64_t binom_mod_prime(int n, int k, const FactTable& ft, uint64_t p) {
    if (k < 0 || k > n) return 0;
    uint64_t res = static_cast<uint64_t>((u128)ft.fact[n] * ft.invfact[k] % p);
    res = static_cast<uint64_t>((u128)res * ft.invfact[n - k] % p);
    return res;
}

uint64_t crt(uint64_t a, uint64_t b) {
    uint64_t t = (b + kP2 - (a % kP2)) % kP2;
    uint64_t inv = mod_inv_prime(kP1 % kP2, kP2);
    uint64_t k = static_cast<uint64_t>((u128)t * inv % kP2);
    return static_cast<uint64_t>((u128)a + (u128)kP1 * k);
}

uint64_t binom_mod_composite(int n, int k, const FactTable& f1, const FactTable& f2) {
    uint64_t a = binom_mod_prime(n, k, f1, kP1);
    uint64_t b = binom_mod_prime(n, k, f2, kP2);
    return crt(a, b);
}

struct WaysEvenMod {
    int max_s = 0;
    std::vector<uint64_t> ways;
};

WaysEvenMod build_ways_even_mod(int k, int r, const SmallComb& comb) {
    WaysEvenMod out;
    out.max_s = k + r;
    out.ways.assign(out.max_s + 1, 0);
    for (int x = 0; x <= k; x += 2) {
        unsigned long long cx = comb.c[k][x];
        for (int y = 0; y <= r; ++y) {
            int s = x + y;
            u128 add = static_cast<u128>(cx) * comb.c[r][y];
            uint64_t val = static_cast<uint64_t>(add % kMod);
            uint64_t sum = out.ways[s] + val;
            if (sum >= kMod) sum -= kMod;
            out.ways[s] = sum;
        }
    }
    return out;
}

int bit_length(long long v) {
    int bits = 0;
    while (v > 0) {
        ++bits;
        v >>= 1;
    }
    return bits == 0 ? 1 : bits;
}

uint64_t count_with_ways_mod(long long S, const WaysEvenMod& ways) {
    int max_s = ways.max_s;
    std::vector<uint64_t> dp(max_s + 1, 0);
    std::vector<uint64_t> next(max_s + 1, 0);
    dp[0] = 1;

    int bits = bit_length(S) + 7; // flush remaining carry.
    for (int b = 0; b < bits; ++b) {
        int bit = static_cast<int>((S >> b) & 1LL);
        std::fill(next.begin(), next.end(), 0);
        for (int carry = 0; carry <= max_s; ++carry) {
            uint64_t base = dp[carry];
            if (base == 0) continue;
            for (int s = 0; s <= max_s; ++s) {
                uint64_t ways_s = ways.ways[s];
                if (ways_s == 0) continue;
                int total = s + carry;
                if ((total & 1) != bit) continue;
                int carry_out = total >> 1;
                uint64_t add = static_cast<uint64_t>((u128)base * ways_s % kMod);
                uint64_t sum = next[carry_out] + add;
                if (sum >= kMod) sum -= kMod;
                next[carry_out] = sum;
            }
        }
        dp.swap(next);
    }
    return dp[0];
}

struct WaysEvenExact {
    int max_s = 0;
    std::vector<u128> ways;
};

WaysEvenExact build_ways_even_exact(int k, int r, const SmallComb& comb) {
    WaysEvenExact out;
    out.max_s = k + r;
    out.ways.assign(out.max_s + 1, 0);
    for (int x = 0; x <= k; x += 2) {
        u128 cx = comb.c[k][x];
        for (int y = 0; y <= r; ++y) {
            int s = x + y;
            out.ways[s] += cx * comb.c[r][y];
        }
    }
    return out;
}

u128 count_with_ways_exact(long long S, const WaysEvenExact& ways) {
    int max_s = ways.max_s;
    std::vector<u128> dp(max_s + 1, 0);
    std::vector<u128> next(max_s + 1, 0);
    dp[0] = 1;

    int bits = bit_length(S) + 7;
    for (int b = 0; b < bits; ++b) {
        int bit = static_cast<int>((S >> b) & 1LL);
        std::fill(next.begin(), next.end(), 0);
        for (int carry = 0; carry <= max_s; ++carry) {
            u128 base = dp[carry];
            if (base == 0) continue;
            for (int s = 0; s <= max_s; ++s) {
                u128 ways_s = ways.ways[s];
                if (ways_s == 0) continue;
                int total = s + carry;
                if ((total & 1) != bit) continue;
                int carry_out = total >> 1;
                next[carry_out] += base * ways_s;
            }
        }
        dp.swap(next);
    }
    return dp[0];
}

u128 binom_exact(int n, int k) {
    if (k < 0 || k > n) return 0;
    if (k > n - k) k = n - k;
    u128 res = 1;
    for (int i = 1; i <= k; ++i) {
        res = res * static_cast<u128>(n - k + i) / static_cast<u128>(i);
    }
    return res;
}

u128 compute_W_exact(int n, int c, const SmallComb& comb) {
    int m = c + 1;
    long long N = n - m;
    u128 total = static_cast<u128>(m) * binom_exact(n, m);
    if (m == 1) return total;

    int k = (m + 1) / 2;
    int r = m - k + 1;
    if (m % 2 == 0) {
        WaysEvenExact ways = build_ways_even_exact(k, r, comb);
        u128 losing = static_cast<u128>(m - 1) * count_with_ways_exact(N, ways);
        return total - losing;
    }

    WaysEvenExact ways_k = build_ways_even_exact(k, r, comb);
    WaysEvenExact ways_km1 = build_ways_even_exact(k - 1, r, comb);
    u128 L0 = count_with_ways_exact(N, ways_k);
    u128 L1 = count_with_ways_exact(N + 1, ways_k) - count_with_ways_exact(N + 1, ways_km1);
    u128 losing = L0 + static_cast<u128>(m - 2) * L1;
    return total - losing;
}

uint64_t compute_W_mod(int n, int c, const SmallComb& comb,
                       const FactTable& f1, const FactTable& f2) {
    int m = c + 1;
    long long N = n - m;
    uint64_t total = static_cast<uint64_t>((u128)m * binom_mod_composite(n, m, f1, f2) % kMod);
    if (m == 1) return total;

    int k = (m + 1) / 2;
    int r = m - k + 1;
    if (m % 2 == 0) {
        WaysEvenMod ways = build_ways_even_mod(k, r, comb);
        uint64_t losing = static_cast<uint64_t>((u128)(m - 1) * count_with_ways_mod(N, ways) % kMod);
        return (total + kMod - losing) % kMod;
    }

    WaysEvenMod ways_k = build_ways_even_mod(k, r, comb);
    WaysEvenMod ways_km1 = build_ways_even_mod(k - 1, r, comb);

    auto fut0 = std::async(std::launch::async, [&]() { return count_with_ways_mod(N, ways_k); });
    auto fut1 = std::async(std::launch::async, [&]() { return count_with_ways_mod(N + 1, ways_k); });
    auto fut2 = std::async(std::launch::async, [&]() { return count_with_ways_mod(N + 1, ways_km1); });

    uint64_t L0 = fut0.get();
    uint64_t L1 = fut1.get();
    uint64_t L1_sub = fut2.get();
    L1 = (L1 + kMod - L1_sub) % kMod;

    uint64_t losing = (L0 + static_cast<uint64_t>((u128)(m - 2) * L1 % kMod)) % kMod;
    return (total + kMod - losing) % kMod;
}

std::string to_string_u128(u128 v) {
    if (v == 0) return "0";
    std::string s;
    while (v > 0) {
        int digit = static_cast<int>(v % 10);
        s.push_back(static_cast<char>('0' + digit));
        v /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

} // namespace

int main() {
    const int max_small = 52;
    SmallComb comb = build_small_comb(max_small);

    std::cout << "--- Validation Checkpoints ---\n";
    u128 v1 = compute_W_exact(10, 2, comb);
    u128 v2 = compute_W_exact(100, 10, comb);
    std::cout << "W(10, 2) = " << to_string_u128(v1)
              << (v1 == 324 ? " [PASS]" : " [FAIL]") << "\n";
    std::cout << "W(100, 10) = " << to_string_u128(v2)
              << (v2 == 1514704946113500ULL ? " [PASS]" : " [FAIL]") << "\n";

    std::cout << "\n--- Final Solution ---\n";
    const int n = 1000000;
    const int c = 100;
    FactTable f1 = build_fact_table(n, kP1);
    FactTable f2 = build_fact_table(n, kP2);

    uint64_t result = compute_W_mod(n, c, comb, f1, f2);
    std::cout << "W(" << n << ", " << c << ") mod " << kMod << " = " << result << "\n";
    return 0;
}

Python

def bit_length(v):
    bits = 0
    while v > 0:
        bits += 1
        v >>= 1
    return 1 if bits == 0 else bits

kMod = 1000036000099
kP1 = 1000003
kP2 = 1000033

class SmallComb:
    def __init__(self, max_n):
        self.c = [[0] * (max_n + 1) for _ in range(max_n + 1)]
        for n in range(max_n + 1):
            self.c[n][0] = self.c[n][n] = 1
            for k in range(1, n):
                self.c[n][k] = self.c[n - 1][k - 1] + self.c[n - 1][k]

def build_small_comb(max_n):
    return SmallComb(max_n)

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

def mod_inv_prime(a, p):
    return pow(a, p - 2, p)

class FactTable:
    def __init__(self, n, p):
        self.fact = [0] * (n + 1)
        self.invfact = [0] * (n + 1)
        self.fact[0] = 1
        for i in range(1, n + 1):
            self.fact[i] = (self.fact[i - 1] * i) % p
        self.invfact[n] = mod_inv_prime(self.fact[n], p)
        for i in range(n, 0, -1):
            self.invfact[i - 1] = (self.invfact[i] * i) % p

def build_fact_table(n, p):
    return FactTable(n, p)

def binom_mod_prime(n, k, ft, p):
    if k < 0 or k > n: return 0
    res = (ft.fact[n] * ft.invfact[k]) % p
    res = (res * ft.invfact[n - k]) % p
    return res

def crt(a, b):
    t = (b + kP2 - (a % kP2)) % kP2
    inv = mod_inv_prime(kP1 % kP2, kP2)
    k = (t * inv) % kP2
    return a + kP1 * k

def binom_mod_composite(n, k, f1, f2):
    a = binom_mod_prime(n, k, f1, kP1)
    b = binom_mod_prime(n, k, f2, kP2)
    return crt(a, b)

class WaysEvenMod:
    def __init__(self, k, r, comb):
        self.max_s = k + r
        self.ways = [0] * (self.max_s + 1)
        for x in range(0, k + 1, 2):
            cx = comb.c[k][x]
            for y in range(0, r + 1):
                s = x + y
                add = (cx * comb.c[r][y]) % kMod
                self.ways[s] = (self.ways[s] + add) % kMod

def build_ways_even_mod(k, r, comb):
    return WaysEvenMod(k, r, comb)

def count_with_ways_mod(S, ways):
    max_s = ways.max_s
    dp = [0] * (max_s + 1)
    dp[0] = 1
    
    bits = bit_length(S) + 7
    for b in range(bits):
        bit = (S >> b) & 1
        nxt = [0] * (max_s + 1)
        for carry in range(max_s + 1):
            base = dp[carry]
            if base == 0: continue
            for s in range(max_s + 1):
                ways_s = ways.ways[s]
                if ways_s == 0: continue
                total = s + carry
                if (total & 1) != bit: continue
                carry_out = total >> 1
                add = (base * ways_s) % kMod
                nxt[carry_out] = (nxt[carry_out] + add) % kMod
        dp = nxt
    return dp[0]

def compute_W_mod(n, c, comb, f1, f2):
    m = c + 1
    N = n - m
    total = (m * binom_mod_composite(n, m, f1, f2)) % kMod
    if m == 1:
        return total
        
    k = (m + 1) // 2
    r = m - k + 1
    if m % 2 == 0:
        ways = build_ways_even_mod(k, r, comb)
        losing = ((m - 1) * count_with_ways_mod(N, ways)) % kMod
        return (total + kMod - losing) % kMod
        
    ways_k = build_ways_even_mod(k, r, comb)
    ways_km1 = build_ways_even_mod(k - 1, r, comb)
    
    L0 = count_with_ways_mod(N, ways_k)
    L1 = count_with_ways_mod(N + 1, ways_k)
    L1_sub = count_with_ways_mod(N + 1, ways_km1)
    L1 = (L1 + kMod - L1_sub) % kMod
    
    losing = (L0 + ((m - 2) * L1) % kMod) % kMod
    return (total + kMod - losing) % kMod

def solve():
    max_small = 52
    comb = build_small_comb(max_small)
    n = 1000000
    c = 100
    f1 = build_fact_table(n, kP1)
    f2 = build_fact_table(n, kP2)
    ans = compute_W_mod(n, c, comb, f1, f2)
    return str(ans)

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

Java

import java.util.*;

public class Euler344 {

    static final long MOD = 1000036000099L;
    static final long P1 = 1000003L;
    static final long P2 = 1000033L;

    static class SmallComb {
        long[][] c;

        SmallComb(int maxN) {
            c = new long[maxN + 1][maxN + 1];
            for (int n = 0; n <= maxN; n++) {
                c[n][0] = c[n][n] = 1;
                for (int k = 1; k < n; k++) {
                    c[n][k] = c[n - 1][k - 1] + c[n - 1][k];
                }
            }
        }
    }

    static long modPow(long base, long exp, long mod) {
        long result = 1;
        while (exp > 0) {
            if ((exp & 1) == 1)
                result = mulMod(result, base, mod);
            base = mulMod(base, base, mod);
            exp >>= 1;
        }
        return result;
    }

    static long mulMod(long a, long b, long mod) {
        return (long) ((((java.math.BigInteger.valueOf(a)).multiply(java.math.BigInteger.valueOf(b)))
                .remainder(java.math.BigInteger.valueOf(mod))).longValue());
    }

    static long modInvPrime(long a, long p) {
        return modPow(a, p - 2, p);
    }

    static class FactTable {
        long[] fact;
        long[] invfact;

        FactTable(int n, long p) {
            fact = new long[n + 1];
            invfact = new long[n + 1];
            fact[0] = 1;
            for (int i = 1; i <= n; i++) {
                fact[i] = (fact[i - 1] * i) % p;
            }
            invfact[n] = modInvPrime(fact[n], p);
            for (int i = n; i >= 1; i--) {
                invfact[i - 1] = (invfact[i] * i) % p;
            }
        }
    }

    static long binomModPrime(int n, int k, FactTable ft, long p) {
        if (k < 0 || k > n)
            return 0;
        long res = (ft.fact[n] * ft.invfact[k]) % p;
        res = (res * ft.invfact[n - k]) % p;
        return res;
    }

    static long crt(long a, long b) {
        long t = (b + P2 - (a % P2)) % P2;
        long inv = modInvPrime(P1 % P2, P2);
        long k = (t * inv) % P2;
        return a + P1 * k;
    }

    static long binomModComposite(int n, int k, FactTable f1, FactTable f2) {
        long a = binomModPrime(n, k, f1, P1);
        long b = binomModPrime(n, k, f2, P2);
        return crt(a, b);
    }

    static class WaysEvenMod {
        int max_s;
        long[] ways;

        WaysEvenMod(int k, int r, SmallComb comb) {
            max_s = k + r;
            ways = new long[max_s + 1];
            for (int x = 0; x <= k; x += 2) {
                long cx = comb.c[k][x];
                for (int y = 0; y <= r; y++) {
                    int s = x + y;
                    long add = mulMod(cx, comb.c[r][y], MOD);
                    ways[s] = (ways[s] + add) % MOD;
                }
            }
        }
    }

    static int bitLength(long v) {
        int bits = 0;
        while (v > 0) {
            bits++;
            v >>= 1;
        }
        return bits == 0 ? 1 : bits;
    }

    static long countWithWaysMod(long S, WaysEvenMod ways) {
        int max_s = ways.max_s;
        long[] dp = new long[max_s + 1];
        long[] next = new long[max_s + 1];
        dp[0] = 1;

        int bits = bitLength(S) + 7;
        for (int b = 0; b < bits; b++) {
            int bit = (int) ((S >> b) & 1L);
            Arrays.fill(next, 0);
            for (int carry = 0; carry <= max_s; carry++) {
                long base = dp[carry];
                if (base == 0)
                    continue;
                for (int s = 0; s <= max_s; s++) {
                    long ways_s = ways.ways[s];
                    if (ways_s == 0)
                        continue;
                    int total = s + carry;
                    if ((total & 1) != bit)
                        continue;
                    int carryOut = total >> 1;
                    long add = mulMod(base, ways_s, MOD);
                    next[carryOut] = (next[carryOut] + add) % MOD;
                }
            }
            long[] temp = dp;
            dp = next;
            next = temp;
        }
        return dp[0];
    }

    static long computeWMod(int n, int c, SmallComb comb, FactTable f1, FactTable f2) {
        int m = c + 1;
        long N = n - m;
        long total = mulMod(m, binomModComposite(n, m, f1, f2), MOD);
        if (m == 1)
            return total;

        int k = (m + 1) / 2;
        int r = m - k + 1;
        if (m % 2 == 0) {
            WaysEvenMod ways = new WaysEvenMod(k, r, comb);
            long losing = mulMod(m - 1, countWithWaysMod(N, ways), MOD);
            return (total + MOD - losing) % MOD;
        }

        WaysEvenMod ways_k = new WaysEvenMod(k, r, comb);
        WaysEvenMod ways_km1 = new WaysEvenMod(k - 1, r, comb);

        long L0 = countWithWaysMod(N, ways_k);
        long L1 = countWithWaysMod(N + 1, ways_k);
        long L1_sub = countWithWaysMod(N + 1, ways_km1);
        L1 = (L1 + MOD - L1_sub) % MOD;

        long losing = (L0 + mulMod(m - 2, L1, MOD)) % MOD;
        return (total + MOD - losing) % MOD;
    }

    public static String solve() {
        int maxSmall = 52;
        SmallComb comb = new SmallComb(maxSmall);
        int n = 1000000;
        int c = 100;
        FactTable f1 = new FactTable(n, P1);
        FactTable f2 = new FactTable(n, P2);

        long ans = computeWMod(n, c, comb, f1, f2);
        return String.valueOf(ans);
    }

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