Probability of no void in a bridge table

Question

52 cards, 13 of four different suits are randomly distributed to 4 people, each person receiving exactly 13 cards. What is the probability that all 4 persons get at least one card of each suit? (That is there is no void anywhere in the table?)

Note: The question is not about one particular person having a void but about no void in the bridge table. It is possible to work out the solution of the above problem from this by considering case by case break-up. But this is very tedious. Is there any elegant solution?

Answer 1

The following solution gives an answer which is relatively concise, but requires a fair amount of computing by a computer algebra system (a few minutes on my aging PC).

There are $\binom{52}{13 \; 13 \; 13 \; 13}=53644737765488792839237440000$ (a multinomial coefficient) ways to deal four bridge hands, all of which we assume are equally likely. We want to count the number of deals in which no player is void in any suit.

To that end, we will find a generating function for the answer to a more general problem: How many ways are there to deal $i$ cards to player 1, $j$ cards to player 2, $k$ cards to player three, and $\ell$ cards to player 4 with no player having a void in any suit? Let's say that number is $a_{ijk \ell}$, and define $f(w,x,y,z)$ to be the four-variable generating function of $a_{ijk \ell}$, i.e. $$f(w,x,y,z) = \sum_{i,j,k,\ell} a_{ijk \ell} w^i x^j y^k z^{\ell}$$ Now define $$S_0(w,x,y,z) = (1+w+x+y+z)^{13} $$ $$S_1(w,x,y,z) = (1+x+y+z)^{13} +(1+w+y+z)^{13}+ (1+w+x+z)^{13}+(1+w+x+y)^{13}$$ $$S_2(w,x,y,z) = (1+w+x)^{13}+(1+w+y)^{13}+(1+w+z)^{13}+(1+x+y)^{13}+(1+x+z)^{13}+(1+y+z)^{13}$$ $$S_3(w,x,y,z) = (1+w)^{13}+(1+x)^{13}+(1+y)^{13}+(1+z)^{13}$$ and $$S_4 = 1$$ The point of these definitions is that $S_0-S_1+S_2-S_3+S_4$ consists of the terms in $(1+w+x+y+z)^{13}$ of the form $w^i x^j y^k z^{\ell}$ with $i,j,k$ and $\ell$ all positive. Then $$f(w,x,y,z) = [S_0(w,x,y,z)-S_1(w,x,y,z)+S_2(w,x,y,z)-S_3(w,x,y,z)+S_4(w,x,y,z)]^4$$ A computer algebra system can expand $f(w,x,y,z)$ and find the coefficient of $w^{13} x^{13} y^{13} z^{13}$, which is the number of bridge deals in which no player has a void in any suit. Courtesy of Mathematica, we have $$[w^{13} x^{13} y^{13} z^{13}] f(x,w,y,z) = 43786641333457372963248979200$$ So the probability that no player has a void in any suit is $$\frac{[w^{13} x^{13} y^{13} z^{13}] f(x,w,y,z)}{\binom{52}{13 \; 13 \; 13 \; 13}} \approx \boxed{0.8162336728}$$

Readers unfamiliar with generating functions may find the answers to this question interesting: How can I learn about generating functions?

Answer 2

The inelegant approach is easy to do with a computer program, so in the hope of spotting some useful patterns which can direct a search for an elegant solution I've generalised it to the case of $p$ players, $s = p$ suits (to keep it simple), and $c \ge p$ cards per suit.

The simplest case is $$\begin{array}{c|c|c} n & c & P(\textrm{no voids}) \\ \hline 2 & 2 & \frac{4}{6} = \frac{2}{3} \\ 2 & 3 & \frac{18}{20} = \frac{9}{10} \\ 2 & 4 & \frac{68}{70} = \frac{34}{35} \\ 2 & 5 & \frac{250}{252} = \frac{125}{126} \\ 2 & 6 & \frac{922}{924} = \frac{461}{462} \\ 2 & 7 & \frac{3430}{3432} = \frac{1715}{1716} \end{array}$$ where probabilities are given in unreduced and reduced form, and I think we can easily prove that for $n=2$ the probability of no voids is $$1 - \frac2{\binom{2c}{c}}$$ The two ways of getting voids are for player 1 to have all of one suit, or for player 1 to have all of the other suit.

---

For $n=3$ it's already quite a bit messier:

$$\begin{array}{c|c|c} n & c & P(\textrm{no voids}) \\ \hline 3 & 3 & \frac{216}{1680} = \frac{9}{70} \\ 3 & 4 & \frac{10368}{34650} = \frac{576}{1925} \\ 3 & 5 & \frac{372000}{756756} = \frac{31000}{63063} \\ 3 & 6 & \frac{11259000}{17153136} = \frac{156375}{238238} \\ 3 & 7 & \frac{310482228}{399072960} = \frac{25873519}{33256080} \\ 3 & 8 & \frac{8146518912}{9465511770} = \frac{452584384}{525861765} \\ 3 & 9 & \frac{208310706288}{227873431500} = \frac{17359225524}{18989452625} \\ 3 & 10 & \frac{5261019174000}{5550996791340} = \frac{1538309700}{1623098477} \\ 3 & 11 & \frac{132223145196648}{136526995463040} = \frac{5509297716527}{5688624810960} \\ 3 & 12 & \frac{3320705738396376}{3384731762521200} = \frac{46120913033283}{47010163368350} \end{array}$$

Is it already necessary to do an inclusion-exclusion in order to get a closed form?

---

If anyone wants to play with my code, here you go:

#!/usr/bin/python3

# https://math.stackexchange.com/q/3443423

from collections import defaultdict
from fractions import Fraction


def binom(n, k):
    nCk = 1
    for i in range(k):
        nCk = nCk * (n - i) // (i + 1)
    return nCk


def distributionShapes(suitQuantities, handSize, playersRemaining):
    if len(suitQuantities) == 0:
        if handSize == 0:
            yield []
        return

    # Draw n of the first suit and recurse.
    #   n >= 1 because we don't have a void in the first suit
    #   handSize - n >= len(suitQuantities) - 1 because we don't have voids in the other suits either
    #   suitQuantities[0] - n >= playersRemaining because the other players don't have voids in the first suit
    maxSuitDraw = min(handSize - len(suitQuantities) + 1, suitQuantities[0] - playersRemaining)
    for n in range(1, maxSuitDraw + 1):
        for tail in distributionShapes(suitQuantities[1:], handSize - n, playersRemaining):
            yield [n] + tail


def distributionWeights(suitQuantities, handSize, playersRemaining):
    """
    Returns a dict from card counts per suit to number of ways of drawing that hand distribution.
    We guarantee to leave at least playersRemaining cards in each suit.
    """
    weights = {}
    for shape in distributionShapes(suitQuantities, handSize, playersRemaining):
        weight = 1
        for n, k in zip(suitQuantities, shape):
            weight *= binom(n, k)
        weights[tuple(shape)] = weight
    return weights


def prob(players, cards):
    if cards < players:
        return 0

    suits = players
    partial = { tuple([0] * suits) : 1 }

    denominator = 1
    for i in range(1, players + 1):
        nextPartial = defaultdict(int)
        for shapei in partial:
            pj = distributionWeights([cards - x for x in shapei], cards, players - i)
            for shapej in pj:
                # Combined shape
                shapeij = [x + y for x, y in zip(shapei, shapej)]
                nextPartial[tuple(shapeij)] += partial[shapei] * pj[shapej]
        partial = nextPartial
        denominator *= binom(i * cards, cards)

    numerator = partial[tuple([cards] * suits)]
    return [numerator / denominator, numerator, denominator, Fraction(numerator, denominator)]

for p in range(2, 5):
    for c in range(0, 10):
        print(p, p + c, prob(p, p + c))