Problem 796: A Grand Shuffle

View on Project Euler

Project Euler Problem 796 Solution

EulerSolve provides an optimized solution for Project Euler Problem 796, A Grand Shuffle, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We study a shoe made from \(D\) decks. Each deck contains \(S \cdot R\) ordinary cards and \(J\) jokers. An ordinary card carries three labels: its deck, its suit, and its rank. A joker contributes only to its deck label, not to any suit or rank. After a uniform random shuffle, cards are revealed from the top until every required deck, every required suit, and every required rank has appeared at least once. For the Project Euler instance, $$D=10,\qquad S=4,\qquad R=13,\qquad J=2,$$ so the shoe contains $$M=D(SR+J)=10(4\cdot 13+2)=540$$ cards in total. The task is to compute the expected stopping time exactly enough to print eight decimal places. Mathematical Approach The implementations evaluate a closed inclusion-exclusion formula. The key is to count, for each family of still-missing categories, how many cards are allowed to appear before that family is broken. Step 1: Define the stopping time Let \(T\) be the draw index at which all required categories have been seen for the first time. Because no more than \(M\) cards can be drawn, the tail-sum identity gives $$\mathbb E[T]=\sum_{n=0}^{M-1}\Pr(T>n).$$ So we only need a formula for the probability that after \(n\) draws at least one required deck, suit, or rank is still missing. Step 2: Choose which categories are still missing Fix three sets of missing categories: \(s\) suits, \(r\) ranks, and \(d\) decks....

Detailed mathematical approach

Problem Summary

We study a shoe made from \(D\) decks. Each deck contains \(S \cdot R\) ordinary cards and \(J\) jokers. An ordinary card carries three labels: its deck, its suit, and its rank. A joker contributes only to its deck label, not to any suit or rank. After a uniform random shuffle, cards are revealed from the top until every required deck, every required suit, and every required rank has appeared at least once.

For the Project Euler instance,

$$D=10,\qquad S=4,\qquad R=13,\qquad J=2,$$

so the shoe contains

$$M=D(SR+J)=10(4\cdot 13+2)=540$$

cards in total. The task is to compute the expected stopping time exactly enough to print eight decimal places.

Mathematical Approach

The implementations evaluate a closed inclusion-exclusion formula. The key is to count, for each family of still-missing categories, how many cards are allowed to appear before that family is broken.

Step 1: Define the stopping time

Let \(T\) be the draw index at which all required categories have been seen for the first time. Because no more than \(M\) cards can be drawn, the tail-sum identity gives

$$\mathbb E[T]=\sum_{n=0}^{M-1}\Pr(T>n).$$

So we only need a formula for the probability that after \(n\) draws at least one required deck, suit, or rank is still missing.

Step 2: Choose which categories are still missing

Fix three sets of missing categories: \(s\) suits, \(r\) ranks, and \(d\) decks. There are

$$\binom{S}{s}\binom{R}{r}\binom{D}{d}$$

ways to choose them.

If those categories are all still unseen after \(n\) draws, then every revealed card must avoid the missing sets. Only the remaining \(D-d\) decks may appear. Inside such a deck, an ordinary card is allowed only if its suit lies among the \(S-s\) allowed suits and its rank lies among the \(R-r\) allowed ranks. Jokers in the surviving decks are always allowed, because they carry no suit or rank label. Therefore the number of allowed cards is

$$A(s,r,d)=(D-d)\bigl((S-s)(R-r)+J\bigr).$$

The complementary number of forbidden cards is

$$B(s,r,d)=M-A(s,r,d).$$

Step 3: Probability that the first \(n\) draws avoid those categories

After a random shuffle, the first \(n\) cards form a uniformly random \(n\)-subset of the whole shoe. Hence, for fixed missing sets, the probability that all first \(n\) cards lie among the \(A(s,r,d)\) allowed cards is

$$\frac{\binom{A(s,r,d)}{n}}{\binom{M}{n}},$$

with the understanding that the expression is \(0\) when \(n>A(s,r,d)\).

Now apply inclusion-exclusion over all nonempty choices of missing suits, ranks, and decks:

$$\Pr(T>n)=\sum_{\substack{0\le s\le S\\0\le r\le R\\0\le d\le D\\s+r+d>0}}(-1)^{s+r+d+1}\binom{S}{s}\binom{R}{r}\binom{D}{d}\frac{\binom{A(s,r,d)}{n}}{\binom{M}{n}}.$$

The sign is positive when an odd number of category families is declared missing and negative when it is even.

Step 4: Collapse the tail sum

Insert the previous formula into the tail-sum expression for \(\mathbb E[T]\) and exchange the order of summation:

$$\mathbb E[T]=\sum_{\substack{0\le s\le S\\0\le r\le R\\0\le d\le D\\s+r+d>0}}(-1)^{s+r+d+1}\binom{S}{s}\binom{R}{r}\binom{D}{d}\sum_{n=0}^{A(s,r,d)}\frac{\binom{A(s,r,d)}{n}}{\binom{M}{n}}.$$

The inner sum is classical. If \(B=M-A\) cards are marked as forbidden, then in a random permutation the expected position of the first marked card is

$$\frac{M+1}{B+1}.$$

The same quantity is equal to

$$\sum_{n=0}^{A}\frac{\binom{A}{n}}{\binom{M}{n}}=\frac{M+1}{M-A+1}=\frac{M+1}{B+1}.$$

So each inclusion-exclusion term collapses to a single rational factor.

Step 5: Final closed form

Substituting \(B(s,r,d)=M-(D-d)\bigl((S-s)(R-r)+J\bigr)\), we obtain

$$\boxed{\mathbb E[T]=\sum_{\substack{0\le s\le S\\0\le r\le R\\0\le d\le D\\s+r+d>0}}(-1)^{s+r+d+1}\binom{S}{s}\binom{R}{r}\binom{D}{d}\frac{M+1}{B(s,r,d)+1}.}$$

This is exactly the quantity evaluated by the implementations. The formula also handles reduced variants: if suits, ranks, or decks are not required, the corresponding summation range collapses to the single value \(0\).

Worked Example: only ranks are required in one deck

One implementation checks the simpler case \(D=1\), \(S=4\), \(R=13\), \(J=2\), where only the \(13\) ranks matter. Then

$$M=1(4\cdot 13+2)=54.$$

There is no summation over suits or decks, so only \(r\) remains. For a chosen set of \(r\) missing ranks, the allowed cards are

$$A(r)=4(13-r)+2=54-4r,$$

hence the forbidden cards are

$$B(r)=54-(54-4r)=4r.$$

The expectation becomes

$$\mathbb E[T]=\sum_{r=1}^{13}(-1)^{r+1}\binom{13}{r}\frac{55}{4r+1}.$$

Evaluating the sum gives

$$\mathbb E[T]\approx 29.05361725,$$

which matches the numerical checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations all evaluate the closed form directly. They first compute the total shoe size \(M\), then iterate over the possible counts of missing suits, missing ranks, and missing decks. For each triple \((s,r,d)\), they compute the binomial multiplicity, the parity sign from \(s+r+d\), the allowed-card count \(A(s,r,d)\), and the forbidden-card count \(B(s,r,d)\), then add the signed contribution

$$\binom{S}{s}\binom{R}{r}\binom{D}{d}\frac{M+1}{B(s,r,d)+1}.$$

The same framework also handles reduced modes by fixing a disabled category count to \(0\). The C++ implementation adds two sanity checks before printing the final answer: the rank-only checkpoint above, and a tiny shoe with \(D=S=R=2\) and \(J=0\), where the closed formula is compared against an exact exhaustive expectation.

Complexity Analysis

Let \(S_0\), \(R_0\), and \(D_0\) denote the active category sizes: each is either the full size or \(0\), depending on whether that category type is required. The closed-form evaluation uses

$$O\bigl((S_0+1)(R_0+1)(D_0+1)\bigr)$$

signed terms and \(O(1)\) memory. The implementations compute binomial coefficients with short multiplicative products, so there is only a tiny extra factor linear in the category size. For the actual Project Euler parameters, the full computation visits just

$$ (4+1)(13+1)(10+1)=770 $$

terms, which is easily fast enough.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=796
  2. Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
  3. Negative hypergeometric distribution: Wikipedia - Negative hypergeometric distribution
  4. Binomial coefficient: Wikipedia - Binomial coefficient

Problem 796 source code

C++

#include <cassert>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>
#include <functional>

using i64 = std::int64_t;

static i64 binom(int n, int k) {
    if (k < 0 || k > n) {
        return 0;
    }
    k = std::min(k, n - k);
    i64 res = 1;
    for (int i = 1; i <= k; ++i) {
        res = res * (n - k + i) / i;
    }
    return res;
}

static long double expected_formula(int decks,
                                    int suits,
                                    int ranks,
                                    int jokers_per_deck,
                                    bool need_suits,
                                    bool need_ranks,
                                    bool need_decks) {
    const i64 total_cards = static_cast<i64>(decks) * (static_cast<i64>(suits) * ranks + jokers_per_deck);

    const int max_s = need_suits ? suits : 0;
    const int max_r = need_ranks ? ranks : 0;
    const int max_d = need_decks ? decks : 0;

    long double ans = 0.0L;

    for (int s = 0; s <= max_s; ++s) {
        const i64 cs = need_suits ? binom(suits, s) : 1;
        for (int r = 0; r <= max_r; ++r) {
            const i64 cr = need_ranks ? binom(ranks, r) : 1;
            for (int d = 0; d <= max_d; ++d) {
                const i64 cd = need_decks ? binom(decks, d) : 1;

                const int missed = s + r + d;
                if (missed == 0) {
                    continue;
                }

                const i64 safe_per_deck = static_cast<i64>(suits - s) * (ranks - r) + jokers_per_deck;
                const i64 safe_cards = static_cast<i64>(decks - d) * safe_per_deck;
                const i64 bad_cards = total_cards - safe_cards;

                const long double weight = static_cast<long double>(total_cards + 1) /
                                           static_cast<long double>(bad_cards + 1);
                const long double mult = static_cast<long double>(cs) *
                                         static_cast<long double>(cr) *
                                         static_cast<long double>(cd);

                if (missed & 1) {
                    ans += mult * weight;
                } else {
                    ans -= mult * weight;
                }
            }
        }
    }

    return ans;
}

struct TinyCard {
    int deck_bit;
    int suit_bit;
    int rank_bit;
};

static long double tiny_bruteforce_expectation() {
    std::vector<TinyCard> cards;
    for (int d = 0; d < 2; ++d) {
        for (int s = 0; s < 2; ++s) {
            for (int r = 0; r < 2; ++r) {
                cards.push_back({1 << d, 1 << s, 1 << r});
            }
        }
    }

    const int n = static_cast<int>(cards.size());
    const int full_deck = (1 << 2) - 1;
    const int full_suit = (1 << 2) - 1;
    const int full_rank = (1 << 2) - 1;

    std::vector<long double> memo(1 << n, -1.0L);

    std::function<long double(int)> dfs = [&](int mask) -> long double {
        long double& ref = memo[mask];
        if (ref >= 0.0L) {
            return ref;
        }

        int dmask = 0;
        int smask = 0;
        int rmask = 0;
        int drawn = 0;

        for (int i = 0; i < n; ++i) {
            if ((mask >> i) & 1) {
                ++drawn;
                dmask |= cards[i].deck_bit;
                smask |= cards[i].suit_bit;
                rmask |= cards[i].rank_bit;
            }
        }

        if (dmask == full_deck && smask == full_suit && rmask == full_rank) {
            ref = 0.0L;
            return ref;
        }

        const int rem = n - drawn;
        long double sum = 0.0L;
        for (int i = 0; i < n; ++i) {
            if (((mask >> i) & 1) == 0) {
                sum += dfs(mask | (1 << i));
            }
        }

        ref = 1.0L + sum / static_cast<long double>(rem);
        return ref;
    };

    return dfs(0);
}

int main() {
    const long double check_rank = expected_formula(1, 4, 13, 2, false, true, false);
    assert(std::fabsl(check_rank - 29.05361725L) < 1e-8L);

    const long double tiny_formula = expected_formula(2, 2, 2, 0, true, true, true);
    const long double tiny_exact = tiny_bruteforce_expectation();
    assert(std::fabsl(tiny_formula - tiny_exact) < 1e-12L);

    const long double ans = expected_formula(10, 4, 13, 2, true, true, true);
    std::cout << std::fixed << std::setprecision(8) << static_cast<double>(ans) << '\n';
    return 0;
}

Python

def binom(n, k):
    if k < 0 or k > n:
        return 0
    k = min(k, n - k)
    res = 1
    for i in range(1, k + 1):
        res = res * (n - k + i) // i
    return res

def expected_formula(decks, suits, ranks, jokers_per_deck, need_suits, need_ranks, need_decks):
    total_cards = decks * (suits * ranks + jokers_per_deck)
    
    max_s = suits if need_suits else 0
    max_r = ranks if need_ranks else 0
    max_d = decks if need_decks else 0
    
    ans = 0.0
    
    for s in range(max_s + 1):
        cs = binom(suits, s) if need_suits else 1
        for r in range(max_r + 1):
            cr = binom(ranks, r) if need_ranks else 1
            for d in range(max_d + 1):
                cd = binom(decks, d) if need_decks else 1
                
                missed = s + r + d
                if missed == 0:
                    continue
                    
                safe_per_deck = (suits - s) * (ranks - r) + jokers_per_deck
                safe_cards = (decks - d) * safe_per_deck
                bad_cards = total_cards - safe_cards
                
                weight = (total_cards + 1.0) / (bad_cards + 1.0)
                mult = float(cs * cr * cd)
                
                if missed % 2 == 1:
                    ans += mult * weight
                else:
                    ans -= mult * weight
                    
    return ans

def solve():
    ans = expected_formula(10, 4, 13, 2, True, True, True)
    # The problem asks for 8 decimal places
    return "{:.8f}".format(ans)

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

Java

import java.text.DecimalFormat;
import java.text.DecimalFormatSymbols;
import java.util.Locale;

public class Euler796 {

    static long binom(int n, int k) {
        if (k < 0 || k > n) {
            return 0;
        }
        k = Math.min(k, n - k);
        long res = 1;
        for (int i = 1; i <= k; ++i) {
            res = res * (n - k + i) / i;
        }
        return res;
    }

    static double expectedFormula(int decks, int suits, int ranks, int jokersPerDeck,
            boolean needSuits, boolean needRanks, boolean needDecks) {
        long totalCards = (long) decks * ((long) suits * ranks + jokersPerDeck);

        int maxS = needSuits ? suits : 0;
        int maxR = needRanks ? ranks : 0;
        int maxD = needDecks ? decks : 0;

        double ans = 0.0; // Using double instead of long double is sufficient for 8 decimal places here
                          // usually, but if precision issues arise, BigDecimal might be needed. For this
                          // logic, Python float (double equivalent) passes test.

        for (int s = 0; s <= maxS; ++s) {
            long cs = needSuits ? binom(suits, s) : 1;
            for (int r = 0; r <= maxR; ++r) {
                long cr = needRanks ? binom(ranks, r) : 1;
                for (int d = 0; d <= maxD; ++d) {
                    long cd = needDecks ? binom(decks, d) : 1;

                    int missed = s + r + d;
                    if (missed == 0)
                        continue;

                    long safePerDeck = (long) (suits - s) * (ranks - r) + jokersPerDeck;
                    long safeCards = (long) (decks - d) * safePerDeck;
                    long badCards = totalCards - safeCards;

                    double weight = (double) (totalCards + 1) / (double) (badCards + 1);
                    double mult = (double) cs * (double) cr * (double) cd;

                    if ((missed & 1) == 1) {
                        ans += mult * weight;
                    } else {
                        ans -= mult * weight;
                    }
                }
            }
        }

        return ans;
    }

    public static String solve() {
        double ans = expectedFormula(10, 4, 13, 2, true, true, true);
        DecimalFormatSymbols symbols = new DecimalFormatSymbols(Locale.US);
        DecimalFormat df = new DecimalFormat("0.00000000", symbols);
        return df.format(ans);
    }

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