, 2025

Any $f$ such that $f(x) \in \{0, 1\}$ for all $x \in \mathbb{R}^{+}$ works. Furthermore, any $f$ such that
$$f(x) = \begin{cases} 0 & \text{or } 1 \\ c & x = x_0 \\ 0 & x \in (x_0, \infty) \end{cases}$$
works as well, where $x_0 > 0$, $c \ge 0$ are arbitrary constants. We now show that these are the only solutions. For $f(x) \ne 0$ easy both-ways induction yields that for all $n \in \mathbb{Z}$ it is true that
$$f(x)^n f(x + y) = f(x + y f(x)^n) \quad (2)$$
Now assume there exist $0 < x_0 < x_1$ such that $f(x_0) \notin \{0, 1\}$ and $f(x_1) \ne 0$ (if such a pair doesn't exist then $f$ must have one of the two forms described above). Then substituting $[x_0, x_1 - x_0]$ into (1) and manipulating $n$ (in particular we consider $n \to -\infty$ if $f(x_0) < 1$, and $n \to +\infty$ if $f(x_0) > 1$) yields that $f$ reaches arbitrarily large values at arbitrarily large arguments. Hence, for every pair of positive reals $c_1, c_2$ there are infinitely many $x$ such that $x > c_1$ and $f(x) > c_2$. Call this fact $(\star)$.
We now multiply the given equation by $f(x + y + z)$, where $z$ is a positive real number, to get
$$f(x + y + z) f(x + y f(x)) = f(x) f(x + y) f(x + y + z) = f(x) f(x + y + z f(x + y)),$$
where we've used the property from the problem statement to obtain the second equality. We now choose $z$ such that $z > y f(x) - y$. Then $x + y + z > x + y f(x)$. Hence, we can apply the problem statement on both the left-most side and the right-most side of the above equation to get
$$\begin{aligned} f(x + y + z) f(x + y f(x)) &= f(x + y f(x) + (z - y f(x) + y) f(x + y f(x))) \\ f(x) f(x + y + z f(x + y)) &= f(x + (y + z f(x + y)) f(x)) \end{aligned}$$
Together with $f(x + y f(x)) = f(x) f(x + y)$, since the LHS's are equal in the above two equations, we get

$$f(x + y f(x) + (z - y f(x) + y) f(x) f(x + y)) = f(x + (y + z f(x + y)) f(x)). \quad (3)$$
If the arguments in the above equation were equal, then by simplification, this would yield the equivalent equality
$$(-y f(x) + y) f(x) f(x + y) = 0. \quad (4)$$
We now choose $x_0, y_0$ such that $f(x_0) \notin \{0, 1\}$ and $f(x_0 + y_0) \neq 0$ and substitute $[x_0, y_0]$ into (2). Note that for this pair, equation (3) does not hold, and hence the arguments in (2) are always distinct. In particular, the arguments on both sides of (2) are linear functions in $z$ with the same positive gradient (namely $f(x_0) f(x_0 + y_0)$), but different $y$-intercept values. Since (2) holds for all large $z$ (namely all $z > y_0 f(x_0) - y_0$), it follows that $f$ is eventually periodic. Hence, there are constants $C, P > 0$ (dependent on $x_0, y_0$), such that $f(x) = f(x + P)$ for all $x > C$.
By $(\star)$ we know that there is an $x_2 > C$ such that $f(x_2) \notin \{0, 1\}$. Then by comparing $[x_2, y]$ with $[x_2, y + P]$ in the original equation we get
$$f(x_2 + y f(x_2)) = f(x_2 + y f(x_2) + P f(x_2)),$$
since the RHS's remains the same (since $x_2 + y > C$). Now let $y = \frac{P}{f(x_2)}$ in the above, to obtain
$$f(x_2 + P) = f(x_2 + P + P f(x_2))$$
and hence $f(x_2) = f(x_2 + P f(x_2)) = f(x_2) f(x_2 + P) = f(x_2)^2$, clear contradiction.