Introduction

Quaternions for Julia

The goal of this package is to provide a simple but flexible and complete implementation of quaternions, without restricting the interpretation of quaternions to being rotations — while also providing extensive support for rotations — along with thorough testing, documentation, and integration with the rest of Julia. Wherever possible, standard functions that work with Complex will also work with Quaternion.

In addition to a basic Quaternion{T} type, we also have Rotor{T} and QuatVec{T} specializations, which can improve the accuracy and efficiency of certain applications. Each of these can be defined over any T<:Number; in addition to the standard primitive types (Float64, etc.), BigFloat, Complex numbers, and Symbolics.Num are tested extensively. Quaternions with complex components also represent Lorentz transformations, as described in Lorentz transformations.

This page is a quick tour of the package. The Manual pages describe each group of functions in more detail, the Functions of time page describes interpolation and angular velocity, and Differentiating by quaternionic arguments describes the support for automatic differentiation. The Theory pages derive the conventions used throughout the package, and the API index at the end lists every documented function and type.

Examples

There are numerous ways to construct a Quaternion — the simplest being to just give the components:

julia> using Quaternionic

julia> q = quaternion(1.0, 2.0, 3.0, 4.0)
1.0 + 2.0𝐢 + 3.0𝐣 + 4.0𝐤

julia> p = quaternion(4, 3, 2, 1)
4 + 3𝐢 + 2𝐣 + 1𝐤

For convenience, you can also pass a vector that will be unpacked for you:

julia> c = [1.0, 2.0, 3.0, 4.0];

julia> q = quaternion(c)
1.0 + 2.0𝐢 + 3.0𝐣 + 4.0𝐤

If the vector has the wrong number of elements, a DimensionMismatch error is raised.

Each quaternion type is parametrized by the types of its components (which are promoted to be all the same type). Any subtype of Number is allowed, and is detected automatically. For example, above q has type Quaternion{Float64}, while p has type Quaternion{Int64}.^[1] The base type may be given explicitly if desired, to override the detected type:

julia> r = Quaternion{Float64}(4, 3, 2, 1)
4.0 + 3.0𝐢 + 2.0𝐣 + 1.0𝐤

The various Float and Int types work well, as do BigFloat and the Num type from Symbolics.jl. In particular, we can use symbolic expressions as components:

julia> using Quaternionic, Symbolics

julia> @variables a b c d e;

julia> quaternion(a-b, b*c, c/d, d+e)
a - b + b*c𝐢 + (c / d)𝐣 + (d + e)𝐤

In analogy with the complex types, the aliases QuaternionF64, QuaternionF32, and QuaternionF16 are provided, as well as the constants imx, imy, and imz, and (for copy-paste convenience) the aliases 𝐢, 𝐣, and 𝐤 (as Unicode bold characters):

julia> QuaternionF64
QuaternionF64 (alias for Quaternion{Float64})
julia> 0.1 + 2.3imx + 4.5imz
0.1 + 2.3𝐢 + 0.0𝐣 + 4.5𝐤
julia> 0.1 + 2.3𝐢 + 0.0𝐣 + 4.5𝐤
0.1 + 2.3𝐢 + 0.0𝐣 + 4.5𝐤

As with the complex im, the result of multiplying imx, etc., with any real number will be a quaternion with the type of the other number.

It is also possible to construct random quaternions using randn with a Quaternion type.

julia> randn(QuaternionF64)
-0.17986445341174084 + 0.5436042462142929𝐢 - 0.20979480846942436𝐣 + 0.3594549687329696𝐤

julia> randn(RotorF32)
rotor(0.18842402 - 0.30743068𝐢 + 0.92128336𝐣 + 0.14567046𝐤)

For a Quaternion, each component is chosen independently from a normal distribution with mean 0 and variance 1/4 (standard deviation 1/2), so that the expected value of abs2(q) is 1. For a QuatVec, each of the three vector components has variance 1/3, and a Rotor is the normalized result. In each case, the resulting quaternion has an equal probability of pointing in any direction — the probability distribution is "isotropic". This is, for example, a good way of choosing a random direction:

julia> normalize(randn(QuatVecF64))
 - 0.3018853063494534𝐢 + 0.4571280910615297𝐣 - 0.8365997670169042𝐤

Note that we have called normalize (originally from LinearAlgebra, but also exported by Quaternionic) to obtain a unit vector in a random direction.

Components of the quaternion are stored as a four-element static array (even for QuatVec):

julia> components(q)
4-element StaticArraysCore.SVector{4, Float64} with indices SOneTo(4):
 1.0
 2.0
 3.0
 4.0

Those components can be indexed directly, just like an ordinary array:

julia> q[1], q[2], q[3], q[4]
(1.0, 2.0, 3.0, 4.0)
julia> q[2:4]
3-element Vector{Float64}:
 2.0
 3.0
 4.0
julia> q[[3, 2]]
2-element Vector{Float64}:
 3.0
 2.0

For convenience, the scalar and vector components can also be accessed in analogy with complex numbers as

julia> real(q)
1.0
julia> imag(q)
3-element view(::StaticArraysCore.SVector{4, Float64}, 2:4) with eltype Float64:
 2.0
 3.0
 4.0

Alternatively, and slightly less efficiently, various parts can be accessed as fields:

julia> q.w, q.x, q.y, q.z
(1.0, 2.0, 3.0, 4.0)
julia> q.vec
3-element Vector{Float64}:
 2.0
 3.0
 4.0
julia> q.re
1.0
julia> q.im
3-element Vector{Float64}:
 2.0
 3.0
 4.0

Again, however, these field accesses incur a slight overhead, so it's more efficient to treat the quaternion as an array and use indexing.

Because quaternions are Numbers, broadcasting treats each quaternion as a single scalar; functions are not broadcast over the components. To apply a function to each component, apply it to the components and construct a new quaternion. For example, this can be particularly helpful when simplifying Symbolics expressions:

julia> @variables q[1:4];  # Defines q[1] through q[4] as symbolic variables

julia> Q = quaternion(q...);

julia> quaternion(simplify.(components(Q * imz * conj(Q)))...)
0 + (2q[1]*q[3] + 2q[2]*q[4])𝐢 + (-2q[1]*q[2] + 2q[3]*q[4])𝐣 + (q[1]^2 - (q[2]^2) - (q[3]^2) + q[4]^2)𝐤

(Though, note that you probably want to use Q(imz) instead of the last expression, when using floating-point numbers, for efficiency reasons.)

The basic algebraic operations work as you would expect:

julia> p + q
5.0 + 5.0𝐢 + 5.0𝐣 + 5.0𝐤
julia> p - q
3.0 + 1.0𝐢 - 1.0𝐣 - 3.0𝐤
julia> p * q
-12.0 + 16.0𝐢 + 4.0𝐣 + 22.0𝐤
julia> q * p  # Note the non-commutativity
-12.0 + 6.0𝐢 + 24.0𝐣 + 12.0𝐤
julia> q / p
0.6666666666666666 + 0.3333333333333333𝐢 + 0.0𝐣 + 0.6666666666666666𝐤

Essential mathematical functions familiar from complex math, such as conj, abs, abs2, log, exp, etc., are also available.

Also note that one of the more useful quaternion operations is group conjugation or "sandwiching". This is the operation that allows a Rotor to actually rotate a QuatVec. For a rotor $R$ and vector $v$, this is expressed mathematically as

\[v′ = R\, v\, R^{-1} = R\, v\, \bar{R}.\]

Here, the first equality is the usual definition of conjugation of $v$ by $R$, while the second equality stems from the fact that for a rotor, $R^{-1} = \bar{R}$. This second operation is sometimes referred to as "sandwiching" or "∗-conjugation" (when working in a ∗-algebra) to distinguish it from the more usual group conjugation involving the inverse. It is this second version that is implemented here by using Quaternions or Rotors as functions. For a Q of either type, and a v<: QuatVec, we have

Q * v * conj(Q) ≈ Q(v)

In particular, Q(v) (the right-hand side) is about twice as efficient as performing the conjugation and two multiplications explicitly (the left-hand side). Note that, for a general Quaternion Q, the result is scaled by abs2(Q).

Rotations

The most common application of quaternions is to represent rotations. A rotation by angle $\vartheta$ about the unit axis $\hat{n}$ is given by the rotor $\exp(\vartheta\, \hat{n} / 2)$, where $\hat{n}$ is a QuatVec. For example, a rotation by $\pi/2$ about the $z$ axis takes the $x$ axis to the $y$ axis:

julia> R = exp(π/4 * imz)
rotor(0.7071067811865476 + 0.0𝐢 + 0.0𝐣 + 0.7071067811865475𝐤)

julia> R(imx) ≈ imy
true

Rotations can be converted to and from other representations, such as rotation matrices, Euler angles, and spherical coordinates (see Rotations and conversions); compared with distance (see Distances and alignment); and interpolated with slerp and squad (see Functions of time):

julia> to_rotation_matrix(R) ≈ [0 -1 0; 1 0 0; 0 0 1]
true

julia> distance(R, one(R)) ≈ π/4
true

julia> slerp(one(R), R, 0.5) ≈ exp(π/8 * imz)
true

Special relativity

Quaternions with complex components represent Lorentz transformations of Minkowski spacetime, through the isomorphism between the complex quaternions and the even subalgebra of the spacetime algebra. The Lorentz type is a Rotor with complex components, and Boost constructs a pure boost. Calling a Lorentz transformation on a 4-vector $[t, x, y, z]$ transforms it:

julia> B = Boost([0.6, 0.0, 0.0]);  # A boost with velocity 0.6 along x

julia> B([1, 0, 0, 0]) ≈ [1.25, 0.75, 0, 0]  # The 4-velocity of a moving particle
true

See Lorentz transformations for more details, and Quaternions and the Spacetime Algebra for the theory.

Sources

Several of the algorithms implemented here are — as far as I know — original to this package, including exp, log, and sqrt. But important contributions came from a number of sources, including the notes by Dantam [1] and Kahan [2]. More advanced functions of time were developed in Boyle et al. [3] and Boyle [4]. The full bibliography is on the References page.