How to solve this candy collectors puzzle? I am trying to solve this Jane Street puzzle with the solution given here.
For convenience I will rewrite the puzzle here:

Five children went trick-or-treating together and decided to randomly split their candy haul at the end of the night.  As it turned out, they got a total of 25 pieces of candy, 5 copies each of 5 different types (they live in a small town).  They distribute the candies by choosing an ordering of the 25 uniformly at random from all shufflings, and then giving the first 5 to the first child, the second 5 to the second, and so on.
What is the probability that each child has one type of candy that they have strictly more of than every other trick-or-treater? Give your (exact!) answer in a lowest terms fraction.

Here is what I have done so far.

Label each candy type as either candy 0, 1, 2, 3, or 4. Label the children as either child A, B, C, D, or E. Suppose each candy is taken randomly out of a bag and placed in order to produce a sequence like this
[[2, 0, 2, 0, 4], [0, 4, 0, 4, 3], [3, 3, 1, 2, 1], [2, 1, 2, 4, 3], [1, 3, 1, 4, 0]],
where the first 5 go to child A, the second 5 go to child B and so on. The number of permutations of this list is
$$\frac{25!}{(5!)^5}.$$
We need to calculate the number of permutations where every child dominates in 1 and only 1 candy type. We can create a table which looks like this
$$\begin{matrix}
  & A & B & C & D & E \\
0 & 2 & 2 & 0 & 0 & 1 \\
1 & 0 & 0 & 2 & 1 & 2 \\
2 & 2 & 0 & 1 & 2 & 0 \\
3 & 0 & 1 & 2 & 1 & 1 \\
4 & 1 & 2 & 0 & 1 & 1 \\
\end{matrix}$$
where each element of the table shows how much of each candy type each child has. For example the sequence above shows that child A has 2 of candy 0, none of candy 1, 2 of candy 2, none of candy 3 and 1 of candy 4.
For a child to dominate a candy type they must have at least 2 of that type. Therefore we can initialise our table like this
$$\begin{matrix}
  & A & B & C & D & E \\
0 & 2 & * & * & * & * \\
1 & * & 2 & * & * & * \\
2 & * & * & 2 & * & * \\
3 & * & * & * & 2 & * \\
4 & * & * & * & * & 2 \\
\end{matrix}$$
Now we need to create some code to solve all the possible ways to fill in the $*$ such that every child domiantes one candy type (and the sum of each row and column equals 5) where the diagonal elements can have values $\in \{2, 3, 4, 5\}$ and the off-diagonal elements can have values $\in \{0, 1, 2\}$. Where we multiply each result by
$$\prod_{i=1}^{5}\frac{5!}{\prod_{j=1}^{5} a_{ij}!},$$
where $a_{ij}$ corresponds to the value in row $i$, column $j$, to account for the different permutations. Finally, we multiply the final result by $5!$ to account for the different permutations associated with different children dominating different candy types.

Is this the right way of approaching the problem? I don't know how I would make such a piece of code I was thinking of using something similar to the backtracking algorithms which are used to solve sodoku problems.
 A: Here is a solution of a very different flavor, which uses exponential generating functions, along with any CAS which can multiply multivariate polynomials, such as Mathematica.
Let $f_1$ be the following polynomial in five variables:
\begin{align}
f_1(a,b,c,d,e) &= \frac{a^5}{5!}+\frac{a^4}{4!}(b+c+d+e)+\\
&\quad\frac{a^3}{3!}\left(\frac{b^2}{2!}+\frac{c^2}{2!}+\frac{d^2}{2!}+\frac{e^2}{2!}+bc+bd+be+cd+ce+de\right)
\\&\quad\frac{a^2}{2!}(bcd+bce+bde+cde)
\end{align}
Roughly the variables represent children, and each monomial represents the ways to distribute candy number $1$ so that child $a$ has more than anyone else. For example, one of the monomials is $\frac{a^3}{3!}cd$, which means child $a$ gets $3$ of candy number one, while $c$ and $d$ each one of candy number one.
Similarly, let $f_2$ represent the number of ways to distribute candy $2$ so that child $b$ has more than anyone else. Concretely, $f_2(a,b,c,d,e)=f_1(b,c,d,e,a)$. Then, do the same for $f_3,f_4,f_5$ for candies $3,4$ and $5$. It turns out that multiplying these polynomials, then multiplying by $5!^5$ represents the number of ways to distribute each of the candies. This is the magic of exponential generating functions; to learn more of why this works, I suggest looking at Chapter 3 of Herbert Wilf's excellent text generatingfunctionology, available for free online from the author's website.
Since we want each child to receive five candies total, the number of valid ways to distribute the candies is the coefficient of $(abcde)^5$ in the product $f_1f_2f_3f_4f_5$. You then multiply by $5!$ one again to account the for other ways to assign which child has the most of which candy, then divide by the total number of candy distributions, $(25!)/(5!)^2$. Here it the Mathematica code which accomplishes this, which you can run online here.
f[a_,b_,c_,d_,e_] := 
a^5/5! + 
a^4/4!*(b + c + d + e) + 
a^3/3!*((b^2+c^2+d^2+e^2)/2!+ b c + b d + b e + c d + c e + d e) + 
a^2/2!*(b c d + b c e + b d e + c d e);

EGF = f[a,b,c,d,e]*f[b,c,d,e,a]*f[c,d,e,a,b]*f[d,e,a,b,c]*f[e,a,b,c,d];

numerator = 5!^6 * Coefficient[EGF, (a*b*c*d*e)^5];

Print[numerator / (25!/ 5!^5)]

