Let $p$ be an odd prime number. Find the number of subsets $\mathbb{A}$ of the set $\{1, 2, \cdots, 2p\}$ such that $\mathbb{A}$ has exactly $p$ elements and the sum of all elements in $\mathbb{A}$ is divisible by $p$.
Consider the polynomial $$f_a(x)=(1+ax)(1+a^2x)(1+a^3x)\cdots(1+a^{2p}x)$$
When the right side is expanded, let $c_{n,k}$ be the number of the terms in the form $(a^{i_1}x)(a^{i_2}x)\cdots(a^{i_k}x)$, where $i_1$, $i_2$, $\cdots$, $i_k$ are integers such that
- $1\le i_1 < i_2 < \cdots < i_k \le 2p$, and
- $i_1 + i_2 +\cdots + i_k =n$
Then, it must hold that $$f_a(x)=1+\sum_{k=1}^{2p}\left(\sum_{n=1}^{\infty}c_{n,k}a^n\right)x^k$$
The coefficient of the term $x^p$ in $f_a(x)$ is simply $$c_{n,p}=\sum_{n=1}^{\infty}c_{n,p}a^n$$
Then the answer to the given question is $$C=c_{p,p} + c_{2p,p} + c_{3p, p} + \cdots$$
This is a sequence of coefficients with interval $p$. Let $\omega \ne 1$ be a complex root to $x^p=1$. Then using the complex number method , it can be derived that $$C=\frac{1}{p}\sum_{j=0}^{p-1}\sum_{n=1}^{\infty}c_{n,p}\omega^{nj}$$
The right side of this relation is the coefficient of $x^p$ $$\frac{1}{p}\sum_{j=0}^{p-1}f_{a=\omega^j}(x)=\frac{1}{p}\sum_{j=0}^{p-1}(1+\omega^jx)(1+\omega^{2j}x)\cdots(1+\omega^{2pj}x)$$
When $j=0$, the product at right side above equals $(1+x)^{2p}$. For $1\le j \le p-1$, let $\theta=\omega^j\ne 1$ and $\theta^p=1$. Then $$\begin{array}{rl} & (1+\omega^jx)(1+\omega^{2j}x)\cdots(1+\omega^{2pj}x) \\ =&\left((1+\theta x)(1+\theta^2 x)\cdots(1+\theta^{p-1} x)(1+x)\right)\\ &\left((1+\theta x)(1+\theta^2 x)\cdots(1+\theta^{p-1}x)(1+x)\right) \\=&(1+x^p)^2 \end{array}$$
Therefore, the previous expression equals $$ \frac{1}{p}\left((1+x)^{2p} +(p-1)(1+x^p)^2\right) $$
Accordingly, the coefficient of $x^p$, i.e. the final answer, is $$\boxed{\frac{1}{p}\left(\binom{2p}{p} +2(p-1)\right)}$$