The quaternions—and type classes

09 Mar 2022

[ examples  quaternions  type classes  ]

The quaternion number system is an extension of the complex numbers to 4 dimensions, introduced by Hamilton in 1843. I translated the HOL Light formalisation of quaternions into Isabelle/HOL some time ago. One notable feature of the formalisation (taken from the Isabelle/HOL formalisation of the complex numbers) is that its definition can be regarded as coinductive. Moreover, continuing the previous post about axiomatic type classes, we have a dramatic demonstration of how quickly a new class of numbers can be made native (so to speak).

Defining the type

Quaternions have the form $a + b \mathbf{i} + c \mathbf{j} + d \mathbf{k}$ where $a$, $b$, $c$ and $d$ are real numbers and $\mathbf{i}$, $\mathbf{j}$, $\mathbf{k}$ are the primitive quaternions, satisfying a number of laws such as

\[\mathbf{i}^2 = \mathbf{j}^2 = \mathbf{k}^2 = \mathbf{i}\mathbf{j}\mathbf{k} = -1.\]

It would be natural to represent quaternions by 4-tuples, but it is even simpler to represent them (like the complex numbers) as a codatatype. A codatatype is dual to a datatype; just as the latter is specified by enumerating its constructor functions, the former is specified by enumerating its selector functions (sometimes called destructors or projections). Category theorists should note that codatatypes are terminal coalgebras. Here we get something close to a 4-tuple.

codatatype quat = Quat (Re: real) (Im1: real) (Im2: real) (Im3: real)

The selectors Im1, Im2, Im3 will return the coefficients of to $\mathbf{i}$, $\mathbf{j}$, $\mathbf{k}$, respectively, while Re will return the real part of a quaternion. It is trivial to prove that two quaternions are equal if and only if their four components all coincide.

lemma quat_eq_iff: "x = y ⟷ Re x = Re y ∧ Im1 x = Im1 y ∧ Im2 x = Im2 y ∧ Im3 x = Im3 y"
  by (auto intro: quat.expand)

Defining the primitive operators

It is now possible to define the primitive quaternions, $\mathbf{i}$, $\mathbf{j}$ and $\mathbf{k}$. With corecursion, new entities are defined in terms of the behaviour of the selectors. For example, the quaternion $\mathbf{i}$ yields zero for all selectors except Im1. The quaternions $\mathbf{j}$ and $\mathbf{k}$ are defined analogously.

primcorec quat_ii :: quat  ("𝗂")
  where "Re 𝗂 = 0" | "Im1 𝗂 = 1" | "Im2 𝗂 = 0" | "Im3 𝗂 = 0"

primcorec quat_jj :: quat  ("𝗃")
  where "Re 𝗃 = 0" | "Im1 𝗃 = 0" | "Im2 𝗃 = 1" | "Im3 𝗃 = 0"

primcorec quat_kk :: quat  ("𝗄")
  where "Re 𝗄 = 0" | "Im1 𝗄 = 0" | "Im2 𝗄 = 0" | "Im3 𝗄 = 1"

Addition and subtraction: an abelian group

The development now proceeds by showing that quaternions are an instance of one type class after another. The first is ab_group_add, the class of abelian groups with the signature $(0, {+}, {-})$. For this instance declaration to be accepted, definitions must be provided for zero, addition, unary minus and subtraction, and they must be shown to satisfy the usual group axioms. The proofs of the latter turn out to take a single line. The definitions again use corecursion, and as we can see below, these operations are all defined componentwise.

instantiation quat :: ab_group_add
begin

primcorec zero_quat
  where "Re 0 = 0" |"Im1 0 = 0" | "Im2 0 = 0" | "Im3 0 = 0"

primcorec plus_quat
  where
    "Re (x + y) = Re x + Re y"
  | "Im1 (x + y) = Im1 x + Im1 y"
  | "Im2 (x + y) = Im2 x + Im2 y"
  | "Im3 (x + y) = Im3 x + Im3 y"

primcorec uminus_quat
  where
    "Re (- x) = - Re x"
  | "Im1 (- x) = - Im1 x"
  | "Im2 (- x) = - Im2 x"
  | "Im3 (- x) = - Im3 x"

primcorec minus_quat
  where
    "Re (x - y) = Re x - Re y"
  | "Im1 (x - y) = Im1 x - Im1 y"
  | "Im2 (x - y) = Im2 x - Im2 y"
  | "Im3 (x - y) = Im3 x - Im3 y"

instance
  by standard (simp_all add: quat_eq_iff)

end

This instantiation declaration not only makes the entire library of facts about ab_group_add available to quaternions, but also any associated operators. This particular lemma (which actually appears later in the theory) states that Re distributes over finite summations.

lemma Re_sum [simp]: "Re(sum f S) = sum (λx.  Re(f x)) S" for f :: "'a ⇒ quat"
  by (induct S rule: infinite_finite_induct) auto

Scalar multiplication: a vector space

The development continues, following the HOL Light original while continually looking for opportunities to instantiate an appropriate type class. Next is to show that quaternions are a vector space, so we instantiate real_vector. The (overloaded) scalar multiplication operator is called scaleR and its corecursive definition is self-explanatory. The required properties (which amount to linearity) have a one line proof.

instantiation quat :: real_vector
begin

primcorec scaleR_quat
  where
    "Re (scaleR r x) = r * Re x"
  | "Im1 (scaleR r x) = r * Im1 x"
  | "Im2 (scaleR r x) = r * Im2 x"
  | "Im3 (scaleR r x) = r * Im3 x"

instance
  by standard (auto simp: quat_eq_iff distrib_left distrib_right scaleR_add_right)

end

Multiplication: real algebra with unit

According to legend, Hamilton struggled for years to define a well-behaved multiplication operation for the three-dimensional space he was working on, suddenly realising that it could be done if he introduced a fourth component. As always, corecursion defines the operations (1 and $\times$) in terms of the four selector functions. A real algebra with unit satisfies the axioms of a ring and the multiplication also commutes with scalar multiplication. However, quaternion multiplication is not commutative. I was pleased to find that all the necessary type classes for this development were already available in Isabelle/HOL.

instantiation quat :: real_algebra_1

begin

primcorec one_quat
  where "Re 1 = 1" | "Im1 1 = 0" | "Im2 1 = 0" | "Im3 1 = 0"

primcorec times_quat
  where
    "Re (x * y) = Re x * Re y - Im1 x * Im1 y - Im2 x * Im2 y - Im3 x * Im3 y"
  | "Im1 (x * y) = Re x * Im1 y + Im1 x *  Re y + Im2 x * Im3 y - Im3 x * Im2 y"
  | "Im2 (x * y) = Re x * Im2 y - Im1 x * Im3 y + Im2 x *  Re y + Im3 x * Im1 y"
  | "Im3 (x * y) = Re x * Im3 y + Im1 x * Im2 y - Im2 x * Im1 y + Im3 x *  Re y"

instance
  by standard
     (auto simp: quat_eq_iff distrib_left distrib_right right_diff_distrib left_diff_distrib)

end

Suddenly by magic we’ve gained the ability to use numerals to denote quaterions. Numerals can be hundreds of digits long, represented internally by a symbolic binary notation. Arithmetic on them is surprisingly efficient. This type class instantiation has also given us the functions of_nat, of_int, of_real, injecting other numeric types (the natural numbers, integers, reals) into the quaternions.

And now division: a real division algebra

The next instantiation, to the type class real_div_algebra—the quaternions do not form a field—requires us to define the multiplicative inverse and then (trivially) division. For the first time, justification of the type class axioms is not trivial, hence the four show commands below.

instantiation quat :: real_div_algebra
begin

primcorec inverse_quat
  where
    "Re (inverse x) = Re x / ((Re x)2 + (Im1 x)2 + (Im2 x)2 + (Im3 x)2)"
  | "Im1 (inverse x) = - (Im1 x) / ((Re x)2 + (Im1 x)2 + (Im2 x)2 + (Im3 x)2)"
  | "Im2 (inverse x) = - (Im2 x) / ((Re x)2 + (Im1 x)2 + (Im2 x)2 + (Im3 x)2)"
  | "Im3 (inverse x) = - (Im3 x) / ((Re x)2 + (Im1 x)2 + (Im2 x)2 + (Im3 x)2)"

definition "x div y = x * inverse y" for x y :: quat

instance
proof
  show "⋀x::quat. x ≠ 0 ⟹ inverse x * x = 1"
    by (auto simp: quat_eq_iff add_nonneg_eq_0_iff
        power2_eq_square add_divide_distrib [symmetric] diff_divide_distrib [symmetric])
  show "⋀x::quat. x ≠ 0 ⟹ x * inverse x = 1"
    by (auto simp: quat_eq_iff add_nonneg_eq_0_iff power2_eq_square add_divide_distrib [symmetric])
  show "⋀x y::quat. x div y = x * inverse y"
    by (simp add: divide_quat_def)
  show "inverse 0 = (0::quat)"
    by (auto simp: quat_eq_iff)
qed

end

Here is a trivial example illustrating the use of both numerals and division for quaternions:

lemma "(x::quat)*1000/1001 = x"

And here’s another thing. If we type in the line above, Isabelle instantly and automatically detects that the claim is false, producing the following message:

Auto Quickcheck found a counterexample:
  x = Quat (real_of_rat (Fract (- 1) 1))
       (real_of_rat (Fract (- 1) 1))
       (real_of_rat (Fract (- 1) 1))
       (real_of_rat (Fract (- 1) 1))
Evaluated terms:
  x * 1000 / 1001 =
    Quat (- (1000 / 1001)) (- (1000 / 1001))
     (- (1000 / 1001)) (- (1000 / 1001))

It also works with much larger numbers. Counterexample detection is not always possible, but it works in much more sophisticated situations than the one shown, and it is a tremendous time saver. Nitpick is another counterexample finder, working on different principles from Quickcheck and also highly effective.

Real normed division algebra

Now we define a norm on quaternions, instantiating the type class real_normed_div_algebra. The norm of $a + b \mathbf{i} + c \mathbf{j} + d \mathbf{k}$ is simply $\sqrt{a^2+b^2+c^2+d^2}$. This type class requires definitions of associated operators such as dist, the metric space distance function. The declaration explicitly proves some required properties of norm, such as the triangle inequality, while the properties of the derived operations are trivial by definition.

instantiation quat :: real_normed_div_algebra
begin

definition "norm z = sqrt ((Re z)2 + (Im1 z)2 + (Im2 z)2 + (Im3 z)2)"

definition "sgn x = x /R norm x" for x :: quat

definition "dist x y = norm (x - y)" for x y :: quat

definition [code del]:
  "(uniformity :: (quat × quat) filter) = (INF e∈{0 <..}. principal {(x, y). dist x y < e})"

definition [code del]:
  "open (U :: quat set) ⟷ (∀x∈U. eventually (λ(x', y). x' = x ⟶ y ∈ U) uniformity)"

lemma norm_eq_L2: "norm z = L2_set (quat_proj z) {..3}"
  by (simp add: norm_quat_def L2_set_def numeral_3_eq_3)

instance
proof
  fix r :: real and x y :: quat and S :: "quat set"
  show "(norm x = 0) ⟷ (x = 0)"
    by (simp add: norm_quat_def quat_eq_iff add_nonneg_eq_0_iff)
  have eq: "L2_set (quat_proj (x + y)) {..3} = L2_set (λi. quat_proj x i + quat_proj y i) {..3}"
    by (rule L2_set_cong) (auto simp: quat_proj_add)
  show "norm (x + y) ≤ norm x + norm y"
    by (simp add: norm_eq_L2 eq L2_set_triangle_ineq)
  show "norm (scaleR r x) = ¦r¦ * norm x"
    by (simp add: norm_quat_def quat_eq_iff power_mult_distrib distrib_left [symmetric] real_sqrt_mult)
  show "norm (x * y) = norm x * norm y"
    by (simp add: norm_quat_def quat_eq_iff real_sqrt_mult [symmetric]
        power2_eq_square algebra_simps)
qed (rule sgn_quat_def dist_quat_def open_quat_def uniformity_quat_def)+

end

Having done this instantiation, Isabelle knows that quaternions are a topological space. Quaternions now enjoy the full vocabulary of open and closed sets, limits and continuity, derivatives and infinite sums. Because HOL Light does not have type classes, its formalisation of quaternions includes whole files of library material, copied and pasted with the obvious type substitutions.

Real inner product space

Our final instantiation is to real_inner, the type class of real-valued inner product spaces. This step unlocks additional library material, for dot products. Five necessary properties of the given definition of inner products, such as the commutative and distributive laws, are shown explicitly.

instantiation quat :: real_inner
begin

definition inner_quat_def:
  "inner x y = Re x * Re y + Im1 x * Im1 y + Im2 x * Im2 y + Im3 x * Im3 y"

instance
proof
  fix x y z :: quat and r :: real
  show "inner x y = inner y x"
    unfolding inner_quat_def by (simp add: mult.commute)
  show "inner (x + y) z = inner x z + inner y z"
    unfolding inner_quat_def by (simp add: distrib_right)
  show "inner (scaleR r x) y = r * inner x y"
    unfolding inner_quat_def by (simp add: distrib_left)
  show "0 ≤ inner x x"
    unfolding inner_quat_def by simp
  show "inner x x = 0 ⟷ x = 0"
    unfolding inner_quat_def by (simp add: add_nonneg_eq_0_iff quat_eq_iff)
  show "norm x = sqrt (inner x x)"
    unfolding inner_quat_def norm_quat_def
    by (simp add: power2_eq_square)
qed

end

Up to this point I have hidden essentially nothing from the full quaternion development. There’s only one small technical lemma and some commands to prevent syntactic ambiguities between our 𝗂 and the version of 𝗂 belonging to the complex numbers. Starting from a few basic definitions, we issued six instantiation declarations to put at our disposal the material from three decades of library construction.

That’s only the beginning. The quaternion development goes on to derive specific properties and applications of quaternions, all “borrowed” from the HOL Light original. There’s also a development of Octonions by Angeliki Koutsoukou-Argyraki.

The point of this example is simply to see how much we can accomplish with type classes alone. What can’t be done with type classes can possibly be done using locales, which are strictly more general, but that’s a topic for another time.