Also, there are only finitely many solutions to that functional equation I believe.
It is clear that if a solution vanishes somewhere it vanishes everywhere (let y be a root to make the RHS vanish and then vary x so the LHS vanishing tells you that f is identically zero).
Let us assume now that f is a nonvanishing solution.
Applying the functional equation tells us:
f(x+y+z)=f(x+(y+z))=f(x)f(y+z)f(xy+xz)=f(x)f(y)f(z)f(yz)f(xy)f(xz)f(x^2yz).
On the other hand, if we group the terms differently, we obtain f(x+y+z)=f(y+(x+z))=f(x)f(y)f(z)f(yz)f(xy)f(xz)f(xy^2z).
Setting z=1 and dividing out the factors (which are nonzero!) tells us f(x^2y)=f(xy^2) for all real x and y.
But for nonzero a,b we can always solve x^2y=a, xy^2=b. (x^3=a^2/b, y^3=b^2/a).
This implies that f must be constant on R \ {0}. Moreover, by setting x = y =/= 0 in the functional equation, we obtain that this constant must satisfy x^3=x, which has only three solutions.
Similarly, setting x=y=0 tells us that the value at the origin must also solve this same cubic.
So we have narrowed down possible solutions to the functional equation to a finite set of functions and we are done.
(It would not take much additional effort to find ALL solutions to the functional equation if one so desired.)