core design

This page explains the mathematics that Luna-Flow/quaternion implements and why its API has the shape it has. Every formula below is the one the code in src/quaternion.mbt evaluates; where the implementation deviates from the mathematics, the deviation is stated.

Design goal

The package gives MoonBit one quaternion type, Quaternion[T], that serves two audiences:

• algebraic code, which wants $\mathbb H$ as a ring over an exact or approximate scalar type and uses it through the luna-generic traits;
• geometric code, which wants unit quaternions as 3D rotations: build them from an axis and an angle or from Euler angles, compose, apply, interpolate and decompose them.

Both are served by the same generic type. Each operation asks only for the traits of T it needs, so ring arithmetic works over Int exactly, while the operations that need square roots or trigonometry go through Double.

Mathematical background

Quaternions as a real algebra

The quaternions $\mathbb H$ are the 4-dimensional real vector space with basis $1, i, j, k$,

$$q = w + x\,i + y\,j + z\,k, \qquad w, x, y, z \in \mathbb R,$$

with a bilinear multiplication fixed by Hamilton’s relations¹1 W. R. Hamilton, “On quaternions; or on a new system of imaginaries in algebra”, 1843. The relations are the ones he carved into Brougham Bridge in Dublin.

$$i^2 = j^2 = k^2 = ijk = -1 .$$

We write $q = (w, \mathbf u)$ with the scalar part $w$ and the vector part $\mathbf u = (x, y, z)$, and identify $\mathbb R^3$ with the pure quaternions $\{(0, \mathbf u)\}$. Quaternion[T] stores exactly this pair: a field r for $w$ and a triple vec for $\mathbf u$. Real numbers $(s, \mathbf 0)$ commute with everything, because the multiplication is $\mathbb R$-bilinear.

Products of the units

Everything about the product follows from the four relations. Multiply $ijk = -1$ on the right by $k$ and use $k^2 = -1$:

$$\begin{aligned} ijk\,k &= -k \\ -ij &= -k \\ ij &= k . \end{aligned}$$

Multiply $ijk = -1$ on the left by $i$ and use $i^2 = -1$: $-jk = -i$, so $jk = i$. The reversed products follow from these two:

$$\begin{aligned} ji &= j\,(jk) = j^2 k = -k, \\ kj &= (ij)\,j = i\,j^2 = -i, \\ ik &= i\,(ij) = i^2 j = -j, \\ ki &= (ij)\,i = i\,(ji) = -ik = j . \end{aligned}$$

So the units multiply like the cross product of the standard basis, with an extra $-1$ on the diagonal:

$\cdot$	$i$	$j$	$k$
$i$	$-1$	$k$	$-j$
$j$	$-k$	$-1$	$i$
$k$	$j$	$-i$	$-1$

The eight elements $\pm 1, \pm i, \pm j, \pm k$ form the quaternion group $Q_8$, which can be realized by $2\times 2$ complex matrices.²2 For example $1 \mapsto I$, $i \mapsto \begin{pmatrix} i & 0 \\ 0 & -i\end{pmatrix}$, $j \mapsto \begin{pmatrix} 0 & 1 \\ -1 & 0\end{pmatrix}$, $k \mapsto \begin{pmatrix} 0 & i \\ i & 0\end{pmatrix}$. These matrices satisfy Hamilton’s relations, and matrix multiplication is associative. The multiplication of $Q_8$ is associative, and by bilinearity so is the product of $\mathbb H$.

The Hamilton product

Expanding $p = w_1 + x_1 i + y_1 j + z_1 k$ times $q = w_2 + x_2 i + y_2 j + z_2 k$ term by term with the table gives

$$\begin{aligned} pq ={} & (w_1 w_2 - x_1 x_2 - y_1 y_2 - z_1 z_2) \\ &+ (w_1 x_2 + x_1 w_2 + y_1 z_2 - z_1 y_2)\, i \\ &+ (w_1 y_2 - x_1 z_2 + y_1 w_2 + z_1 x_2)\, j \\ &+ (w_1 z_2 + x_1 y_2 - y_1 x_2 + z_1 w_2)\, k . \end{aligned}$$

The same product is clearer in scalar–vector form. For two pure quaternions $\mathbf a$ and $\mathbf b$, the squares of the units contribute $-a_1 b_1 - a_2 b_2 - a_3 b_3$, and the mixed terms pair up, for example $a_1 b_2\, ij + a_2 b_1\, ji = (a_1 b_2 - a_2 b_1)k$. Hence

$$\mathbf a\,\mathbf b = -\,\mathbf a\cdot\mathbf b + \mathbf a\times\mathbf b .$$

With bilinearity and the fact that scalars commute,

$$\begin{aligned} (s_1, \mathbf v_1)(s_2, \mathbf v_2) &= s_1 s_2 + s_1\mathbf v_2 + s_2\mathbf v_1 + \mathbf v_1\mathbf v_2 \\ &= \big(s_1 s_2 - \mathbf v_1\cdot\mathbf v_2,\ \ s_1\mathbf v_2 + s_2\mathbf v_1 + \mathbf v_1\times\mathbf v_2\big). \end{aligned}$$

This is literally the Mul implementation: the scalar part is self.r * other.r - dot(self.vec, other.vec) and the vector part is cross(self.vec, other.vec) plus the two scaled vectors.

Swapping the factors only flips the sign of the cross product, so

$$pq - qp = 2\,\mathbf v_1 \times \mathbf v_2 ,$$

and two quaternions commute exactly when their vector parts are parallel. In particular $\mathbb H$ is not commutative: $ij = k$ but $ji = -k$.

The derivation used that the components commute with each other ($s_1\mathbf v_2 = \mathbf v_2 s_1$ and so on). The generic implementation therefore assumes that T has a commutative multiplication, which holds for every numeric type luna-generic provides.

Conjugate and norm

The conjugate is $\bar q = (w, -\mathbf u)$. Multiplying a quaternion by its conjugate leaves a real number, because the cross product of parallel vectors vanishes:

$$q\bar q = \big(w^2 + \mathbf u\cdot\mathbf u,\ -w\mathbf u + w\mathbf u - \mathbf u\times\mathbf u\big) = \big(w^2 + x^2 + y^2 + z^2,\ \mathbf 0\big) = |q|^2 ,$$

and in the same way $\bar q q = |q|^2$. Quaternion::square_len computes this $|q|^2$, and Quaternion::dot is the inner product of $\mathbb R^4$ whose norm it is.

Conjugation reverses products. With $p = (s_1, \mathbf v_1)$ and $q = (s_2, \mathbf v_2)$,

$$\begin{aligned} \overline{pq} &= \big(s_1 s_2 - \mathbf v_1\cdot\mathbf v_2,\ -s_1\mathbf v_2 - s_2\mathbf v_1 - \mathbf v_1\times\mathbf v_2\big), \\ \bar q\,\bar p &= \big(s_2 s_1 - \mathbf v_2\cdot\mathbf v_1,\ -s_2\mathbf v_1 - s_1\mathbf v_2 + \mathbf v_2\times\mathbf v_1\big), \end{aligned}$$

and the two agree because $\mathbf v_2\times\mathbf v_1 = -\mathbf v_1\times\mathbf v_2$. So $\overline{pq} = \bar q\,\bar p$.

The norm is multiplicative. Using associativity, $\overline{pq} = \bar q\bar p$, and the fact that the real number $|q|^2$ commutes with $p$:

$$\begin{aligned} |pq|^2 &= (pq)\,\overline{(pq)} = p\,(q\bar q)\,\bar p \\ &= p\,|q|^2\,\bar p = |q|^2\,p\bar p = |p|^2\,|q|^2 . \end{aligned}$$

Over the integers this is Euler’s four-square identity, which the tutorial checks with Quaternion[Int].

Inverse and the two divisions

If $q$ is not zero, then $|q|^2 > 0$, and $q\bar q = \bar q q = |q|^2$ shows that

$$q^{-1} = \frac{\bar q}{|q|^2}$$

is a two-sided inverse. This is the Inverse implementation: the conjugate scaled by one() / square_len(). Every non-zero quaternion is invertible, so $\mathbb H$ is a division ring (a skew field). It is not a field, because it is not commutative.

Without commutativity, ”$q$ divided by $r$” has two meanings, one for each side on which the unknown sits:

$$\begin{aligned} x\,r = q &\iff x = q\,r^{-1} && \text{(right division, } q / r\text{)}, \\ r\,x = q &\iff x = r^{-1} q && \text{(left division, } q\text{.left\_div}(r)\text{)}. \end{aligned}$$

The implementation avoids forming $r^{-1}$ and divides once at the end. With $q = (s_q, \mathbf v_q)$, $r = (s_r, \mathbf v_r)$ and $\bar r = (s_r, -\mathbf v_r)$, the scalar–vector product gives

$$\begin{aligned} q\,r^{-1} &= \frac{q\,\bar r}{|r|^2} = \frac{\big(s_q s_r + \mathbf v_q\cdot\mathbf v_r,\ \ s_r\mathbf v_q - s_q\mathbf v_r - \mathbf v_q\times\mathbf v_r\big)}{|r|^2}, \\ r^{-1} q &= \frac{\bar r\,q}{|r|^2} = \frac{\big(s_q s_r + \mathbf v_q\cdot\mathbf v_r,\ \ s_r\mathbf v_q - s_q\mathbf v_r - \mathbf v_r\times\mathbf v_q\big)}{|r|^2}, \end{aligned}$$

which are the Div implementation (with c = cross(q.vec, r.vec)) and Quaternion::left_div (with c = cross(r.vec, self.vec)). The quotients have the same scalar part and differ by

$$q\,r^{-1} - r^{-1} q = -\frac{2\,\mathbf v_q\times\mathbf v_r}{|r|^2}.$$

Example. Take $q = 1 + 2i + 3j + 4k$ and $r = 5 + 6i + 7j + 8k$. Then $|r|^2 = 174$, $s_q s_r + \mathbf v_q\cdot\mathbf v_r = 5 + 12 + 21 + 32 = 70$, $s_r\mathbf v_q - s_q\mathbf v_r = (4, 8, 12)$ and $\mathbf v_q\times\mathbf v_r = (-4, 8, -4)$, so

$$\begin{aligned} q / r &= \tfrac{1}{174}\big(70 + 8i + 0j + 16k\big) \approx 0.4023 + 0.0460\,i + 0.0920\,k, \\ q.\mathrm{left\_div}(r) &= \tfrac{1}{174}\big(70 + 0i + 16j + 8k\big) \approx 0.4023 + 0.0920\,j + 0.0460\,k . \end{aligned}$$

The API page checks these values and both defining equations.

Unit quaternions and rotations

A unit quaternion ($|q| = 1$) can be written

$$q = \Big(\cos\tfrac{\theta}{2},\ \sin\tfrac{\theta}{2}\,\mathbf n\Big), \qquad |\mathbf n| = 1 ,$$

which is exactly what from_axis_angle(axis, angle) returns, with $\mathbf n$ the normalized axis. For unit $q$ we have $q^{-1} = \bar q$. Consider the map $\mathbf v \mapsto q\,\mathbf v\,\bar q$ on pure quaternions. Write $q = (w, \mathbf u)$. First,

$$q\,\mathbf v = \big(-\mathbf u\cdot\mathbf v,\ w\mathbf v + \mathbf u\times\mathbf v\big).$$

Multiplying on the right by $\bar q = (w, -\mathbf u)$, the scalar part is

$$-w\,\mathbf u\cdot\mathbf v + (w\mathbf v + \mathbf u\times\mathbf v)\cdot\mathbf u = -w\,\mathbf u\cdot\mathbf v + w\,\mathbf v\cdot\mathbf u + 0 = 0 ,$$

so the image is again a pure quaternion, and the vector part is

$$\begin{aligned} \mathbf v' &= (\mathbf u\cdot\mathbf v)\,\mathbf u + w\,(w\mathbf v + \mathbf u\times\mathbf v) - (w\mathbf v + \mathbf u\times\mathbf v)\times\mathbf u \\ &= (\mathbf u\cdot\mathbf v)\,\mathbf u + w^2\mathbf v + 2w\,\mathbf u\times\mathbf v + \mathbf u\times(\mathbf u\times\mathbf v) \\ &= (w^2 - |\mathbf u|^2)\,\mathbf v + 2(\mathbf u\cdot\mathbf v)\,\mathbf u + 2w\,\mathbf u\times\mathbf v , \end{aligned}$$

using $\mathbf a\times\mathbf u = -\mathbf u\times\mathbf a$ and $\mathbf u\times(\mathbf u\times\mathbf v) = (\mathbf u\cdot\mathbf v)\mathbf u - |\mathbf u|^2\mathbf v$. Substituting $w = \cos\frac\theta2$ and $\mathbf u = \sin\frac\theta2\,\mathbf n$, and the half-angle identities $\cos^2\frac\theta2 - \sin^2\frac\theta2 = \cos\theta$, $2\sin^2\frac\theta2 = 1 - \cos\theta$, $2\sin\frac\theta2\cos\frac\theta2 = \sin\theta$:

$$\mathbf v' = \cos\theta\,\mathbf v + (1 - \cos\theta)(\mathbf n\cdot\mathbf v)\,\mathbf n + \sin\theta\,(\mathbf n\times\mathbf v).$$

This is Rodrigues’ rotation formula: $\mathbf v \mapsto q\mathbf v q^{-1}$ is the rotation by $\theta$ about $\mathbf n$, counter-clockwise by the right-hand rule. Since $|q\mathbf v\bar q| = |q||\mathbf v||\bar q| = |\mathbf v|$, the map is an isometry, as a rotation must be.

Three consequences shape the API:

• Composition is multiplication. $p\,(q\mathbf v q^{-1})\,p^{-1} = (pq)\,\mathbf v\,(pq)^{-1}$, so the rotation “first $q$, then $p$” is p * q.
• Double cover. $(-q)\,\mathbf v\,(-q)^{-1} = q\mathbf v q^{-1}$, so $q$ and $-q$ are the same rotation. Indeed the angle $\theta + 2\pi$ gives $-q$. Every rotation has exactly two unit quaternions, which is why == does not compare rotations and why slerp flips the sign of one input.
• Cheap evaluation. Quaternion::rotate evaluates $\mathbf t = 2\,\mathbf u\times\mathbf v$ and $\mathbf v' = \mathbf v + w\,\mathbf t + \mathbf u\times\mathbf t$. Expanding, $\mathbf u\times\mathbf t = 2(\mathbf u\cdot\mathbf v)\mathbf u - 2|\mathbf u|^2\mathbf v$, so $$\mathbf v + w\,\mathbf t + \mathbf u\times\mathbf t = (1 - 2|\mathbf u|^2)\,\mathbf v + 2(\mathbf u\cdot\mathbf v)\,\mathbf u + 2w\,\mathbf u\times\mathbf v ,$$ which equals the formula above exactly when $1 - 2|\mathbf u|^2 = w^2 - |\mathbf u|^2$, that is when $w^2 + |\mathbf u|^2 = 1$. This is why rotate requires a unit quaternion: it needs no division, but it relies on the norm being 1.

The same formula, applied to the basis vectors, gives the rotation matrix of a unit quaternion:

$$R(q) = \begin{pmatrix} 1 - 2(y^2 + z^2) & 2(xy - wz) & 2(xz + wy) \\ 2(xy + wz) & 1 - 2(x^2 + z^2) & 2(yz - wx) \\ 2(xz - wy) & 2(yz + wx) & 1 - 2(x^2 + y^2) \end{pmatrix}.$$

The package does not return this matrix, but the Euler-angle extraction below reads its entries.

Euler angles

Write $q_X(\alpha)$, $q_Y(\alpha)$, $q_Z(\alpha)$ for the rotations by $\alpha$ about the coordinate axes, and $R_X, R_Y, R_Z$ for their matrices. A sequence of three rotations about axes $A, B, C$ with angles $a, b, c$ can be read in two ways:

• extrinsic (external=true): about the fixed axes, first $A$, then $B$, then $C$. Later rotations multiply on the left: $q = q_C(c)\,q_B(b)\,q_A(a)$.
• intrinsic (external=false): about the axes of the moving body, first $A$, then the moved $B$, then the twice-moved $C$: $q = q_A(a)\,q_B(b)\,q_C(c)$.

The same product read in both ways shows that intrinsic $ABC$ with angles $(a, b, c)$ is extrinsic $CBA$ with angles $(c, b, a)$. The implementation uses exactly this: every to_euler_internal_* function calls the extrinsic function of the reversed order and reverses the triple.

Extraction, extrinsic XYZ. For $R = R_Z(c)\,R_Y(b)\,R_X(a)$, multiplying out the elementary matrices gives

$$R = \begin{pmatrix} \cos b\cos c & \cdots & \cdots \\ \cos b\sin c & \cdots & \cdots \\ -\sin b & \cos b\sin a & \cos b\cos a \end{pmatrix}.$$

Comparing with $R(q)$ entry by entry:

$$\begin{aligned} \sin b &= -R_{31} = 2(wy - xz), \\ a &= \operatorname{atan2}(R_{32}, R_{33}) = \operatorname{atan2}\big(2(wx + yz),\ 1 - 2(x^2 + y^2)\big), \\ c &= \operatorname{atan2}(R_{21}, R_{11}) = \operatorname{atan2}\big(2(wz + xy),\ 1 - 2(y^2 + z^2)\big), \end{aligned}$$

valid while $\cos b > 0$ (dividing both atan2 arguments by $\cos b$ does not change the angle). These are the expressions in to_euler_external_XYZ, with $b = \arcsin(\cdot) \in [-\pi/2, \pi/2]$. The functions to_euler_external_YZX and to_euler_external_ZYX are the same derivation with the axes permuted; a permutation that is not cyclic flips the sign inside the arcsine and the atan2 numerators. The module normalizes q before extracting, so the identities $w^2 + x^2 + y^2 + z^2 = 1$ used in the diagonal of $R(q)$ hold up to rounding.

Gimbal lock. When $\cos b = 0$, the entries $R_{32}, R_{33}, R_{21}, R_{11}$ all vanish and the first and third rotations act about the same physical axis. For $b = \pi/2$ in the XYZ case,

$$R_Z(c)\,R_Y(\tfrac\pi2)\,R_X(a) = R_Y(\tfrac\pi2)\,R_X(a - c),$$

so only $a - c$ is determined (and $a + c$ for $b = -\pi/2$). A conventional fix sets $c = 0$ and recovers $a$ from the entries that remain, here $a = \operatorname{atan2}(R_{12}, R_{22})$. The implementation switches to a special branch when $|\sin b| \ge 0.9998$ (about $88.85^\circ$): it prints a warning, sets $b = \pm\pi/2$ and $c = 0$, but takes the first angle from $\operatorname{atan2}(R_{13}, R_{22})$ for XYZ, and from the regular-branch formula (whose arguments are near zero) for the other orders. Neither isolates $a \mp c$, so the branch currently returns angles that do not reproduce the input; see known deviations.

Construction. from_euler(roll, pitch, yaw) with the default "XYZ" multiplies out $q_X(\text{roll})\,q_Y(\text{pitch})\,q_Z(\text{yaw})$ in closed form. With $c_r = \cos\frac{\text{roll}}{2}$, $s_r = \sin\frac{\text{roll}}{2}$ and so on, two applications of the Hamilton product give

$$q_X q_Y q_Z = \big(c_r c_p c_y - s_r s_p s_y,\ \ s_r c_p c_y + c_r s_p s_y,\ \ c_r s_p c_y - s_r c_p s_y,\ \ c_r c_p s_y + s_r s_p c_y\big),$$

the expression in the source. to_euler_internal_XYZ inverts it.

Spherical linear interpolation

Unit quaternions form the 3-sphere $S^3 \subset \mathbb R^4$, and the shortest path between two orientations is a great-circle arc. Let $q_1, q_2$ be unit with $\cos\Omega = q_1\cdot q_2$, $0 < \Omega < \pi$. The unit vector

$$q_\perp = \frac{q_2 - \cos\Omega\, q_1}{\sin\Omega}$$

is orthogonal to $q_1$, and the arc is $\gamma(\varphi) = \cos\varphi\, q_1 + \sin\varphi\, q_\perp$. At $\varphi = t\Omega$,

$$\begin{aligned} \gamma(t\Omega) &= \frac{\sin\Omega\cos t\Omega - \cos\Omega\sin t\Omega}{\sin\Omega}\,q_1 + \frac{\sin t\Omega}{\sin\Omega}\,q_2 \\ &= \frac{\sin\big((1-t)\Omega\big)}{\sin\Omega}\,q_1 + \frac{\sin(t\Omega)}{\sin\Omega}\,q_2 , \end{aligned}$$

which is the formula slerp evaluates. It moves at constant angular speed, and the rotation angle grows linearly in $t$. Because of the double cover, slerp first replaces $q_2$ by $-q_2$ when $q_1\cdot q_2 < 0$, so that $\Omega \le \pi/2$ and the rotation takes the short way around. It clamps $\cos\Omega$ to $[-1, 1]$ before acos, since rounding can push the dot product of two unit quaternions slightly past 1.

Powers

pow_by_int uses $q^{2m} = (q^m)^2$ and $q^{2m+1} = q\,(q^m)^2$, which needs only associativity; all factors are powers of the same $q$, so they commute with each other. Negative exponents use $q^{-n} = (q^{-1})^n$.

pow_by_T uses the polar form. Every non-zero $q$ is $q = |q|\,(\cos\varphi + \hat{\mathbf n}\sin\varphi)$ with $\varphi \in [0, \pi]$ and a unit vector $\hat{\mathbf n}$, and the span of $1$ and $\hat{\mathbf n}$ is a copy of $\mathbb C$ (because $\hat{\mathbf n}^2 = -1$). De Moivre’s formula then defines

$$q^t = |q|^t\,\big(\cos t\varphi + \hat{\mathbf n}\sin t\varphi\big).$$

The implementation computes $\varphi = \arcsin\lvert\mathbf u\rvert/|q|$, which is correct only for $\varphi \le \pi/2$, that is for a non-negative scalar part; see known deviations.

Design decisions

Generic components with per-function bounds

Problem. Quaternions are useful over exact integers (number theory, testing identities) and over Double (geometry), and the luna-generic ecosystem wants one type that fits both.

Options. A Double-only type; a type parameterized by a single “real number” trait; or a type with no bound on T whose functions each require only what they use.

Choice. The last. Quaternion[T] has no bound; + needs Add, * needs Mul + Sub + Add, square_len needs Mul + Add, and only functions with square roots or trigonometry ask for DoubleConvert. This follows the Luna-Flow principle that code depends on the smallest trait composition that states its needs, and it lets Quaternion[Int] check algebraic identities exactly.

Scalar–vector storage

The type stores r : T and vec : (T, T, T) rather than four separate fields. The scalar–vector form is how the product, the conjugate, the inverse and the rotation are derived above, so the code reads like the derivations (dot, cross, scale on vec). The type is abstract, which keeps the representation free to change; the cost is that there are no component accessors yet.

Why Ring and not Field

luna-generic’s Field is Ring + Inverse + Div, and generic code written against it is entitled to the laws of a commutative field, such as $a b = b a$ or $a / b \cdot c = a c / b$. Quaternions satisfy every ring law and have inverses, but not commutativity, so claiming Field would let such code compute wrong answers silently. The package therefore implements Zero, One, AddMonoid, MulMonoid, Semiring and Ring (all lawful for $\mathbb H$ over a commutative T), plus the operation traits Inverse and Conjugate, and stops there. Division is still available through Div and left_div, but generic code has to ask for it explicitly.

/ is right division

Problem. With non-commuting multiplication, q / r must pick a side. Before 0.2.0 it computed $r^{-1} q$.

Choice. Since 0.2.0, q / r is $q\,r^{-1}$, the solution of $x r = q$. This is the usual convention for a division ring and keeps Div consistent with Inverse in the way it is for the luna-generic scalar types: a / b == a * b.inv(). It also reads naturally for rotations: to / from is the rotation that, applied after from, gives to. The other quotient stays available as left_div, with a name that says which side the inverse is on. The change is breaking and is recorded in the changelog.

The Double bridge

Square roots and trigonometric functions exist for Double only, and the package depends on nothing but luna-generic, which has no analytic traits. DoubleConvert is a two-method trait that converts a component to Double and back, so magnitude, normalize, slerp, pow_by_T and the Euler conversions stay generic in signature while computing in Double. It is implemented for Double (as the identity) and Int (truncating from_double). The trait is declared pub rather than pub(open), so other packages cannot implement it for their own types.

Exact structural equality

Eq and Hash are derived and compare components exactly. Equality “as rotations” ($q \sim -q$) or “within a tolerance” depends on the application, so it is left to the caller (the tutorial shows one way). Exact equality also keeps Eq consistent with Hash.

Explicit method promotion

MoonBit 0.10 no longer turns trait implementations into methods implicitly. src/extends.mbt promotes the operators, equal, hash, to_string, conjugate, zero, one and inv, which are natural operations on a quaternion. The method forms not_equal, hash_combine, output, default and to_repr are kept, deprecated and hidden, for compatibility.

Correctness and invariants

Algebraic laws

Over a commutative ring T with exact arithmetic (Int, BigInt), the implementation satisfies:

• the ring laws: $(\mathbb H, +, 0)$ is an abelian group, $(\mathbb H, \cdot, 1)$ a monoid, and $\cdot$ distributes over $+$ on both sides;
• $\overline{pq} = \bar q\,\bar p$, $\bar{\bar q} = q$, $q\bar q = \bar q q = |q|^2$ and $|pq|^2 = |p|^2|q|^2$;
• for exact division (a field T such as rationals): $q\,q^{-1} = q^{-1}q = 1$, $(q / r)\,r = q$ and $r\,(q.\mathrm{left\_div}(r)) = q$.

The package tests check the ring laws on integer samples, non-commutativity, and the two division identities.

Floating-point error

Over Double the laws hold up to rounding. Each component of $pq$ is a sum of four products of a component of $p$ with a component of $q$, evaluated with four roundings of products and three of additions. The standard bound for such an inner product³3 N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM 2002, §3.1: $\lvert\mathrm{fl}(x^{\mathsf T}y) - x^{\mathsf T}y\rvert \le \gamma_n \lvert x\rvert^{\mathsf T}\lvert y\rvert$ for any order of summation. gives, componentwise,

$$\lvert \mathrm{fl}(pq)_m - (pq)_m \rvert \le \gamma_4 \sum_{l=1}^{4} \lvert p_l\rvert\,\lvert q_{\sigma_m(l)}\rvert \le \gamma_4\, |p|\,|q|, \qquad \gamma_n = \frac{n u}{1 - n u},$$

where $\sigma_m$ is the permutation of $q$‘s components that appears in component $m$ and the second inequality is Cauchy–Schwarz. Summing over the four components, $\lVert\mathrm{fl}(pq) - pq\rVert \le 2\gamma_4 |p||q|$, so by the multiplicativity of the norm

$$(1 - 2\gamma_4)\,|p|\,|q| \le |\mathrm{fl}(pq)| \le (1 + 2\gamma_4)\,|p|\,|q| .$$

A product of $n$ unit quaternions therefore has norm $1 + \delta$ with $|\delta| \lesssim 8nu$ to first order, where $u = 2^{-53} \approx 1.1 \times 10^{-16}$; in practice the errors partly cancel and the drift is smaller. Because rotate assumes $|q| = 1$, a drifted quaternion scales vectors by about $1 + 2\delta$. Long chains should call normalize periodically; its result has norm $1$ within a few units of $u$ (for example $1 + 2i + 3j + 4k$ normalizes to a quaternion of computed norm $0.9999999999999999$).

Division by zero and degenerate inputs

The package does not check its inputs; degenerate cases follow the arithmetic of T:

Input	Double	Int
inv, /, left_div with a zero divisor	NaN components	runtime trap (integer division by zero)
normalize of zero	zero returned unchanged	zero returned unchanged
pow_by_T of zero	zero returned unchanged	zero returned unchanged
from_axis_angle with a zero axis	NaN vector part	runtime trap
slerp of a zero input	the zero input is used unnormalized	not meaningful

Over Int, inv is zero unless $|q|^2 = 1$, and every function that goes through DoubleConvert truncates.

Thresholds

slerp falls back to normalized linear interpolation when $\cos\Omega > 0.9995$, that is $\Omega < 0.0316$ rad. There $\sin\Omega < 0.032$, and dividing by it would amplify rounding errors in the weights; the linear path, on the other hand, deviates from the arc by less than $5.1 \times 10^{-7}$ rad on $S^3$ (about $10^{-6}$ rad of rotation angle) at the threshold, and less below it.

The Euler conversions treat $|\sin b| \ge 0.9998$ as gimbal lock. Rounding $b$ to $\pm\pi/2$ there costs up to $\pi/2 - \arcsin 0.9998 \approx 0.020$ rad in the middle angle.

Complexity

Every operation is $O(1)$ on four components: a Hamilton product costs 16 multiplications and 12 additions, rotate 18 multiplications and 12 additions. pow_by_int(n) uses $O(\log |n|)$ products and recursion depth.

Known deviations

These are deviations of the current implementation from the mathematics above. They are documented here rather than hidden, and are reported for fixing:

• from_euler with "ZYX" or "YZX" does not evaluate a product of axis rotations; the results are not unit quaternions. "ZXY" takes its angles by position while "YXZ" takes them by axis.
• to_euler_external_XZY (and therefore to_euler_internal_YZX) evaluates the intrinsic X-Y-Z extraction instead of the named sequence.
• The gimbal-lock branch of the Euler extraction does not isolate the free angle (see above), and it reports the warning with println.
• pow_by_T computes $\varphi$ with asin, which cannot exceed $\pi/2$; $\operatorname{atan2}(\lvert\mathbf u\rvert, w)$ would cover $[0, \pi]$.

Alternatives rejected

• Implementing Field. Rejected because generic Field code may rely on commutativity (see above).
• Keeping / as left division. It contradicted a / b == a * b.inv() and the usual reading of to / from for rotations.
• A Double-only type. It would lose exact arithmetic over Int and the luna-generic ring structure for other scalar types.
• Storing four named fields or a rotation matrix. Four fields hide the scalar–vector structure the formulas use; a matrix has nine entries, needs re-orthogonalization instead of a cheap normalization, and cannot be interpolated as simply.
• Comparing with a tolerance in Eq. A tolerance-based == is not transitive and cannot agree with Hash.

Boundaries

The package deliberately does not:

• provide rotation matrices, exponential and logarithm maps, or quaternion calculus (dual quaternions, derivatives of orientations);
• check its inputs or return Result: degenerate inputs follow the arithmetic of T, and rotate trusts that its receiver is a unit quaternion;
• implement luna-generic Field or Num, or make quaternions an ordered type;
• depend on arithmetic, linear-algebra or luna-complex; its only Luna-Flow dependency is luna-generic;
• choose a tolerance for comparing rotations; that is left to the caller.

Footnotes

1. W. R. Hamilton, “On quaternions; or on a new system of imaginaries in algebra”, 1843. The relations are the ones he carved into Brougham Bridge in Dublin. ↩

2. For example $1 \mapsto I$, $i \mapsto \begin{pmatrix} i & 0 \\ 0 & -i\end{pmatrix}$, $j \mapsto \begin{pmatrix} 0 & 1 \\ -1 & 0\end{pmatrix}$, $k \mapsto \begin{pmatrix} 0 & i \\ i & 0\end{pmatrix}$. These matrices satisfy Hamilton’s relations, and matrix multiplication is associative. ↩

3. N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM 2002, §3.1: $\lvert\mathrm{fl}(x^{\mathsf T}y) - x^{\mathsf T}y\rvert \le \gamma_n \lvert x\rvert^{\mathsf T}\lvert y\rvert$ for any order of summation. ↩