Algebra and mathematical functions

conj(q)

Return the quaternion conjugate, which flips the sign of each "vector" component.

Examples

julia> conj(quaternion(1,2,3,4))
1 - 2𝐢 - 3𝐣 - 4𝐤

a × b
cross(a, b)

Return the cross product of two pure-vector quaternions. Equal to ½ of the commutator product a*b-b*a.

This extends LinearAlgebra.cross, of which × is an alias, so the same function is used whether × comes from Quaternionic or from LinearAlgebra.

This can be typed as \times<tab> in a Julia-aware editor.

p ⋅ q
dot(p, q)

Return the quaternion conj(p) * q, or, if p and q are both QuatVecs, the usual dot product p[2]*q[2] + p[3]*q[3] + p[4]*q[4] of their vector parts.

For general quaternions, this is the definition that LinearAlgebra assumes for the dot of two numbers, as for Complex numbers, so generic linear algebra works for arrays of quaternions. For example, x' * y for vectors x and y of quaternions is sum(conj.(x) .* y), and qr, svd, least-squares solutions with \, and eigen of a Hermitian matrix all give the quaternionic results. A few functions of LinearAlgebra assume that the entries of a matrix commute, and throw an error for matrices of quaternions; see Matrices of quaternions.

For two QuatVecs, which represent vectors in three dimensions, the dot product of vectors is more useful. It is the scalar part of conj(p) * q. To compare the directions of two general quaternions as vectors in four dimensions — for example, to decide whether two rotors are more anti-parallel than parallel — use real(p ⋅ q), which is the scalar part of p ⋅ q.

This method extends LinearAlgebra.dot, of which ⋅ is an alias. It can be typed as \cdot<tab> in a Julia-aware editor.

normalize(q)

Return a copy of this quaternion, normalized.

Note that this returns the same type as the input quaternion. If you want to convert to a Rotor, just call rotor(q), which includes a normalization step.

This extends LinearAlgebra.normalize for quaternion types.

a ×̂ b

Return the direction of the cross product between a and b; the normalized vector along a×b — unless the magnitude is zero, in which case the zero vector is returned. Both cases give the same element type, which is a floating-point type for integer or rational inputs.

For complex components, the magnitude is the spinor norm computed by absvec, which can vanish for a nonzero null vector. Such a vector cannot be normalized, so it is returned unchanged.

This can be typed as \times<tab>\hat<tab> in a Julia-aware editor.

q ^ s
^(q, s)

Exponentiation operator, equivalent to $\exp(s \log(q))$.

When s is a real number, this is useful for natural "linear" interpolation/extrapolation through quaternion space, starting from q^0 = 1 and going directly to q^1 = q. Specifically, when q is a Rotor, this provides a geodesic on the unit 3-sphere going between those two points. For more general interpolations, see the slerp function.

The type of the result depends on the types of the arguments:

• A Rotor raised to a real or complex power is a Rotor. The norm of the input is assumed to be 1, so the result always has unit norm. This includes Lorentz rotors, with complex components.
• A Quaternion or QuatVec raised to a non-integer power is a Quaternion.
• When s is itself a quaternion, the result is the Quaternion $\exp(s \log(q))$, with the product taken in that order.
• When s is an Integer, the power is computed by repeated multiplication. Then a Rotor gives a Rotor, and a Quaternion or QuatVec gives a Quaternion. With integer components, a nonnegative power keeps integer components, and a negative power gives floating-point components.

The branch cut is the same as for log: a quaternion whose vector part is exactly zero and whose scalar part is negative is treated as $|q| \exp(π𝐤)$. Just as for log, derivatives are not defined there for non-integer s, and will be returned as zero with respect to q.

Non-integer powers choose their branch by comparing component values, so with symbolic component types such as Symbolics.Num they are evaluated as exp(s * log(q)), with the general formula for log. Integer powers work with them directly.

abs(q)

Square root of the sum of the squares of all four components of the quaternion.

For real components, this function uses Julia's built-in hypot function to avoid overflow and underflow.

Note that the result for a Rotor is identically 1, even if that is not true numerically.

For complex components, this is the principal square root of the complex spinor norm abs2, computed with a rescaling that likewise avoids overflow and underflow. It is a complex number in general, and may vanish for a nonzero quaternion. Use norm to obtain the real Euclidean norm.

Examples

julia> abs(quaternion(1,2,4,10))
11.0

abs2(q)

Sum of the squares of all four components of the quaternion.

Note that the result for a Rotor is identically 1, even if that is not true numerically.

For complex components, as used to represent elements of the spacetime algebra (see Lorentz), this is the complex "spinor norm" $\sum_i z_i^2$, which is a complex number that may even vanish for a nonzero quaternion; it is not the Euclidean sum $\sum_i |z_i|^2$. Use norm to obtain the Euclidean norm.

Examples

julia> abs2(quaternion(1,2,4,10))
121

angle(q)

Phase angle in radians of the (spinorial) rotation represented by this quaternion.

Note that this may be different from your interpretation of the angle of a complex number in an important way. Because quaternions act on vectors by conjugation — as in q*v*conj(q) — there are two copies of q involved in that expression; in some sense, a quaternion acts "twice". Therefore, this angle may be twice what you expect from an analogy with complex numbers — depending on how you interpret the correspondence between complex numbers and quaternions.

Also, because quaternions are spinors, rotation by 2π is not the identity operation; only a full rotation through 4π is. Rotation by 2π changes the sign of the quaternion, which has no effect on the rotation of vectors, but does have an effect on rotations of more general objects (spinors). This issue of sign is important when considering the continuity of quaternionic functions (as in interpolations, differentiation, and integration, for example). Therefore, the angle returned by this function will be in the range [0, 2π], rather than the range [0, π] used for complex numbers. If you only care about the effects on vectors, and want to map the result of this function to the range [0, π], you can use θ -> min(θ, 2π-θ).

For real components, this is computed as 2atan(absvec(q), q[1]), which is equivalent to 2absvec(log(q)) but is accurate, along with its derivatives, even near -1. The angle has cone-shaped kinks where the vector part vanishes, at which derivatives are returned as zero.

For complex components, as for a Lorentz transformation, the result is the complex number 2absvec(log(q)). When the rotation and the boost share an axis, this is $θ ± iη$, where $θ$ is the rotation angle, $η$ is the rapidity of the boost, and the sign is positive when the boost is parallel to the rotation axis. In general, the result does not separate so simply into a rotation angle and a rapidity, and the imaginary part may be negative.

Examples

julia> θ=1.2;

julia> R=exp(θ * imz / 2);  # R*v*conj(R) rotates v by θ about the z axis

julia> angle(R)
1.2

julia> angle(exp(6.2 * imz / 2))  # Note that this is greater than π
6.2

exp(q)

Exponential of a quaternion.

The exponential of a quaternion is defined as usual by its power series, which converges for all finite quaternions:

\[\exp(q) = \sum_{k=0}^\infty \frac{q^k}{k!}.\]

The exponential of a QuatVec is a Rotor. The exponential of a Quaternion or of a Rotor is a Quaternion.

Examples

julia> R = exp(imx*π/4)  # Rotation through π/2 (note the extra 1/2) about the x axis
rotor(0.7071067811865476 + 0.7071067811865475𝐢 + 0.0𝐣 + 0.0𝐤)

julia> R * imx * conj(R)
0.0 + 1.0𝐢 + 0.0𝐣 + 0.0𝐤

julia> R * imy * conj(R)
0.0 + 0.0𝐢 + (2.220446049250313e-16)𝐣 + 1.0𝐤

julia> R * imz * conj(R)
0.0 + 0.0𝐢 - 1.0𝐣 + (2.220446049250313e-16)𝐤

Notes

The quaternionic exponential is very easy to calculate by analogy with the complex exponential. We can write a quaternion as $q = s + v\hat{v}$, where $s$ is the scalar part, $v$ is the norm of the pure-vector part, and $\hat{v}$ is a unit vector. Then, $\hat{v}^2 = -1$, so it acts exactly like the imaginary unit $i$ in complex numbers. Obviously, $s$, $v$, and $\hat{v}$ all commute with each other, so the math is simple and we can immediately calculate in analogy with Euler's formula for complex numbers:

\[\exp(q) = \exp(s) \left(\cos(v) + \hat{v}\sin(v)\right).\]

The only special case is when $v=0$, in which case there is no unique choice of $\hat{v}$, but then exp(q) = exp(s). Unfortunately, this means that we need to use a separate branch for this isolated case, which means that automatic differentiation with respect to the vector components will incorrectly produce a derivative of zero at that point if we do nothing but return the value. Moreover, when $v$ is small but nonzero, $v$ is the square root of a small number, and higher derivatives computed through that square root and through $\sin(v)/v$ lose accuracy to cancellation. Therefore, whenever $v^2$ is small, we evaluate $\cos(v)$ and $\sin(v)/v$ as Taylor series in $v^2$, which are accurate to machine precision there, and give accurate derivatives.

With symbolic component types such as Symbolics.Num, the vector part cannot be tested for being small, so the general formula is returned. That expression has a removable singularity where the vector part vanishes, so it should not be evaluated there by direct substitution.

log(q)

Logarithm of a quaternion.

Examples

julia> log(exp(1.2imy))
 + 0.0𝐢 + 1.2𝐣 + 0.0𝐤

julia> log(quaternion(exp(7)))
7.0 + 0.0𝐢 + 0.0𝐣 + 0.0𝐤

julia> log(quaternion(-exp(7)))
7.0 + 0.0𝐢 + 0.0𝐣 + 3.141592653589793𝐤

Notes

The log function is very analogous to the sqrt function [which is essentially exp(log(q)/2)], in that both functions have discontinuous behavior along the negative real axis, and derivatives that blow up as you approach that axis. Therefore, the same caveats about computing derivatives using automatic differentiation apply to both functions.

Integer and rational components are converted to floating-point numbers first. Because the algorithm chooses its branch by comparing component values, symbolic component types such as Symbolics.Num instead get the general formula $\log|q| + \arctan(|\vec{v}|, w)\, \vec{v}/|\vec{v}|$ (from the Symbolics extension), which should not be evaluated on the real axis by direct substitution.

If we decompose the result of this function as $\log(q) = s + \vec{v}$, where $s$ is the scalar part and $\vec{v}$ is the pure-vector part, then clearly $s$ and $\vec{v}$ commute, so their exponentials also commute, and we have

\[q = \exp\left(\log(q)\right) = \exp\left(s\right) \exp\left(\vec{v}\right).\]

Note that the exponential of a pure-vector quaternion is a unit quaternion, so we have decomposed $q$ into a product of a positive real number and a unit quaternion. We already know how to take the logarithm of a positive real number; when $|q|$ is close to 1, we compute it as log1p(abs2(q) - 1)/2, with abs2(q) - 1 evaluated without cancellation, for accuracy. Therefore, in the following we assume that $q$ is a unit quaternion, so that $s=0$ and we only need to find $\vec{v}$.

We now write $\log(q)$ as $v\hat{v}$, where $v$ is just the scalar norm. Note that, because of the periodicity of the exp function, we can assume that $v \in [0, \pi]$ and, in particular, $\sin(v) \geq 0$. Now, expand the exponential as

\[\exp\left(\vec{v}\right) = \exp\left(v \hat{v}\right) = \cos(v) + \hat{v} \sin(v).\]

The input to this function is the right-hand side, but we do not yet know its decomposition into $v$ and $\hat{v}$. But we can find $\cos(v)$ as the scalar part of the input, and $\sin(v)$ as the absvec (since we know that $\sin(v) \geq 0$). Then, we can compute $v = \mathrm{atan}(\sin(v), \cos(v))$. And finally, we simply multiply the vector part of the input by $v / \sin(v)$ to obtain the logarithm. This factor is given accurately by invsinc(v) whenever $|v| \leq \pi/2$. Near the identity, where the vector part is small compared to the scalar part, $\sin(v)$ is the square root of a small number, whose derivatives lose accuracy; there, we instead evaluate $v/\sin(v)$ as a Taylor series in $x = |\vec{q}|^2 / q[1]^2 = \tan^2(v)$, so that automatic differentiation remains accurate.

When $|v| > \pi/2$ (which also implies $\cos(v)<0$), we instead define $v' = \pi - v$, so that $\cos(v') = -\cos(v)$ and $\sin(v') = \sin(v)$, and we can rewrite the problem as

\[\exp\left(v \hat{v}\right) = \cos(v) + \hat{v} \sin(v) = -\cos(v') + \hat{v} \sin(v').\]

We can compute $v' = \mathrm{atan}(\sin(v), -\cos(v))$ accurately, and then we need to multiply the unit vector $\hat{v}$ by $v = \pi - v'$. This algorithm is surprisingly accurate, even when $v$ is extremely close to $\pi$, which implies that the vector part of the input is extremely small.

The only special case remaining to handle is when $\cos(v) < 0$ but $\sin(v)$ is identically zero. In this case, we could throw an error, but this is not usually helpful. Instead, we arbitrarily choose to return $\pi 𝐤$.

For complex components, the same formulas are used with the complex spinor norm abs in place of the Euclidean norm. The branch is then chosen by the sign of the real part of $\cos(v) = q[1]/|q|$. When it is nonnegative, the angle is computed as $v = -i \log\left(\cos(v) + i \sin(v)\right)$, where $i$ is the imaginary unit of the complex components and the logarithm is the principal complex logarithm; this is smooth even where $q[1] = 0$. Otherwise, $v'$ is computed with the principal complex branch of atan.

If q is a Rotor, we return a QuatVec; if q is a general Quaternion, we return a general Quaternion — though if q has unit norm, the scalar part of the result will be zero up to roundoff. Note that, because there is no geometric reason to take the logarithm of a QuatVec, that case is not implemented; if you really need to compute it, you can convert the QuatVec to a Quaternion first.

sqrt(q)

Square root of a quaternion.

The square root of a Rotor is a Rotor, and the square root of a Quaternion or of a QuatVec is a Quaternion. (The square root of a pure vector is not a pure vector.) Integer and rational components are converted to floating-point numbers first.

Examples

julia> q = quaternion(1.2, 3.4, 5.6, 7.8);

julia> sqrtq = √q;

julia> sqrtq^2 ≈ q
true

julia> √quaternion(4.0)
2.0 + 0.0𝐢 + 0.0𝐣 + 0.0𝐤

julia> √quaternion(-4.0)
0.0 + 0.0𝐢 + 0.0𝐣 + 2.0𝐤

Notes

The general formula whenever the denominator is nonzero is

\[\sqrt{q} = \frac{|q| + q} {\sqrt{2(|q| + q[1])}}\]

This can be proven by squaring the numerator, using q = q[1] + q⃗ and |q|^2 = q[1]^2 - q⃗^2 = q[1]^2 + |q⃗|^2. When the denominator is zero — or quite simply whenever q[1] < 0 so that the denominator may be subject to cancellation — we can use the fact that

\[|q| + q[1] = \frac{|q⃗|^2} {|q| − q[1]}\]

to evaluate the expression in a more stable form:

\[\sqrt{q} = \frac{|q⃗|}{\sqrt{2(|q| − q[1])}} + \frac{q⃗}{|q⃗|} \sqrt{\frac{|q| − q[1]}{2}}.\]

For real floating-point components, a quaternion whose norm is so large or so small that these expressions would overflow or underflow is first rescaled by a power of 4, and the result is rescaled by the corresponding power of 2. In the second expression, the unit vector $q⃗/|q⃗|$ is computed from the unscaled vector part, after dividing it by its largest component, so that it keeps its precision even when the vector part is tiny compared to the scalar part. For complex components, the first expression is used whenever the real part of q[1] is nonnegative, and the second with $|q⃗|^2$ evaluated as the complex sum abs2vec(q) otherwise.

Note that whenever the vector part is zero and the scalar part is negative, the solution is not unique (and the denominator above is zero), because it necessarily involves the square root of -1, of which there are infinitely many in the space of quaternions. In this case, we arbitrarily choose the vector part of the result to be proportional to 𝐤, as mentioned above. A reasonable alternative would be to throw an error; instead it is left to the user to check for that condition if it would be a problem.

Analytically, any derivative of this function will blow up as you approach the negative real axis, because the function is discontinuous there. Therefore, you should not expect to be able to accurately compute derivatives of this function at points near the negative real axis (including q=0) using automatic differentiation. Specifically on the negative real axis the derivative is not defined at all, but because of our fixed choice of result here, automatic differentiation will typically (and incorrectly) return the derivative as zero. Ultimately, the reason for this is geometrical, so it should be avoided by the user in any case. Therefore, for the sake of efficiency and accuracy when the problem is geometrically well conditioned, we do not attempt to special-case the derivative at these points.

Because the algorithm chooses its branch by comparing component values, it does not work with symbolic component types such as Symbolics.Num.

norm(q)

Euclidean norm of the quaternion q.

For real components, this equals abs(q); in particular, it is identically 1 for a Rotor. For complex components, it is instead the real, nonnegative Euclidean norm $\sqrt{\sum_i |z_i|^2}$ of the eight real numbers making up the four components, whereas abs and abs2 return the complex spinor norm. This is the norm with which ≈ (isapprox) compares quaternions with complex components.

Examples

julia> norm(quaternion(1,2,4,10))
11.0

julia> norm(Quaternion{ComplexF64}(1, 1im, 0, 0))
1.4142135623730951

abs2vec(q)

Sum of the squares of the three "vector" components of the quaternion.

For complex components, this is the complex sum $\sum_i z_i^2$ over the vector components, as with abs2.

Examples

julia> abs2vec(quaternion(1,2,3,6))
49

absvec(q)

Square root of the sum of the squares of the three "vector" components of the quaternion.

For real components, this function uses Julia's built-in hypot function to avoid overflow and underflow. For complex components, this is the principal square root of the complex sum abs2vec, as with abs.

Examples

julia> absvec(quaternion(1,2,3,6))
7.0