Explicitly representing an isometry as a composition of circle inversions Let $\mathcal C$ be a circle $|z - z_0| = a$ in the complex plane.  Then inversion in $\mathcal C$ is the map $z \mapsto z_0 + \frac{a^2}{\overline{z} - \overline{z_0}}$, which is easily seen to be the complex conjugate of a Möbius transformation.  Composing inversions in different circle we find that composing an even number of circle inversions produces a Möbius transformation, and composing an odd number of circle inversions produces the reflection of one.
This Math Overflow answer asserts that all Möbius transformations and their reflections can be generated by circle inversions.  If so, then any isometry of $\mathbb R^2$ (which can be represented as a Möbius transformation if we identify $\mathbb R^2$ with $\mathbb C$ in the usual way) should be expressible as a composition of inversions.  Since the isometries are generated by the reflections, it is enough to show that reflection across any line $\ell$ can be expressed as a composition of some odd number (presumably 3?) of inversions.
How do we do this explicitly?  Since the line of reflection is determined by two points, presumably we should take two of the circles to be centered on $\ell$, but where do we place the center of the third circle?  Do all three (?) circles have the same radius?  How exactly does this work?
I'd like a description of the general case, but I'll settle for an explanation of how the conjugation map $z \mapsto \overline{z}$ can be explicitly represented as a composition of an odd number of circle inversions.  I can probably work out the rest from that.
 A: Let $\ell$ be the axis we want to reflect through, and let $c$ be any circle whose centre is not in $\ell$.
Then the reflection $R_\ell$ through $\ell$ is sent to the inversion $I_{\ell'}$ when the whole plane is transformed by the inversion $I_c$ through $c$, where $\ell'=I_c(\ell)$.
Algebraically it can be expressed as conjugation by $I_c$.
In other words, since $I_c^{-1}=I_c$, we have
$$R_\ell\ =\ I_c\circ I_{\ell'}\circ I_c\,.$$
We can directly see that the right hand side is involutive (its square is the identity), and that exactly the points of $\ell$ are fixed.
Ultimately, you can also verify it by calculation for the case $R_\ell=z\mapsto\bar z$ e.g. with $c:\ |z-2i|=2$.
(Or, it's somewhat simpler if we move it all down so that the centre of $c$ is the origin, i.e. conjugate all by $z\mapsto z-2i$.)
A: When we stereographically project the Riemann sphere to the complex plane, circles through the pole (the 'point at $\infty$') become straight lines. Inversions include both inversions across circles and inversions across lines, which are just reflections. Thus, any reflection across a line is a product of an odd number of inversions, since the reflection is itself an inversion (and $1$ is an odd number).
The fact that all rigid motions are the product of two reflections is a special case of the Cartan-Dieudonne theorem, which says any rotation (special orthogonal matrix) in $n$ dimensions is expressible as a product of an even number (less than $n$) of reflections (across hyperplanes). Indeed, this follows from the spectral theorem, which tells us any rotation is a product of 2D rotations in orthogonal planes.
I'll discuss how to visualize creating Mobius transformations from inversions.

It helps to have an idea of the four following orthogonal coordinate systems:

*

*Rectangular coordinates. Involves a family of lines parallel to the $x$-axis, and an orthogonal family of lines parallel to the $y$-axis.

*Polar coordinates. Involves a family of lines through the origin and an orthogonal family of circles centered at the origin.

*Bipolar coordinates. Involves the family of circles through two given points (the "poles") and an orthogonal family of so-called Apollonian circles.

*Smith chart. Involves a family of circles incident to the origin symmetric across the $x$-axis, orthogonal to the family of circles incident to the origin symmetric across the $y$-axis. (I don't think this system has an official name, so I've decided to use "Smith chart.")

It also helps to understand what these systems look like stereographically projected to the sphere.
First, imagine a plane tangent to the sphere (at the pole we project from). If we rotate this plane around an axis through it, we carve out a family of circles on the sphere (illustrating the sphere is a so-called smash product of lower-dimensional spheres), and if we rotate it around a perpendicular tangent axis we get an orthogonal family of circles. This is the rectangular coordinate system projected. If we project from the sphere to the plane from a different point, this turns into a Smith chart. Equivalently, we get the Smith chart if we apply a circular inversion to rectangular coordinates.
If we project spherical coordinates from one of its poles, we get polar coordinates in the plane. If we project from a different point, we get bipolar coordinates. Thus, polar coordinates can be considered a limiting case of bipolar coordinates with one of the poles reaching infinity.

All Mobius transformations are conjugate to an element of one of three subgroups:

*

*$A$ ("abelian") consists of hyperbolic transformations aka homotheties or dilations. In spherical or bipolar coordinates, think of one pole as a "source" and one pole as a "sink." Represented by real diagonal matrices.

*$N$ ("nilpotent") consists of parabolic transformations. In rectangular coordinates, these are translations. On the sphere, I like to visualize this as rubber bands being pulled around. Represented by unitriangular matrices.

*$K$ ("kompact") consists of elliptic transformations. In spherical and polar coordinates these are just rotations. In bipolar coordinates, restricting to a half-plane these model rotations around a point in hyperbolic geometry. Represented by special unitary matrices (aka $\mathrm{SU}(2)\cong\mathrm{Sp}(1)=S^3$).

The letters $ANK$ are traditional from the Iwasawa decomposition (in Lie theory) of the Mobius group $\mathrm{SL}_2\mathbb{C}$. Hyperbolic and parabolic transformations together model the isometries of hyperbolic geometry (corresponding to $\mathrm{SL}_2\mathbb{R}$ for the half-plane model or $\mathrm{SU}(1,1)$ for the disk model, which are conjugate via the Cayley transform). Parabolic and elliptic transformations together model the isometries of Euclidean geometry (corresponding to triangular matrices in the Mobius group, its so-called Borel subgroup, I think conjugate to $\mathrm{SO}_2\mathbb{C}$).
Parabolic transformations have a fixed point (with double multiplicity), while hyperbolic and elliptic transformations each have two fixed points.

There are three possible configurations of two circles in the plane, each corresponding (by composing their inversions) to the three types of transformations: (a) nonintersecting circles, corresponding to hyperbolic transformations; (b) tangent (or incident) circles, corresponding to parabolic transformations; (c) intersecting circles, corresponding to elliptic transformations.
Two nontangent circles lie in a unique bipolar coordinate system, unless they are concentric circles or intersecting lines in which case they lie in a unique polar coordinate system. Two tangent circles lie in a unique Smith chart system, unless they are parallel lines in which case they live in a unique rectangular coordinate system. These coordinate systems illustrate what the transformations look like when we compose the inversions.
Indeed, for any pair of circles, if we use the appropriate one-parameter subgroup $\phi$, then if $\phi(s)$ is the transformation that sends the first circle to the second, then the composition of their inversions is $\phi(2s)$; there is a "doubling effect." For example, (a) the product of reflections across two parallel lines is a translation in the perpendicular direction by twice the distance between them, (b) the product of reflections across two intersecting lines $\theta$ apart is a rotation (around their intersection) by $2\theta$, and (c) the product of inversions across two intersecting circles acts locally as a rotation by $2\theta$ in neighborhoods of the two poles, where $\theta$ is the acute angle between the circles.
