Solution (source code)

= Solution

Ignore the zero-area boundary sets where a reciprocal or its digit is undefined, and put $a=\lfloor1/x\rfloor$. On this digit branch the <natural extension of the continued-fraction Gauss map> has inverse
$$
x=\frac1{a+x'},\qquad y=\frac1{y'}-a,
\qquad a=\left\lfloor\frac1{y'}\right\rfloor.
$$
The images are the disjoint horizontal strips $1/(a+1)<y'<1/a$, covering the square almost everywhere. Thus the branch calculation proves global invariance without an extra sum over overlapping branches.

Its <Jacobian determinant> has magnitude $J=[x^2(a+y)^2]^{-1}$. Since
$$
1+x'y'=\frac{a+y+1/x-a}{a+y}=\frac{1+xy}{x(a+y)},
$$
the transformed density satisfies
$$
P(x',y')J=\frac1{(\log2)(1+x'y')^2}\frac1{x^2(a+y)^2}=P(x,y).
$$
The <change of variables formula> therefore proves invariance. Normalization follows from
$$
\int_0^1\frac{dy}{(1+xy)^2}=\frac1{1+x},\qquad
\int_0^1\frac{dx}{(\log2)(1+x)}=1.
$$
Projecting onto the first coordinate commutes with the <Gauss continued-fraction map>, so its invariant marginal is the <Gauss measure>:
$$
\boxed{p_G(x)=\frac1{(\log2)(1+x)},\quad 0<x<1.}
$$
One may also check it directly. The inverse branches $h_a(x)=(a+x)^{-1}$ give
$$
\sum_{a\geq1}p_G(h_a(x))|h_a'(x)|
=\frac1{\log2}\sum_{a\geq1}\frac1{(a+x)(a+x+1)}=p_G(x),
$$
where the last sum telescopes. Values assigned at the exceptional boundary points do not affect the absolutely continuous <invariant measure>.