Find all functions $f\colon \mathbb{R}^* \to \mathbb{R}^* $ satisfying $f(xyz)=f(xy+yz+xz)\big(f(x)+f(y)+f(z)\big)$ when $xy+yz+xz\ne 0$ 
Find all functions $f\colon \mathbb{R}^* \to \mathbb{R}^* $ from the non-zero reals to the non-zero reals, such that
$$f(xyz)=f(xy+yz+xz)\big(f(x)+f(y)+f(z)\big)$$
for all non-zero reals $x, y, z$ such that $xy+yz+xz\ne 0$.

I think that only two solutions are: $f(x)=\frac{1}{3}$ and $f(x)=\frac{1}{x}$.
I would appreciate any suggestions.
 A: Plugging in $(x,y,z)=(1,1,-1)$ and $(x,y,z)=(1,-1,-1)$, we get
$$\begin{eqnarray*}
  f(-1)&=&f(-1)(2f(1)+f(-1))\\
  f(1)&=&f(-1)(2f(-1)+f(1))
\end{eqnarray*}$$
Solving, $f(1)=1$ or $\frac13$.
Case 1: $f(1)=1$
For any $a\neq0$, plugging in $(x,y,z)=(a,-a,1)$ gives
$$
  f(-a^2)=f(-a^2)(f(a)+f(-a)+1).
$$
Hence $f(-a)=-f(a)$. For any $a,b>0$, putting $(x,y,z)=(\sqrt{a},-\sqrt{a},b)$ gives
$$
  f(-ab)=f(-a)f(b)
$$
Hence $f(ab)=f(a)f(b)$.
Let $g(x)=\frac1{f(x)}$, so $g(ab)=g(a)g(b)$ for $a,b>0$. In particular if $a>0$ then $g(a)=g(\sqrt{a})^2>0$. The original equation becomes
$$
  g(xy+yz+zx)=g(xyz)\left(\frac1{g(x)}+\frac1{g(y)}+\frac1{g(z)}\right)
    =g(xy)+g(yz)+g(zx).
$$
For any $a,b,c>0$, setting $(x,y,z)=(\sqrt{ac/b},\sqrt{ab/c},\sqrt{bc/a})$ gives
$$
  g(a+b+c)=g(a)+g(b)+g(c).
$$
In particular if $a>b>0$ then $g(a)=g(b)+2g(\frac{a-b}2)>g(b)$. Note that $g(3)=3g(1)=3$, so $g(3^k)=3^k$ for all integers $k$. Thus for any $a,b>0$,
$$
  g(a+b)\leq g\left(a+b+3^k\right)=g(a)+g(b)+3^k.
$$
Sending $k\rightarrow-\infty$, $g(a+b)\leq g(a)+g(b)$. On the other hand, pick $k$ such that $3^k<b$. Then
$$
  g(a)+g(b)\leq g(a)+g(b-3^k)+g(3^k)=g(a+b).
$$
Hence $g(a+b)=g(a)+g(b)$. Since $g(1)=1$, $g(x)=x$ on positive rationals. Also $g$ is increasing, so $g(x)=x$ on positive reals. Since $f$ is odd, $f(x)=\frac1x$ for all $x\in\mathbb R^*$.
Case 2: $f(1)=\frac13$
For any $a\neq0$, plugging in $(x,y,z)=(a,-a,1)$ gives
$$
  f(-a^2)=f(-a^2)(f(a)+f(-a)+\frac13).
$$
Hence $f(a)+f(-a)=\frac23$. Putting $(x,y,z)=(a,-a,b^2)$ gives
$$
  f(-a^2b^2)=f(-a^2)(\frac23+f(b^2))=f(-a^2)(\frac43-f(-b^2)).
$$
Thus
$$
  f(-a^2b^2)+f(-a^2)f(-b^2)=\frac43f(-a^2).
$$
Swapping $a$ and $b$, $f(-a^2)=f(-b^2)$. In particular, putting $b=1$, $f(-a^2)=\frac13$. Hence $f(a)=\frac13$ for $a<0$, and for $a>0$ we have $f(a)=\frac23-f(-a)=\frac13$ also.
