Here is a heuristic.

Observe that $f$ is decreasing on any interval of the form $I_n:=(1/(n+1),1/n)$. If $x\in I_n$ and $f(x)=a$, then $x=\frac{1}{a+n}$, and so $f^{-1}(a,1)$ is the union of the intervals $(1/(n+1),1/(n+a))$. Supposing there were an $f$-invariant measure $\mu=gdx$, you can see that

$\displaystyle{\sum_{n=1}^\infty \int_{1/n+1}^{1/n+a}g(x)dx=\int_a^1g(x)dx}$.

Supposing that $g$ is continuous, take the derivative of both sides:

(**) $\displaystyle{\sum_{n=1}^\infty g\left(1/(n+a)\right)/(n+a)^2=g(a)}$.

Approximating the sum by an integral, we get

$\displaystyle{g(a)\approx\int_1^\infty g\left(\frac{1}{x+a}\right)\frac{1}{(x+a)^2}dx}$.

Supposing that $g$ is never too big or too small, this gives

$g(a)\approx \frac{C}{1+a}$.

One can then plug this kind of function into (**) to see that such a function works for any C. Then just pick C to normalize.

Probably someone from the era of Gauss (especially with his acuity at mathematics) did not need to do anything past (**) since people back then seem like magicians when it comes to expressions with infinite sums.