hp41programs

Quaternions Quaternions for the HP-41

Overview

1°) Addition & Subtraction
2°) Multiplication
3°) Reciprocal
4°) Rectangular-Polar Conversion
5°) Polar-Rectangular Conversion
6°) Exponential
7°) Logarithm
8°) Raising a Quaternion to a Power

a) Real Exponent
b) Quaternion Exponent

9°) Hyperbolic Functions

a)   Sinh q
b)   Cosh q
c)   Tanh q

10°) Inverse Hyperbolic Functions

a)   ArcSinh q
b)   ArcCosh q
c)   ArcTanh q

11°) Trigonometric Functions

a)    Sin q
b)    Cos q
c)    Tan q

12°) Inverse Trigonometric Functions

a)    ArcSin q
b)    ArcCos q
c)    ArcTan q

13°) Quaternion Polynomials

14°) Quaternion Equations

a) Program#1 f(q) = q
b) Program#2 f(q) = 0

15°) Quaternions & Rotations

16°) 2 Programs

a) TANQ - Q*Q - POL - 1/Q - REC - E^Q - LNQ - Q^Q - SINQ - COSQ
b) = a) + ATANQ - ACOSQ - ASINQ - ATHQ - THQ - Q^R - ASHQ - ACHQ - SHQ - CHQ

Latest Update: Paragraph 16°)a) & 16°)b)

-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 H = R⁴
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     i² = j² = k² = -1         Therefore, the product of 2 quaternions is not commutative. But it has the associative and distributive properties.
   k.i = -i.k = j

-So, ( H , + , * ) is a non-commutative field.

>>>> Some of the following programs use Micro-code routines listed in "A few M-Code Routines for the HP-41"
but alternative "focal" programs are also given.

1°) Addition and Subtraction

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
Subroutines: /

01 LBL "Q+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

( 41 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+Q" yields    3                  and         SF 01   XEQ "Q+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

2°) 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'
Flags: /
Subroutines: /

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

01 LBL "Q*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

( 80 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*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 components of the imaginary part of q.q' are the 3 components 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*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

-You can save several bytes if your HP-41 has the M-code routines ST<>A CROSS Z*Z

01 LBL "Q*Q"
02 RCL 08
03 RCL 07
04 RCL 06
05 ST<>A
06 RCL 04
07 RCL 03
08 RCL 02
09 CROSS
10 ST<>A

11 RCL 06
12 RCL 05
13 RCL 02
14 RCL 01
15 Z*Z
16 X<>Y
17 ST+ M
18 CLX
19 RCL 03
20 RCL 07

21 *
22 -
23 RCL 04
24 RCL 08
25 *
26 -
27 RCL 01
28 RCL 07
29 *
30 RCL 03

31 RCL 05
32 *
33 +
34 ST+ N
35 CLX
36 RCL 01
37 RCL 08
38 *
39 RCL 04
40 RCL 05

41 *
42 +
43 ST+ O
44 CLX
45 ENTER^
46 ENTER^
47 ST<>A
48 R^
49 END

( 65 bytes / SIZE 009 )

3°) Reciprocal of a Quaternion

Data Registers: /
Flags: /
Subroutines: /

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

( 48 bytes / SIZE 000 )


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

with q = a + b.i + c.j + d.k ; 1/q = x + y.i + z.j + t.k and µ = ( a² + b² + c² + d² )^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 = µ² = 78 )

Note: Here is a shorter program to compute 1/q via 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 )

-The results may be, however, slightly less accurate

4°) Rectangular-Polar conversion

-A quaternion q = x + y.i + z.j + t.k can also be defined by its modulus µ = ( x² + y² + z² + t² )^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 of the quaternion
                       z = µ sin A sin B cos C                     C = the longitude of the quaternion
                       t = µ sin A sin B sin C

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

Data Registers: /
Flags: /
Subroutines: /

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 78^1/2 )
                                RDN    76.9115 = A
                                RDN 110.4104 = B
                                RDN -60.2551 = C

5°) Polar-Rectangular conversion

Data Registers: /
Flags: /
Subroutines: /

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°

XEQ "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 )

6°) Natural Exponential of a Quaternion

-The exponential of a quaternion q is defined by the series exp ( q ) = 1 + q + q²/2! + q³/3! + ....... + qⁿ/n! + ........
but the polar form provides a much faster way to compute e^q

Data Registers: /
Flags: /
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 ; e^q = x + y.i + z.j + t.k

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

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

N.B:

ln ( e^q ) is not always equal to q ( like ordinary complex numbers )
e^q.e^q' is not always equal to e^q+q'

-We can also use the following formula to compute exp(q)

exp( x + y i + z j + t k ) = e^x [ cos m + ( sin m ).( y i + z j + t k )/m ] where m = ( y² + z² + t² )^1/2

-Lines 06 to 19 may be replaced by RDN NORM LASTX X<>Y
where NORM is listed in "A few M-code routines for the HP-41"

01 LBL "E^Q"
02 RAD
03 E^X
04 STO M
05 STO N
06 X<> T
07 ENTER^

08 X^2
09 STO O
10 X<> T
11 ENTER^
12 X^2
13 ST+ O
14 X<> T

15 ENTER^
16 X^2
17 ST+ O
18 X<> O
19 SQRT
20 COS
21 ST* N

22 X<> L
23 X#0?
24 ST/ M
25 SIN
26 ST* M
27 X<> M
28 ST* T

29 ST* Z
30 *
31 RCL N
32 CLA
33 DEG
34 END

( 59 bytes / SIZE 000 )

-The results may be more accurate than those given by the first version.

7°) Natural Logarithm of a Quaternion

Data Registers: /
Flags: /
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 )

The variant hereunder uses the formula:

Ln( x + y i + z j + t k ) = Ln ( x² + y² + z² + t² )^1/2 + Atan2(m,x).( y i + z j + t k )/m where m = ( y² + z² + t² )^1/2

If m = 0 , ( y i + z j + t k )/m is replaced by i

-Lines 05 to 07 may be replaced by X^2 RCL Y X^2 + R^ X^2 + SQRT STO O RCL M

01 LBL "LNQ"
02 STO M
03 RDN
04 STO N
05 NORM
06 STO O
07 R^

08 ST* M
09 RAD
10 R-P
11 DEG
12 CLX
13 RCL O
14 X#0?

15 GTO 00
16 SIGN
17 STO N
18 LBL 00
19 /
20 RCL N
21 X<>Y

22 ST* T
23 ST* Z
24 *
25 RCL O
26 X^2
27 ST+ M
28 X<> M

29 SQRT
30 LN
31 CLA
32 END

( 54 bytes / SIZE 000 )

8°) Raising a Quaternion to a power

a) Real exponent

Data Register: • R00 = r is to be initialized before executing Q^R
Flags: /
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 ; q^r = x + y.i + z.j + t.k

Example1: q = 2 - 3i + 4j - 7k evaluate q^3.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   q^3.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 q^1/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     q^1/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 b² + c² + d² = 1 then ( b.i + c.j + d.k )² = -1
-This program gives ( -1 )^1/2 = i

b) Quaternion exponent

-This program calculates q^q' with the definition q^q' = e^q'.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'
Flags:   /
Subroutines: "LNQ" "E^Q" "Q*Q"

-When this program stops, registers R01 thru R04 contain q' and registers R05 thru R08 contain 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*Q"
19 XEQ "E^Q"
20 END

( 44 bytes / SIZE 009 )


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

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

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

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 q^q' = - 45.8455 + 68.7134 i + 8.2012 j - 39.0781 k ( rounded to 4D )

Note:

-One can also employ the definition q^q' = e^(lnq).q' which is used by the routine hereafter.

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

( 36 bytes / SIZE 009 )

-With the same example: q = 2 - 3i + 4j - 7k and q' = 1 - 4i + 2j + 5k

2   STO 01     1 STO 05          XEQ "Q^Q"        yields    -45.8455
-3 STO 02    -4 STO 06                                      RDN     18.3841
4 STO 03      2 STO 07                                      RDN    -55.4506
-7 STO 04     5 STO 08                                      RDN    -53.8808

-Thus q^q' = - 45.8455 + 18.3841 i - 55.4506 j - 53.8808 k ( rounded to 4D )

9°) Hyperbolic Functions

a) Sinh q

Formula: Sinh ( x + y i + z j + t k ) = Sinh x Cos m + ( Cosh x Sin m )( y i + z j + t k )/m where m = ( y² + z² + t² )^1/2

Data Registers: /
Flags: /
Subroutines: /

-Lines 03 to 11 may be replaced by

    E^X-1        X<> L     ST- M      2              ST/ N          X^2          ENTER^      X<> T         ST+ O
    STO M      CHS        ST+ N      ST+ N     X<> T         STO O     X^2              ENTER^     X<> O
    STO N       E^X-1     CLX         ST/ M      ENTER^     X<> T      ST+ O          X^2             SQRT

01 LBL "SHQ"
02 DEG
03 SINH
04 STO M
05 X<> L
06 COSH

07 STO N
08 RDN
09 NORM
10 LASTX
11 X<>Y
12 X#0?

13 ST/ N
14 R-D
15 COS
16 ST* M
17 X<> L
18 SIN

19 ST* N
20 X<> N
21 ST* T
22 ST* Z
23 *
24 RCL M

25 CLA
26 END

( 48 bytes / SIZE 000 )

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

With Sinh ( x + y i + z j + t k ) = x' + y' i + z' j + t' k

Example:

   5 ENTER^
   4 ENTER^
   3 ENTER^
   2 XEQ "SHQ" >>>>   2.5582                --- In 3 seconds ---
                            RDN   1.1315
                             RDN   1.5086
                             RDN   1.8858

Sinh ( 2 + 3 i + 4 j + 5 k ) = 2.5582 + 1.1315 i + 1.5086 j + 1.8858 k ( rounded to 4D )

b) Cosh q

Formula: Cosh ( x + y i + z j + t k ) = Cosh x Cos m + ( Sinh x Sin m )( y i + z j + t k )/m where m = ( y² + z² + t² )^1/2

Data Registers: /
Flags: /
Subroutines: /

-Lines 03 to 11 may be replaced by

   E^X-1        X<> L     ST+ M      2              ST/ N          X^2          ENTER^      X<> T         ST+ O
   STO M      CHS        ST- N       ST+ M     X<> T         STO O     X^2              ENTER^     X<> O
   STO N       E^X-1     CLX         ST/ M      ENTER^     X<> T      ST+ O          X^2             SQRT

01 LBL "CHQ"
02 DEG
03 COSH
04 STO M
05 X<> L
06 SINH

07 STO N
08 RDN
09 NORM
10 LASTX
11 X<>Y
12 X#0?

13 ST/ N
14 R-D
15 COS
16 ST* M
17 X<> L
18 SIN

19 ST* N
20 X<> N
21 ST* T
22 ST* Z
23 *
24 RCL M

25 CLA
26 END

( 48 bytes / SIZE 000 )

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

With Cosh ( x + y i + z j + t k ) = x' + y' i + z' j + t' k

Example:

   5 ENTER^
   4 ENTER^
   3 ENTER^
   2 XEQ "CHQ" >>>>   2.6537                --- In 3 seconds ---
                             RDN   1.0908
                             RDN   1.4543
                             RDN   1.8179

Cosh ( 2 + 3 i + 4 j + 5 k ) = 2.6537 + 1.0908 i + 1.4543 j + 1.8179 k ( rounded to 4D )

c) Tanh q

Tanh q is defined by Tanh q = ( Sinh q ) ( Cosh q )^-1

Data Registers: R00 unused R01 to R08: temp
Flags: /
Subroutines: "SHQ" "CHQ" "1/Q" "Q*Q"

01 LBL "THQ"
02 STO 01
03 R^
04 STO 04
05 R^
06 STO 03

07 R^
08 STO 02
09 R^
10 XEQ "SHQ"
11 X<> 01
12 R^

13 X<> 04
14 R^
15 X<> 03
16 R^
17 X<> 02
18 R^

19 XEQ "CHQ"
20 XEQ "1/Q"
21 STO 05
22 RDN
23 STO 06
24 RDN

25 STO 07
26 X<>Y
27 STO 08
28 XEQ "Q*Q"
29 END

( 57 bytes / SIZE 009 )

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

With Tanh ( x + y i + z j + t k ) = x' + y' i + z' j + t' k

Example:

   4 ENTER^
   3 ENTER^
   2 ENTER^
   1 XEQ "THQ" >>>>    1.0249                --- In 10 seconds ---
                             RDN   -0.1023
                             RDN   -0.1534
                             RDN   -0.2046

Tanh ( 1 + 2 i + 3 j + 4 k ) = 1.0249 - 0.1023 i - 0.1534 j - 0.2046 k ( rounded to 4D )

-The variant hereafter is faster ( 7 seconds instead of 10 ). It uses the formula Tanh q = ( e^2q - 1 ) ( e^2q + 1 )^-1

Data Registers: R00 unused R01 to R08: temp
Flags: /
Subroutines: "E^Q" "1/Q" "Q*Q"

01 LBL "THQ"
02 ST+ X
03 R^
04 ST+ X
05 R^
06 ST+ X
07 R^

08 ST+ X
09 R^
10 XEQ "E^Q"
11 STO 01
12 R^
13 STO 04
14 R^

15 STO 03
16 R^
17 STO 02
18 R^
19 SIGN
20 ST* X
21 ST- 01

22 ST+ L
23 X<> L
24 XEQ "1/Q"
25 STO 05
26 RDN
27 STO 06
28 RDN

29 STO 07
30 X<>Y
31 STO 08
32 XEQ "Q*Q"
33 END

( 61 bytes / SIZE 009 )

10°) Inverse Hyperbolic Functions

a) ArcSinh q

Formula: ArcSinh q = Ln [ q + ( q² + 1 )^1/2 ]

Data Registers: R00 = 2 and then 1/2
Flags: /
Subroutines: "Q^R" "LNQ"

-Lines 15-16 may be replaced by X<> 00 SIGN ST+ 00 X<> L

01 LBL "ASHQ"
02 STO 00
03 CLX
04 2
05 X<> 00
06 STO M
07 R^

08 STO P
09 R^
10 STO O
11 R^
12 STO N
13 R^
14 XEQ "Q^R"

15 X+1
16 X<> 00
17 1/X
18 X<> 00
19 XEQ "Q^R"
20 ST+ M
21 X<> T

22 RCL P
23 +
24 R^
25 RCL O
26 +
27 R^
28 RCL N

29 +
30 RCL M
31 XEQ "LNQ"
32 CLA
33 END

( 68 bytes / SIZE 001 )

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

With ArcSinh ( x + y i + z j + t k ) = x' + y' i + z' j + t' k

Example:

   5 ENTER^
   4 ENTER^
   3 ENTER^
   2 XEQ "ASHQ" >>>>   2.6837                --- In 14 seconds ---
                                RDN   0.5484
                                RDN   0.7313
                                RDN   0.9141

ArcSinh ( 2 + 3 i + 4 j + 5 k ) = 2.6837 + 0.5484 i + 0.7313 j + 0.9141 k ( rounded to 4D )

b) ArcCosh q

Formula: ArcCosh q = Ln [ q + ( q + 1 )^1/2 ( q - 1 )^1/2 ]

Data Registers: R00 to R12: temp
Flags: /
Subroutines: "Q^R" "Q*Q" "LNQ"

01 LBL "ACHQ"
02 STO 00
03 STO 05
04 STO 09
05 CLX
06 SIGN
07 ST+ 00
08 ST- 05
09 ST+ X
10 1/X
11 X<> 00

12 R^
13 STO 08
14 STO 12
15 R^
16 STO 07
17 STO 11
18 R^
19 STO 06
20 STO 10
21 R^
22 XEQ "Q^R"

23 STO 01
24 R^
25 STO 04
26 R^
27 STO 03
28 R^
29 STO 02
30 RCL 08
31 RCL 07
32 RCL 06
33 RCL 05

34 XEQ "Q^R"
35 STO 05
36 R^
37 STO 08
38 R^
39 STO 07
40 R^
41 STO 06
42 XEQ "Q*Q"
43 ST+ 09
44 X<> T

45 RCL 12
46 +
47 R^
48 RCL 11
49 +
50 R^
51 RCL 10
52 +
53 RCL 09
54 XEQ "LNQ"
55 END

( 86 bytes / SIZE 013 )

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

With ArcCosh ( x + y i + z j + t k ) = x' + y' i + z' j + t' k

Example:

   5 ENTER^
   4 ENTER^
   3 ENTER^
   2 XEQ "ACHQ" >>>>   2.6916                --- In 16 seconds ---
                                RDN   0.5505
                                RDN   0.7340
                                RDN   0.9175

ArcCosh ( 2 + 3 i + 4 j + 5 k ) = 2.6916 + 0.5505 i + 0.7340 j + 0.9175 k ( rounded to 4D )

Note:

-This program returns the principal value of arccosh, but other definitions are possible, for instance: arccosh q = Ln [ q + ( q² - 1 )^1/2 ]
-The following routine uses this formula.
-It's shorter and uses 1 data register instead of 13
-Lines 01 & 15 are the unique differences with "ASHQ"

-Lines 15-16 may be replaced by X<> 00 SIGN ST- 00 X<> L

01 LBL "ACHQ2"
02 STO 00
03 CLX
04 2
05 X<> 00
06 STO M
07 R^

08 STO P
09 R^
10 STO O
11 R^
12 STO N
13 R^
14 XEQ "Q^R"

15 X-1
16 X<> 00
17 1/X
18 X<> 00
19 XEQ "Q^R"
20 ST+ M
21 X<> T

22 RCL P
23 +
24 R^
25 RCL O
26 +
27 R^
28 RCL N

29 +
30 RCL M
31 XEQ "LNQ"
32 CLA
33 END

( 69 bytes / SIZE 001 )

-The sign of the results will be different sometimes.

c) ArcTanh q

Formula: ArcTanh q = (1/2) [ Ln ( 1 + q ) - Ln ( 1 - q ) ]

Data Registers: /
Flags: /
Subroutines: "LNQ" ( first version )

-Lines 20 to 23 may be replaced by R^ X<> P R^ X<> O R^ X<> N R^ X<> M

01 LBL "ATHQ"
02 R^
03 STO P
04 CHS
05 R^
06 STO O
07 CHS
08 R^
09 STO N

10 CHS
11 R^
12 STO M
13 CHS
14 SIGN
15 ST*X
16 ST+ L
17 ST+ M
18 X<> L

19 XEQ "LNQ"
20 ST<>A
21 R^
22 X<> P
23 RDN
24 XEQ "LNQ"
25 ST- M
26 X<> T
27 RCL P

28 -
29 R^
30 RCL O
31 -
32 R^
33 RCL N
34 -
35 2
36 ST/ T

37 ST/ Z
38 ST/ M
39 /
40 RCL M
41 CHS
42 CLA
43 END

( 79 bytes / SIZE 000 )

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

With ArcTanh ( x + y i + z j + t k ) = x' + y' i + z' j + t' k

Example:

   4 ENTER^
   3 ENTER^
   2 ENTER^
   1 XEQ "ATHQ" >>>>   0.0323                --- In 7 seconds ---
                                RDN   0.5173
                                RDN   0.7760
                                RDN   1.0347

ArcTanh ( 1 + 2 i + 3 j + 4 k ) = 0.0323 + 0.5173 i + 0.7760 j + 1.0347 k ( rounded to 4D )

11°) Trigonometric Functions

a) Sin q

Formula: Sin ( x + y i + z j + t k ) = Sin x Cosh m + ( Cos x Sinh m )( y i + z j + t k )/m where m = ( y² + z² + t² )^1/2

Data Registers: /
Flags: /
Subroutines: /

-Lines 09 to 12 may be replaced by

     X<> T         STO O       X^2          ENTER^      X<> O
     ENTER^     X<> T        ST+ O      X^2              SQRT
     X^2             ENTER^    X<> T      ST+ O

-Lines 15 to 18 may be replaced by

     E^X-1           CHS           CLX           ST/ O           X<> P
     STO O          E^X-1        SIGN         ST/ P
     STO P          ST+ O        ST+ X        X<> O
     X<> L          ST- P          ST+ O        ST* M

01 LBL "SINQ"
02 DEG
03 R-D
04 SIN
05 STO M
06 X<> L

07 COS
08 STO N
09 RDN
10 NORM
11 LASTX
12 X<>Y

13 X#0?
14 ST/ N
15 COSH
16 ST* M
17 X<> L
18 SINH

19 ST* N
20 X<> N
21 ST* T
22 ST* Z
23 *
24 RCL M

25 CLA
26 END

( 49 bytes / SIZE 000 )

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

With Sin ( x + y i + z j + t k ) = x' + y' i + z' j + t' k

Example:

   1.2 ENTER^
    1    ENTER^
   0.8 ENTER^
   0.6 XEQ "SINQ"   >>>>   1.6816              --- Execution time = 3 seconds ---
                                  RDN    1.0554
                                  RDN    1.3192
                                  RDN    1.5831

Sin ( 0.6 + 0.8 i + j + 1.2 k ) = 1.6816 + 1.0554 i + 1.3192 j + 1.5831 k ( rounded to 4D )

b) Cos q

Formula: Cos ( x + y i + z j + t k ) = Cos x Cosh m - ( Sin x Sinh m )( y i + z j + t k )/m where m = ( y² + z² + t² )^1/2

Data Registers: /
Flags: /
Subroutines: /

-Lines 10 to 13 may be replaced by

     X<> T         STO O       X^2          ENTER^      X<> O
     ENTER^     X<> T        ST+ O      X^2              SQRT
     X^2             ENTER^    X<> T      ST+ O

-Lines 16 to 19 may be replaced by

01 LBL "COSQ"
02 DEG
03 R-D
04 COS
05 STO M
06 X<> L

07 SIN
08 CHS
09 STO N
10 RDN
11 NORM
12 LASTX

13 X<>Y
14 X#0?
15 ST/ N
16 COSH
17 ST* M
18 X<> L

19 SINH
20 ST* N
21 X<> N
22 ST* T
23 ST* Z
24 *

25 RCL M
26 CLA
27 END

( 50 bytes / SIZE 000 )

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

With Cos ( x + y i + z j + t k ) = x' + y' i + z' j + t' k

Example:

   1.2 ENTER^
    1    ENTER^
   0.8 ENTER^
   0.6 XEQ "COSQ" >>>>    2.4580                  --- Execution time = 3 seconds ---
                                  RDN    -0.7220
                                  RDN    -0.9025
                                  RDN    -1.0831

Cos ( 0.6 + 0.8 i + j + 1.2 k ) = 2.4580 - 0.7220 i - 0.9025 j - 1.0831 k ( rounded to 4D )

c) Tan q

Formula: Tan q = ( Sin q ) ( Cos q )^-1

Data Registers: /
Flags: /
Subroutines: "SINQ" "COSQ" "1/Q" "Q*Q"

01 LBL "TANQ"
02 STO 01
03 R^
04 STO 04
05 R^
06 STO 03

07 R^
08 STO 02
09 R^
10 XEQ "SINQ"
11 X<> 01
12 R^

13 X<> 04
14 R^
15 X<> 03
16 R^
17 X<> 02
18 R^

19 XEQ "COSQ"
20 XEQ "1/Q"
21 STO 05
22 RDN
23 STO 06
24 RDN

25 STO 07
26 X<>Y
27 STO 08
28 XEQ "Q*Q"
29 END

( 60 bytes / SIZE 009 )

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

With Tan ( x + y i + z j + t k ) = x' + y' i + z' j + t' k

Example:

    1    ENTER^
    2    ENTER^
    3    ENTER^
    4    XEQ "TANQ"   >>>     0.001113                     --- Execution time = 10 seconds ---
                                  RDN     0.8019
                                  RDN     0.5346
                                  RDN     0.2673

Tan ( 4 + 3 i + 2 j + k ) = 0.001113 + 0.8019 i + 0.5346 j + 0.2673 k

12°) Inverse Trigonometric Functions

a) ArcSin q

Formula: if q = x + y i + z j + t k
ArcSin q = - [ ( y i + z j + t k ) / m ] ArcSinh [ q ( y i + z j + t k ) / m ] where m = ( y² + z² + t² )^1/2

Data Registers: R00 to R08: temp
Flags: /
Subroutines: "Q*Q" "ASHQ"

-Lines 09 to 12 may be replaced by R-P R^ STO 02 STO 06 R-P

01 LBL "ASINQ"
02 STO 01
03 R^
04 STO 04
05 STO 08
06 R^
07 STO 03
08 STO 07

09 R^
10 STO 02
11 STO 06
12 NORM
13 X#0?
14 GTO 00
15 SIGN
16 STO 06

17 LBL 00
18 ST/ 06
19 ST/ 07
20 ST/ 08
21 CLX
22 STO 05
23 XEQ "Q*Q"
24 XEQ "ASHQ"

25 X<> 05
26 STO 01
27 RDN
28 X<> 06
29 CHS
30 STO 02
31 RDN
32 X<> 07

33 CHS
34 STO 03
35 X<>Y
36 X<> 08
37 CHS
38 STO 04
39 XEQ "Q*Q"
40 END

( 72 bytes / SIZE 009 )

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

With ArcSin ( x + y i + z j + t k ) = x' + y' i + z' j + t' k

Example:

   5 ENTER^
   4 ENTER^
   3 ENTER^
   2 XEQ "ASINQ" >>>>   0.2732                --- In 19 seconds ---
                                 RDN   1.1419
                                 RDN   1.5226
                                 RDN   1.9032

ArcSin ( 2 + 3 i + 4 j + 5 k ) = 0.2732 + 1.1419 i + 1.5226 j + 1.9032 k ( rounded to 4D )

b) ArcCos q

Formula: ArcCos q = PI/2 - ArcSin q

Data Registers: R00 to R08: temp
Flags: /
Subroutine: "ASINQ"

-Line 06 is an M-code function ( not X^2 )
-Lines 05-06 may be replaced by RAD SIGN ASIN DEG

01 LBL "ACOSQ"
02 XEQ "ASINQ"
03 STO 00
04 CLX
05 PI
06 X/2
07 ST- 00
08 R^
09 CHS
10 R^
11 CHS
12 R^
13 CHS
14 RCL 00
15 CHS
16 END

( 34 bytes / SIZE 009 )

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

With ArcCos ( x + y i + z j + t k ) = x' + y' i + z' j + t' k

Example:

   5 ENTER^
   4 ENTER^
   3 ENTER^
   2 XEQ "ACOSQ" >>>>   1.2976                --- In 19 seconds ---
                                  RDN   -1.1419
                                  RDN   -1.5226
                                  RDN   -1.9032

ArcCos ( 2 + 3 i + 4 j + 5 k ) = 1.2976 - 1.1419 i - 1.5226 j - 1.9032 k ( rounded to 4D )

c) ArcTan q

Formula: if q = x + y i + z j + t k
ArcTan q = - [ ( y i + z j + t k ) / m ] ArcTanh [ q ( y i + z j + t k ) / m ] where m = ( y² + z² + t² )^1/2

Data Registers: R00 to R08: temp
Flags: /
Subroutines: "Q*Q" "ATHQ"

-Lines 09 to 12 may be replaced by R-P R^ STO 02 STO 06 R-P

01 LBL "ATANQ"
02 STO 01
03 R^
04 STO 04
05 STO 08
06 R^
07 STO 03
08 STO 07

09 R^
10 STO 02
11 STO 06
12 NORM
13 X#0?
14 GTO 00
15 SIGN
16 STO 06

17 LBL 00
18 ST/ 06
19 ST/ 07
20 ST/ 08
21 CLX
22 STO 05
23 XEQ "Q*Q"
24 XEQ "ATHQ"

25 X<> 05
26 STO 01
27 RDN
28 X<> 06
29 CHS
30 STO 02
31 RDN
32 X<> 07

33 CHS
34 STO 03
35 X<>Y
36 X<> 08
37 CHS
38 STO 04
39 XEQ "Q*Q"
40 END

( 72 bytes / SIZE 009 )

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

With ArcTan ( x + y i + z j + t k ) = x' + y' i + z' j + t' k

Example:

   5 ENTER^
   4 ENTER^
   3 ENTER^
   2 XEQ "ATANQ" >>>> 1.5331                --- In 14 seconds ---
                                  RDN   0.0558
                                  RDN   0.0744
                                  RDN   0.0930

ArcTan ( 2 + 3 i + 4 j + 5 k ) = 1.5331 + 0.0558 i + 0.0744 j + 0.0930 k ( rounded to 4D )

13°) Quaternion Polynomials

-The routine hereunder evaluates: P(q) = A_n qⁿ + A_n-1 q^n-1 + ......... + A₁ q + A₀

where A_k = a_k + b_k i + c_k j + d_k t and q = x + y i + z j + t k are quaternions.

Data Registers: R00 thru R09: temp ( Registers Rbb thru Ree & Rbb' thru Ree' are to be initialized before executing "PLQ" )

                                      • Rbb   = a_n       • Rb+4 = a_n-1       ..................................     • Re-3 = a₀
                                      • Rb+1 = b_n       • Rb+5 = b_n-1      ..................................     • Re-2 = b₀
                                      • Rb+2 = c_n       • Rb+6 = c_n-1       .................................      • Re-1 = c₀
                                      • Rb+3 = d_n       • Rb+7 = d_n-1       .................................      • Ree   = d₀            bb and bb' > 09

• Rbb' = x • Rbb'+1 = y • Rbb'+2 = z • Ree' = t ( ee' = bb' + 3 )

>>>> When the program stops, P(q) is in registers R01 R02 R03 R04 and q is also in registers R05 R06 R07 R08

Flags: /
Subroutine: "Q*Q"

-Lines 02 to 06 may be replaced by X<>Y STO 00 X<>Y STO 09 5 XEQ "LCO" ( cf "Miscellaneous Short Routines" )

01 LBL "PLQ"
02 STO 09
03 5
04 LCO
05 X<>Y
06 STO 00
07 CLX
08 STO 01

09 STO 02
10 STO 03
11 STO 04
12 LBL 01
13 XEQ "Q*Q"
14 STO 01
15 CLX
16 RCL IND 00

17 ST+ 01
18 ISG 00
19 CLX
20 RCL IND 00
21 +
22 STO 02
23 ISG 00
24 CLX

25 RCL IND 00
26 +
27 STO 03
28 ISG 00
29 CLX
30 RCL IND 00
31 +
32 STO 04

33 ISG 00
34 GTO 01
35 RCL 03
36 RCL 02
37 RCL 01
38 END

( 61 bytes )

STACK	INPUTS	OUTPUTS
T	/	t'
Z	/	z'
Y	bbb.eee	y'
X	bbb'.eee'	x'

Where bbb.eee = control number of the polynomial P All bbb's > 009 ( or bbb' = 005 )
bbb'.eee' = control number of the quaternion q = x + y i + z j + t k

and P ( x + y i + z j + t k ) = x' + y' i + z' j + t' k

Example: P(q) = ( 2 + 3 i + 4 j + 5 k ) q³ + ( -1 + 2 i - 3 j + 4 k ) q² + ( i + j + k ) q + 4 - 3 i - 5 j - 7 k Evaluate P( 6 - i - 2 j - 3 k )

If we store [ 2 3 4 5 -1 ..................... 4 -3 -5 -7 ] into [ R10 R11 R12 .................. R25 ]
and [ 6 -1 -2 -3 ] into [ R31 R32 R33 R34 ]

10.025 ENTER^
31.034 XEQ "PLQ" >>>>   2456                                  --- Execution time = 10s ---
                                    RDN    -222
                                    RDN    -159
                                    RDN    -894

-Thus, P( 6 - i - 2 j - 3 k ) = 2456 - 222 i -159 j -894 k

Note:

-The program may be simplified if you always use registers R05-R06-R07-R08 to store q
-In this case, delete lines 02 to 05 and place bbb.eee in X-register.

14°) Quaternion equation

a) Program#1 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') | < h | q - q' | with h < 1 , provided q and q' are close to the solution,
then the sequence q_n+1 = f ( q_n ) converges to a root.

Data Registers: • R00 = function name ( Registers R00 thru R04 are to be initialized before executing "SOLVE" )

      • R01 = a                        R09 to R12 also contain   a ; b ; c ; d    when the subroutine is executed.
      • R02 = b                        R13 = function name too.
      • R03 = c
      • R04 = d                        with    q₀ = a + b.i + c.j + d.k = guess

Flags: /
Subroutine: the function f ( assuming q = x + y.i + z.j + t.k is in registers R01 thru R04 , R09 thru R12 and in X- Y- Z- T-registers 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.

-Line 05 , the real parts of the successive approximations are displayed.
-To avoid a possible infinite loop, replace line 37 by the 2 instructions:

E-8 ( or another "small" number )
X<Y?

-Or you can add RND after line 36
and the accuracy will be controlled by the display setting.

01 LBL "SOLVE"
02 RCL 00
03 STO 13
04 LBL 01
05 VIEW 01
06 RCL 04
07 STO 12
08 RCL 03
09 STO 11

10 RCL 02
11 STO 10
12 RCL 01
13 STO 09
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 X<>Y
24 STO 04
25 ST- 12
26 RCL 09
27 ABS

28 RCL 10
29 ABS
30 RCL 11
31 ABS
32 RCL 12
33 ABS
34 +
35 +
36 +

37 X#0?
38 GTO 01
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: q³ + ( 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 ( q³ + 2 - 3.i + 4.j - 7.k ) = f(q) we load the subroutine:

   LBL "T"
   STO 00
   CLX
   3
   X<> 00
   XEQ "Q^R"
   STO 05
   RDN
   STO 06
   RDN
   STO 07
   X<>Y
   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
   X<>Y
   STO 04
   XEQ "Q*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                                      ---Execution time = 3mn34s---
                        RDN -1.290564600
                        RDN -0.290931374
                        RDN   0.490521464       whence   q = 0.808540975 - 1.290564600 i - 0.290931374 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 equation g(q) = 0 may be rewritten f(q) = q with f(q) = q - µ g(q) where µ is constant. This may sometimes lead to convergence.
-The last decimals may vary according to the first guess and/or the method for computing f(q).

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

b) Program#2 f(q) = 0

-The program listed hereunder uses the secant method to solve directly f(q) = 0
-I'm not sure it's really sound theoretically...

-It requires 2 initial guesses q = a + b i + c j + d k and q' = a' + b' i + c' j + d' k which are to be stored in registers R01 thru R08

Data Registers: • R00 = function name ( Registers R00 thru R08 are to be initialized before executing "SOLVE2" )

           • R01 = a          • R05 = a'
           • R02 = b          • R06 = b'                             q is also stored in R14 to R17
           • R03 = c          • R07 = c'                             q' is also stored in R18 to R21             R13 to R25: temp
           • R04 = d          • R08 = d'

Flags: /
Subroutine:

-The function f must return f(q) = x' + y' i + z' j + t' k in the stack registers X , Y , Z , T
assuming q = x + y i + z j + t k is registers X , Y , Z , T and R01 thru R04 upon entry.

-Registers R00 thru R12 and R26 or greater can be used ( and altered ) to compute f(q) but without changing registers R13 to R25

-Lines 02-03 may be replaced by .008 ENTER^ 13 XEQ "LCO" ( cf "Miscellaneous Short Routines" )
-Lines 92-93 may be replaced by 14.017 ENTER^ 18 XEQ "LCO"
-Line 114 is a three-byte GTO 01
-Lines 116 to 128 are only useful to return | f(q) | in L-register

01 LBL "SOLVE2"
02 .013009
03 REGMOVE
04 RCL 08
05 STO 04
06 RCL 07
07 STO 03
08 RCL 06
09 STO 02
10 RCL 05
11 STO 01
12 XEQ IND 00
13 STO 22
14 RDN
15 STO 23
16 RDN
17 STO 24
18 X<>Y
19 STO 25
20 LBL 01
21 RCL 17
22 STO 04
23 RCL 16
24 STO 03
25 RCL 15
26 STO 02
27 RCL 14
28 STO 01

29 VIEW 01
30 XEQ IND 13
31 STO 09
32 RDN
33 STO 10
34 RDN
35 STO 11
36 X<>Y
37 STO 12
38 RCL 18
39 RCL 14
40 -
41 STO 01
42 RCL 19
43 RCL 15
44 -
45 STO 02
46 RCL 20
47 RCL 16
48 -
49 STO 03
50 RCL 21
51 RCL 17
52 -
53 STO 04
54 E-16
55 RCL 09
56 RCL 22

57 -
58 STO 05
59 X^2
60 RCL 23
61 RCL 10
62 -
63 STO 06
64 X^2
65 +
66 RCL 24
67 RCL 11
68 -
69 STO 07
70 X^2
71 +
72 RCL 25
73 RCL 12
74 -
75 STO 08
76 X^2
77 +
78 X<Y?
79 GTO 02
80 ST/ 05
81 ST/ 06
82 ST/ 07
83 ST/ 08
84 XEQ "Q*Q"

85 STO 01
86 RDN
87 STO 02
88 RDN
89 STO 03
90 X<>Y
91 STO 04
92 14.018004
93 REGMOVE
94 RCL 09
95 STO 05
96 STO 22
97 RCL 10
98 STO 06
99 STO 23
100 RCL 11
101 STO 07
102 STO 24
103 RCL 12
104 STO 08
105 STO 25
106 XEQ "Q*Q"
107 ST+ 14
108 RDN
109 ST+ 15
110 RDN
111 ST+ 16
112 X<>Y

113 ST+ 17
114 GTO 01
115 LBL 02
116 RCL 09
117 X^2
118 RCL 10
119 X^2
120 RCL 11
121 X^2
122 RCL 12
123 X^2
124 +
125 +
126 +
127 SQRT
128 SIGN
129 RCL 13
130 STO 00
131 RCL 17
132 STO 04
133 RCL 16
134 STO 03
135 RCL 15
136 STO 02
137 RCL 14
138 STO 01
139 END

( 213 bytes / SIZE 026 )

STACK	INPUTS	OUTPUTS
T	/	t
Z	/	z
Y	/	y
X	/	x
L	/	\| f(q) \|

With q = x + y i + z j + t k is a solution of f(q) = 0 provided | f(q) | in L-register is "small"

Example: find a solution to q³ + ( 1 + 2 i + 3 j + 4 k ) q + ( 2 - 3 i + 4 j - 7 k ) = 0

-Load the following subroutine:

LBL "T"       XEQ "Q^R"       STO 11           RCL 02          STO 08           3                       ST+ 09          CLX          ST+ 09            RCL 10
STO 00       STO 09              X<>Y             STO 06          1                      STO 03            RDN              7                RCL 12           3
CLX            RDN                  STO 12           RCL 03          STO 01           4                       ST+ 10          -                RCL 11            -
3                 STO 10              RCL 01           STO 07           2                     STO 04            RDN              ST+12       4                      RCL 09
X<> 00       RDN                  STO 05           RCL 04           STO 02          XEQ"Q*Q"       ST+ 11          2                +                     RTN

"T" ASTO 00

0 STO 01 STO 02 STO 03 STO 04 if you choose 0
0.1 STO 05 STO 06 STO 07 STO 08 and 0.1 + 0.1 i + 0.1 j + 0.1 k as initial guesses.

XEQ "SOLVE2" >>>> yields after ( too ) many iterations: ---Execution time = 25mn10s---

0.808540975 RDN -1.290564601 RDN -0.290931376 RDN 0.490521463 LASTX gives | f(q) | = 0.000000005

-So, a solution is q = 0.808540975 - 1.290564601 i - 0.290931376 j + 0.490521463 k

Note:

-In fact, when the subroutine is called, q is also in registers R14-R15-R16-R17 except at line 12.
-Add 14.018004 REGSWAP after line 19 and after line 03 ( this is not absolutely necessary! )

-As you can see, the convergence is very slow and irregular.
-If you know a better method...

15°) Rotations and Quaternions

-Quaternions can describe the spin of elementary particles, and they are useful in geometry too:

-A vectorial rotation ( in space ) can be defined by an angle µ and a 3D-vector u(a;b;c)
-This rotation transforms a 3D-vector V(x;y;z) into V'(x';y';z')

Formula: x'.i + y'.j + z'.k = q^-1 ( x.i + y.j + z.k ) q where q = cos(µ/2) - sin(µ/2) ( a.i + b.j + c.k ) / ( a² + b² + c² )^1/2

-The following program deals with a rotation r defined by an angle µ , a point A(x_A,y_A,z_A) and a vector u(a,b,c)

-It takes a point M(x,y,z) and returns M'(x',y',z') where M' = r(M)

Data Registers: • R00 = µ ( Registers R00 thru R06 are to be initialized before executing "ROT" )

• R01 = x_A • R02 = y_A • R03 = z_A
• R04 = a • R05 = b • R06 = c

R07 to R15 are used for temporary data storage
Flags: /
Subroutine: "Q*Q"

01 LBL "ROT"
02 1.009006
03 STO 15
04 REGMOVE
05 CLX
06 RCL 01
07 -
08 STO 06
09 CLX
10 RCL 02
11 -
12 STO 07
13 CLX
14 RCL 03

15 -
16 STO 08
17 CLX
18 STO 05
19 RCL 00
20 2
21 /
22 1
23 P-R
24 STO 01
25 CLX
26 RCL 12
27 STO 02
28 X^2

29 RCL 13
30 STO 03
31 X^2
32 RCL 14
33 STO 04
34 X^2
35 +
36 +
37 SQRT
38 /
39 ST* 02
40 ST* 03
41 ST* 04
42 XEQ "Q*Q"

43 X<> 01
44 STO 05
45 RDN
46 X<> 02
47 CHS
48 STO 06
49 RDN
50 X<> 03
51 CHS
52 STO 07
53 X<>Y
54 X<> 04
55 CHS
56 STO 08

57 XEQ "Q*Q"
58 CLX
59 RCL 15
60 REGSWAP
61 CLX
62 RCL 03
63 ST+ T
64 CLX
65 RCL 02
66 ST+ Z
67 CLX
68 RCL 01
69 +
70 END

( 104 bytes / SIZE 016 )

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

Example: µ = 24° A(2,3,7) U(4,6,9) Find M' = r(M) if M(1,4,8)

24 STO 00

      2 STO 01      4 STO 04
      3 STO 02      6 STO 05
      7 STO 03      9 STO 06

-Set the HP-41 in DEG mode

     8 ENTER^
     4 ENTER^
     1 XEQ "ROT" >>>> 1.009250425   RDN   3.497956694   RDN   8.330584238         ---Execution time = 6s ---

-Whence M'( 1.009250425 , 3.497956694 , 8.330584238 )

Note:

-The sign of the rotation-angle µ is determined by the right-hand rule.

16°) 2 Programs

a) TANQ - Q*Q - POL - 1/Q - REC - E^Q - LNQ - Q^Q - SINQ - COSQ

-HP41 in mode DEG

01 LBL "TANQ"
02 STO 01
03 RDN
04 STO 02
05 RDN
06 STO 03
07 RDN
08 STO 04
09 RDN
10 XEQ 07
11 X<> 01
12 RDN
13 X<> 02
14 RDN
15 X<> 03
16 RDN
17 X<> 04
18 RDN
19 XEQ 08
20 XEQ 00
21 STO 05
22 RDN
23 STO 06
24 RDN
25 STO 07
26 X<>Y
27 STO 08
28 LBL "Q*Q"
29 LBL 01
30 RCL 01
31 RCL 05
32 *
33 RCL 02
34 RCL 06
35 *

36 -
37 RCL 03
38 RCL 07
39 *
40 -
41 RCL 04
42 RCL 08
43 *
44 -
45 STO 09
46 RCL 01
47 RCL 06
48 *
49 RCL 02
50 RCL 05
51 *
52 +
53 RCL 03
54 RCL 08
55 *
56 +
57 RCL 04
58 RCL 07
59 *
60 -
61 STO 10
62 RCL 01
63 RCL 08
64 *
65 RCL 04
66 RCL 05
67 *
68 +
69 RCL 02
70 RCL 07

71 *
72 +
73 RCL 03
74 RCL 06
75 *
76 -
77 RCL 01
78 RCL 07
79 *
80 RCL 03
81 RCL 05
82 *
83 +
84 RCL 02
85 RCL 08
86 *
87 -
88 RCL 04
89 RCL 06
90 *
91 +
92 RCL 10
93 RCL 09
94 RTN
95 LBL "POL"
96 LBL 02
97 RDN
98 RDN
99 R-P
100 R^
101 R-P
102 R^
103 R-P
104 RTN
105 LBL "1/Q"

106 LBL 00
107 R^
108 X^2
109 STO 00
110 X<> L
111 CHS
112 R^
113 X^2
114 ST+ 00
115 X<> L
116 CHS
117 R^
118 X^2
119 ST+ 00
120 X<> L
121 CHS
122 R^
123 X^2
124 ST+ 00
125 CLX
126 RCL 00
127 ST/ L
128 ST/ Y
129 ST/ Z
130 ST/ T
131 X<> L
132 RTN
133 LBL "REC"
134 LBL 03
135 P-R
136 RDN
137 P-R
138 RDN
139 P-R
140 RDN

141 RDN
142 RTN
143 LBL "E^Q"
144 LBL 04
145 RDN
146 RDN
147 R-P
148 R^
149 R-P
150 R-D
151 R^
152 E^X
153 GTO 03
154 LBL "LNQ"
155 LBL 05
156 XEQ 02
157 LN
158 RDN
159 D-R
160 P-R
161 RDN
162 P-R
163 RDN
164 RDN
165 RTN
166 LBL "Q^Q"
167 RCL 04
168 RCL 03
169 RCL 02
170 RCL 01
171 XEQ 05
172 X<> 05
173 STO 01
174 RDN

175 X<> 06
176 STO 02
177 RDN
178 X<> 07
179 STO 03
180 RDN
181 X<> 08
182 STO 04
183 XEQ 01
184 GTO 04
185 LBL "SINQ"
186 LBL 07
187 R-D
188 SIN
189 STO 09
190 X<> L
191 COS
192 STO 10
193 GTO 06
194 LBL "COSQ"
195 LBL 08
196 R-D
197 COS
198 STO 09
199 X<> L
200 SIN
201 CHS
202 STO 10
203 LBL 06
204 X<> T
205 ENTER
206 X^2
207 STO 11
208 X<> T

209 ENTER
210 X^2
211 ST+ 11
212 X<> T
213 ENTER
214 X^2
215 ST+ 11
216 X<> 11
217 SQRT
218 X#0?
219 ST/ 10
220 E^X-1
221 STO 11
222 STO 12
223 X<> L
224 CHS
225 E^X-1
226 ST+ 11
227 ST- 12
228 CLX
229 2
230 ST+ 11
231 ST/ 11
232 ST/ 12
233 X<> 11
234 ST* 09
235 X<> 12
236 ST* 10
237 X<> 10
238 ST* T
239 ST* Z
240 *
241 RCL 09
242 END

( 363 bytes / SIZE 013 )

b) = a) + ATANQ - ACOSQ - ASINQ - ATHQ - THQ - Q^R - ASHQ - ACHQ - SHQ - CHQ

-HP41 in mode DEG

01 LBL "ATANQ"
02 XEQ 10
03 XEQ 04
04 GTO 12
05 LBL "ACOSQ"
06 XEQ 14
07 STO 00
08 CLX
09 ACOS
10 D-R
11 ST- 00
12 R^
13 CHS
14 R^
15 CHS
16 R^
17 CHS
18 RCL 00
19 CHS
20 RTN
21 LBL "ASINQ"
22 LBL 14
23 XEQ 10
24 XEQ 13
25 LBL 12
26 X<> 05
27 STO 01
28 RDN
29 X<> 06
30 CHS
31 STO 02
32 RDN
33 X<> 07
34 CHS
35 STO 03
36 X<>Y
37 X<> 08
38 CHS
39 STO 04
40 XEQ 01
41 RTN
42 LBL 10
43 STO 01
44 R^
45 STO 04
46 STO 08
47 R^
48 STO 03
49 STO 07
50 R-P
51 R^
52 STO 02
53 STO 06
54 R-P
55 X#0?
56 GTO 14
57 SIGN
58 STO 06
59 LBL 14
60 ST/ 06
61 ST/ 07
62 ST/ 08
63 CLX
64 STO 05
65 XEQ 01
66 RTN
67 LBL "ATHQ"
68 LBL 04
69 R^
70 STO 12
71 CHS
72 R^
73 STO 11
74 CHS
75 R^

76 STO 10
77 CHS
78 R^
79 STO 09
80 CHS
81 SIGN
82 ST* X
83 ST+ L
84 ST+ 09
85 X<> L
86 XEQ 05
87 R^
88 X<> 12
89 R^
90 X<> 11
91 R^
92 X<> 10
93 R^
94 X<> 09
95 XEQ 05
96 ST- 09
97 X<> T
98 RCL 12
99 -
100 R^
101 RCL 11
102 -
103 R^
104 RCL 10
105 -
106 2
107 ST/ T
108 ST/ Z
109 ST/ 09
110 /
111 RCL 09
112 CHS
113 RTN
114 LBL "THQ"
115 STO 01
116 RDN
117 STO 02
118 RDN
119 STO 03
120 RDN
121 STO 04
122 RDN
123 XEQ 11
124 XEQ 06
125 XEQ 12
126 GTO 04
127 LBL "TANQ"
128 STO 01
129 RDN
130 STO 02
131 RDN
132 STO 03
133 RDN
134 STO 04
135 RDN
136 XEQ 07
137 XEQ 06
138 XEQ 08
139 LBL 04
140 XEQ 00
141 STO 05
142 RDN
143 STO 06
144 RDN
145 STO 07
146 X<>Y
147 STO 08
148 LBL "Q*Q"
149 LBL 01
150 RCL 01

151 RCL 05
152 *
153 RCL 02
154 RCL 06
155 *
156 -
157 RCL 03
158 RCL 07
159 *
160 -
161 RCL 04
162 RCL 08
163 *
164 -
165 STO 09
166 RCL 01
167 RCL 06
168 *
169 RCL 02
170 RCL 05
171 *
172 +
173 RCL 03
174 RCL 08
175 *
176 +
177 RCL 04
178 RCL 07
179 *
180 -
181 STO 10
182 RCL 01
183 RCL 08
184 *
185 RCL 04
186 RCL 05
187 *
188 +
189 RCL 02
190 RCL 07
191 *
192 +
193 RCL 03
194 RCL 06
195 *
196 -
197 RCL 01
198 RCL 07
199 *
200 RCL 03
201 RCL 05
202 *
203 +
204 RCL 02
205 RCL 08
206 *
207 -
208 RCL 04
209 RCL 06
210 *
211 +
212 RCL 10
213 RCL 09
214 RTN
215 LBL 06
216 X<> 01
217 RDN
218 X<> 02
219 RDN
220 X<> 03
221 RDN
222 X<> 04
223 RDN
224 RTN
225 LBL "POL"

226 LBL 02
227 RDN
228 RDN
229 R-P
230 R^
231 R-P
232 R^
233 R-P
234 RTN
235 LBL "1/Q"
236 LBL 00
237 R^
238 X^2
239 STO 00
240 X<> L
241 CHS
242 R^
243 X^2
244 ST+ 00
245 X<> L
246 CHS
247 R^
248 X^2
249 ST+ 00
250 X<> L
251 CHS
252 R^
253 X^2
254 ST+ 00
255 CLX
256 RCL 00
257 ST/ L
258 ST/ Y
259 ST/ Z
260 ST/ T
261 X<> L
262 RTN
263 LBL "Q^R"
264 LBL 09
265 XEQ 02
266 X<>Y
267 X<> 00
268 ST* 00
269 Y^X
270 LASTX
271 X<> 00
272 X<>Y
273 LBL "REC"
274 LBL 03
275 P-R
276 RDN
277 P-R
278 RDN
279 P-R
280 RDN
281 RDN
282 RTN
283 LBL "Q^Q"
284 RCL 04
285 RCL 03
286 RCL 02
287 RCL 01
288 XEQ 05
289 X<> 05
290 STO 01
291 RDN
292 X<> 06
293 STO 02
294 RDN
295 X<> 07
296 STO 03
297 RDN
298 X<> 08
299 STO 04
300 XEQ 01

301 LBL "E^Q"
302 RDN
303 RDN
304 R-P
305 R^
306 R-P
307 R-D
308 R^
309 E^X
310 GTO 03
311 LBL "ASHQ"
312 LBL 13
313 STO 00
314 CLX
315 2
316 X<> 00
317 STO 09
318 RDN
319 STO 10
320 RDN
321 STO 11
322 RDN
323 STO 12
324 RDN
325 XEQ 09
326 X<> 00
327 SIGN
328 ST+ 00
329 X<> L
330 1/X
331 X<> 00
332 XEQ 09
333 ST+ 09
334 X<> T
335 RCL 12
336 +
337 R^
338 RCL 11
339 +
340 R^
341 RCL 10
342 +
343 RCL 09
344 GTO 05
345 LBL "ACHQ"
346 STO 00
347 STO 05
348 STO 13
349 CLX
350 SIGN
351 ST+ 00
352 ST- 05
353 ST+ X
354 1/X
355 X<> 00
356 R^
357 STO 08
358 STO 12
359 R^
360 STO 07
361 STO 11
362 R^
363 STO 06
364 STO 14
365 R^
366 XEQ 09
367 STO 01
368 R^
369 STO 04
370 R^
371 STO 03
372 R^
373 STO 02
374 RCL 08
375 RCL 07

376 RCL 06
377 RCL 05
378 XEQ 09
379 STO 05
380 R^
381 STO 08
382 R^
383 STO 07
384 R^
385 STO 06
386 XEQ 01
387 ST+ 13
388 X<> T
389 RCL 12
390 +
391 R^
392 RCL 11
393 +
394 R^
395 RCL 14
396 +
397 RCL 13
398 LBL "LNQ"
399 LBL 05
400 XEQ 02
401 LN
402 RDN
403 D-R
404 P-R
405 RDN
406 P-R
407 RDN
408 RDN
409 RTN
410 LBL "SINQ"
411 LBL 07
412 R-D
413 SIN
414 STO 09
415 X<> L
416 COS
417 STO 10
418 GTO 06
419 LBL "COSQ"
420 LBL 08
421 R-D
422 COS
423 STO 09
424 X<> L
425 SIN
426 CHS
427 STO 10
428 LBL 06
429 X<> T
430 ENTER
431 X^2
432 STO 11
433 X<> T
434 ENTER
435 X^2
436 ST+ 11
437 X<> T
438 ENTER
439 X^2
440 ST+ 11
441 X<> 11
442 SQRT
443 X#0?
444 ST/ 10
445 E^X-1
446 STO 11
447 STO 12
448 X<> L
449 CHS
450 E^X-1

451 ST+ 11
452 ST- 12
453 CLX
454 2
455 ST+ 11
456 ST/ 11
457 ST/ 12
458 X<> 11
459 ST* 09
460 X<> 12
461 ST* 10
462 X<> 10
463 ST* T
464 ST* Z
465 *
466 RCL 09
467 RTN
468 LBL "SHQ"
469 LBL 11
470 E^X-1
471 STO 09
472 STO 10
473 X<> L
474 CHS
475 E^X-1
476 ST- 09
477 ST+ 10
478 CLX
479 2
480 ST+ 10
481 GTO 10
482 LBL "CHQ"
483 LBL 12
484 E^X-1
485 STO 09
486 STO 10
487 X<> L
488 CHS
489 E^X-1
490 ST+ 09
491 ST- 10
492 CLX
493 2
494 ST+ 09
495 LBL 10
496 ST/ 09
497 ST/ 10
498 X<> T
499 ENTER
500 X^2
501 STO 11
502 X<> T
503 ENTER
504 X^2
505 ST+ 11
506 X<> T
507 ENTER
508 X^2
509 ST+ 11
510 X<> 11
511 SQRT
512 X#0?
513 ST/ 10
514 R-D
515 COS
516 ST* 09
517 X<> L
518 SIN
519 ST* 10
520 X<> 10
521 ST* T
522 ST* Z
523 *
524 RCL 09
525 END

( 817 bytes / SIZE 015 )

-We save several bytes !