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.
I will answer a much more general problem:
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,
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.
${}^1$ Stanley, Enumerative Combinatorics, volume 2, Chapter 7, Section 5.