Counting bridge deals if we only distinguish face cards Bridge is a game of four players in which each player is dealt 13 cards from a standard 52 card deck. Bridge players (such as myself) are interested in the number of possible deals, where each player is distinct. This can be counted by
$$\binom{52}{13}\binom{39}{13}\binom{26}{13}\binom{13}{13}=5.364\times10^{28}$$
However, this number is misleadingly large, since bridge players usually only care about the face cards (jack, queen, king, and ace) in each suit. We often consider the cards with denominations 2-10 as indistinguishable. Supposing we distinguish only face cards, what is the number of possible deals? 
This source puts the figure at $8.110\times10^{15}$ based on a computer program. I am curious if there is a more elegant mathematical solution.
 A: (not enough rep to comment)
FWIW, Wikipedia cites another site which gives the same number as the OP's link. Both this site and the OP's link say that there doesn't seem to be any simple formula to answer this question.
A: There are $16$ cards of interest, so distributing them would give you $4^{16}$ possibilities.  The problem is that no hand can contain more than $13$ cards, so we subtract off the ones that have a hand of more than that.  There are $4$ hands where one hand has all $16$ of the interesting cards.  For hands with $15$ cards, there are $4$ ways to choose the hand that gets them, $16$ ways to choose the card taken out, and $3$ ways to choose the hand that gets the odd one.  For hands with $14$ cards, there are $4$ ways to choose the hand that gets them, $16 \choose 2$ ways to choose the two other cards, and $3^2$ ways to choose which hand gets them.  The grand total is 
$$4^{16}-4-4\cdot16\cdot 3-4{16 \choose 2}3^2=4294962780\approx 4.3\cdot 10^9$$
The source you link considers $10$s to be distinct as well, which is why the number is so much higher.  A similar analysis works if you considered $10$s.  I get
$$4^{20}-4\sum_{k=0}^6{20 \choose k}3^k=1099381833744\approx 1.1\cdot 10^{12}$$
A: I will answer a much more general problem:

Suppose you have a collection of playing cards, each falling into one of $k$ types, so that cards of the same type are indistinguishable. For each $i\in \{1,\dots,k\}$, there are $n_i$ cards of type $i$. How many ways are there to deal these cards to $p$ players so that for each $j\in \{1,\dots,p\}$, the $j^{th}$ player receives $m_j$ cards?

This can be solved using generating functions. Specifically, the enumerator for dealing $n_i$ indistinguishable cards to $p$ players is
$$
\sum_{a_1+\dots+a_p=n_i}x_1^{a_1}\dots x_p^{a_p}=h_{n_1}(x_1,x_2,\dots,x_p),
$$
where $h_{n_i}$ is the $(n_i)^{th}$ homogenous symmetric polynomial in $p$ variables. Each summand corresponds to a particular way to deal the $n_i$ indistinguishable cards to the $p$ players; the summand $x_1^{a_1}\dots x_p^{a_p}$ corresponds to giving $a_j$ cards of type $i$ to the $j^{th}$ player, for $j\in \{1,\dots,p\}$. 
Furthermore, the action of dealing out all of the cards is simply accomplished by multiplying the enumerators for each type of card together. Each summand in this product with powers $x_1^{b_1}\cdots x_p^{b_p}$ comes with a coefficient equal to the number of ways to deal the cards so that player $j$ receives $b_j$ cards for $j\in \{1,\dots,p\}$. Therefore, 

The number of dealings is equal to the coefficient of $x_1^{m_1}\cdots x_p^{m_p}$ in $\prod_{i=1}^k h_{n_i}(x_1,\dots,x_p)$.

In your case, you want the coefficient of $x_1^{13}x_2^{13}x_3^{13}x_4^{13}$ in $h_{36}\cdot h_1^{16}$. If you considered the lower ranks within a suit to be indistinguishable, yet distinct from other suits (as the author of the linked webpage did), then you would want the coefficient of the same monomial in $h_{9}^4\cdot h_1^{16}$.
We can leverage the knowledge of symmetric functions to attack this computationally. Specifically, let $\lambda$ be the decreasing sorted list of numbers of cards in each type, $(n_1,n_2,\dots,n_k)$, and let $\mu$ be the sorted list $(m_1,m_2,\dots,m_p)$. It can be shown${}^1$ that the coefficient of $x_1^{m_1}\cdots x_p^{m_p}$ in $\prod_{i=1}^k h_{n_i}(x_1,\dots,x_p)$ is equal to $N_{\lambda,\mu}$, defined to be the number of $k\times p$ matrices with entries in the nonnegative integers whose vector of row sums is equal to $\lambda$ and whose vector of column sums is equal to $\mu$. 
I am not sure what the best way to compute $N_{\lambda,\mu}$ is. This procedure is built into the symmetric functions package of Sage. The following program computes the number of deals  when the ranks $2$ to $9$ are indistinguishable within a suit but different suits are distinct. It gives the count of $8110864720503360$, taking about $8$ minutes to run on the online interpreter CoCalc. This agrees with your source. Furthermore, the program is easily configurable to work for any number of suits, players, and number of ranks considered indistinguishable. 
from time import time

Sym = SymmetricFunctions(QQ)
m = Sym.monomial()
h = Sym.homogeneous()

suits      = 4
low_ranks  = 9
high_ranks = 4
players    = 4

lamb = [low_ranks]*suits + [1]*(high_ranks*suits)
targ = [(low_ranks + high_ranks) * suits // players] * (players)

t0 = time()
print('Number of hands:    ', m(h(lamb)).coefficient(targ) )
print('Seconds to compute: ', time() - t0 )

${}^1$ Stanley, Enumerative Combinatorics, volume 2, Chapter 7, Section 5.
