We know that in two dimensions, two Möbius transformations performed one after another can be replaced with a single Möbius transformation that is the composite of those two, which can be easily found if the Möbius transformations are given on projective matrix representation. And in three dimesnions or higher, we know that all conformal transformations are Möbius transformations, which can be given on the form
$$ f(\vec{x}) = \vec{b} + \frac{\alpha \textbf{A}(\vec{x} - \vec{a})}{|\vec{x} - \vec{a}|^\varepsilon} \tag{1}\label{1} $$
where $\vec{a},\vec{b} \in \mathbb R^n$, $\alpha\in\mathbb R$, $A$ is an orthogonal $n$-by-$n$ matrix, and $\varepsilon$ is 0 or 2. And since all Möbius transformations are conformal, two Möbius transformations performed one after another should also correspond to a single Möbius transformation. But how can we find the representation of the composite of two Möbius transformations that is given on the form given by $\eqref{1}$ in dimensions higher than 2?
What I have tried
If $\varepsilon = 0$, the Möbius transformation reduces to an affine transformation, which it is easy to composite with another Möbius transformation, so let's focus on the case where $\varepsilon = 2$, and let's for simplicity say that the composite transformation isn't affine either. Say we have to Möbvius transformations,
$$ f_1(\vec{x}) = \vec{b}_1 + \frac{\textbf{M}_1(\vec{x} - \vec{a}_1)}{|\vec{x} - \vec{a}_1|^2} \qquad\text{and}\qquad f_2(\vec{x}) = \vec{b}_2 + \frac{\textbf{M}_2(\vec{x} - \vec{a}_2)}{|\vec{x} - \vec{a}_2|^2}, $$
where $\textbf{M}_1$ and $\textbf{M}_2$ each is an orthogonal matrix times a real-valued scalar. What I want to find is parameters $\vec{a}$, $\vec{b}$ and $\textbf{M}$ such that
$$ f(\vec{x}) = \vec{b} + \frac{\textbf{M}(\vec{x} - \vec{a})}{|\vec{x} - \vec{a}|^2} = f_1(f_2(\vec{x})) = \vec{b}_1 + \frac{\displaystyle\textbf{M}_1\left(\vec{b}_2 + \frac{\textbf{M}_2(\vec{x} - \vec{a}_2)}{|\vec{x} - \vec{a}_2|^2} - \vec{a}_1\right)} {\displaystyle\left|\vec{b}_2 + \frac{\textbf{M}_2(\vec{x} - \vec{a}_2)}{|\vec{x} - \vec{a}_2|^2} - \vec{a}_1\right|^2}. $$
Maybe I'm abusing notation here, but for the limit when $\vec{x}\to\infty$, we get
$$ \vec{b} = \vec{b}_1 + \frac{\displaystyle\textbf{M}_1\left(\vec{b}_2 - \vec{a}_1\right)} {\displaystyle\left|\vec{b}_2 - \vec{a}_1\right|^2} $$
so we get an expression for $\vec{b}$, which is finite if we assume that $\vec{b}_2 \neq \vec{a}_1$ (I suppose that if $\vec{b}_2$ and $\vec{a}_1$ are equal, we are looking at an affine transformation, since all other transformations should have finite values for $\vec{a}$, $\vec{b}$ and $\textbf{M}$). Also, for the limit where $\vec{x}\to\vec{a}$, we get
$$ \infty = \vec{b}_1 + \frac{\displaystyle\textbf{M}_1\left(\vec{b}_2 + \frac{\textbf{M}_2(\vec{a} - \vec{a}_2)}{|\vec{a} - \vec{a}_2|^2} - \vec{a}_1\right)} {\displaystyle\left|\vec{b}_2 + \frac{\textbf{M}_2(\vec{a} - \vec{a}_2)}{|\vec{a} - \vec{a}_2|^2} - \vec{a}_1\right|^2}, $$
and the only way for the RHS to also be $\infty$ is for $\displaystyle\vec{b}_2 + \frac{\textbf{M}_2(\vec{a} - \vec{a}_2)}{|\vec{a} - \vec{a}_2|^2} - \vec{a}_1$ to be zero. So we have
$$ \frac{\textbf{M}_2(\vec{a} - \vec{a}_2)}{|\vec{a} - \vec{a}_2|^2} = \vec{a}_1 - \vec{b}_2, $$
which can be rewritten as
$$ \textbf{M}_2(\vec{a} - \vec{a}_2) = (\vec{a}_1 - \vec{b}_2)(a^2 - 2\vec{a}\cdot\vec{a}_2 + a_2^2), $$
where $a=|\vec{a}|$ and $a_2=|\vec{a}_2|$. So the only unknown in this equation is $\vec{a}$; however, I'm not sure how isolate $\vec{a}$ so I can get an expression for it.
Edit: I just figured out how to solve this equation, but I haven't solved for $\textbf M$ yet. However, if there is some other representation of the Möbius transform for dimension three or higher that was simpler to work with, that would be great, because this seems to involve a lot of computation. I'm going to use this to transform things in a computer graphics application and speed is imperative.
(The solution to the equation is $ \vec{a} = \vec{a}_2 + \displaystyle \frac{\textbf{M}_2^{-1}(\vec{a}_1 - \vec{b}_2)}{\left|\textbf{M}_2^{-1}(\vec{a}_1 - \vec{b}_2)\right|^2}$ )
There are three or four kinds of transformations used:
Every Mobius transformation has a unique standard form $T_b\circ M\circ C\circ T_{-a}$ or $T_b\circ M\circ T_{-a}$, where we combine the dilation and isometry as $M=A\circ D_\alpha=D_\alpha\circ A$. To simplify a composition of two Mobius transformations, it suffices to know how to "slide" these four basic transformations past each other. For example, dilations commute with isometries and isometries commute with inversion, but with dilation and inversion we have the sliding rule $C\circ D_\alpha=D_{1/\alpha}\circ C $. For translations sliding past isometries and dilations we have the rule $T_v\circ M=M\circ T_{M^{-1}v}$ All of these kinds of rules can be written out explicitly with algebraic expressions, of course.
As you've said, it should be clear how to simplify a composition of two Mobius transformations to this "standard form" unless both of them have an inversion $C$ in their standard forms. Then, the only mysterious case is out to simplify $C\circ T_{-v}\circ C$. In other words,
$$ \frac{\displaystyle\frac{x}{\|x\|^2}-v~~}{\displaystyle \Big\|\frac{x}{\|x\|^2}-v\,\Big\|^2} \quad= \quad ? $$
If we plug in $x=\infty$ we get the outer translation vector is $-v/\|v\|^2$. Then we can add $v/\|v\|^2$ to our expression, and multiply numerator/denominator by $\|x\|^2$ and get:
$$ \frac{\displaystyle x-\|x\|^2v~~}{\displaystyle \Big\|\frac{x}{\|x\|}-\|x\|v\,\Big\|^2} + \frac{v}{\|v\|^2} $$
$$ = \frac{x-\|x\|^2v}{1-2x\cdot v+\|x\|^2\|v\|^2}+\frac{v}{\|v\|^2} $$
$$ = \frac{\displaystyle x-\|x\|^2v + \frac{v}{\|v\|^2}\big(1-2x\cdot v+\|x\|^2\|v\|^2\big)}{1-2x\cdot v+\|x\|^2\|v\|^2} $$
$$ = \frac{\displaystyle x-\|x\|^2v+\frac{v}{\|v\|^2}-2x\cdot v\frac{v}{\|v\|^2}+v\|x\|^2}{1-2x\cdot v+\|x\|^2\|v\|^2} $$
$$ = \frac{\displaystyle\left(I-2\hat{v}\hat{v}^T\right)x+\frac{v}{\|v\|^2}}{1-2x\cdot v+\|x\|^2\|v\|^2}. $$
Note, with $\hat{v}=v/\|v\|$, this uses the Householder reflection $R_v(x)=(I-\hat{v}\hat{v}^T)x$, an isometry which reflects vectors across the hyperplane perpendicular to $v$ (in particular, $R_v(v)=-v$). Notice the denominator is symmetric in $x$ and $v$, so we can replace it with what it once was but with $x,v$ swapped:
$$ \frac{\displaystyle R_v(x)+\frac{v}{\|v\|^2}}{\displaystyle \Big\|\|v\|x-\frac{v}{\|v\|} \Big\|^2} $$
$$ =\frac{1}{\|v\|^2} \frac{\displaystyle R_v\Big(x-\frac{v}{\|v\|^2}\Big)}{\displaystyle\Big\|x-\frac{v}{\|v\|^2} \Big\|^2} $$
Therefore we may finally conclude
$$ \frac{\displaystyle\frac{x}{\|x\|^2}-v~~}{\displaystyle \Big\|\frac{x}{\|x\|^2}-v\,\Big\|^2} = -\frac{v}{\|v\|^2} + \frac{R_v}{\|v\|^2}\frac{\displaystyle x-\frac{v}{\|v\|^2}}{\displaystyle\Big\|x-\frac{v}{\|v\|^2}\Big\|^2}. $$
In other words, we can say
$$ CT_{-v}C=T_{-C(v)}D_{1/\|v\|^2}R_vCT_{-C(v)}. $$