Finding the covariance of a mixed pair of r.v.'s given one's distribution and a conditional distribution $
\newcommand{\E}{\mathbb{E}}
\newcommand{\dd}{\mathrm{d}}
\newcommand{\al}[1]{\begin{align*} #1 \end{align*}}
\newcommand{\set}[1]{\left\{ #1 \right\}}
\newcommand{\para}[1]{\left( #1 \right)}
\newcommand{\br}[1]{\left[ #1 \right]}
\newcommand{\B}[1]{\Big[ #1 \Big]}
\newcommand{\var}{\mathrm{var}}
\newcommand{\Beta}{\mathrm{B}}
\newcommand{\b}[1]{\big[ #1 \big]}
\newcommand{\cov}{\mathrm{cov}}
$
The Problem: Let $X,Y$ be jointly distributed r.v.'s such that
\begin{align*}
Y \mid X = x &\sim \mathrm{Binomial}(n,x) \\
X &\sim \mathrm{Uniform}(0,1)
\end{align*}
Among other things, I seek to calculate $\mathrm{cov}(X,Y)$, without using Bayesian techniques.
(This comes to me in the form of a homework assignment, exercise $35$ in Chapter $2$ of Allan Gut's An Intermediate Course in Probability, second edition.)

My Work So Far: Note that, with the given info, we have
$$\al{
f_{Y \mid X = x}(k) &= \binom n k x^k (1-x)^{n-k}, \; k \in \set{0,1,\cdots,n} \\
f_X(x) &= 1, \; x \in (0,1)\\
\E\br{Y \mid X = x} &= nx \\
\E\br{X} &= \frac 1 2\\
\var \para{Y \mid X} &= nx(1-x) \\
\var \para{X} &= \frac{1}{12}
}$$
We'll need to calculate $Y$'s expectation as well:
$$\al{
\E \b{Y}
&= \E \B{ \E \b{Y \mid X} } \tag{known result} \\
&= \E \br{ nX } \tag{known value} \\
&= n \cdot \E \br X \tag{linearity}\\
&= \frac n 2 \tag{known value} 
}$$
We can calculate the covariance as so:
$$\al{
\cov(X,Y)
&= \E \B{ \big( X - \E[X] \big) \big( Y - \E[Y] \big) } \tag{definition} \\
&= \E \br{ \para{ X - \frac 1 2 } \para{ Y - \frac n 2 } } \tag{known values} \\
&= \E \br{ XY - \frac n 2 X - \frac 1 2 Y +\frac n 4 }  \tag{expand} \\
&= \E \br{XY} - \frac n 2 \E \br{X} - \frac 1 2 \E \br{Y} + \E \br{ \frac n 4 } \tag{linearity}
}$$
Everything except one thing works to me: how would one find $\E \br{XY}$? We could note that...
$$
f_{Y \mid X = x}(y) = \frac{f_{X,Y}(x,y)}{f_X(x)} \implies f_{X,Y}(x,y) = f_{Y \mid X = x}(y) \cdot f_X(x)
$$
Then I think we have this:
$$\al{
\E \br{XY}
&= \sum_{y=0}^n \int_0^1 xy \binom n y x^y (1-x)^{n-y} \, \dd x\\
&= \sum_{y=0}^n y \binom n y \int_0^1 x^{y+1} (1-x)^{n-y} \, \dd x\\
&= \sum_{y=0}^n y \binom n y \Beta(y+2,n-y+1)
}$$
where $\Beta(\cdot,\cdot)$ denotes the usual beta function. And I think this form follows since, essentially, covariance is like
$$
\cov(X,Y) = \kern-30pt \sum \limits_{(x,y) \,\in\, \mathrm{supp}(X) \times \mathrm{supp}(Y)} \kern-30pt  xy \, f_{X,Y}(x,y)
$$
(though with a continuous variable like $X$ we would interpret the respective sum as an integral, though we use sum for discrete r.v.'s like $Y \mid X$).

My Question: Simply put, is this correct?
It doesn't feel like it is at all. I'm struggling with the fact that $Y \mid X$ and $X$ are different "kinds" (discrete versus continuous) of distributions, and how to handle them in a clean way. And the final answer just feels too complicated to be right.
Does anyone have any insights they can give me on this problem? Particularly with respect to finding $\E[XY]$ within the given parameters.

Edit: (3/24/2021)
The suggestion of using that
$$
\E \br{XY} = \E \B{\E \br{XY \mid X}} = \E \B{X \cdot \E \br{Y \mid Y}}
$$
(per StubbornAtom in the comments) has helped in achieving the right answer, namely
$$
\cov(X,Y) = \frac{n}{12}
$$
However, I'm not totally sure if it's kosher to use this in the class, at least without proof; I'm pretty sure this identity hasn't come up at any event. Of course, one can prove it, and some links have been provided hinting at how to do so. But of course this begs the question -- how might one achieve the desired answer, without this identity?
Of course, perhaps the end result is "you just get that nasty expression you derived," but I figured it'd be worth asking to see some alternative approaches. Or maybe my answer can be somehow easily shown to get the same answer anyways?
 A: You can in fact evaluate the sum directly; it is just a bit tedious compared to using iterated expectations.
$$\begin{align}
\operatorname{E}[XY]
&= \sum_{y=0}^n \int_{x=0}^1 xy \binom{n}{y} x^y (1-x)^{n-y} \, dx \\
&= \sum_{y=0}^n y \binom{n}{y} \int_{x=0}^1 x^{y+1} (1-x)^{n-y} \, dx \\
&= \sum_{y=0}^n y \binom{n}{y} \frac{\Gamma(y+2)\Gamma(n-y+1)}{\Gamma(n+3)} \\
&= \sum_{y=0}^n y \frac{n!}{y!(n-y)!} \frac{(y+1)!(n-y)!}{(n+2)!} \\
&= \sum_{y=0}^n \frac{y(y+1)}{(n+1)(n+2)}. \\
\end{align}$$
At this point, we can evaluate the sum directly, using the familiar formulas $\sum y = n(n+1)/2$ and $\sum y^2 = n(n+1)(2n+1)/6$, or we can use the hockey stick identity
$$\sum_{y=0}^n \frac{y(y+1)}{(n+1)(n+2)} = \frac{2}{(n+1)(n+2)} \sum_{y=1}^n \binom{y+1}{2} = \frac{2}{(n+1)(n+2)} \binom{n+2}{3} = \frac{n}{3}.$$
Combining this with the other terms in the covariance expression, we get $$\operatorname{Cov}[X,Y] = \frac{n}{3} - \frac{n}{4} = \frac{n}{12}$$ as claimed.
As an exercise, can you evaluate $\operatorname{E}[XY]$ with the summation and integration order exchanged; i.e., how would we proceed with $$\int_{x=0}^1 \sum_{y=0}^n xy \binom{n}{y} x^y (1-x)^{n-y} \, dx?$$  In fact, is this a more natural way to evaluate this?  Does it mimic the iterated expectation evaluation $$\operatorname{E}[XY] = \operatorname{E}[\operatorname{E}[XY \mid X]]?$$
