How can I find all functions $ f : \mathbb R \setminus \{ 0 \} \to \mathbb R \setminus \{ 0 \} $ such that the following two functional equations are fulfilled:
$$ f ( x y ) = f ( x ) f ( y ) \quad \text {and} \quad f \big( f ( x ) \big) = f ( x ) $$
There are at least four such functions:
Also, if some $ f ( x ) $ is a solution then $ g ( x ) = \frac x { f ( x ) } $ is a solution as well.
It can be shown that the only regular functions $ f : \mathbb R \setminus \{ 0 \} \to \mathbb R \setminus \{ 0 \} $ satisfying $$ f ( x y ) = f ( x ) f ( y ) \tag 0 \label 0 $$ for all $ x , y \in \mathbb R \setminus \{ 0 \} $ and $$ f \big( f ( x ) \big) = f ( x ) \tag 1 \label 1 $$ for all $ x \in \mathbb R \setminus \{ 0 \} $, are exactly those listed by yourself (and there are also wild solutions). To see this, it's easier to first characterize all functions $ f : ( 0 , + \infty ) \to ( 0 , + \infty ) $ satisfying \eqref{0} and \eqref{1}, which I do below.
The functions $ f : ( 0 , + \infty ) \to ( 0 , + \infty ) $ satisfying \eqref{0} and \eqref{1} for all $ x , y \in ( 0 , + \infty ) $, are exactly those of the form $ f ( x ) = \exp \big( A ( \log x ) \big) $ for some additive and compositionally idempotent $ A $; i.e. those $ A : \mathbb R \to \mathbb R $ that satisfy $$ A ( x + y ) = A ( x ) + A ( y ) \tag 2 \label 2 $$ for all $ x , y \in \mathbb R $ and $$ A \big( A ( x ) \big) = A ( x ) \tag 3 \label 3 $$ for all $ x \in \mathbb R $. It's straightforward to check that a function $ f $ of this form will satisfy \eqref{0} and \eqref{1}. To see that these are the only ones, note that for any $ f $ which is a solution, you can define $ A : \mathbb R \to \mathbb R $ with $ A ( x ) = \log \big( f ( \exp x ) \big) $, and then \eqref{0} will imply \eqref{2}, and \eqref{1} will imply \eqref{3}.
It seems that we've only rephrased the problem. But this new formulation has the additional advantage that we can use linear algebra to characterize all the solutions. Here, we look at $ \mathbb R $ as a vector space over the field $ \mathbb Q $.
First, note that assuming some regularity conditions on $ f $, it can be shown that the only solutions are $ f ( x ) = 1 $ and $ f ( x ) = x $. Any of these regularity conditions will work: differentiability at some point, continuity at some point, boundedness from either above or below on some interval, or Lebesgue measurability. These conditions will also hold for $ A $, and as $ A $ satisfies \eqref{2}, there must be a constant $ a \in \mathbb R $ such that $ A ( x ) = a x $ for all $ x \in \mathbb R $. \eqref{3} shows that $ a \in \{ 0 , 1 \} $, and that implies what was desired. See this post for more information about functions satisfying Cauchy's functional equation.
Without regularity conditions, there are other solutions as well. Characterizing all functions $ A $ satisfying \eqref{2} and \eqref{3} can be done using Hamel bases of $ \mathbb R $, which can be shown to exist using the axiom of choice. Choose a Hamel basis $ \mathcal H $ and a nonempty subset $ \mathcal I $ of $ \mathcal H $, map every member of $ \mathcal I $ to itself and every member of $ \mathcal H \setminus \mathcal I $ to an arbitrary member of the linear span of $ \mathcal I $, and finally extend this mapping linearly to a function $ A : \mathbb R \to \mathbb R $. It's then straightforward to check that \eqref{2} and \eqref{3} hold. Conversely, every $ A $ with those properties comes from a Hamel basis following the method we described. To see this, first note that using \eqref{2}, we can see that range of $ A $ is linear subspace of $ \mathbb R $. By \eqref{3}, $ A $ maps every point of this particular subspace to itself. Now, choose a basis $ \mathcal I $ for this subspace, and extend it to a basis $ \mathcal H $ for the whole $ \mathbb R $. Then range of $ A $ is exactly the linear span of $ \mathcal I $, and everything, including members of $ \mathcal H \setminus \mathcal I $, must be mapped into it.
Suppose $ f : \mathbb R \setminus \{ 0 \} \to \mathbb R \setminus \{ 0 \} $ satisfies \eqref{0} and \eqref{1}. For any $ x \in ( 0 , + \infty ) $, we can substitute $ \sqrt x $ for both $ x $ and $ y $ in \eqref{0} to see that $ f ( x ) = f \left( \sqrt x \right) ^ 2 > 0 $. Thus if $ g $ is the restriction of $ f $ to $ ( 0 , + \infty ) $, we can see that $ g : ( 0 , + \infty ) \to ( 0 , + \infty ) $, and it satisfies equations similar to \eqref{0} and \eqref{1}, which means that it is one of the functions discussed above. Now, letting $ x = y = - 1 $ in \eqref{0} we get $ f ( - 1 ) ^ 2 = f ( 1 ) = 1 $, or equivalently $ f ( - 1 ) \in \{ - 1 , 1 \} $. Putting $ y = - 1 $ in \eqref{0}, we see that $ f ( - x ) = f ( - 1 ) f ( x ) $, which means that either $ f ( x ) = g ( | x | ) $ for all $ x \in \mathbb R \setminus \{ 0 \} $, or $ f ( x ) = ( \operatorname {sgn} x ) g ( | x | ) $ for all $ x \in \mathbb R \setminus \{ 0 \} $. Finally, it's straightforward to check that in any of these cases, with $ g $ being any of the functions discussed before, $ f $ will satisfy \eqref{0} and \eqref{1}, and therefore we've found all the solutions. Assuming regularity conditions on $ f $, $ g $ will be regular, too, and thus it must either be the constant one function or the identity function, which proves that $ f $ must be one of the functions listed by you.