Symbolic Quaternion Operations with SymPy

Representing a quaternion as a scalar-first list [w, x, y, z] keeps every component visible to SymPy. This is useful when deriving, expanding, or simplifying formulas for special rotation sequences. For ordinary quaternion arithmetic, SymPy also provides the native sympy.algebras.quaternion.Quaternion class.

The axis-angle constructor below normalizes its axis before creating a unit quaternion. Without this step, a non-unit input axis would scale the vector part and the result would no longer represent a pure rotation.

from sympy import Matrix, cos, sin, sqrt


def axis_angle(axis, angle):
    """Return [w, x, y, z] for a rotation around a nonzero 3D axis."""
    if len(axis) != 3:
        raise ValueError("The rotation axis must have three components.")

    vector = Matrix(axis)
    length = sqrt(vector.dot(vector))
    if length.is_zero is True:
        raise ValueError("The rotation axis must be nonzero.")

    unit_axis = vector / length
    half_angle = angle / 2

    return [
        cos(half_angle),
        *(component * sin(half_angle) for component in unit_axis),
    ]


def multiply(left, right):
    """Return the Hamilton product left * right in scalar-first order."""
    w1, x1, y1, z1 = left
    w2, x2, y2, z2 = right

    return [
        w1*w2 - x1*x2 - y1*y2 - z1*z2,
        w1*x2 + x1*w2 + y1*z2 - z1*y2,
        w1*y2 - x1*z2 + y1*w2 + z1*x2,
        w1*z2 + x1*y2 - y1*x2 + z1*w2,
    ]


def rotation_x(angle):
    return axis_angle([1, 0, 0], angle)


def rotation_y(angle):
    return axis_angle([0, 1, 0], angle)


def rotation_z(angle):
    return axis_angle([0, 0, 1], angle)

Symbolic Rotation Composition

With the usual active-rotation convention, the product below applies the X rotation first, followed by Y and then Z. Quaternion multiplication is not commutative, so changing the product order changes the resulting orientation.

from sympy import simplify, symbols

roll, pitch, yaw = symbols("roll pitch yaw", real=True)

rotation = multiply(
    multiply(rotation_z(yaw), rotation_y(pitch)),
    rotation_x(roll),
)

for component in rotation:
    print(simplify(component))

The component equations above are the explicit Hamilton product. Keeping them as lists makes it straightforward to pass individual expressions to expand(), simplify(), factor(), or SymPy's code-generation functions.

Checking Against SymPy's Quaternion Class

The list implementation can be checked directly against SymPy's built-in class:

from sympy.algebras.quaternion import Quaternion

left = [1, 2, 3, 4]
right = [5, 6, 7, 8]

explicit = multiply(left, right)
native = Quaternion(*left) * Quaternion(*right)

assert explicit == [native.a, native.b, native.c, native.d]
print(explicit)  # [-60, 12, 30, 24]

For numerical JavaScript applications after the symbolic formulas have been reduced, Quaternion.js uses the same scalar-first component convention.