 The Museum of HP Calculators

# Quaternions for the HP-41

This program is Copyright © 2004 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

## Overview

-A quaternion ( or hypercomplex number ) is an expression of the form  q = x + y.i + z.j + t.k
where x , y , z , t  are 4 real numbers and  i ; j ; k  are 3 imaginary units without any linear relation between them.
-More precisely, a quaternion can be regarded as an element ( x ; y ; z ; t ) of  R4
and 1 ; i ; j ; k as the standard basis ( 1 ; 0 ; 0 ; 0 ) ( 0 ; 1 ; 0 ; 0 ) ( 0 ; 0 ; 1 ; 0 ) ( 0 ; 0 ; 0 ; 1 )
-Addition and subtraction are defined like ordinary vectors, but the multiplication is defined by the rules:
i.j = -j.i = k
j.k = -k.j = i       and     i2 = j2 = k2 = -1          Therefore, the product of 2 quaternions is not commutative.But it has the associative and distributive properties.
k.i =  -i.k = j
-The following programs perform a few basic calculations, namely: addition, subtraction, multiplication, reciprocal, polar-rectangular conversions,
natural exponential and logarithm, and raising a quaternion to a power.
-Then, we solve equations in one unknown.
-And the last program illustrates a link between quaternions and 3-dimensional rotations.

Addition and Subtraction of 2 Quaternions

Data Registers:             R00:  unused
•   R01 = a        •  R05 = a'
•   R02 = b        •  R06 = b'       with    q  =  a + b.i + c.j + d.k            R01 thru R08 are to be initialized
•   R03 = c        •  R07 = c'        and    q' = a' + b'.i + c'.j + d'.k
•   R04 = d        •  R08 = d'

Flag:  F01    CF 01 = addition
SF 01 = subtraction
Subroutine:  none

01  LBL "Q+"
02  RCL 04
03  RCL 08
04  FS? 01
05  CHS
06  +
07  RCL 03
08  RCL 07
09  FS? 01
10  CHS
11  +
12  RCL 02
13  RCL 06
14  FS? 01
15  CHS
16  +
17  RCL 01
18  X<> 05
19  FS? 01
20  CHS
21  ST+ 05
22  FS? 01
23  CHS
24  X<> 05
25  END

( 40 bytes / SIZE 009 )

 STACK INPUTS OUTPUTS T / t Z / z Y / y X / x

with      q + q'   =  x + y.i + z.j + t.k    ( if CF 01 )
or       q - q'    =  x + y.i + z.j + t.k    ( if SF 01 )

Example:   q = 2 - 3i + 4j - 7k  and  q' = 1 - 4i + 2j + 5k       calculate  q + q'  and   q - q'

2   STO 01     1  STO 05       CF 01   XEQ "Q+"  yields    3                  and         SF 01   XEQ "Q+"  ( or R/S )    produces     1
-3  STO 02    -4  STO 06                                       RDN    -7                                                                                    RDN        1
4  STO 03      2  STO 07                                       RDN    6                                                                                     RDN         2
-7  STO 04     5  STO 08                                       RDN    -2                                                                                    RDN      -12

Thus      q + q'  = 3 -7i + 6j - 2k      and     q - q' =  1 + i + 2j - 12k

Multiplication of 2 Quaternions

Data Registers:           R00:  unused
•   R01 = a        •  R05 = a'
•   R02 = b        •  R06 = b'       with    q  =  a + b.i + c.j + d.k            R01 thru R08 are to be initialized
•   R03 = c        •  R07 = c'        and    q' = a' + b'.i + c'.j + d'.k
•   R04 = d        •  R08 = d'
Flag: none
Subroutine:  none

N.B:  Synthetic registers M and N can of course be replaced by any other unused registers,
but avoid registers usage conflict between "Q*" and other programs when "Q*" is used as a subroutine.

01  LBL "Q*"
02  RCL 01
03  RCL 05
04  *
05  RCL 02
06  RCL 06
07  *
08  -
09  RCL 03
10  RCL 07
11  *
12  -
13  RCL 04
14  RCL 08
15  *
16  -
17  STO M
18  RCL 01
19  RCL 06
20  *
21  RCL 02
22  RCL 05
23  *
24  +
25  RCL 03
26  RCL 08
27  *
28  +
29  RCL 04
30  RCL 07
31  *
32  -
33  STO N
34  RCL 01
35  RCL 08
36  *
37  RCL 04
38  RCL 05
39  *
40  +
41  RCL 02
42  RCL 07
43  *
44  +
45  RCL 03
46  RCL 06
47  *
48  -
49  RCL 01
50  RCL 07
51  *
52  RCL 03
53  RCL 05
54  *
55  +
56  RCL 02
57  RCL 08
58  *
59  -
60  RCL 04
61  RCL 06
62  *
63  +
64  0
65  X<> N
66  0
67  X<> M
68  END

( 79 bytes / SIZE 009 )

 STACK INPUTS OUTPUTS T / t Z / z Y / y X / x

with      q.q'  =  x + y.i + z.j + t.k

Example1:   q = 2 - 3i + 4j - 7k  and  q' = 1 - 4i + 2j + 5k       calculate  q.q'

2   STO 01     1  STO 05              XEQ "Q*"    yields    17
-3  STO 02    -4  STO 06                                   RDN    23
4  STO 03      2  STO 07                                   RDN    51
-7  STO 04     5  STO 08                                   RDN    13

Thus    q.q' = 17 + 23i + 51j + 13k                      ( remark:    q'.q = 17 - 45i -35j - 7k )

Example2:   A theorem states that  if  q = b.i + c.j + d.k   and   q' = b'.i + c'.j +d'.k  are 2 purely imaginary quaternions,
then the real part of  q.q'  is -U.U'  where  U.U'  is the dot product of the 3D-vectors   U(b;c;d)   and   U'(b';c';d')
and the constituents of the imaginary part of  q.q'  are the 3 constituents of the cross product   UxU'

-For instance, calculate the dot product and the cross product of the 2 vectors:    U(2;-4;9)    U'(3;8;-6)

0  STO 01          0  STO 05         XEQ "Q*"    gives     80
2  STO 02          3  STO 06                             RDN   -48          therefore    U.U'  =   - 80
-4  STO 03          8  STO 07                            RDN     39              and        UxU' =  ( -48 ; 39 ; 28 )
9  STO 04         -6  STO 08                            RDN     28

Reciprocal of a Quaternion

01  LBL "1/Q"
02  STO N
03  X^2
04  STO M
05  X<> T
06  CHS
07  ENTER^
08  X^2
09  ST+ M
10  X<> T
11  CHS
12  ENTER^
13  X^2
14  ST+ M
15  X<> T
16  CHS
17  ENTER^
18  X^2
19  ST+ M
20  CLX
21  X<> M
22  ST/ N
23  ST/ T
24  ST/ Z
25  /
26  0
27  X<> N
28  END

( 49 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T d t Z c z Y b y X a x L / µ2

with      q =  a + b.i + c.j + d.k     ;    1/q =  x + y.i + z.j + t.k    and    µ = ( a2 + b2 + c2 + d2 )1/2 = the modulus of the quaternion q

Example:   q = 2 - 3i + 4j - 7k  evaluate 1/q

-7   ENTER^
4   ENTER^
-3   ENTER^
2   XEQ "1/Q"   yields    0.0256
RDN    0.0385
RDN  -0.0513
RDN    0.0897             thus   1/q = 0.0256 + 0.0385 i  -0.0513 j + 0.0897 k   ( rounded to 4D )
( L = µ2 = 78 )

Note:  Here is a shorter program to compute 1/q using the rectangular-polar conversions:

01  LBL "1/Q"
02  XEQ "POL"
03  1/X
04  X<>Y
05  CHS
06  X<>Y
07  XEQ "REC"
08  END

( 24 bytes / SIZE 000 )

Rectangular-Polar conversions

-A quaternion   q =  x + y.i + z.j + t.k   can also be defined by its modulus  µ = ( x2 + y2 + z2 + t2 )1/2  and 3 angles A , B , C such that:

x = µ cos A                                      A = the amplitude of the quaternion
y = µ sin A cos B                              B = the co-latitude ---------------
z = µ sin A sin B cos C                     C = the longitude   ---------------
t = µ sin A sin B sin C

-The 2 following programs perform these conversions and work in all angular mode.

1°) Rectangular-Polar

01  LBL "POL"
02  R^
03  R^
04  R-P
05  R^
06  R-P
07  R^
08  R-P
09  END

( 17 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T t C Z z B Y y A X x µ L / x

Example:   q = 2 - 3i + 4j - 7k      ( in DEG mode )

-7   ENTER^
4   ENTER^
-3   ENTER^
2   XEQ "POL"     gives       8.8318  = µ      ( actually  781/2 )
RDN    76.9115  = A
RDN  110.4104  = B
RDN  -60.2551  = C

2°) Polar-Rectangular

01  LBL "REC"
02  P-R
03  RDN
04  P-R
05  RDN
06  P-R
07  R^
08  R^
09  END

( 17 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T C t Z B z Y A y X µ x

Example:   Find the rectangular form of the quaternion q defined by:       µ = 7 ;  A = 41° ; B = 37°  ;  C = -168°

[DEG]    -168   ENTER^
37  ENTER^
41  ENTER^
7   XEQ "REC"      yields      5.2830
RDN      3.6677
RDN     -2.7034
RDN     -0.5746        thus      q = 5.2830 + 3.6677 i  -2.7034 j -0.5746 k     ( rounded to 4D )

Natural Exponential and Logarithm

1°) Exponential

-The exponential of a quaternion q is defined by the serie  exp ( q ) = 1 + q + q2/2! + q3/3! + ....... + qn/n! + ........
but the polar form provides a much faster way to compute eq

Flag: none
Subroutine:  "REC"

01  LBL "E^Q"
02  DEG
03  R^
04  R^
05  R-P
06  R^
07  R-P
08  R-D
09  R^
10  E^X
11  XEQ "REC"
12  END

( 24 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T d t Z c z Y b y X a x

with      q =  a + b.i + c.j + d.k     ;    eq =  x + y.i + z.j + t.k

Example:   q = 2 - 3i + 4j - 7k  evaluate eq

-7   ENTER^
4   ENTER^
-3   ENTER^
2   XEQ "E^Q"     produces    -5.0277
RDN       -1.8884
RDN        2.5178
RDN       -4.4062       and    eq  =   -5.0277 -1.8884 i + 2.5178 j - 4.4062 k   ( rounded to 4D )

N.B:  a)  ln (  eq ) is not always equal to q       ( like ordinary complex numbers )
b)  eq.eq'    ---------------------  eq+q'
c)  Hyperbolic functions can easily be defined too.

2°) Logarithm

Flag: none
Subroutine:  "POL"

01  LBL "LNQ"
02  DEG
03  XEQ "POL"
04  LN
05  RDN
06  D-R
07  P-R
08  RDN
09  P-R
10  R^
11  R^
12  END

( 24 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T d t Z c z Y b y X a x

with      q =  a + b.i + c.j + d.k     ;    ln q =  x + y.i + z.j + t.k

Example:   q = 2 - 3i + 4j - 7k  evaluate  ln q

-7   ENTER^
4   ENTER^
-3   ENTER^
2   XEQ "LNQ"    yields       2.1784
RDN      -0.4681
RDN       0.6242
RDN     -1.0923         ln q = 2.1784 -0.4681 i + 0.6242 j -1.0923 k         ( once again rounded to 4D )

Raising a Quaternion to a power

1°) Real exponent

Data Register:
• R00 = r  is to be initialized before executing Q^R
Flag: none
Subroutines:  "POL" "REC"

01  LBL "Q^R"
02  XEQ "POL"
03  X<>Y
04  X<> 00
05  ST* 00
06  Y^X
07  LASTX
08  X<> 00
09  X<>Y
10  XEQ "REC"
11  END

( 30 bytes / SIZE 001 )

 STACK INPUTS OUTPUTS T d t Z c z Y b y X a x

with      q =  a + b.i + c.j + d.k     ;    qr =  x + y.i + z.j + t.k

Example1:   q = 2 - 3i + 4j - 7k  evaluate  q3.14

3.14   STO 00
-7   ENTER^
4   ENTER^
-3   ENTER^
2   XEQ "Q^R"   yields    -445.8830
RDN     286.4187
RDN    -381.8916
RDN     668.3104              therefore   q3.14 =   - 445.8830 + 286.4187 i  - 381.8916 j + 668.3104 k  ( to 4D )

Example2:   q = 2 - 3i + 4j - 7k  calculate one cube root of q  i-e  q1/3

3  1/X   STO 00
-7   ENTER^
4   ENTER^
-3   ENTER^
2   XEQ "Q^R"   yields   1.8635
RDN  -0.3119
RDN   0.4159
RDN  -0.7278        thus     q1/3 =  1.8635 - 0.3119 i + 0.4159 j - 0.7278 k  ( 4D )

Notes:   -A non-zero quaternion has in general n  n-th roots.
-However, -1 has an infinity of square roots!   Namely, if  b2 + c2 + d2 = 1  then   ( b.i + c.j + d.k )2 = -1
-This program gives  ( -1 )1/2 = i

2°) Quaternion exponent

-This program calculates  qq' with the definition   qq' =   eq'.lnq

Data Registers:           R00:  unused
•   R01 = a        •  R05 = a'
•   R02 = b        •  R06 = b'       with    q  =  a + b.i + c.j + d.k            R01 thru R08 are to be initialized
•   R03 = c        •  R07 = c'        and    q' = a' + b'.i + c'.j + d'.k
•   R04 = d        •  R08 = d'
Flag: none
Subroutines:  "LNQ"  "E^Q"  "Q*"

Note:  When this program stops, registers R01 thru R04 contain q'
and registers R05 thru R08  ------- ln q

01  LBL "Q^Q"
02  RCL 04
03  RCL 03
04  RCL 02
05  RCL 01
06  XEQ "LNQ"
07  X<> 05
08  STO 01
09  RDN
10  X<> 06
11  STO 02
12  RDN
13  X<> 07
14  STO 03
15  RDN
16  X<> 08
17  STO 04
18  XEQ "Q*"
19  XEQ "E^Q"
20  END

( 43 bytes / SIZE 009 )

 STACK INPUTS OUTPUTS T / t Z / z Y / y X / x

with      qq'  =  x + y.i + z.j + t.k

Example:   q = 2 - 3i + 4j - 7k  and  q' = 1 - 4i + 2j + 5k       calculate  qq'

2   STO 01     1  STO 05              XEQ "Q^Q"    yields    -45.8455
-3  STO 02    -4  STO 06                                      RDN     68.7134
4  STO 03      2  STO 07                                      RDN       8.2012
-7  STO 04     5  STO 08                                      RDN    -39.0781
Thus    qq' = - 45.8455 +  68.7134 i + 8.2012 j  - 39.0781 k   ( rounded to 4D )

Solving a quaternion equation:       f ( q ) = q

-The following program uses an iterative method:
-The equation must be rewritten in the form:    f ( q ) = q
and if  f  satisfies a Lipschitz condition   | f(q) - f(q') | < m | q - q' |    ,    m<1   ,   provided q and q' are close to the solution,
then the sequence  qn+1 = f ( qn )  converges to a root.

Data Registers:

•  R00 = function name            R09 to R12 also contain   a ; b ; c ; d    when the subroutine is executed.
•  R01 = a                               R13 = function name too.
•  R02 = b
•  R03 = c
•  R04 = d                               with    q0  =  a + b.i + c.j + d.k  = guess      ( R00 thru R04 are to be initialized )

Flag:  none
Subroutine:  the function f  ( assuming  q = x + y.i + z.j + t.k  is registers R01 thru R04 upon entry )
>>>      f(q) = X + Y.i + Z.j + T.k  must end up in the stack registers X ; Y ; Z ; T
-Registers R00 thru R08 and R14 or greater can be used ( and altered ) to compute f(q)
-Registers R09 thru R13 can't be modified but can be used too.

01  LBL "SOLVE"
02  RCL 00
03  STO 13
04  LBL 01
05  VIEW 01                the real parts of the successive approximations are displayed
06  RCL 01
07  STO 09
08  RCL 02
09  STO 10
10  RCL 03
11  STO 11
12  RCL 04
13  STO 12
14  XEQ IND 13
15  STO 01
16  ST- 09
17  RDN
18  STO 02
19  ST- 10
20  RDN
21  STO 03
22  ST- 11
23  RDN
24  STO 04
25  ST- 12
26  RCL 09
27  X^2
28  RCL 10
29  X^2
30  RCL 11
31  X^2
32  RCL 12
33  X^2
34  +
35  +
36  +
37  X#0?                     to avoid a possible infinite loop, replace this line by the 2 instructions:    E-16        ( or another "small" number  )
38  GTO 01                                                                                                                           X<Y?
39  RCL 13
40  STO 00
41  RCL 04
42  RCL 03
43  RCL 02
44  RCL 01
45  END

( 62 bytes / SIZE 014 )

 STACK INPUTS OUTPUTS T / t Z / z Y / y X / x

with   q =  x + y.i + z.j + t.k  = one solution     ( x ; y ; z ; t are also in registers R01 to R04 )

Example:    Find a root of:    q3 + ( 1 + 2.i + 3.j + 4.k ) q + ( 2 - 3.i + 4.j - 7.k ) = 0
-This equation can be put in the form  q = f ( q ) in many ways but, unfortunately,  f doesn't necessarily satisfy the required Lipschitz condition!
-Let's try:      q = ( -1-2.i -3.j -4.k )-1  ( q3 + 2 - 3.i + 4.j - 7.k ) = f(q)      and let's key in the subroutine:

LBL "T"
3
STO 00
RCL 04
RCL 03
RCL 02
RCL 01
XEQ "Q^R"
STO 05
RDN
STO 06
RDN
STO 07
RDN
STO 08
2
ST+ 05
3
ST- 06
4
ST+ 07
7
ST- 08
4
CHS
3
CHS
2
CHS
1
CHS
XEQ "1/Q"
STO 01
RDN
STO 02
RDN
STO 03
RDN
STO 04
XEQ "Q*"
RTN

Then,   T  ASTO 00
CLX   STO 01  STO 02  STO 03  STO 04       ( if we choose q = 0 as initial value )
XEQ "SOLVE"

The successive approximations are displayed:
0.000000000
0.666666666
0.873580246
.....................
.....................
and eventually:             0.808540975
RDN  -1.290564599
RDN  -0.290931376
RDN   0.490521464       whence   q =  0.808540975 - 1.290564599 i  - 0.290931376 j + 0.490521464 k  is a root of this polynomial.

Notes:

-Convergence may be slow.
-If  f doesn't satisfy the required Lipschitz condition or if we choose a bad initial guess, the algorithm may be divergent.
-The last decimals may vary according to the first guess and/or the method for computing f(q).
For instance, the short version of "1/Q" is slightly less accurate than the long one.
-I've tried to generalize the secant method to solve directly  f(q) = 0  but the results are somewhat disappointing and convergence is very slow and chaotic.

Exercise:     Find a root of the equation:         k.ln q + q2 + i + j + k = 0
answer:       A solution is          q = 0.791739122  - 0.754317889 i - 0.231635888 j - 0.844420665 k

Rotations and Quaternions

-Quaternions can describe the spin of elementary particles, and they are useful in geometry too:
-A rotation ( in space ) can be defined by an angle A and a 3D-vector u(a;b;c)
-This rotation transforms a 3D-vector  V(x;y;z)  into  V'(x';y';z') and the following program computes  x' ; y' ; z'  from  x ; y ; z
-The formula is:
x'.i + y'.j + z'.k = q-1 ( x.i + y.j + z.k ) q     where  q = cos(A/2) -  sin(A/2) ( a.i + b.j + c.k ) / ( a2 + b2 + c2 )1/2

Data Registers:    • R00 = A      • R01 = a     • R02 = b    • R03 = c         are to be initialized before executing  "ROT"
R04 thru R11 are used for temporary data storage
Flag: none
Subroutines:   "Q*"

Notes:  -Registers R01-R02-R03  are altered during the calculations but they are restored at the end.
-The sign of the rotation angle A is determined by the right-hand rule.

01  LBL "ROT"
02  STO 06
03  RDN
04  STO 07
05  X<>Y
06  STO 08
07  CLX
08  STO 05
09  RCL 00
10  2
11  /
12  1
13  P-R
14  X<> 01
15  STO 09
16  X^2
17  LASTX
18  X<> 02
19  STO 10
20  X^2
21  LASTX
22  X<> 03
23  STO 04
24  STO 11
25  X^2
26  +
27  +
28  SQRT
29  /
30  ST* 02
31  ST* 03
32  ST* 04
33  XEQ "Q*"
34  X<> 01
35  STO 05
36  RDN
37  X<> 02
38  CHS
39  STO 06
40  RDN
41  X<> 03
42  CHS
43  STO 07
44  X<>Y
45  X<> 04
46  CHS
47  STO 08
48  XEQ "Q*"
49  CLX
50  RCL 09
51  STO 01
52  CLX
53  RCL 10
54  STO 02
55  CLX
56  RCL 11
57  STO 03
58  RDN
59  END

( 83 bytes / SIZE 012 )

 STACK INPUTS OUTPUTS Z z z' Y y y' X x x'

Example:     If r is the rotation defined by  A = 37°     and    u(1;2;4)      evaluate  V' = r(V)  where  V( 2 ; 3 ; 7 )

-Set your HP-41 in DEG mode and    37  STO 00     1  STO 01   2  STO 02   4  STO 03

-Then       7   ENTER^
3   ENTER^
2   XEQ "ROT"   produces    2.2051
RDN       3.2176
RDN       6.8399      and   V' = (  2.2051 ;  3.2176 ; 6.8399 )         ( rounded to 4D )