Find all functions $f:\mathbb{R}\rightarrow\mathbb{R}$ where $f(xf(y)+f(x)+y)=xy+f(x)+f(y)$ 
Find all functions $f:\mathbb{R}\rightarrow\mathbb{R}$ for two real numbers $x$ and $y$ where $f(xf(y)+f(x)+y)=xy+f(x)+f(y)$

For $x=0$ and $y=-f(0)$ then $f(-f(0))=0$. So, there is a real root $r_0$ for function $f$.
For $x=r_0$ and $y=r_0$ we have $r_0^2=0$, so $f(0)=0$ and zero root is unique.
Please help me to complete the proof.
 A: Serveral good ideas have been presented in the solutions and comments. But I was surprised that seemingly noone tried to use calculus. Therefore I have tried to fill the gap.
However, as @PierreCarre correctly observed in a comment this is an incomplete solution to the problem of finding all functions for which the equation
$$f(xf(y)+f(x)+y)=xy+f(x)+f(y)\tag{1}$$
holds, because it assumes that all derivatives of $f$ with repect to its argument exist.
But it might be still of interest.
Main idea
Under this supposition taking successive derivatives with respect to $x$ at the point $(x,y) = (-1,-1)$ gives
$$f(-1) = -1$$
$$f'(-1) = +1$$
$$f''(-1) = 0$$
and all higher derivatives of $f$ vanish at this point.
Hence $$f(x) = x$$
Derivation in detail
The main non trivial point is to show that for $n\ge2$ the n-th derivate of $f$ at $x=-1$ vanishes.
Let
$$h(x,y) = f(g(x,y))- r(x,y)$$
where
$$g(x,y) = x f(y) + f(x)+y$$
$$r(x,y) = x y + f(x) + f(y)$$
then the defining equation becomes
$$h(x,y) = 0$$
For the partial derivatives with respect to $x$ we have
$$\frac{\partial^n}{\partial x^n}h(x,y) = 0$$
Let us write down the derivatives of $g$ and $r$
$$\frac{\partial}{\partial x}g(x,y) = f(y)+f'(x)$$
$$\frac{\partial^n}{\partial x^n}g(x,y) = f''(x), \text{for }n\ge2$$
$$\frac{\partial}{\partial x}r(x,y) = y+f'(x)$$
$$\frac{\partial^n}{\partial x^n}r(x,y) = f''(x), \text{for }n\ge2$$
Now the derivatives of $h$
$$\frac{\partial}{\partial x}h(g(x,y)) = f'(g)g_x-r_x$$
$$\frac{\partial^2}{\partial x^2}h(g(x,y)) = f''(g)(g_x)^2+ f'(g)g_{xx}-r_{xx}\\=
f''(g)(g_x)^2+ f'(g)f''(x)-f''(x)$$
$$\frac{\partial^3}{\partial x^3}h(g(x,y)) = 
f^{(3)}(g)(g_x)^3+ 2f''(g)g_x g_{xx} + f''(g)g_x f''(x) + f'(g) f^{(3)}(x)-f^{(3)}(x)\\
=\Theta(g_x) +f^{(3)}(x)(f'(g)-1)\\$$
here $\Theta(g_x)$ comprises terms containing at least one factor $g_x\text{, }g_{xx} \text{ and higher}$
$$\frac{\partial^4}{\partial x^4}h(g(x,y)) = 
\Theta(g_x)+\frac{\partial}{\partial x} f^{(3)}(x)(f'(g)-1)\\=
\Theta(g_x) + f^{(4)}(x)(f'(g)-1)+f^{(3)}(x)(f'(g)g_x)\\=
\Theta(g_x) + f^{(4)}(x)(f'(g)-1)
$$
and generally for $n\ge2$
$$\frac{\partial^n}{\partial x^n}h(g(x,y)) = 
\Theta(g_x) + f^{(n)}(x)(f'(g)-1)
$$
Now, as suggested first in a comment, we look at the point $P=(x=-1, y=-1)$
We have
$$g(p) = x f(y) + f(x)+y\to - f(-1) + f(-1)-1 = -1$$
$$r(p) = x y + f(x) + f(y)\to 1 + 2f(-1)$$
Now for $h$ we have at $P$
$$h(P)=f(g)-r \to f(-1) - 1 - 2f(-1)= - 1- f(-1) $$
Letting $h=0$ leads to
$$f(-1) = -1$$
Next, for $h_x$ we get
$$h_x(P)=f'(g)g_x-r_x=f'(g)(f(y)+f'(x))-(y+f'(x))=f'(-1)(f(-1)+f'(-1))-(-1+f'(-1))=f'(-1)(-1+f'(-1))-(-1+f'(-1))=-f'(-1)^2-2f'(-1)+1=(f'(-1)-1)^2$$
Letting $h_x=0$ gives
$$f'(-1)=+1$$
Next
$$h_{xx}=f''(g)(g_x)^2+ f''(x)(f'(g)-1)$$
observing that $g_x(P)=0$ and $f'(P)=-1$ we have
$$0=h_{xx}|_P=-2 f''(-1)$$
And all derivatives with $n\ge2$ vanish at $P$ as seen from
$$0=h_{x^{n}}= \Theta(g_x) + f^{(n)}(x)(f'(g)-1)|_P=-2 f^{(n)}(-1)
$$
A: Introduction
I'm not really an expert at these functional equations, but it seems that I'm on a bit of a streak (two of my previous three answers come from this tag), and as everybody knows: being in flow is a wonderful thing.
The first thing I did was check Approach0 and the likes for similar questions. I didn't find the exact question anywhere, but I must admit that I found some really similar questions. However, similarity can easily count for nothing, especially with olympiad problems. That's why I'm not going to bother citing those similar problems.
I will definitely try and explain my thought process every step of the way so that others can also understand how to attack similar problems. I will also, briefly, step down a blind alley to illustrate a useful-looking method that I believe will not work for this scenario.

First thoughts and notation
When I hit this problem for the first time, I made some observations that I have here on record.
Let $P(x,y)$ be the assertion $$f(xf(y)+f(x)+y) = xy+f(x)+f(y)$$ for $x,y \in \mathbb R$. I observed that

*

*The right-hand side is symmetric in $x,y$.


*The left hand side is $f(\cdot)$ for some $\cdot$. Therefore, every RHS produced by $P(x,y)$ must be in the range of $f$.

On $f(0)$ and the range of $f$
Then I began with the obvious substitutions
$$
P(x,0) : f(xf(0)+f(x)) = f(x)+f(0)
$$
$$
P(0,y) : f(f(0)+y) = f(0)+f(y)
$$
$$
P(x,-f(x)) : f(xf(-f(x))) = -xf(x) + f(x) - f(-f(x))
$$
The third one is useful : using $x=0$, $P(0,-f(0))$ gives $f(-f(0)) = 0$ following cancellations.

On the uniqueness of $f(y)=0$, and the range of $f$
As the OP mentioned, $0$ has the unique preimage $0$.
Indeed, if $f(y)=0$ then $P(y,y)$ gives $y^2=0$ so that $y=0$.
We derived earlier that $f(-f(0)) = 0$. Therefore, it follows that $-f(0)=0$, hence $f(0)=0$ and $0$ is the unique such value by the previous paragraph.
Once we get this, look back at $P(x,0)$ and observe that $f(0)=0$ to see that $f(f(x))=f(x)$ for all $x$.

The range of $f$, and a blind alley
We saw that $f(f(x)) = f(x)$ is true for all $x$. We can be tempted to show that $f$ is surjective now. In this direction, I'll explain what I did, because the commenters above, who very eagerly attempted to solve this question, will be happy to know that I followed them and treated them like my gurus.
Let $x,y$ be arbitrary. Start with $P(f(x),f(y))$ and use $f(f(x))$ to get $$
f(f(x)f(y)+f(x)+f(y)) = f(x)f(y)+f(x)+f(y)
$$
That is, we've proven that under the function $g(x,y)=xy+x+y$, the range of $f$ is closed. This can be used to derive various corollaries from the comments.
Another observation can be made by plugging in $P(-1,-1)$ which gives $f(-1)=-1$. (THIS will be truly crucial)
One can also try to prove the following : $f(f(x)^2)= f(x)^2$ for all $x$ , and $f(-f(x)) = -f(x)$ for all $x$. Thus, we've shown that the range is closed under various operations.
However, this is actually a blind alley. While I cannot assert it with utmost confidence, observe that

*

*One can't attempt substitutions with rational numbers in general: try $f(\frac 12)$, for example, and see where you get. Or maybe try setting $y = \frac 1x$, and you'll have a problem. It seems that you can't get out of the integers.


*Even if you get to the rational numbers, there's no way to place any monotonicity conditions on $f$, or continuity conditions, which allow you to move from the rational numbers to the real numbers.
Therefore, attempts to make this work are likely to be futile, though I absolutely, absolutely invite everyone to try.
EDIT: It seems that Sil, in the comments, has found a way out using a clever substitution! I'll have a go and change everything I've said above if it works out.

Hitting a right approach
A right approach , in this case, comes from observing that if we expect, in light of everything we've said before, that $f(x)=x$ for $x \in \mathbb R$, then the function $g(x) = f(x)-x$ should be identically zero for $x \in \mathbb R$.
However, when you form the functional equation for $g(x)$, something that one may not observe in the original equation manifests itself, and we see niceties occur.
Let's do it : let $g(x) = f(x)-x$ for $x \in \mathbb R$. Start with $P(x,y)$ and begin to substitute $g$ for $f$ everywhere, knowing that $f(z) = g(z)+z$.
\begin{align}
f(xf(y)+f(x)+y) &= xy+f(x)+f(y) \\
\implies g(xf(y)+f(x)+y) + xf(y)+f(x)+y &= xy+g(x)+x+g(y)+y \\
\implies g(xg(y)+xy+g(x)+x+y)+xg(y)+xy+g(x)+x+y &= xy+g(x)+x+g(y)+y \\
\implies g(xg(y)+xy+g(x)+x+y)+xg(y) &= g(y) \\
\implies g(xg(y)+xy+g(x)+x+y) &= (1-x)g(y)
\end{align}
We make a crucial observation about $g$ now. Let's call the last identity in the chain above $P^*(x,y)$.

Claim : if there is a $y_0$ such that $g(y_0) \neq 0$, then $g$ is surjective.

Proof : If $g(y_0) \neq 0$, then $g(y_0)$ belongs in the range of $g$ obviously. Let $T \in \mathbb R$ be arbitrary and let $x_0 = 1-\frac{T}{g(y_0)}$. Then, $(1-x_0)g(y_0)=T$. By observing $P^*(x_0,y_0)$, one sees that $T$ is in the range of $g$, as desired.
We will show that $g$ cannot be surjective in the next section, creating the contradiction.

An ideal choice, and a contradiction to $g$ being surjective
From surjectivity, we will now pick the most suitable candidate for the next substitution, leading to a vast simplification.
The key idea is the following : look at the left hand side of $P^*$, and you have $g(xg(y)+g(x)+xy+x+y)$. We want a value of $x$ or $y$ that can perhaps cancel a great number of terms. One such substitution is easily observed: $x=-1$. Recall that $f(-1)=-1$ (Yes, this was important!) so $g(-1) = 0$.
Now $P^*(-1,y)$ gives $$
g(-g(y)+0-y-1+y) = 2g(y) \implies g(-g(y)-1) = 2g(y)
$$
However, we mentioned that $g$ was surjective! Therefore, the above equation actually holds for all $z$ by merely substituting $y = g^{-1}(z)$ above where $y$ is any preimage of $z$.
Thus, we obtain $$
g(-z-1) = 2z
$$
for all $z$. This is a contradiction : take $z=-1$, then $g(0)=-2$ is obtained, which is a contradiction because $f(0)=0$ therefore $g(0)=0$.
Finally, we obtain that $g \equiv 0$, and $f(x) =x$ everywhere.
