On adding quaternions to Scheme

The numbers supported by the Scheme programming language form a tower, where the set of numbers at each level is a subset of the set at the level above. The levels are, respectively: integers, rationals, reals,¹ and complex numbers. This article investigates the potential extension of Scheme's numeric tower upward to include quaternions [3], a set of hypercomplex numbers that properly includes the complex numbers, and by extension, all the numbers below.

1  Quaternions

A quaternion q is a 4-tuple (w,x,y,z), where w, x, y, and z are real. To give it the aspect of a number, it is written w + ix + jy + kz, where i, j, and k are imaginary units that satisfy the conditions ii = jj = kk = - 1, ij = - ji = k, jk = - kj = i, and ki = - ik = j.² Clearly, multiplication on quaternions does not commute. w is q's real (or scalar) part, and x, y, and z are its i-, j-, and k-part respectively. ix + jy + kz is q's vector part. A quaternion with zero j- and k-parts is an ordinary complex number, and if, further, its i-part is also zero, it is just a real number. A quaternion with zero real part is called a vector quaternion or, simply, vector.³

Notation: We will use q_n (where n is either an integer or absent) interchangeably with w_n + ix_n + jy_n + kz_n. q_n's vector part ix_n + jy_n + kz_n is also written v_n and its unit vector (see below) is written u_n. E => e means that the Scheme expression E evaluates to a Scheme value e. In order to prevent confusion, we add a subscript q to the name of a Scheme procedure that has been made quaternion-aware. For example, exp is the exponential procedure described in R5RS [4]; whereas expq is the same procedure that has been extended to accept a quaternion argument.

2  Lexical representation of quaternions

Recall that Scheme complex numbers that are not also real are represented in rectangular (or cartesian) notation as a+bi or a-bi, where a is optional (for numbers that are purely imaginary) or real, and b is a non-negative real. Examples:

+i
-i
+2i
-2i
3+4i

(The sign is never optional, eg, i must be represented as +i rather than i, which latter is parsed as an identifier.) We can naturally extend this to quaternions: thus q above is written w+xi+yj+zk. The pluses can be minuses, and any part that is zero may be omitted, with the proviso that a vector quaternion must always start with a sign or an explicit zero real part, in order to prevent being mistaken for an identifier. Examples of some Scheme quaternions that are not also complex numbers:

+j
+3k
-2+3i-3k
-6i-6k
3-4i+5j+6k

3  Polar notation

Scheme also permits users to enter complex numbers in a notation that uses polar coordinates, viz, µ@, where µ and are real. A corresponding notation for quaternions is also possible. For example, q = w + ix + jy + kz has the polar (or spherical or curvilinear) coordinates µ, , , , which are related to the rectangular coordinates as follows:⁴

µ, , , are respectively the magnitude, angle, colatitude, and longitude of q.⁵ Extending Scheme's notation for complex numbers, q may be represented as µ @ % & . Any zero angle may be omitted along with its prefix. If the colatitude and longitude are both zero, the number is complex.

As for complex numbers, polar notation is used only to enter quaternions. Scheme outputs quaternion values exclusively in rectangular. Thus,

-3@4%-5&6
=> 
1.960930862590836+0.6440287493492101i+2.0904336369493484-0.6083291310311992k

4  Building quaternions and taking them apart

The predicate quaternion? checks if its argument is a quaternion. Because quaternions now form the pinnacle of Scheme's numeric tower, quaternion? returns true (#t) for any number, and thus agrees with Scheme's number? predicate.

The Scheme complex-number procedures real-part, imag-part, magnitude, and angle are extended to accept quaternion arguments. They return, respectively, the real part, the i-part, the magnitude, and the angle, of their argument. To these are added four new related procedures: jmag-part, kmag-part, colatitude, and longitude return, respectively, the j-part, k-part, colatitude, and longitude of their argument.

The procedure vector-part returns the vector part of its argument.

The procedures make-rectangular and make-polar are extended to take two extra arguments in order to produce quaternions that aren't just complex numbers.

(make-rectangularq 3 4 5 6)
=>  3+4i+5j+6k

(make-polarq -3 4 -5 6)
=> 
1.960930862590836+0.6440287493492101i+2.0904336369493484-0.6083291310311992k

Note that the result printed by Scheme is the rectangular equivalent of the polar coordinates specified.

Every quaternion q has a special vector quaternion u called its unit vector, which is defined as

In Scheme,

(define unit-vector
  (lambda (q)
    (*q (vector-part q)
       (/ (magnitudeq q)))))

The magnitude of a unit vector is of course unity.

5  Arithmetic

Scheme numerical procedures that work on complex numbers can be naturally extended to accept quaternion arguments.⁶

The procedures =, +, - extend in the expected manner for quaternion arguments -- they all operate on the parts severally.

(=q 3+4j 3+4j)      => #t
(=q 3+4j 3+4j 3+4k) => #f
(+q 3+4j 3+4j)      => 6+8j
(+q 3+4j 3+4j 3+4k) => 9+8j+4k
(-q 4+5j 2-3k)      => 2+8k
(-q 2+3j-4k)        => -2-3j+4k

The multiplication procedure *q (using the identities ii = jj = kk = ijk = - 1) can be implemented as follows:

Example:

(*q 4-3i+2j+6k 4+5i-7j+9k)
=>  -9+68i+37j+71k

Generally, multiplication is not commutative. However if two quaternions have the same unit vector disregarding sign, they do commute on multiplication, and furthermore, their product has the same unit vector. The square of a unit vector (its product with itself) is - 1.

If q₁, q₂ are vector quaternions, ie, w₁ = w₂ = 0, then the real part of q₁ q₂ is the negative of the dot product q₁ · q₂ and the vector part of q₁ q₂ is the cross product q₁ q₂.

(define dot-product
  (lambda (q1 q2)
    (if (not (and (= (real-partq q1) 0)
                  (= (real-partq q2) 0)))
        (error 'dot-product "arguments must be vector quaternions")
        (- (real-partq (*q q1 q2))))))

(define cross-product
  (lambda (q1 q2)
    (if (not (and (= (real-partq q1) 0)
                  (= (real-partq q2) 0)))
        (error 'cross-product "arguments must be vector quaternions")
        (vector-part (*q q1 q2)))))

If the imaginary parts of q₁ and q₂ are severally the negatives of each other, then the product q₁ q₂ is real. Such q₁, q₂ are termed conjugates of each other. They have the same magnitude, and their product is the square of that magnitude.

(define conjugate
  (lambda (q)
    (make-rectangularq
      (real-partq q)
      (- (imag-partq q))
      (- (jmag-part q))
      (- (kmag-part q)))))

We can use fact (1) to implement the reciprocal of a quaternion. Multiply both numerator and denominator of 1/q by q's conjugate. Unary / is the Scheme procedure for reciprocals. Example:

(/q 8+3i+4j+5k)
=>  
0.07017543859649122-0.02631578947368421i-0.03508771929824561j-0.043859649122807015k

Multiplication by the reciprocal leads to a definition of division. Dividing q₁ by q₂ is viewed as multiplying q₁ by q₂^-1. The noncommutativity of multiplication allows a choice between putting the reciprocal before or after q₁. Thus division from the left and from the right are different.

In keeping with the visual layout of Scheme's / procedure, we can define (/q q1 q2 ...) with (*q q1 (/q q2) ...), but to prevent confusion, it is best to always use explicit reciprocals (unary /q) and multiplication (*q).

6  Exponential and logarithmic functions

As for complex numbers, we use the Maclaurin series to define quaternion extensions of the transcendental functions. Thus, for quaternion q = w + uv, where uq is its unit vector and w and v are real, we can derive

u functions much the same as i for complex numbers. The difference from the complex case is that u depends on (ie, varies with) q, whereas i is the same for all complex z.

In Scheme,

(define expq
  (lambda (q)
    (let ((w (real-partq q))
          (u (unit-vector q))
          (v (magnitudeq (vector-part q))))
      (*q (exp w)
         (+q (cos v) (*q u (sn v)))))))

Note that the unit vectors of q and e^q are the same, disregarding sign. The inverse of of the exponential function is the natural logarithm,⁷ and is given by:

(where the log on the right only need be defined for positive reals). To make the log function single-valued, we choose tan^-1 (v/w) between - adn . In Scheme,

(define logq
  (lambda (q)
    (let ((w (real-partq q))
          (u (unit-vector q))
          (v (magnitudeq (vector-part q))))
      (+q (*q 1/2 (log (magnitudeq q)))
         (*q u (atan v w))))))

We can use logq to define exptq, since:

In Scheme,

(define exptq
  (lambda (q1 q2)
    (expq (*q (logq q1) q2))))

The sqrtq function then is:

(define sqrtq
  (lambda (q)
    (exptq q 1/2)))

7  Trigonometric functions

For sin q and cos q, we have (using the Maclaurin's series for sinh and cosh):

In Scheme,⁸

(define sinq
  (lambda (q)
    (let ((w (real-partq q)) 
          (u (unit-vector q))
          (v (magnitudeq (vector-part q))))
      (+q (*q (sin w) (cosh v))
         (*q u (cos w) (sinh v))))))

(define cosq
  (lambda (q)
    (let ((w (real-partq q)) 
          (u (unit-vector q))
          (v (magnitudeq (vector-part q))))
      (-q (*q (cos w) (cosh v))
         (*q u (sin w) (sinh v))))))

tan q can be defined as (sin q)/(cos q). Since sin q, and cos q have the same unit vector as q, the direction of the division doesn't matter.

8  Inverse trigonometric functions

The following are derived using Maclaurin's and straightforward quadratic solving:

In Scheme,

(define asinq
  (lambda (q)
    (let ((u (unit-vector q)))
      (-q (*q u (logq (+q (*q u q)
                              (sqrtq (-q 1 (*q q q))))))))))

(define acosq
  (lambda (q)
    (let ((u (unit-vector q)))
      (-q (*q u (logq (+q q 
                              (sqrtq (-q (*q q q) 1)))))))))

(define atanq
  (lambda (q)
    (let ((u (unit-vector q)))
      (*q 1/2 u (logq (*q (+q u q) (/q (-q u q))))))))

9  Signatures

The following collects the signatures of all the quaternion procedures proposed above. The types number and quaternion are synonymous. I shall use the name quaternion merely to emphasize that a function result can be a quaternion. Arguments enclosed in brackets are optional.

quaternion? : any boolean

real-partq : number real
imag-partq : number real
jmag-part : number real
kmag-part : number real

vector-part : number quaternion
unit-vector : number quaternion

magnitudeq : number real
angleq : number real
colatitude : number real
longitude : number real

make-rectangularq : real real [ real real] quaternion
make-polarq : real real [ real real] quaternion

+q : number* quaternion
-q : number+ quaternion
*q : number* quaternion
/q : number+ quaternion

dot-product : number number real
cross-product : number number quaternion
conjugate : number quaternion

expq : number quaternion
logq : number quaternion
exptq : number number quaternion
sqrtq : number quaternion

sinq : number quaternion
cosq : number quaternion
tanq : number quaternion

asinq : number quaternion
acosq : number quaternion
atanq : number quaternion
atan : real real real (two arguments)

10  References

¹ In practice, an implementation does not differentiate between rationals and reals.

² Invoking the associativity of multiplication, these conditions can be more succinctly stated as ii = jj = kk = ijk = - 1.

³ These are indeed the vectors familiarly used in physics to represent 3-dimensional quantities. They are not be confused with Scheme's one-dimensional array data type called vector.

⁴ The rectangular coordinates can be deduced from the polar as follows:

⁵ The more usual terms for magnitude and angle are modulus and amplitude respectively, but we wish to remain compatible with Scheme's already established names for complex numbers.

⁶ Procedures like <, >, <=, >=, positive?, negative?, min, max, floor, ceiling, truncate, round, rationalize two-argument atan, make-rectangular and make-polar are not defined for quaternions, but they are not defined for complex numbers either.

⁷ Following Scheme practice, we use log to mean log_e rather than log₁₀.

⁸ sinh and cosh for reals are not part of the Scheme standard [4] but they can be readily defined using the identities