Square root of an intrinsic rotation in 2D inner product space

Question

I have a "rotation" in a 2D inner product space, i.e. a matrix satisfying $$R^TMR = M,\quad \det R>0$$ for positive-definite inner product $M$. I want to compute the square root "half-rotation," i.e. a matrix $S$ with $$S^T M S = M, \quad \det S>0, \quad S^2 = R.$$

What's the best way to do this? My current approach is:

1. Compute $A = M^{1/2} R M^{-1/2}$, where $M^{1/2}$ is the unique positive-definite square root of $M$. The matrix $A$ is now orthogonal.

2. Compute the angle of rotation of $A$, $\theta = \arctan(A_{21}/A_{11})$.

3. Build a new rotation matrix $B = \begin{bmatrix}\cos(\theta/2) &-\sin(\theta/2)\\\sin(\theta/2) & \cos(\theta/2)\end{bmatrix}.$
4. Finally, we have $S = M^{-1/2}BM^{1/2}.$

This seems to work in practice but surely there's a more elegant approach (based on some kind of direct decomposition of the original $R$ matrix)?