``` *     ┌ (x^2 * y + z^2) / (x^2 + z^2),   x,   x * z * (y - 1) / (x^2 + z^2) ┐
* M = │  x * (y^2 - 1)  / (x^2 + z^2),   y,    z * (y^2 - 1)  / (x^2 + z^2) │
*     └ x * z * (y - 1) / (x^2 + z^2),   z,   (x^2 + z^2 * y) / (x^2 + z^2) ┘
* ```
* * This is stable as long as v (the bone) is not too much aligned with +/-Y * (i.e. x and z components are not too close to 0). * * Since v is normalized, we have `x^2 + y^2 + z^2 = 1`, * hence `x^2 + z^2 = 1 - y^2 = (1 - y)(1 + y)`. * * This allows to simplifies M like this: *
``` *     ┌ 1 - x^2 / (1 + y),   x,     -x * z / (1 + y) ┐
* M = │                -x,   y,                   -z │
*     └  -x * z / (1 + y),   z,    1 - z^2 / (1 + y) ┘
* ```
* * Written this way, we see the case v = +Y is no more a singularity. * The only one * remaining is the bone being aligned with -Y. * * Let's handle * the asymptotic behavior when bone vector is reaching the limit of y = -1. * Each of the four corner elements can vary from -1 to 1, * depending on the axis a chosen for doing the rotation. * And the "rotation" here is in fact established by mirroring XZ plane by that given axis, * then inversing the Y-axis. * For sufficiently small x and z, and with y approaching -1, * all elements but the four corner ones of M will degenerate. * So let's now focus on these corner elements. * * We rewrite M so that it only contains its four corner elements, * and combine the `1 / (1 + y)` factor: *
``` *                    ┌ 1 + y - x^2,        -x * z ┐
* M* = 1 / (1 + y) * │                            │
*                    └      -x * z,   1 + y - z^2 ┘
* ```
* * When y is close to -1, computing 1 / (1 + y) will cause severe numerical instability, * so we ignore it and normalize M instead. * We know `y^2 = 1 - (x^2 + z^2)`, and `y < 0`, hence `y = -sqrt(1 - (x^2 + z^2))`. * * Since x and z are both close to 0, we apply the binomial expansion to the first order: * `y = -sqrt(1 - (x^2 + z^2)) = -1 + (x^2 + z^2) / 2`. Which gives: *
``` *                        ┌  z^2 - x^2,  -2 * x * z ┐
* M* = 1 / (x^2 + z^2) * │                         │
*                        └ -2 * x * z,   x^2 - z^2 ┘
* ```