First, we show that $f(ab) = f(a) + f(b)$ for $a, b \in (0, \infty)$ such that $a^2b \geq 4$. In fact, if $a^2b \geq 4$ and $a, b > 0$, then we can find $x, y > 0$ such that $x + y = ab$ and $xy = b$, namely
$$x = \frac{a + \sqrt{a^2 - (4/b)}}{(2/b)} > 0 \quad \text{and} \quad y = \frac{a - \sqrt{a^2 - (4/b)}}{(2/b)} > 0.$$

Thus, for $a, b \in (0, \infty)$ such that $a^2b \ge 4$, we plug $x, y$ as above in the functional equation to get $f(ab) = f(x+y) = f(\frac{x+y}{xy}) + f(xy) = f((ab)/b) + f(b) = f(a) + f(b)$. Now let $x, y > 0$. Let
$$z = \max \left\{ \frac{4}{x^2 y^2}, \frac{4}{x^2 y}, \frac{4}{y^2} \right\} > 0.$$
Then $(xy)^2z \ge 4$, $x^2(yz) \ge 4$ and $y^2z \ge 4$ which imply that $f(xyz) = f(xy) + f(z)$, $f(xyz) = f(x) + f(yz)$ and $f(yz) = f(y) + f(z)$ respectively. Therefore,
$$$$
$$\begin{aligned} f(xy) &= f(xyz) - f(z) \\ &= (f(x) + f(yz)) - f(z) \\ &= f(x) + f(yz) - f(z) \\ &= f(x) + (f(y) + f(z)) - f(z) = f(x) + f(y). \end{aligned}$$