Solution (source code)

= Solution

For real $x$ at which the denominator is nonzero, mapping to the <unit circle> gives $|ax+b|^2=|cx+d|^2$. There are infinitely many such real $x$, so the resulting quadratic <polynomial> identity has all coefficients zero:
$$
|a|^2=|c|^2,\qquad
\operatorname{Re}(a\overline b)=\operatorname{Re}(c\overline d),\qquad
|b|^2=|d|^2.
$$
If $c=0$, the first equality gives $a=0$, contradicting the nonzero <determinant> $ad-bc$. Thus $a,c\ne0$. Put $\lambda=a/c$, $k=b/a$ and $\ell=d/c$. The equalities become
$$
|\lambda|=1,\qquad\operatorname{Re}k=\operatorname{Re}\ell,\qquad |k|=|\ell|.
$$
Two <complex numbers> with equal real parts and equal moduli have imaginary parts equal in magnitude. Consequently $\ell=k$ or $\ell=\overline k$. The first option would give $ad-bc=ac(\ell-k)=0$, so it is excluded. The second gives
$$
\boxed{w=\lambda\frac{z+k}{z+\overline k},\qquad |\lambda|=1.}
$$
Moreover $k$ must be nonreal, since otherwise the two options coincide and the map would be constant. This proves the <Möbius transformations from the real axis to the unit circle> representation directly from coefficients.

Conversely, for real $x$, the numbers $x+k$ and $x+\overline k$ are <complex conjugates>, and their moduli agree. Nonreal $k$ ensures the denominator is nonzero on the <real line>; multiplication by $\lambda$ preserves modulus. At the point at infinity on the <Riemann sphere>, the map takes the value $\lambda$. Thus the representation also respects the extended real axis.