Interpolating Moebius transformations

In a previous post on stretching cusps in the Mandelbrot set I used Moebius transformations to map between generalized circles (circles plus straight lines), in particular mapping three points on a circle to a straight line through \(0\), \(1\) and \(\infty\). Yesterday I was wondering how to animate the transition, which requires interpolating between Moebius transformations.

Some research online led me to David Speyer's answer on mathoverflow, which suggested using a matrix representation and interpolating that:

\[h(t) = f \exp(t \log (f^{-1} g)) \]

The matrix representation for a Moebius transformation is:

\[ \left( \begin{array}{cc} a & b \\ c & d \end{array} \right) \sim \frac{a z + b}{c z + d} \]

and the trace \(\mathrm{tr}\), determinant \(\det\) and inverse \(.^{-1}\) of a 2x2 matrix are:

\[ \begin{aligned} \mathrm{tr} \left( \begin{array}{cc} a & b \\ c & d \end{array} \right) &= a + d \\ \det \left( \begin{array}{cc} a & b \\ c & d \end{array} \right) &= a d - b c \\ \left( \begin{array}{cc} a & b \\ c & d \end{array} \right)^{-1} &= \frac{1}{a d - b c} \left( \begin{array}{cc} d & -b \\ -c & a \end{array} \right) \end{aligned} \]

Matrix \(\exp\) and \(\log\) can be defined by power series, but it's also possible to compute them using diagonalization: the diagonalization of a matrix \(M\) is a pair of matrices \(D\), \(P\) such that \(D\) is diagonal (all elements not on the diagonal are \(0\)) and \(M = P D P^{-1}\). Then \(\exp\) and \(\log\) on \(M\) simplify to element-wise on the diagonal elements of \(D\):

\[ \begin{aligned} \exp M &= P \exp(D) P^{-1} \\ \log M &= P \log(D) P^{-1} \end{aligned} \]

Diagonalizing a matrix involves computing eigenvalues and eigenvectors, which can be quite an involved process. But for 2x2 matrices there is a simple closed form solution:

\[ \begin{aligned} M &= \left( \begin{array}{cc} a & b \\ c & d \end{array} \right) \\ \lambda_\pm &= \frac{\mathrm{tr}{M}}{2} \pm \sqrt{\frac{(\mathrm{tr}M)^2}{4} - \det{M}} \\ D &= \left( \begin{array}{cc} \lambda_+ & 0 \\ 0 & \lambda_- \end{array} \right) \\ P &= \left\{ \begin{array}{l l} \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \right) & \quad b = c = 0 \\ \left( \begin{array}{cc} b & b \\ \lambda_+ - a & \lambda_- - a \end{array} \right) & \quad b \ne 0 \\ \left( \begin{array}{cc} \lambda_+ - d & \lambda_- - d \\ c & c \end{array} \right) & \quad c \ne 0 \end{array} \right. \end{aligned} \]

I implemented this method implemented this method, but sometimes the transitions go a long way round instead of a more direct route. This is hinted at in the other answers on the mathoverflow page, but I didn't find a solution for complex-valued matrices yet (negating when the eigenvalues are negative is tricky, because how do you define negative for a complex number..).

I uploaded a short video demoing this Moebius transformation interpolation: Mandelbrot Moebius Experiments. You can download the C99 source code that generated the video frames. It depends on my pre-release mandlebrot mandlebrot libraries mandelbrot-numerics (git HEAD at 3d55cfc99cc97decb0ba0d9c2fb271a504b8e504) and mandelbrot-graphics (git HEAD at ac34f197c1fc7ce13da2d3a056f4e6d588a82f1f).