Solution:
It can be easily checked that the functions $f(x)=x$, $f(x)=-x$ and $f(x)=\frac{1}{2}-x$ satisfy the given condition. We will show that these are the only functions doing so. Let $y=-f(x)$ in the original equation to obtain
$$f\left(2 f(x)^{2}+2 x f(-f(x))\right)=0$$
for all $x$. In particular, $0$ is a value of $f$. Suppose that $u$ and $v$ are such that $f(u)=0=f(v)$. Plugging $x=u$ or $v$ and $y=u$ or $v$ in the original equations we get $f\left(u^{2}\right)=u^{2}$, $f\left(u^{2}\right)=uv$, $f\left(v^{2}\right)=uv$ and $f\left(v^{2}\right)=v^{2}$. We conclude that $u^{2}=uv=v^{2}$ and hence $u=v$. So there is exactly one $a$ mapped to $0$, and
$$f(x)^{2}+x f(-f(x))=\frac{a}{2}$$
for all $x$.
Suppose that $f\left(x_{1}\right)=f\left(x_{2}\right) \neq 0$ for some $x_{1}$ and $x_{2}$. Using (*) we obtain $x_{1} f\left(-f\left(x_{1}\right)\right)=x_{2} f\left(-f\left(x_{2}\right)\right)=x_{2} f\left(-f\left(x_{1}\right)\right)$ and hence either $x_{1}=x_{2}$ or $f\left(x_{1}\right)=f\left(x_{2}\right)=-a$. In the second case, letting $x=a$ and $y=x_{1}$ in the original equation we get $f\left(x_{1}^{2}-2 a^{2}\right)=0$, hence $x_{1}^{2}-2 a^{2}=a$. Similarly, $x_{2}^{2}-2 a^{2}=a$, and it follows that $x_{1}=x_{2}$ or $x_{1}=-x_{2}$ in this case.
Using the symmetry of the original equation we have
$$f\left(f(x)^{2}+y^{2}+2 x f(y)\right)=(x+f(y))(y+f(x))=f\left(f(y)^{2}+x^{2}+2 y f(x)\right)$$
for all $x$ and $y$. Suppose $f(x)^{2}+y^{2}+2 x f(y) \neq f(y)^{2}+x^{2}+2 y f(x)$ for some $x$ and $y$. Then by the observations above, $(x+f(y))(y+f(x)) \neq 0$ and $f(x)^{2}+y^{2}+2 x f(y)=-\left(f(y)^{2}+x^{2}+2 y f(x)\right)$. But these conditions are contradictory as the second one can be rewritten as $(f(x)+y)^{2}+(f(y)+x)^{2}=0$.
Therefore from $(**)$ now it follows that
$$f(x)^{2}+y^{2}+2 x f(y)=f(y)^{2}+x^{2}+2 y f(x)$$
for all $x$ and $y$. In particular, letting $y=0$ we obtain $f(x)^{2}=(f(0)-x)^{2}$ for all $x$. Let $f(x)=s(x)(f(0)-x)$ where $s: \mathbf{R} \rightarrow \{1,-1\}$. Plugging this in $(***)$ gives
$$x(y s(y)+f(0)(1-s(y)))=y(x s(x)+f(0)(1-s(x)))$$
for all $x$ and $y$. So $s(x)+f(0)(1-s(x))/x$ must be constant for $x \neq 0$.
If $f(0)=0$ it follows that $s(x)$ is constant for $x \neq 0$, and therefore either $f(x)=x$ for all $x$ or $f(x)=-x$ for all $x$. Suppose that $f(0) \neq 0$. If $s(x)$ is $-1$ for all $x \neq 0$, then $-1+2 f(0)/x$ must be constant for all $x \neq 0$, which is not possible. On the other hand, if there exist nonzero $x$ and $y$ such that $s(x)=-1$ and $s(y)=1$, then $-1+2 f(0)/x=1$. That is, there can be only one such $x$, that $x$ is $f(0)$, and hence $f(x)=f(0)-x$ for all $x$. Putting this back in the original equation gives $2 f(0)^{2}=f(0)$ and hence $f(0)=1/2$. We are done.