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

-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.