# hp41programs

Quaternions Quaternions for the HP-41

Overview

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 = R4
and 1 , i , j , k as the standard basis ( 1 , 0 , 0 , 0 ) ( 0 , 1 , 0 , 0 ) ( 0 , 0 , 1 , 0 ) ( 0 , 0 , 0 , 1 )
-Addition and subtraction are defined like ordinary vectors, but the multiplication is defined by the rules:

i.j = -j.i = k
j.k = -k.j = i       and     i2 = j2 = k2 = -1         Therefore, the product of 2 quaternions is not commutative. But it has the associative and distributive properties.
k.i =  -i.k = j

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

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 / µ2

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

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

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

Note:     Here is a shorter program to compute 1/q  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  µ = ( x2 + y2 + z2 + t2 )1/2  and 3 angles A , B , C such that:

x = µ cos A                                      A = the amplitude of the quaternion
y = µ sin A cos B                              B = the co-latitude 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  781/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 + q2/2! + q3/3! + ....... + qn/n! + ........
but the polar form provides a much faster way to compute eq

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     ;    eq =  x + y.i + z.j + t.k

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

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

N.B:

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

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

exp( x + y i + z j + t k ) = ex [ cos m + ( sin m ).( y i + z j + t k )/m ]     where  m = ( y2 + z2 + t2 )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 ( x2 + y2 + z2 + t2 )1/2 + Atan2(m,x).( y i + z j + t k )/m      where  m = ( y2 + z2 + t2 )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     ;    qr =  x + y.i + z.j + t.k

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

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

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

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

Notes:

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

b) Quaternion exponent

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

Data Registers:           R00:  unused

•   R01 = a        •  R05 = a'
•   R02 = b        •  R06 = b'       with    q  =  a + b.i + c.j + d.k            R01 thru R08 are to be initialized
•   R03 = c        •  R07 = c'        and    q' = a' + b'.i + c'.j + d'.k
•   R04 = d        •  R08 = d'
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      qq'  =  x + y.i + z.j + t.k

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

2   STO 01     1  STO 05          XEQ "Q^Q"        yields    -45.8455
-3  STO 02    -4  STO 06                                      RDN     68.7134
4  STO 03      2  STO 07                                      RDN       8.2012
-7  STO 04     5  STO 08                                      RDN    -39.0781

-Thus    qq' = - 45.8455 +  68.7134 i + 8.2012 j  - 39.0781 k   ( rounded to 4D )

Note:

-One can also employ the definition   qq' =   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    qq' = - 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 = ( y2 + z2 + t2 )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 = ( y2 + z2 + t2 )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 = ( e2q - 1 ) ( e2q + 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 + ( q2 + 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 + ( q2 - 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 = ( y2 + z2 + t2 )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 = ( y2 + z2 + t2 )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

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 "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 = ( y2 + z2 + t2 )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 = ( y2 + z2 + t2 )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) = An qn + An-1 qn-1 + ......... + A1 q + A0

where   Ak = ak + bk i + ck j + dk 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   = an       •  Rb+4  = an-1       ..................................     •  Re-3  = a0
•  Rb+1 = bn       •  Rb+5  = bn-1      ..................................     •  Re-2  = b0
•  Rb+2 = cn       •  Rb+6  = cn-1       .................................      •  Re-1  = c0
•  Rb+3 = dn       •  Rb+7  = dn-1       .................................      •  Ree   = d0            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 ) q3 + ( -1 + 2 i - 3 j + 4 k ) q2 + ( 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  qn+1 = f ( qn )  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    q0  =  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:    q3 + ( 1 + 2.i + 3.j + 4.k ) q + ( 2 - 3.i + 4.j - 7.k ) = 0

-This equation can be put in the form  q = f ( q ) in many ways but, unfortunately,  f doesn't necessarily satisfy the required Lipschitz condition !
-Let's try:      q = ( -1-2.i -3.j -4.k )-1  ( q3 + 2 - 3.i + 4.j - 7.k ) = f(q)      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 + q2 + 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  XY   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      q3 + ( 1 + 2 i + 3 j + 4 k ) q + ( 2 - 3 i + 4 j - 7 k ) = 0

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 ) / ( a2 + b2 + c2 )1/2

-The following program deals with a rotation r defined by an angle µ , a point A(xA,yA,zA) 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 = xA     • R02 = yA    • R03 = zA
• 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.