They are indeed not independent as you stated. We know that two random variables $X, Y$ (for which the probability densities $f_X$, $f_Y$ and the joint probability density $f_{X,Y}$ exist) are independent iff:
$f_{X,Y}(x, y) = f_X(x)f_Y(y)$ for all $x, y$
An here it is easy to find a counterexample:
$f_{X,Y}(1/3, 1/3) = 0$
$f_X(1/3) = 1$, $f_Y(1/3) = 1$
So the condition for the variables to be independent does not hold.
However when we are talking about correlation we are usually referring to the Pearson Correlation Coefficient which measures $linear$ relationship between two random variables. Here there is no linear relationship between the two variables, so the correlation is $0$. You can see that in the bottom of the following picture some examples of obviously dependent random variables (like in your example), where the correlation is $0$.
In real-world examples the Pearson Coefficient works well in many cases. For when it doesn't, there are some other correlation measures such as Spearman's rank coefficient, Kendall rank coefficient and others that can be used (although there is no guarantee that they will always manage to point out the existence of a relationship between your variables).
Edit:
As user asked, I am explaining how the Pearson Coefficient can be calculated here. The formula is:
$ρ_{X, Y} = \frac{cov(X, Y)}{σ_Xσ_Y} = \frac{E[XY] - E[X]E[Y]}{σ_Xσ_Y}$
Here the most parts of the formula can be easily calculated:
$E[X] = E[Y] = 1/2$
$σ_X = σ_Y = \sqrt{1/12}$
The only difficult part is the calculation of $E[XY]$:
$E[XY] = \int_0^1 \int_0^1 xy f_{XY}(x, y) dx dy$
The problem here is that $f_{XY}$ isn't continuous in $[0, 1] \times [0, 1]$, so to calculate this integral we need to calculate the eight integrals for the eight areas where the probability density function is non-zero, and add them:
$\int_0^1 \int_0^1 xy f_{XY}(x, y) dx dy = \\
\int_{1/4}^{1/2} \int_0^{1/4} xy f_{XY}(x, y) dx dy +
\int_0^{1/4} \int_{1/4}^{1/2} xy f_{XY}(x, y) dx dy + \\
\int_{1/4}^{1/2} \int_{1/2}^{3/4} xy f_{XY}(x, y) dx dy +
\int_{1/2}^{3/4} \int_{1/4}^{1/2} xy f_{XY}(x, y) dx dy + \\
\int_{0}^{1/4} \int_{3/4}^{1} xy f_{XY}(x, y) dx dy +
\int_{3/4}^{1} \int_{0}^{1/4} xy f_{XY}(x, y) dx dy + \\
\int_{1/2}^{3/4} \int_{3/4}^{1} xy f_{XY}(x, y) dx dy +
\int_{3/4}^{1} \int_{1/2}^{3/4} xy f_{XY}(x, y) dx dy
$
Here you can notice that due to the symmetry in your PDF, the two integrals in each of the lines have the same value, so you only need to calculate four integrals instead of eight. Perhaps there are smarter things you can do to calculate the original integral with less sweat, but I can't think of anything else.
Best Answer
Answer for the first part: it is assumed, though not stated explicitly, that $EX^{2} $ and $EY^{2}$ are finite. To show that covariance is $0$ you need $EX^{2}-EY^{2}=(EX)^{2}-(EY)^{2}$ (as you have already observed). Just write this as $EX^{2}-(EX)^{2}=EY^{2}-(EY)^{2}$ or $var(X)=var(Y)$. This equation holds becasue $X$ and $Y$ have the same distribution.
Let $X$ and $Y$ be i.i.d. with Cauchy distribution. Then it is easy to verify that $Ee^{i(tU+sV)}=Ee^{itU} Ee^{isV}$ is not satisfied. [The characteristic function of Cauchy distribution is $e^{-|t|}$]. This gives an example where $U$ and $V$ are not independent.
The last part is a well known result which involves solving the functional equation $\phi (t+s)\phi (t-s)=\phi(t)^{2} |\phi (s)|^{2}$.