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 )
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall