QUATERNION  Manual


OVERVIEW
 

-The following programs compute elementary functions of quaternions.
-They also generalize several special functions to quaternionic variables, but the order and/or degrees must remain real.
-Identical series expansions are employed for real and hypercomplex variables.
-Though the product of 2 quaternions is not commutative, the results are unambiguous.
-But it would not be the case if the orders or degrees were complex or hypercomplex...
 

Usage:

-A quaternion variable  x + y i + z j + t k  is to be stored in the stack registers X Y Z T  ( t ENTER^ z ENTER^ y ENTER^ x  XEQ "...." )
-And the resulting quaternion  x' + y' i + z' j + t' k  is logically returned in  X-Y-Z-T-registers too !

-If there is a unique order or degree, store it into R00.
-If there are several orders and/or degrees, store them into R09 R10 .....
 ( the idea is that R01 to R08 are always used to perform products of quaternions )

-These registers are sometimes altered during the calculations, but their contents are restored before the end.

-This "rule" about degree/order has 1 exception: the generalized hypergeometric function "HGFQ+"

 >>>  In some cases, the functions depends on 2 or 3 quaternions ( Q*Q , Q^Q , Sylvester equation )
          They are to be stored into R01 to R04 and R05 to R08 ( and R09 to R12 if 3 quaternions are involved )
-Take care that the product is not commutative.

-The French names are used for the hyperbolic functions:   "SH"  "CH"  "TH"  "ASH"  "ACH"  "ATH"

Notes:

-Routines names between "..." are focal programs.
-The other ones are M-Code routines.
 
 
XROM  Function  Desciption
 15,00
 15,01
 15,02
 15,03
 15,04
 15,05
 15,06
 15,07
 15,08
 15,09
 15,10
 15,11
 15,12
 15,13
 15,14
 15,15
 15,16
 15,17
 15,18
 15,19
 15,20
 15,21
 15,22
 15,23
 15,24
 15,25
 15,26
 15,27
 15,28
 15,29
 15,30
 15,31
 15,32
 15,33
 15,34
 15,35
 15,36
 15,37
 15,38
 15,39
 15,40
 15,41
 15,42
 15,43
 15,44
 15,45
 15,46
 15,47
 15,48
 15,49
 15,50
 15,51
 15,52
 15,53
 15,54
 15,55
 15,56
 15,57
 15,58
 15,59
 15,60
 15,61
 15,62
 -QUATERNION
 1/Q
 ABSQ
 ASTOQ
 CHSQ
 NRM2Q
 Q^2
 Q*Q
 Q<>
 Q<>A
 QCHS
 QRCL
 QSTO
 QSTOA
 QSTX
 QSTX*
 QSTX/
 QSUMAB
 QVIEW
 QZERO
 SQRTQ
 SUMABQ
 VIEWQ
 -ELEM FNS
 "Q^Q"
 "Q^R"
 "E^Q"
 "LNQ"
 "SHQ"
 "CHQ"
 "THQ"
 "ASHQ"
 "ACHQ"
 "ATHQ"
 "SINQ"
 "COSQ"
 "TANQ"
 "ASINQ"
 "ACOSQ"
 "ATANQ"
 -EQUATION
 "3DROT"
 "FNQ=Q"
 "FQ=0"
 "FQ=Q"
 "PEVALQ"
 "POL"
 "QLS"
 "REC"
 "SYLV"
 -ORTHOPOL
  1/GM
  DSQ
  GEU
  QSUM
  SUMQ
 "CHBQ"
 "HMTQ"
 "JCPQ"
 "LAGQ"
 "LEGQ"
 "PMNQ"
 "USPQ"
 Section Header
 Inverse of a Quaternion
 Euclidean Norm
 Stores M N O P into X Y Z T
 Changes the sign of the quaternion in X Y Z T
 Square of the Euclidean Norm
 Square of a Quaternion
 Product of 2 Quaternions in R01-R04 & R05-R08
 Exchanges Quaternion in the Stack with Specif. Reg.
 Exchanges X Y Z T with M N O P
 Changes the sign of a Quaternion in Specified Registers
 Recalls a Quaternion from Specified Registers
 Stores a Quaternion in Specified Registers
 Stores X Y Z T into M N O P
 Stores register X into 4 Specified Registers
 Multiplies by X the content of 4 Specified Registers
 Divides by X the content of 4 Specified Registers
 Returns the Sum of the Magnitudes of 4 Specif. Reg.
 Displays the 4 components of a Quaternion in Spec. Reg.
 Stores 0 in a block of 4 Specified Registers
 Square-Root of a Quaternion
 Returns | X | + | Y | + | Z | + | T |
 Displays the 4 components of a Quaternion in the Stack
 Section Header
 Power of 2 Quaternions
 Real Power of a Quaternion
 Exponential
 Natural Logarithm
 Hyperbolic Sine
 Hyperbolic Cosine
 Hyperbolic Tangent
 Inverse Hyperbolic Sine
 Inverse Hyperbolic Cosine
 Inverse Hyperbolic Tangent
 Sine
 Cosine
 Tangent
 Inverse Sine
 Inverse Cosine
 Inverse Tangent
 Section Header
 3-Dimensional Rotations
 Solves Systems of nxn Equations written Fi(qj) = qi
 Solves F(q) = 0
 Solves F(q) = q
 Evaluate a Polynomial P(q)
 Rectangular-Polar Conversion
 Quaternionic Unilateral Linear Systems
 Polar-Rectangular Conversion
 Sylvester Equation a.q + q.b = c
 Section Header
 Inverse of the Gamma Function ( real arguments )
 Difference between 2 consecutive Sums
 Euler's Gamma Constant
 Sum of the 4 components of a Quaternion in Specif. Reg.
 Sum of the 4 components of a Quaternion in the Stack
 Chebyshev Polynomials
 Hermite Polynomials
 Jacobi Polynomials
 Generalized Laguerre Polynomials
 Legendre Polynomials
 Associated Legendre Functions ( Integer Order & Deg. )
 UltraSpherical Polynomials
 16,00
 16,01
 16,02
 16,03
 16,04
 16,05
 16,06
 16,07
 16,08
 16,09
 16,10
 16,11
 16,12
 16,13
 16,14
 16,15
 16,16
 16,17
 16,18
 16,19
 16,20
 16,21
 16,22
 16,23
 16,24
 16,25
 16,26
 16,27
 16,28
 16,29
 -SPEC FCNS
 "AIRYQ"
 "ALFQ"
 "ANWEBQ"
 "CATQ"
 "CHIQ"
 "CIQ"
 "DNQ"
 "EIQ"
 "ENQ"
 "ERFQ"
 "GAMQ"
 "HGFQ+"
 "INQ"
 "JCFQ"
 "JEFQ"
 "JNQ"
 "KNQ"
 "LANQ"
 "PSIQ"
 "QDFT"
 "QLW"
 "RCWFQ"
 "SHIQ"
 "SIQ"
 "STRVHQ"
 "STRVLQ"
 "SWFQ"
 "USFQ"
 "YNQ"
 Section Header
 Airy Functions  Ai & Bi
 Associated Legendre Functions
 Anger & Weber Functions 
 Catalan Numbers
 Hyperbolic Cosine Integral 
 Cosine Integral 
 Parabolic Cylinder Functions
 Exponential Integral  Ei
 Exponential Integral  En
 Error Function
 Gamma Function
 Generalized Hypergeometric Functions
 Modified Bessel Functions 1st Kind
 Jacobi Functions 
 Jacobian Elliptic Functions: sn , cn , dn
 Bessel Functions 1st Kind
 Modified Bessel Functions 2nd Kind
 Generalized Laguerre's Functions
 Digamma Function
 Quaternionic Discrete Fourier Transform
 Lambert-W Function
 Regular Coulomb Wave Functions
 Hyperbolic Sine Integral 
 Sine Integral 
 Struve Functions H
 Struve Functions L
 Spheroidal Wave Functions
 UltraSpherical ( Gegenbauer ) Functions
 Bessel Functions 2nd Kind

 

Euler's Gamma Constant
 
 

-GEU is an M-code routine that simply returns  0.5772156649  in register X.
 
 

Inverse of the Gamma Function ( real arguments )
 

-The following asymptotic expansion is used:

    Gam(x) = [ x^(x-1/2) ] sqrt(2.Pi) exp [ -x + ( 1/12 )/( x + ( 1/30 )/( x + ( 53/210)/( x + (195/371)/( x + ...  )))) ]

  with the relation:  Gam(x+n) = (x+n-1) (x+n-2) ...... (x+1) x Gam(x)   until  x+n > 5
 
 
      STACK        INPUT      OUTPUT
           X             x       1/Gam(x)
           L             /            x

  PI   XEQ "1/GM"  >>>  1/Gam(Pi) = 0.437055717
 

Difference between 2 consecutive Quaternionic Sums
 

-When computing an "infinte" series of hypercomplex numbers, we usually stop the loop when 2 consecutive partial sums are equal.

-DSQ  takes the quaternion  uk in registers  X , Y , Z , T  and the adress of the first register Rnn that contains the partial sum from the next program-line.
-The partial sum is incremented by uk and it returns 0 if and only if 2 consecutive sums are equal.

-For instance, the loop    LBL 01  ................

  STO 01       STO 14     STO 03      -                  RCL 16         ABS           STO 13         X#0?
  RDN           LASTX     RCL 15      ABS            +                    +                LASTX         GTO 01
  STO 02       -                +                +                 STO 16         RCL 01       -
  RCL 14       ABS          STO 15      X<>Y         LASTX         RCL 13       ABS
  +                 X<>Y       LASTX       STO 04       -                    +                +

  may advantageously be replaced by

  LBL 01  ...................   QSTO  1   DSQ   13   X#0?   GTO 01

Example2:    May be another example is not unuseful...

-Compute  S(q) = 1 + q + q2 + q3 + ................ + qn + ...................   for  | q |  <  1
 
 
 01  LBL "TEST"
 02  QSTO
 03  5
 04  CLST
 05  SIGN
 06  QSTO
 07  1
 08  QSTO
 09  9
 10  LBL 01
 11  Q*Q
 12  QSTO
 13  1
 14  DSQ
 15  9
 16  X#0?
 17  GTO 01
 18  QRCL
 19  9
 20  END

 
-With  q = 0.1 + 0.2 i + 0.3 j + 0.4 k

  0.4  ENTER^
  0.3  ENTER^
  0.2  ENTER^
  0.1  XEQ "TEST"  >>>>   0.818181818
                                RDN   0.181818182
                                RDN   0.272727273
                                RDN   0.363636364

-Thus,  S( 0.1 + 0.2 i + 0.3 j + 0.4 k ) = ( 9 + 2 i + 3 j + 4 k ) / 11    which could be obtained directly since  S(q) = 1/(1-q)
 

Notes:

-In this example, the result is also in registers R09 thru R12.
-R01 to R04 = qn , R05 to R08 = q

-This M-code routine may also be used to simply add the quaternion in the stack to a block of 4 consecutive  registers.
 

Storing M N O P into X Y Z T
 

-QSTOA is equivalent to STO M  RDN  STO N  RDN  STO O  RDN  STO P  RDN
-Since register P is altered if a number is displayed or by any digit entry line, use it only in a program.
 

Storing X Y Z T into M N O P
 

-ASTOQ is equivalent to RCL P   RCL O   RCL N   RCL M
-Same remarks...
 

Exchanging X Y Z T with M N O P
 

-Q<>A is equivalent to  X<> M  RDN  X<> N  RDN  X<> O  RDN  X<> P  RDN
-Same remarks
 

Storing X into 4 Specified Registers
 

-QSTX nn  stores the real in X-register into registers Rnn , Rn+1 , Rn+2 , Rn+3.
-Indirect addressing is available too.
-In a program, simply add 128 to nn: For instance,

  QSTX
  128

  will perform  STO (IND 00)  STO (IND 00)+1  STO (IND 00)+2  STO (IND 00)+3.
 

Storing a null Quaternion into a block of 4 Specified Registers
 

-QZERO nn  stores 0 into registers Rnn , Rn+1 , Rn+2 , Rn+3.
-Same remarks as above.
 

Recalling a Quaternion from Specified Registers
 

-QRCL nn  is equivalent to RCLn+3  RCLn+2  RCLn+1  RCLnn
-Same remarks as above.
 

Storing a Quaternion in Specified Registers
 

-QSTO nn  is equivalent to STOnn  RDN  STOn+1  RDN  STOn+2  RDN  STOn+3  RDN
-Same remarks as above.
 

Exchanging a Quaternion in the Stack with Specified Registers
 

-Q<> nn  is equivalent to  X<> nn   RDN   X<> n+1   RDN   X<> n+2   RDN   X<> n+3   RDN
-Same remarks about indirect addressing...
 

Displaying the 4 components of a Quaternion in the Stack
 

-VIEWQ  successively displays  X , Y , Z , T.
-The display stops while a key is held, resuming when it's released.
 

Displaying the 4 components of a Quaternion in Specified Registers
 

-QVIEW nn  successively displays  Rnn , Rn+1 , Rn+2 , Rn+3
-The display stops while a key is held, resuming when it's released.
-Same remarks about indirect addressing...
 

Changing the sign of the Quaternion in the Stack
 

-CHSQ  is equivalent to  CHS  RDN  CHS  RDN  CHS  RDN  CHS  RDN
 

Changing the sign of the Quaternion in 4 Specified Registers
 

-QCHS nn  multiplies the content of registers  Rnn , Rn+1 , Rn+2 , Rn+3 by (-1)
-Same remarks about indirect addressing...
 

Multiplying by X the content of 4 Specified Registers
 

-QSTX* nn   performs  ST* nn   ST* n+1  ST* n+2  ST* n+3
-Same remarks about indirect addressing...
 

Dividing by X the content of 4 Specified Registers
 

-QSTX/ nn   performs  ST/ nn   ST/ n+1  ST/ n+2  ST/ n+3
-Same remarks about indirect addressing...
 

Sum of the 4 components of a Quaternion in the Stack
 

-SUMQ  returns  X+Y+Z+T  in X-register.
-X is saved in L-register, Y- Z- T-registers are unchanged.
 

Sum of the 4 components of a Quaternion in Specified Registers
 

-QSUM nn  returns  (Rnn)+(Rn+1)+(Rn+2)+(Rn+3)  in X-register.
-X is saved in L-register, Y- Z- T-registers are unchanged.
 

Sum of the Magnitudes of a Quaternion in the Stack
 

-SUMABQ  returns  | X | + | Y | + | Z | + | T |  in X-register.
-X is saved in L-register, Y- Z- T-registers are unchanged.
 

Sum of the Magnitudes of 4 Specified Registers
 

-QSUMAB nn  returns  | Rnn | + | Rn+1 | + | Rn+2 | + | Rn+3 |  in X-register.
-X is saved in L-register, Y- Z- T-registers are unchanged.
 

Inverse of a Quaternion
 
 
 
 
      STACK        INPUTS      OUTPUTS
           T             d             t
           Z             c             z
           Y             b             y
           X             a             x
           L             /             a

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

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 )
 

Square of a Quaternion
 
 
 
      STACK        INPUTS      OUTPUTS
           T             d             t
           Z             c             z
           Y             b             y
           X             a             x
           L             /             a

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

Example:       q = 1 + 2 i + 3 j + 4 k

   4  ENTER^
   3  ENTER^
   2  ENTER^
   1  XEQ "Q^2"  gives   -28
                           RDN     4
                           RDN     6
                           RDN     8

-So,  ( 1 + 2 i + 3 j + 4 k )2 =  -28 + 4 i + 6 j + 8 k
 

Square-Root of a Quaternion
 
 
 
      STACK        INPUTS      OUTPUTS
           T             d             t
           Z             c             z
           Y             b             y
           X             a             x
           L             /             a

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

Example:       q = 1 + 2 i + 3 j + 4 k

   4  ENTER^
   3  ENTER^
   2  ENTER^
   1  XEQ "SQRTQ"  yields   1.7996
                                  RDN   0.5557
                                  RDN   0.8335
                                  RDN   1.1113

   ( 1 + 2 i + 3 j + 4 k ) 1/2 =  1.7996 + 0.5557 i + 0.8335 j + 1.1113 k   ( Rounded to 4D )
 

Euclidean Norm
 
 
 
      STACK        INPUTS      OUTPUTS
           T            d             d
           Z             c             c
           Y             b             b
           X             a             µ
           L             /             a

 with      q =  a + b.i + c.j + d.k     ;    µ = ( a2 + b2 + c2 + d2 )1/2

Example:       q = 1 + 2 i + 3 j + 4 k

   4  ENTER^
   3  ENTER^
   2  ENTER^
   1  XEQ "ABSQ"  yields   5.477225575
 

Square of the Euclidean Norm
 
 
 
      STACK        INPUTS      OUTPUTS
           T             d             d
           Z             c             c
           Y             b             b
           X             a             µ2
           L             /             a

 with      q =  a + b.i + c.j + d.k     ;    µ2 = ( a2 + b2 + c2 + d2 )

Example:       q = 1 + 2 i + 3 j + 4 k

   4  ENTER^
   3  ENTER^
   2  ENTER^
   1  XEQ "NRM2Q"  yields   µ2 = 30
 

Product of 2 Quaternions
 

Formula:
 

  ( x + y i + z j + t k ) ( x' + y' i + z' j + t' k ) = ( x.x' - y.y' - z.z' - t.t' ) + ( x.y' + x'.y + z.t' - z'.t ) i + ( x.z' + x'.z + y'.t - y.t' ) j + ( x.t' + x'.t + y.z' - y'.z ) k
 

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

Remarks:

  q'.q = 17 - 45i -35j - 7k

-The alpha "register" is cleared.
-To store q & q' into R01 to R08, you can of course use QSTO _ _

  -7  ENTER^  4  ENTER^  -3  ENTER^  2  QSTO 01  ( or even simpler  QSTO A )
   5  ENTER^  2  ENTER^  -4  ENTER^  1  QSTO 05  ( or even simpler  QSTO E )
 

Power of 2 Quaternions
 

 Q^Q  calculates  qq' with the definition   qq' =   e(lnq).q'
 

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

  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 )

Note:

-To store q & q' into R01 to R08, you can of course use QSTO _ _

  -7  ENTER^  4  ENTER^  -3  ENTER^  2  QSTO 01  ( or even simpler  QSTO A )
   5  ENTER^  2  ENTER^  -4  ENTER^  1  QSTO 05  ( or even simpler  QSTO E )
 

Real Power of a Quaternion
 

Data Register:      R00 = r  is to be initialized before executing Q^R
Flags:  /
 
 
 
      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 )

Note:

-This program returns  0^0 = 1
 

Exponential of a Quaternion
 

 Exp(q)  could be computed by   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

 but "E^Q" uses the rectangular-polar conversion instead.
 

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

   eq  =   -5.0277 -1.8884 i + 2.5178 j - 4.4062 k   ( rounded to 4D )
 

Natural Logarithm of a Quaternion
 

 Ln q  may be calculated by  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

 but "LNQ" uses the rectangular-polar conversion instead.
 

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

Hyperbolic Sine
 

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

Hyperbolic Cosine
 

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

Hyperbolic Tangent
 

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

Data Registers:  R00 unused , R01 to R08: temp
Flag:  F09 is cleared
 
 
      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
                             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 )
 

Inverse Hyperbolic Sine
 

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

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

Inverse Hyperbolic Cosine
 

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

Data Registers:  R00 to R12: temp
Flags: /
 
 
      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
                                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 )
 

Inverse Hyperbolic Tangent
 

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

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

Sine
 

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

Cosine
 

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

Tangent
 

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

Data Registers:    R01 thru R08
Flag:  F09 is cleared
 
 
      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
                                  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
 

Inverse Sine
 

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:  R01 to R08: temp
Flag:  F09 is cleared
 
 
      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
                                 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 )
 

Inverse Cosine
 

Formula:               ArcCos q = PI/2 - ArcSin q
 

Data Registers:  R00 to R08: temp
Flag: F09 is cleared
 
 
      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
                                  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 )

Note:

 4 subroutine-levels are used by this function.
 

Inverse Tangent
 

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:  R01 to R08: temp
Flag:  F09 is clear
 
 
      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
                                  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 )
 

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

Polar-Rectangular Conversion
 

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

Evaluating a Quaternionic Polynomial
 

-"PEVALQ" computes:    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 = bb.eee = control number of the polynomial   ( Registers R00 & Rbb thru Ree are to be initialized before executing "PEVALQ" )

                                        Rbb   = an         Rb+4  = an-1       ..................................       Re-3  = a0
                                        Rb+1 = bn         Rb+5  = bn-1      ..................................       Re-2  = b0            R01 to R10: temp
                                        Rb+2 = cn         Rb+6  = cn-1       .................................        Re-1  = c0
                                        Rb+3 = dn         Rb+7  = dn-1       .................................        Ree   = d0            bb and bb' > 10

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

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

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

  Provided  R00 =  bbb.eee = control number of the polynomial P                                            bbb > 010

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 [ R11  R12  R13  ..................  R26 ]   ====>   Control number = 11.026

   11.026   STO 00

    -3   ENTER^
    -2   ENTER^
    -1   ENTER^
     6   XEQ "PEVALQ"  >>>>   2456
                                       RDN    -222
                                       RDN    -159
                                       RDN    -894

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

3-Dimensional Rotations
 

-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
 

-"3DROT" 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     R04 = a
                               R02 = yA     R05 = b                 R07 to R15 are used for temporary data storage
                               R03 = zA     R06 = c
Flags:  /
 
 
      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 "3DROT"  >>>>  1.009250425   RDN   3.497956694   RDN   8.330584238

-Whence    M'( 1.009250425 , 3.497956694 , 8.330584238 )
 

Note:

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

Chebyshev Polynomials
 

Formulae:

      1st kind          Tn(q) = 2q.Tn-1(q) - Tn-2(q)   ;  T0(q) = 1  ;   T1(q) = q
      2nd kind       Un(q) = 2q.Un-1(q) - Un-2(q)  ;  U0(q) = 1  ;  U1(q) = 2q
 

Data Registers:        R00 = n                                  ( Register R00 is to be initialized before executing "CHBQ" )

                              When the program stops,  Tn(q) or Un(q)  are in registers R05 thru R08
                                                                    Tn-1(q) or Un-1(q)  are in registers R09 thru R12

Flag: F02 ( F09 is cleared )

  if F02 is clear, "CHBQ" calculates the Chebyshev polynomials of the first kind, and the HP-41 displays "TNQ" at the end
  if F02 is set,    "CHBQ" calculates the Chebyshev polynomials of the second kind, and the HP-41 displays "UNQ" at the end
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Tn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k        if   CF 02       with       R00 = n    integer    ( 0 <= n < 1000 )
      or         Un(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k       if    SF 02

Example:         n = 7     q = 1 + 2 i + 3 j + 4 k

   7   STO 00

      CF 02      Chebyshev Polynomials 1st kind

   4   ENTER^
   3   ENTER^
   2   ENTER^
   1   XEQ "CHBQ"   >>>>       -9524759      after displaying  "TNQ"
                                  RDN       -1117678
                                  RDN       -1676517
                                  RDN       -2235356

        T7( 1 + 2 i + 3 j + 4 k ) =  -9524759 - 1117678 i - 1676517 j - 2235356 k

      SF 02      Chebyshev Polynomials 2nd kind

   4   ENTER^
   3   ENTER^
   2   ENTER^
   1   XEQ "CHBQ"   >>>>     -18921448      after displaying  "UNQ"
                                  RDN       -2198096
                                  RDN       -3297144
                                  RDN       -4396192

        U7( 1 + 2 i + 3 j + 4 k ) =  -18921448 - 2198096 i - 3297144 j - 4396192 k
 

Hermite Polynomials
 
 

Formulae:        H0(q) = 1  ,  H1(q) = 2 q
                         Hn(q) = 2.q Hn-1(q) - 2 (n-1) Hn-2(q)
 

Data Registers:        R00 = n                                 ( Register R00 is to be initialized before executing "HMTQ" )

                 At the end,  R05 to R08 = Hn(q)   &   R09 to R12 = Hn-1(q)
Flags: /
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Hn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k         with       R00 = n    integer    ( 0 <= n < 1000 )

Example:         n = 7     q = 1 + 2 i + 3 j + 4 k

   7   STO 00

   4   ENTER^
   3   ENTER^
   2   ENTER^
   1   XEQ "HMTQ"   >>>>      -23716432 = R05
                                  RDN         -3653024 = R06
                                  RDN         -5479536 = R07
                                  RDN         -7306048 = R08

        P7( 1 + 2 i + 3 j + 4 k ) =  -23716432 - 3653024 i - 5479536 j - 7306048 k            and in registers R09 thru R12:

       P6( 1 + 2 i + 3 j + 4 k ) =  -1122232 + 682816 i + 1024224 j + 1365632 k
 

Jacobi Polynomials
 

Formulae:      P0(a;b) (q) = 1  ;   P1(a;b) (q) = (a-b)/2 + q (a+b+2)/2

        2n(n+a+b)(2n+a+b-2) Pn(a;b) (q) = [ (2n+a+b-1).(a2-b2) + q (2n+a+b-2)(2n+a+b-1)(2n+a+b) ] Pn-1(a;b) (q)
                                                              - 2(n+a-1)(n+b-1)(2n+a+b) Pn-2(a;b) (q)
 

Data Registers:          R00 thru R22: temp              ( Registers R09-R10-R11 are to be initialized before executing "JCPQ" )

                                   R09 = a
                                   R10 = b
                                   R11 = n

                 At the end,  R05 to R08 = Pn(a;b)(q)               R16 to R19 = q
                             &   R12 to R15 = Pn-1(a;b)(q)
Flags:  /
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Pn(a;b)(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k         with       R09 = a  ,  R10 = b  ,  R11 = n         ( n integer  0 <= n < 1000 )

Example:        a = sqrt(2)     b = sqrt(3)    n = 7      q = 1 + i/2 + j/3 + k/4

   2   SQRT   STO 09
   3   SQRT   STO 10
   7   STO 11

   4   1/X
   3   1/X
   2   1/X
   1   XEQ "JCPQ"   >>>>   143.5304513
                                RDN   -310.8681606
                                RDN   -207.2454403
                                RDN   -155.4340799

        P7sqrt(2);sqrt(3) ( 1 + i/2 + j/3 + k/4 ) =  143.5304513 - 310.8681606 i - 207.2454403 j - 155.4340799 k
 

Jacobi Functions
 

Formula:           Pn(a;b) (q) = [ Gam(a+n+1) / Gam(n+1) ]  2F~1 ( -n , a+b+n+1 , a+1 , (1-q)/2 )

             where  2F1 tilde is the regularized hypergeometric function
 

Data Registers:        R09 = a                             ( Registers R09-R10-R11 are to be initialized before executing "JCFQ" )
                                   R10 = b
                                   R11 = n                                     R00 thru R17: temp

                 At the end,  R01 to R04 = Pn(a;b) (q)

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

   where      Pn(a;b)(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k         with       R09 = a  ,  R10 = b  ,  R11 = n

Examples:        a = sqrt(2)     b = sqrt(3)    n = PI      q = 1 + i/2 + j/3 + k/4

   2   SQRT   STO 09
   3   SQRT   STO 10
  PI  STO 11

   4   1/X
   3   1/X
   2   1/X
   1   XEQ "JCFQ"   >>>>   -10.14397527
                                RDN     11.74139998
                                RDN       7.827599988
                                RDN       5.870699987

        PPIsqrt(2);sqrt(3) ( 1 + i/2 + j/3 + k/4 ) =  -10.14397527 + 11.74139998 i + 7.827599988 j + 5.870699987 k
 

-Likewise,   PPI(-4;+4) ( 1 + i/2 + j/3 + k/4 ) =  0.360213558 - 0.125970281 i - 0.083980187 j - 0.062985140 k

Note:

-Unless the function is a polynomial, the series is convergent if  | 1 - q | < 2
 

Generalized Laguerre Polynomials
 

Formulae:      L0(a) (q) = 1  ,   L1(a) (q) = a+1-q   ,    n Ln(a) (q) = (2.n+a-1-q) Ln-1(a) (x) - (n+a-1) Ln-2(a) (q)
 

Data Registers:        R09 = a     R10 = n              ( Registers R09 & R10 are to be initialized before executing "LAGQ" )

                 At the end,  R05 to R08 = Ln(a)(q)   &   R11 to R14 = Ln-1(a)(q)
Flags: /
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Ln(a)(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k         with       R09 = a  &  R10 = n   integer    ( 0 <= n < 1000 )

Example:        a = sqrt(2)     n = 7     q = 1 + 2 i + 3 j + 4 k

   2   SQRT   STO 09
   7   STO 10

   4   ENTER^
   3   ENTER^
   2   ENTER^
   1   XEQ "LAGQ"   >>>>  872.2661287  = R05
                                 RDN    47.12527437 = R06
                                 RDN    70.68791171 = R07
                                 RDN    94.25054886 = R08

        L7sqrt(2) ( 1 + 2 i + 3 j + 4 k ) =  872.2661287 + 47.12527437 i + 70.68791171 j + 94.25054886 k                           and

       L6sqrt(2) ( 1 + 2 i + 3 j + 4 k ) =   335.6848018 + 123.7427961 i + 185.6141942 j + 247.4855923 k           ( in registers R11 thru R14 )
 

Laguerre's Functions
 

Formula:      Ln(a)(q) = [ Gam(n+a+1) / Gam(n+1) ]  1F~1( -n , a+1 , q )

   where  Gam = Gamma function
     and     1F~1 = Regularized Kummer's function
 

Data Registers:        R09 = a     R10 = n           ( Registers R09 & R10 are to be initialized before executing "LANQ" )

                                     R00 thru R18: temp
Flags: /
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Ln(a)(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k         with       R09 = a  &  R10 = n , n # -1 , -2 , .............

Example:        a = sqrt(3)     n = PI     q = 1 + 2 i + 3 j + 4 k

   3   SQRT   STO 09
  PI   STO 10

   4   ENTER^
   3   ENTER^
   2   ENTER^
   1   XEQ "LANQ"   >>>>  -59.37663940
                                  RDN     3.135355693
                                  RDN     4.703033568
                                  RDN     6.270711359

        Lpisqrt(3) ( 1 + 2 i + 3 j + 4 k ) =  -59.37663940 + 3.135355693 i + 4.703033568 j + 6.270711359 k
 

Legendre Polynomials
 

Formulae:      n.Pn(q) = (2n-1).q.Pn-1(q) - (n-1).Pn-2(q)  ,  P0(q) = 1  ,  P1(q) = q
 

Data Registers:        R00 = n                 ( Register R00 is to be initialized before executing "LEGQ" )

                                    R01 to R04 = q                   At the end,  R05 to R08 = Pn(q)   &   R09 to R12 = Pn-1(q)
Flags: /
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Pn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k         with       R00 = n    integer    ( 0 <= n < 1000 )

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

   7   STO 00

   4   1/X
   3   1/X
   2   1/X
   1   XEQ "LEGQ"   >>>>    36.20819589 = R05
                                 RDN   -51.58337290 = R06
                                 RDN   -34.38891527 = R07
                                 RDN   -25.79168646 = R08

        P7( 1 + i/2 + j/3 + k/4 ) =  36.20819589 - 51.58337290 i - 34.38891527 j - 25.79168646 k            and in registers R09 thru R12:

       P6( 1 + i/2 + j/3 + k/4 ) =  -9.248135367 - 26.20689563 i - 17.47126375 j - 13.10344782 k
 

Associated Legendre Functions 1st kind ( Integer Order & Degree )
 
 

Formulae:      (n-m) Pnm(q) = q (2n-1) Pn-1m(q) - (n+m-1) Pn-2m(q)

                               Pmm(q) = (-1)m (2m-1)!!  ( 1-q2 )m/2                    where (2m-1)!! = (2m-1)(2m-3)(2m-5).......5.3.1
 
 

Data Registers:      R09 = m            ( Registers R09 &  R10 are to be initialized before executing "PMNQ" )
                                  R10 = n

                      When the program stops, the result is in R05 to 508

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

   where      Pmn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k                  with       R09 = m   &  R10 = n         m & n integers   0 <= m <= n

Example:      m = 3    n = 7       q = 1 + i/2 + j/3 + k/4

       3   STO 09    7   STO 10

      4   1/X
      3   1/X
      2   1/X
      1   XEQ "PMNQ"   >>>>   -12188.53940
                                     RDN    -8670.127573
                                     RDN    -5780.085045
                                     RDN    -4335.063808

     P37 ( 1 + i/2 + j/3 + k/4 ) = -12188.53940 - 8670.127573 i - 5780.085045 j - 4335.063808 k
 

Associated Legendre Functions
 

Formulae:
 

   1st Kind:      Pnm(q) = [ (q+1)/(1-q) ]m/2  2F~1(-n , n+1 ; 1-m ; (1-q)/2 )        (  q # 1 )         F~ = regularized hypergeometric function

  2nd Kind:      Qnm(q) = 2m pi1/2  (1-q2)-m/2 [ -Gam((1+m+n)/2)/(2.Gam((2-m+n)/2)) . sin pi(m+n)/2 .  2F1(-n/2-m/2 ; 1/2+n/2-m/2 ; 1/2 ; q2 )
                                                                     + q Gam((2+n+m)/2) / Gam((1+n-m)/2) . cos pi(m+n)/2 . 2F1((1-m-n)/2 ; (2+n-m)/2 ; 3/2 ; q2 )  ]
 

Data Registers:      R09 = m            ( Registers R09 &  R10 are to be initialized before executing "ALF1Q" )
                                  R10 = n

Flag:  F02

   CF 02   "ALFQ"   will return  ALF 1st kind      The HP-41 displays  "ALF1" at the end in this case
   SF 02    "ALFQ"  will return   ALF 2nd kind     The HP-41 displays  "ALF2"  at the end in this case
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Pmn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k           if  CF 02               (   with       R09 = m   &   R10 = n   )
    and        Qmn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k           if  SF 02

Example1:      m = sqrt(2)    n = sqrt(3)       q = 0.4 + 0.5 i + 0.6 j + 0.7 k

      2  SQRT   STO 09
      3  SQRT   STO 10     CF 02

      0.7   ENTER^
      0.6   ENTER^
      0.5   ENTER^
      0.4   XEQ "ALFQ"   >>>>    -0.866899244           after displaying  "ALF1"
                                       RDN    -1.808960474
                                       RDN    -2.170752569
                                       RDN    -2.532544663

     Psqrt(2)sqrt(3) ( 0.4 + 0.5 i + 0.6 j + 0.7 k ) = -0.866899244 - 1.808960474 i - 2.170752569 j - 2.532544663 k

Example2:       m = 0.4    n = 1.3       q = 0.1 + 0.2 i + 0.3 j + 0.4 k

  0.4   STO 09
  1.3   STO 10          SF 02

  0.4    ENTER^
  0.3    ENTER^
  0.2    ENTER^
  0.1    XEQ "ALFQ"  >>>>   -0.944581517           after displaying  "ALF2"
                                   RDN   -0.368022564
                                   RDN   -0.552033846
                                   RDN   -0.736045128

    Q1.30.4( 0.1 + 0.2 i + 0.3 j + 0.4 k ) = -0.944581517 - 0.368022564 i - 0.552033846 j - 0.736045128 k

-Likewise,   Q43( 0.1 + 0.2 i + 0.3 j + 0.4 k ) = 131.9022421 - 13.96653103 i - 20.94979655 j - 27.93306208 k

Notes:

-For the Associated Legendre Function of the 1st kind, the series usually converges  if   | 1 - q | < 2 , q # 1
-For the Associated Legendre Function of the 2nd kind, the series usually converge  if   | q | < 1
 

UltraSpherical Polynomials
 

Formulae:      C0(a) (q) = 1  ;   C1(a) (q) = 2.a.q   ;    (n+1).Cn+1(a) (q) = 2.(n+a).x.Cn(a) (q) - (n+2a-1).Cn-1(a) (q)      if  a # 0

                      Cn(0) (q) = (2/n).Tn(q)
 

Data Registers:        R09 = a     R10 = n > 0           ( Registers R09 & R10 are to be initialized before executing "USPQ" )

                 At the end,  R05 to R08 = Cn(a)(q)
                             &   R11 to R14 = Cn-1(a)(q)  if  a # 0

Flag:   F02 & F09  if a = 0
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Cn(a)(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k         with       R09 = a  &  R10 = n   integer    ( 0 <= n < 1000 )

Example:        a = sqrt(2)     n = 7     q = 1 + i/2 + j/3 + k/4

   2   SQRT   STO 09
   7   STO 10

   4   1/X
   3   1/X
   2   1/X
   1   XEQ "USPQ"   >>>>    324.5443960
                                 RDN   -689.5883620
                                 RDN   -459.7255746
                                 RDN   -344.7941809

        C7sqrt(2) ( 1 + i/2 + j/3 + k/4 ) =  324.5443960 - 689.5883620 i - 459.7255746 j - 344.7941809 k

-Likewise,   C7(0) ( 1 + i/2 + j/3 + k/4 ) =  29.24554794 - 33.27910357 i - 22.18606907 j - 16.63955179 k
 

UltraSpherical ( Gegenbauer ) Functions
 

 Formula:    Assuming  a # 0

          Cn(a)(q) = [ Gam(n+2a) / Gam(n+1) / Gam(2a) ]  2F1( -n , n+2a , n+1/2 , (1-q)/2 )

                     where  2F1 is the hypergeometric function and Gam = Gamma function
 

Data Registers:        R09 = a # 0     R10 = n          ( Registers R09 & R10 are to be initialized before executing "USFQ" )

    >>>      At the end,  R01 to R04 = Cn(a)(q)

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

   where      Cn(a)(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k         with       R09 = a # 0  &  R10 = n

Example:        a = sqrt(2)     n = sqrt(3)     q = 1 + i/2 + j/3 + k/4

   2   SQRT   STO 09
   3   SQRT   STO 10

   4   1/X
   3   1/X
   2   1/X
   1   XEQ "USFQ"   >>>>    3.241115514
                                 RDN    4.848277901
                                 RDN    3.232185266
                                 RDN    2.424138950

        Csqrt(3)sqrt(2) ( 1 + i/2 + j/3 + k/4 ) =  3.241115514 + 4.848277901 i + 3.232185266 j + 2.424138950 k
 

Airy Functions
 

Formulae:

    Ai(q) =   f(q) - g(q)                             with            f(q) = [ 3 -2/3 / Gamma(2/3) ]  0F1( 2/3 ; q3/9 )
    Bi(q) = [ f(q) + g(q) ] sqrt(3)               and            g(q) = [ 3 -1/3 / Gamma(1/3) ]  0F1( 4/3 ; q3/9 )   q
 

Data Registers:    R00 thru R24:  temp

                   At the end,  R01 to R04 = Ai(q)
                         and        R05 to R08 = Bi(q)
Flags: /
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Ai (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k   = R01 thru R04
     and       Bi (x+y.i+z.j+t.k)  = x"+y".i+z".j+t".k  = R05 thru R08

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

   4   1/X
   3   1/X
   2   1/X
   1   XEQ "AIRYQ"    >>>>    0.105299445 = R01           the HP-41 displays  "BI=R05-R08"
                                    RDN   -0.078442074 = R02
                                    RDN   -0.052294716 = R03
                                    RDN   -0.039221037 = R04

Thus,    Ai ( 1 + i/2 + j/3 + k/4 ) = 0.105299445 - 0.078442074 i - 0.052294716 j - 0.039221037 k

 and     Bi ( 1 + i/2 + j/3 + k/4 ) = 0.973463977 + 0.394832415 i + 0.263221610 j + 0.197416207 k    in   R05 to R08
 

Note:

-This program returns accurate results for "small" arguments only.
 

Anger & Weber Functions
 

Formulae:

  Jn(q) = + (q/2) sin ( 90°n ) 1F~2( 1 ; (3-n)/2 , (3+n)/2 ; -q2/4 )            Anger's functions

                  + cos ( 90°n ) 1F~2( 1 ; (2-n)/2 , (2+n)/2 ; -q2/4 )

   En(q) = - (q/2) cos ( 90°n )  1F~2( 1 ; (3-n)/2 , (3+n)/2 ; -q2/4 )         Weber's functions

                    + sin ( 90°n ) 1F~2( 1 ; (2-n)/2 , (2+n)/2 ; -q2/4 )
 

Data Registers:           R00 = n                                   ( Registers R00 is to be initialized before executing "ANWEBQ" )

        When the program stops,   R01 to R04 = Jn(q)        ( Anger functions )
                                                 R05 to R08 = En(q)       ( Weber functions )              R19 to R22 = q
Flags:  /
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Jn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k   = R01 thru R04            provided        R00 = n
     and       En(x+y.i+z.j+t.k)  = x"+y".i+z".j+t".k  = R05 thru R08

Example:         n = PI     q = 1 + 2 i + 3 j + 4 k

  PI  STO 00

   4   ENTER^
   3   ENTER^
   2   ENTER^
   1   XEQ "ANWEBQ"   >>>>  -11.24234895 = R01                    HP-41 displays  "WEB=R05-R08"
                                       RDN    -3.564968430 = R02
                                       RDN    -5.347452640 = R03
                                       RDN    -7.129936859 = R04

            JPi( 1 + 2 i + 3 j + 4 k ) = -11.24234895 - 3.564968430 i - 5.347452640 j - 7.129936859 k

 and      EPi( 1 + 2 i + 3 j + 4 k ) = -9.566817351 + 4.177204888 i + 6.265807332 j + 8.354409774 k             in  R05 to R08

Notes:

-If n is an integer,  Jn(x) = Jn(x) = Bessel function of the first kind.
-To recall En(q) in the stack,  QRCL 05  ( or QRCL E ) is the simplest way, especially if QRCL _ _ is assigned to a key )
 

Sine Integral
 

-The following series is used

   Si(q)  = Sumn=0,1,2,..... (-1)n q2n+1/((2n+1).(2n+1)!)
 

Data Registers:  R00 thru R12: temp

                             When the program stops, the result is in R13 to R16
Flag:   F01
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Si (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k

Example:              q = 1 + 2 i + 3 j + 4 k

  4   ENTER^
  3   ENTER^
  2   ENTER^
  1   XEQ "SIQ"   >>>>    17.98954626
                             RDN      7.075096293
                             RDN    10.61264444
                             RDN    14.15019259

   So,   Si ( 1 + 2 i + 3 j + 4 k ) =  17.98954626 + 7.075096293 i + 10.61264444 j + 14.15019259 k
 

Hyperbolic Sine Integral
 

-The following series is used

  Shi(q) = Sumn=0,1,2,.....  q2n+1/((2n+1).(2n+1)!)
 

Data Registers:  R00 thru R12: temp

                             When the program stops, the result is in R13 to R16
Flag:   F01
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where    Shi (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k

Example:              q = 1 + 2 i + 3 j + 4 k

  4   ENTER^
  3   ENTER^
  2   ENTER^
  1   XEQ "SHIQ"       >>>>    -0.160779307
                                   RDN      0.522153863
                                   RDN      0.783230794
                                   RDN      1.044307726

  And   Shi ( 1 + 2 i + 3 j + 4 k ) =  -0.160779307 + 0.522153863 i + 0.783230794 j + 1.044307726 k
 

Cosine Integral
 

   Ci(q) is computed by

   Ci(q)  = C + ln(q) + Sumn=1,2,..... (-1)n q2n/(2n.(2n)!)                where  C = 0.5772156649... =  Euler's constant.
 

Data Registers:  R00 thru R12: temp

                             When the program stops, the result is in R13 to R16

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

   where      Ci (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k

Example:              q = 1 + 2 i + 3 j + 4 k

  4   ENTER^
  3   ENTER^
  2   ENTER^
  1   XEQ "CIQ"   >>>>    19.04999179
                             RDN     -6.098016773
                             RDN     -9.147025159
                             RDN   -12.19603355

   So,   Ci ( 1 + 2 i + 3 j + 4 k ) =  19.04999179 - 6.098016773 i - 9.147025159 j - 12.19603355 k
 

Hyperbolic Cosine Integral
 

    Chi(q) is computed by

   Chi(q)= C + ln(q) + Sumn=1,2,.....  q2n/(2n.(2n)!)                 where  C = 0.5772156649... =  Euler's constant.
 

Data Registers:  R00 thru R12: temp

                             When the program stops, the result is in R13 to R16

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

   where    Chi (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k

Example:              q = 1 + 2 i + 3 j + 4 k

  4   ENTER^
  3   ENTER^
  2   ENTER^
  1   XEQ "CHIQ"      >>>>    -0.220285753
                                  RDN      0.529901613
                                  RDN      0.794852418
                                  RDN      1.059803225

   So,   Chi ( 1 + 2 i + 3 j + 4 k ) =  -0.220285753 + 0.529901613 i + 0.794852418 j + 1.059803225 k
 

Exponential Integral  Ei
 

Formula:      Ei(q)  = C + ln(q) + Sumn=1,2,.....  qn/(n.n!)                where  C = 0.5772156649... = Euler's constant.
 

Data Registers:  R00 thru R12: temp

                             When the program stops, the result is in R13 to R16
Flags:  /
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X            x            x'

   where      Ei (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k

Example:              q = 1 + 2 i + 3 j + 4 k

  4   ENTER^
  3   ENTER^
  2   ENTER^
  1   XEQ "EIQ"   >>>>     -0.381065062
                             RDN      1.052055476
                             RDN      1.578083213
                             RDN      2.104110953

Whence,   Ei ( 1 + 2 i + 3 j + 4 k ) =  -0.381065062 + 1.052055476 i + 1.578083213 j + 2.104110953 k
 

Exponential Integral  En
 

Formulae:

        if  n = 0                                   En(q) =  [ exp(-q) ]  / q
        if  n # 1 , 2 , 3 , .......               En(q) = qn-1 Gam(1-n) - [1/(1-n)] M( 1-n , 2-n ; -q )      where    M = Kummer's function
        if n = 1 , 2 , 3 , ........              En(q) = (-q)n-1 ( -Ln q - gamma + Sumk=1,...,n-1 1/k ) / (n-1)!  - Sumk#n-1 (-q)k / (k-n+1) / k!

                                                                                where  gamma = Euler's Constant = 0.5772156649...

Asymptotic Expansion:                En(q)  ~  (1/q) exp(-q)  2F0( 1 , n ; -1/q )
 

Data Registers:             R00 = n                             ( Register R00 is to be initialized before executing "ENQ" )

                                          R01 to R18: temp
Flag:  F04

   CF 04  to use ascending series
   SF 04 to use asymptotic expansion ( large arguments )
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      En(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k                 with       R00 = n

Example1:         n = PI     q = 1 + i/2 + j/3 + k/4

  PI   STO 00    CF 04

   4   1/X
   3   1/X
   2   1/X
   1   XEQ "ENQ"   >>>>    0.066444332
                               RDN   -0.059916128
                               RDN   -0.039944086
                               RDN   -0.029958064

   EPI( 1 + i/2 + j/3 + k/4 ) = 0.066444332 - 0.059916128 i - 0.039944086 j - 0.029958064 k

Example2:         n = 2     q = 1 + i/2 + j/3 + k/4

   2   STO 00    CF 04

   4   1/X
   3   1/X
   2   1/X
   1   XEQ "ENQ"   >>>>    0.082262051
                               RDN   -0.087391297
                               RDN   -0.058260865
                               RDN   -0.043695649

   E2( 1 + i/2 + j/3 + k/4 ) = 0.082262051 - 0.087391297 i - 0.058260865 j - 0.043695649 k

Example3:         n = PI     q = 41 + 42 i + 43 j + 44 k

  PI   STO 00    SF 04

   44   ENTER^
   43   ENTER^
   42   ENTER^
   41   XEQ "ENQ"   >>>>    1.789481147 E-20
                                 RDN   -1.315072133 E-21
                                 RDN   -1.346383379 E-21
                                 RDN   -1.337694607 E-21

   EPI( 41 + 42 i + 43 j + 44 k ) = 1.789481147 E-20 - 1.315072133 E-21 i - 1.346383379 E-21 j - 1.337694607 E-21 k
 

Error Function
 

Formulae:

     erf q  = (2/pi1/2)  SUMn=0,1,2,.....  (-1)n q2n+1 / (n! (2n+1))

     erf q ~  [ q / sqrt(q2) ] [ exp(-q2) / ( q sqrt(PI) ] 2F0( 1 , 1/2 ; -1/q2 )
 

Data Registers:  R00 thru R17: temp
Flag:  F04

   CF 04  to use ascending series
   SF 04  to use asymptotic expansion
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Erf (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k

Example1:              q = 1 + 1.1 i + 1.2 j + 1.3 k

   CF 04

  1.3   ENTER^
  1.2   ENTER^
  1.1   ENTER^
   1     XEQ "ERFQ"   >>>>     -2.355991408
                                  RDN      -3.358261142
                                  RDN      -3.663557605
                                  RDN      -3.968854074

Whence,   Erf ( 1 + 1.1 i + 1.2 j + 1.3 k ) =  -2.355991408 - 3.358261142 i - 3.663557605 j - 3.968854074 k
 

Example2:             q = 6 + 5 i + 4 j + 3 k

   SF 04

   3   ENTER^
   4   ENTER^
   5   ENTER^
   6   XEQ "ERFQ"   >>>>    46022.51173
                                RDN    -40275.87530
                                RDN    -32220.70024
                                RDN    -24165.52518

And   Erf ( 6 + 5 i + 4 j + 3 k ) = 46022.51173 - 40275.87530 i - 32220.70024 j - 24165.52518 k
 

Gamma Function
 

Formula:    Gam(q) =  exp [ (q-1/2) Ln q + Ln (2.PI)1/2  - q + ( 1/12 )/( q + ( 1/30 )/( q + ( 53/210)/( q + (195/371)/( q + ...  )))) ]    for Re(q) > 5
 

Data Registers:   R00 unused     R01 to R16: temp
Flags: /
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Gam (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k

Example:              q = 4 + 3 i + 2 j + k

  1   ENTER^
  2   ENTER^
  3   ENTER^
  4   XEQ "GAMQ"   >>>>     0.541968821
                                  RDN    -0.740334194
                                  RDN    -0.493556130
                                  RDN    -0.246778065

     Gam ( 4 + 3 i + 2 j + k ) =  0.541968821 - 0.740334194 i - 0.493556130 j - 0.246778065 k
 

Catalan Numbers
 

-The Catalan numbers are defined by:            C(q) = 4q Gam(q+1/2) / [ sqrt(PI) Gam(q+2) ]
 

Data Registers:   R00 unused     R01 to R24: temp
Flags: /
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

        where     C (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k

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

  4   1/X
  3   1/X
  2   1/X
  1   XEQ "CATQ"     >>>>    0.844036225
                                  RDN     0.239299403
                                  RDN     0.159532935
                                  RDN     0.119649701

   Whence     C ( 1 + i/2 + j/3 + k/4 ) =  0.844036225 + 0.239299403 i + 0.159532935 j + 0.119649701 k
 

Digamma Function
 

Formula:             Psi(q) ~ Ln q - 1/(2q) -1/(12q2) + 1/(120q4) - 1/(252q6) + 1/(240q8)    is used if  Re(q) > 8

                        Psi(q+1) = Psi(q) + 1/q    is used recursively until  Re(q+1+....+1) > 8
 

Data Registers:   R00 unused     R01 to R16: temp   When the program stops,  Psi(q) is in registers R01 thru R04
Flags: /
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Psi (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k

Example:              q = 1 + 2 i + 3 j + 4 k

  4   ENTER^
  3   ENTER^
  2   ENTER^
  1   XEQ "PSIQ"   >>>>     1.686531556
                               RDN     0.548896352
                               RDN     0.823344528
                               RDN     1.097792703

     Psi ( 1 + 2 i + 3 j + 4 k ) =  1.686531556 + 0.548896352 i + 0.823344528 j + 1.097792703 k
 

Generalized Hypergeometric Functions
 

 "HGFQ+" computes   mFp( a1,a2,....,am ; b1,b2,....,bp ; q ) =  SUMk=0,1,2,.....    [(a1)k(a2)k.....(am)k] / [(b1)k(b2)k.....(bp)k] . qk/k!         if   R00 > 0

           where (ai)k = ai(ai+1)(ai+2) ...... (ai+k-1)   &  (ai)0 = 1    ,    likewise for  (bj)k    ( Pochhammer's symbol )

 or the regularized function F tilde:

           mF~p( a1,a2,....,am ; b1,b2,....,bp ; q ) =  SUMk=0,1,2,.....    [ (a1)k(a2)k.....(am)k ] / [Gam(k+b1) Gam(k+b2).....Gam(k+bp)] . qk/k!      if  R00 < 0

    where Gam = Euler's Gamma function.
 

Data Registers:            R00 =  +/- m.ppp       ( Registers R00 & R16 thru R15+m+p are to be initialized before executing "HGFQ+"  )

                                        R01 to R04 = q
                                        R05 to R08 = uk
                         >>>>      R09 to R12 = the result when the program stops.                       R13-R14-R15: temp

                                      R16 = a1        R17 = a2   ........................     Rm+15 = am
                                      Rm+16 = b1   Rm+17 = b2  ....................     Rm+p+15 = bp
Flags: /
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where     x'+y'.i+z'.j+t'.k = mFp ( a1,a2,....,am ; b1,b2,....,bp ;  x+y.i+z.j+t.k )           if   R00 = + m.ppp     and      R16 = a1  ,  R17 = a2  , ............

      or        x'+y'.i+z'.j+t'.k = mF~p ( a1,a2,....,am ; b1,b2,....,bp ;  x+y.i+z.j+t.k )          if    R00 = - m.ppp

Example1:             a1 = PI  ,   a2 = e  ;  b1 = 1 ,  b2 = 2 ,   b3 = 4  ;   q = 1 + 2 i + 3 j + 4 k

   PI   STO 16    1  E^X  STO 17  ,   1  STO 18   2  STO 19   4  STO 20

    2.003  STO 00

        4   ENTER^
        3   ENTER^
        2   ENTER^
        1   XEQ "HGFQ+"  >>>>   -6.691126905
                                       RDN      1.302530584
                                       RDN      1.953795874
                                       RDN      2.605061165

           2F3( PI , e ; 1 , 2 , 4 ;  1 + 2 i + 3 j + 4 k ) = -6.691126905 + 1.302530584 i + 1.953795874 j + 2.605061165 k

    2.003  CHS  STO 00

        4   ENTER^
        3   ENTER^
        2   ENTER^
        1   XEQ "HGFQ+"  >>>>   -1.115187818
                                       RDN      0.217088431
                                       RDN      0.325632646
                                       RDN      0.434176861

           2F~3( PI , e ; 1 , 2 , 4 ;  1 + 2 i + 3 j + 4 k ) = -1.115187818 + 0.217088431 i + 0.325632646 j + 0.434176861 k

-The first result has been simply divided by  Gam(1) Gam(2) Gam(4) = 6
 

Example2:             a1 = PI  ,   a2 = e  ;  b1 = 1 ,  b2 = -2 ,   b3 = -4  ;   q = 1 + 2 i + 3 j + 4 k

   PI   STO 16    1  E^X  STO 17  ,   1  STO 18   -2  STO 19   -4  STO 20

    2.003  CHS  STO 00        ( the non-regularized hypergeometric function does not exist for these parameters )

        4   ENTER^
        3   ENTER^
        2   ENTER^
        1   XEQ "HGFQ+"  >>>>   -2910223.713
                                       RDN      192140.2265
                                      RDN      288210.3384
                                       RDN      384280.4532

           2F~3( PI , e ; 1 , -2 , -4 ;  1 + 2 i + 3 j + 4 k ) = -2910223.713 + 192140.2265 i + 288210.3384 j + 384280.4532 k
 

Jacobian Elliptic Functions: sn , cn , dn
 

 "JEFQ" computes sn ( q | m ) , cn ( q | m ) , dn ( q | m )  where m is a real number

-The ascending Landen transformation is used if  0 < m < 1 ( cf reference  [1] or [2] )

-If m > 1, the program also uses the relations:

      sn ( q | m ) = sn ( q m1/2 | 1/m ) / m1/2
      cn ( q | m ) = dn ( q m1/2 | 1/m )
      dn ( q | m ) = cn ( q m1/2 | 1/m )

-If m = 1:      sn ( q | m ) = tanh q ; cn ( q | m ) = dn ( q | m ) = sech q

-If m < 0:

      sn ( q | m ) = sd ( q (1-m)1/2 | -m/(1-m) ) / (1-m)1/2
      cn ( q | m ) = cd ( q (1-m)1/2 | -m/(1-m) )
      dn ( q | m ) = nd ( q (1-m)1/2 | -m/(1-m) )

-If m = 0 ,   sn ( q | 0 ) = Sin q  ;   cn ( q | 0 ) = Cos q  ;  dn ( q | 0 ) = 1  are obtained directly.
 

Data Registers:             R00 = m                               ( Register R00 is to be initialized before executing "JEFQ" )

                                         R01 to R04 = sn ( q | m )                      R13 thru R26...... : temp
                                         R05 to R08 = cn ( q | m )
                                         R09 to R12 = dn ( q | m )

Flags:  F06-F07

 >>>  The number of data registers may reach 40 or more
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where     x'+y'.i+z'.j+t'.k = sn ( x+y.i+z.j+t.k | m )      provided  R00 = m

Example:      m = 1/PI   q = 1 + 2 i + 3 j + 4 k

    PI   1/X   STO 00

     4    ENTER^
     3    ENTER^
     2    ENTER^
     1    XEQ "JEFQ"   >>>>    1.498262250            The HP-41 displays   "SNCNDN=R01-R12"
                                  RDN     0.205407205
                                  RDN     0.308110808
                                  RDN     0.410814411

     sn ( 1 + 2 i + 3 j + 4 k | 1/PI ) = 1.498262250 + 0.205407205 i + 0.308110808 j + 0.410814411 k    and

     cn ( 1 + 2 i + 3 j + 4 k | 1/PI ) = -0.694940084 + 0.442849489 i + 0.664274234 j + 0.885698978 k   in  R05 to 508

     dn ( 1 + 2 i + 3 j + 4 k | 1/PI ) = -0.719248903 + 0.136199160 i +0.204298740 j + 0.272398320 k    in R09 to R12

Note:

-You can execute  QRCL 05  or  QRCL E  to get cn ( q | m )  &  QRCL 09 or QRCL I  to get  dn ( q | m ).
 

Bessel Functions 1st Kind
 

Formula:         Jn(q) = (q/2)n0F~1( n+1 ; -q2/4 )
 

Data Registers:             R00 = n               ( Register R00 is to be initialized before executing "JNQ" )

                                         R01 thru R23: temp

Flags:   F01 & F09
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Jn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k                    with       R00 = n

Example:         n = PI     q = 1 + 2 i + 3 j + 4 k

   PI   STO 00

   4   ENTER^
   3   ENTER^
   2   ENTER^
   1   XEQ "JNQ"   >>>>   -11.22987515
                              RDN    -3.570801505
                              RDN    -5.356202248
                              RDN    -7.141603011

        JPI( 1 + 2 i + 3 j + 4 k ) = -11.22987515 - 3.570801505 i - 5.356202248 j - 7.141603011 k
 

Modified Bessel Functions 1st Kind
 

Formula:         In(q) = (q/2)n0F~1( n+1 ; q2/4 )
 

Data Registers:             R00 = n               ( Registers R00 is to be initialized before executing "INQ" )

                                         R01 thru R23: temp

Flags:   F01 & F09
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      In(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k                    with       R00 = n

Example:         n = PI     q = 1 + 2 i + 3 j + 4 k

   PI   STO 00

   4   ENTER^
   3   ENTER^
   2   ENTER^
   1   XEQ "INQ"   >>>>     0.315197912
                              RDN    -0.129988659
                              RDN    -0.194982989
                              RDN    -0.259977319

         IPI( 1 + 2 i + 3 j + 4 k ) = 0.315197912 - 0.129988659 i - 0.194982989 j - 0.259977319 k
 

Bessel Functions 2nd Kind
 

Formulae:

      Yn(q) = [ Jn(q) cos(n(pi)) - J-n(q) ] / sin(n(pi))      if n is not an integer.

-Otherwise:

  Yn(q) = -(1/pi) (q/2)-n SUMk=0,1,....,n-1  (n-k-1)!/(k!) (q2/4)k + (2/pi) ln(q/2) Jn(x) - (1/pi) (q/2)n SUMk=0,1,.....  ( psi(k+1) + psi(n+k+1) ) (-q2/4)k / (k!(n+k)!)

   where  psi = digamma function
 

Data Registers:        R00 = n                              ( Register R00 is to be initialized before executing "YNQ" )

                                     R01 thru R33: temp

Flags:  F01-F09
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Yn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k                    with       R00 = n

Example1:         n = PI     q = 1 + 2 i + 3 j + 4 k

   PI   STO 00

   4   ENTER^
   3   ENTER^
   2   ENTER^
   1   XEQ "YNQ"   >>>>     9.617564067
                               RDN    -4.171358813
                               RDN    -6.257038265
                               RDN    -8.342717668

   So,   YPI( 1 + 2 i + 3 j + 4 k ) = 9.617564067 - 4.171358813 i - 6.257038265 j - 8.342717668 k

Example2:         n = 3     q = 1 + 2 i + 3 j + 4 k

    3   STO 00

   4   ENTER^
   3   ENTER^
   2   ENTER^
   1   XEQ "YNQ"   >>>>    7.690027238
                              RDN    -5.224812676
                              RDN    -7.837219010
                              RDN   -10.44962535

       Y3( 1 + 2 i + 3 j + 4 k ) = 7.690027238 - 5.224812676 i - 7.837219010 j - 10.44962535 k
 

Modified Bessel Functions 2nd Kind
 

Formulae:

      Kn(q) = (pi/2) [ I-n(q) - In(q) ] / sin(n(pi))            if n is not an integer.

-Otherwise:

Kn(q) =  (1/2) (q/2)-n SUMk=0,1,....,n-1  (n-k-1)!/(k!) (-q2/4)k - (-1)n ln(q/2) In(x) + (1/2) (-1)n (q/2)n SUMk=0,1,.....  ( psi(k+1) + psi(n+k+1) ) (q2/4)k / (k!(n+k)!)

   where  psi = digamma function
 

Data Registers:        R00 = n                              ( Register R00 is to be initialized before executing "KNQ" )

                                     R01 thru R31: temp

Flags:  F01-F09
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Kn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k                    with       R00 = n

Example1:         n = PI     q = 1 + 2 i + 3 j + 4 k

   PI   STO 00

   4   ENTER^
   3   ENTER^
   2   ENTER^
   1   XEQ "KNQ"   >>>>      0.203067067
                                RDN    -0.055030892
                                RDN    -0.082546334
                                RDN    -0.110061785

  Whence,   KPI( 1 + 2 i + 3 j + 4 k ) = 0.203067067 - 0.055030892 i - 0.082546334 j - 0.110061785 k

Example2:         n = 3     q = 1 + 2 i + 3 j + 4 k

   3   STO 00

   4   ENTER^
   3   ENTER^
   2   ENTER^
   1   XEQ "KNQ"   >>>      0.208767798
                              RDN    -0.047983196
                              RDN    -0.071974793
                              RDN    -0.095966393

       K3( 1 + 2 i + 3 j + 4 k ) = 0.208767798 - 0.047983196 i - 0.071974793 j - 0.095966393 k
 

Parabolic Cylinder Functions
 

Formulae:            Dn(q) = Hn(q/sqrt(2))  2 -n/2  exp(-q2/4)    where   Hn(q) = Hermite function

Asymtotic Series:         Dn(q) ~ qn exp(-q2/4) 2F0( (1-n)/2 , -n/2 ;; -2/q2 )
 

Data Registers:        R00 = n                                 ( Register R00 is to be initialized before executing "DNQ" )

                                    R01 to R26:  temp
Flag:   F04

   CF 04 = Ascending series
   SF 04 = Asymptotic series.
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Dn(x+y.i+z.j+t.k) = x'+y'.i+z'.j+t'.k                  with       R00 = n

Example1:         n = 0.4      q = 1 + 1.1 i + 1.2 j + 1.3 k

  0.4   STO 00   CF 04

   1.3   ENTER^
   1.2   ENTER^
   1.1   ENTER^
    1     XEQ "DNQ"   >>>>    2.606105105
                                  RDN    -0.969116887
                                  RDN    -1.057218423
                                  RDN    -1.145319956

   D0.4( 1 + 1.1 i + 1.2 j + 1.3 k ) = 2.606105105 - 0.969116887 i - 1.057218423 j - 1.145319956 k

Example2:         n = PI       q = 1 + 2 i + 3 j + 4 k

  PI   STO 00    SF 04

  4    ENTER^
  3    ENTER^
  2    ENTER^
  1    XEQ "DNQ"   >>>>   -33138.87205
                               RDN      93318.15129
                               RDN      139977.2271
                               RDN      186636.3026

  DPI( 1 + 2 i + 3 j + 4 k ) =  -33138.87205 + 93318.15129 i + 139977.2271 j + 186636.3026 k

Notes:

-As usual with asymptotic series ( SF 04 ), infinte loops will occur if q is too small.
-However, they may also be used if q is relatively small when n is a positive integer.
 

Lambert-W Function
 

 QLW  solves  q = w(q) exp [ w(q) ]  by Newton's method.
 

Data Registers:  R00 thru R16: temp
Flags: /
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

     with      W(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k

Example:    q = 1 + 2 i + 3 j + 4 k

  4  ENTER^
  3  ENTER^
  2  ENTER^
  1  XEQ "QLW"   >>>>>    1.281459811
                               RDN      0.304051226
                               RDN      0.456076840
                               RDN      0.608102453

 Whence,    W(1 + 2 i + 3 j + 4 k) =  1.281459811 + 0.304051226 i + 0.456076840 j + 0.608102453 k
 

Regular Coulomb Wave Functions
 

Formulae:         FL(n,r) = CL(n) r L+1 Sum k>L AkL (n) r k-L-1

    with   CL(n) = (1/Gam(2L+2)) 2L e -pi.n/2 | Gam(L+1+i.n) |
    and    AL+1L = 1 ;  AL+2L = n/(L+1)  ;  (k+L)(k-L-1) AkL = 2n Ak-1L - Ak-2L   ( k > L+2 )

   | Gam( 1+i y ) |2 = (Pi.y) / Sinh (Pi y)   and

   | Gam( 1+L+i y ) |2 = [ L2 + y2 ] [ (L-1)2 + y2 ] .................. [ 1 + y2 ]  (Pi.y) / Sinh (Pi y)
 

Data Registers:             R09 = L            R10 = n           ( Registers R09 & R10 are to be initialized before executing "RCWFQ" )

                                         R00 to R17: temp
Flag:  F10 is cleared
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   with      FL( h ; x+y.i+z.j+t.k )  = x'+y'.i+z'.j+t'.k
 where   L is a non-negative integer and  h is a real number.

Example:         L = 2  ,   h = 0.75  ;   q = 1 + 1.2 i + 1.4 j + 1.6 k

     2     STO 09
  0.75   STO 10

  1.6     ENTER^
  1.4     ENTER^
  1.2     ENTER^
    1      XEQ "RCWFQ"   >>>>    -0.478767859
                                         RDN    -0.179694772
                                         RDN    -0.209643900
                                         RDN    -0.239593029

-So,   F2( 0.75 ;  1 + 1.2 i + 1.4 j + 1.6 k ) = -0.478767859 - 0.179694772 i - 0.209643900 j - 0.239593029 k
 

Struve Functions H
 

Formula:         Hn(q) = (q/2)n+1 1F~2( 1 ; 3/2 , n + 3/2 ; - q2/4 )           where  1F~2  is a regularized hypergeometric function.
 

Data Registers:             R00 = n          ( Register R00 is to be initialized before executing "STRVHQ" )

                                         R01 thru R23: temp
Flags:   F01 & F09
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Hn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k                    with       R00 = n

Example:         n = PI     q = 1 + 2 i + 3 j + 4 k

   PI   STO 00

   4   ENTER^
   3   ENTER^
   2   ENTER^
   1   XEQ "STRVHQ"   >>>>    8.546633217
                                      RDN    -4.085079932
                                      RDN    -6.127619898
                                      RDN    -8.170159859

   So,   Hpi( 1 + 2 i + 3 j + 4 k ) = 8.546633217 - 4.085079932 i - 6.127619898 j - 8.170159859 k
 

Struve Functions L
 

Formula:        Ln(q) = (q/2)n+1 1F~2( 1 ; 3/2 , n + 3/2 ;  q2/4 )           where  1F~2  is a regularized hypergeometric function.
 

Data Registers:             R00 = n          ( Register R00 is to be initialized before executing "STRVLQ" )

                                         R01 thru R23: temp
Flags:   F01 & F09
 
 
      STACK        INPUTS      OUTPUTS
           T             t            t'
           Z             z            z'
           Y             y            y'
           X             x            x'

   where      Ln(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k                    with       R00 = n

Example:         n = PI     q = 1 + 2 i + 3 j + 4 k

   PI   STO 00

   4   ENTER^
   3   ENTER^
   2   ENTER^
   1   XEQ "STRVLQ"   >>>>     1.864582785
                                      RDN    -0.116220198
                                      RDN    -0.174330298
                                      RDN    -0.232440397

  Whence,   Lpi( 1 + 2 i + 3 j + 4 k ) = 1.864582785 - 0.116220198 i - 0.174330298 j - 0.232440397 k
 

Spheroidal Wave Functions
 

-We assume that  | q | <= 1  and  Smn(q) is computed by      Smn(q) = ( 1 - q2 )m/2   Sumk=0,1,.... ak qk

     with     (k+1)(k+2) ak+2 - [ k ( k + 2m + 1 ) - Lmn + m ( m + 1 ) ] ak  - c2  ak-2 = 0
 

Data Registers:            R00 thru R08 & R13 thru R18: temp           ( Registers R09 thru R12 are to be initialized before executing "SWFQ" )

                                     R09 = m        R11 = c2
                                     R10 = n         R12 = Lmn
Flag:  F00

   CF 00       Flammer's Scheme      Smn(0) = Pmn(0)  &  S'mn(0) = P' mn(0)   where  Pmn(x)  =  Associated Legendre Functions of the first kind

      a0 = [ (-1)(n+m)/2 (n+m)! ] / [ 2n ((n-m)/2)! ((n+m)/2)! ]                   a1 = 0           if  n - m is even
      a1 = [ (-1)(n+m-1)/2 (n+m+1)! ] / [ 2n ((n-m-1)/2)! ((n+m+1)/2)! ]    a0 = 0            if  n - m is odd

   SF 00       Non-Normalized Functions

      a0 = 1          a1 = 0         if  n - m is even
      a0 = 0          a1 = 1         if  n - m is odd
 
 
      STACK        INPUTS      OUTPUTS
           T             t             t'
           Z             z             z'
           Y             y             y'
           X             x             x'

  with   Smn( x + y i + z j + t k ) = x' + y' i + z' j + t' k  ,   |  x + y i + z j + t k | <= 1

Example1:     m = 0 ,  n = 1 ,  c2 = 2 ,   Lmn = 3.172127920       These 4 numbers are to be stored in  R09 thru R12

   CF 00

   0.4   ENTER^
   0.3   ENTER^
   0.2   ENTER^
   0.1   XEQ "SWFQ"  >>>>    0.117347028                    HP-41 displays "FLAMMER"
                                    RDN     0.210312365
                                    RDN     0.315468547
                                    RDN     0.420624729

   S( 0.1 + 0.2 i + 0.3 j + 0.4 k ) = 0.117347028 + 0.210312365 i + 0.315468547 j + 0.420624729 k

Example2:     m = n = 2  ,  c2 = 3  ,  Lmn =  6.412006610           These 4 numbers are to be stored in  R09 thru R12

   CF 00

   0.4   ENTER^
   0.3   ENTER^
   0.2   ENTER^
   0.1   XEQ "SWFQ"  >>>>    4.058321831                    HP-41 displays "FLAMMER"
                                    RDN   -0.160164520
                                    RDN   -0.240246780
                                    RDN   -0.320329040

   S( 0.1 + 0.2 i + 0.3 j + 0.4 k ) = 4.058321831 - 0.160164520 i - 0.240246780 j - 0.320329040 k
 

Notes:

-If you set F00, you get the same result in the 1st example, and the same result divided by 3 in the 2nd example.
-The HP-41 displays "NON-NORMALZD" in this case.

 Lmn may be computed - for example - by "LMN" that is available in the JMB_MATH  module,
 but you can also use different values for this parameter.
 

Sylvester Equation
 

-"SYLV" solves the linear equation:   a q + q b = c   where   a , b , c  are 3 given quaternions and  q  is unknown.
 

Formulae:

     q = ( 2 Ré(b) + a + | b |2 a -1 ) -1  ( c + a -1 c Conj(b) )

     q = ( c + Conj(a) c b -1 ) ( 2 Re(a) + b + | a |2 b -1 ) -1
 

-To get a better accuracy, the 1st formula is used if  | b | <= | a |
  and the second one if  | b | > | a |
 

Data Registers:               R00: temp                       ( Registers R01 thru R12 are to be initialized before executing "SYLV" )

                                        R01 = a1                R05 = b1                R09 = c1              R13 thru R21:  temp
                                        R02 = a2                R06 = b2                R10 = c2
                                        R03 = a3                R07 = b3                R11 = c3
                                        R04 = a4                R08 = b4                R12 = c4
Flag:  F00 is cleared
 
 
      STACK        INPUTS      OUTPUTS
           T             /             t
           Z             /             z
           Y             /             y
           X             /             x

  With   q = x + y i + z j + t k

Example:      Solve     ( 2 - 3 i + 4 j - 7 k ) q + q ( 3 + 4 i - 5 j + 6 k ) = 1 + 2 i - 3 j + 4 k

   2   STO 01     3   STO 05     1   STO 09
  -3  STO 02     4   STO 06     2   STO 10
   4   STO 03   -5   STO 07    -3   STO 11
  -7  STO 04     6   STO 08     4   STO 12

  or   -7  ENTER^   4  ENTER^  -3  ENTER^  2   QSTO A
          6  ENTER^  -5  ENTER^  4  ENTER^  3   QSTO E
          4  ENTER^  -3  ENTER^  2  ENTER^  1   QSTO I

  XEQ "SYLV"    >>>>   0.239980450
                           RDN    0.418866080
                           RDN   -0.576246334
                           RDN    0.795210166

-Thus,    q = 0.239980450 + 0.418866080 i  -  0.576246334 j  +  0.795210166  k
 

Unilateral Linear Systems
 

-QLS solves unilateral systems of the form  B = A Q  or  B = Q A
-Gaussian elimination with partial pivoting is used.
-The square of the norm of the product of the pivots is also computed and returned in X-register and R00.
-There is a unique solution if and only if it is different from zero.

-Several unilateral systems with the same matrix A may also be solved at the same time.
-Underdetermined & overdetermined systems too.
-Actually, the program is structured like "LS3" in JMB_MATH module.

-For 1 system of n equations in n unknowns, re-write it as:
 

            b1 = a1,1 q1 + ................ + a1,n qn

            ...................................................                          where  bi , qi , ai,j  are quaternions,  qi are the unknowns

            bn = an,1 q1 + ................ + an,n qn
 

-More generally, the coefficients form a nxm matrix P = [ pi,j ]  1 <= i <= n , 1 <= j <= m , which are to be stored in column order from register R01.
 

Data Registers:              R00 = | ProdPiv |2                   ( Registers R01 thru R4nm are to be initialized before executing "QLS" )

                                         R01 to R04 = p1,1        ...................................       R4n(m-1)+1 to R4n(m-1)+4 = p1,m
                                          ........................................................................................................................

                                        R4n-3 to R4n = pn,1      ....................................       R4n.m-3 to R4n.m = pn,m

Flags:   F00-F01-F02

  CF 00   A pivot is regarded as 0 if it's smaller than   E-7
  SF 00   Only 0 is regarded as 0

  CF 01   Partial pivoting
  SF 01   No Pivoting ( not recommended in general, with a few exceptions however )

  CF 02   B = A Q   is solved
  SF 02   B = Q A   is solved
 
 
      STACK        INPUTS      OUTPUTS
           Y             n             /
           X             m      | ProdPiv |2

  where   n = number of rows , m = number of columns
     &      | ProdPiv |2  is the square of the magnitude of the product of all the pivots.

Example1:       Find 3 quaternions  u , v , w  such that

    1 + 2 i + 3 j + 4 k = ( 2 + 3 i + j + 2 k ) u + ( 1 + 3 i + 4 j + 2 k ) v + ( 3 + 4 i + 2 j + k ) w
    2 +  i  + 4 j + 3 k  = ( 1 + 2 i + j + 3 k ) u + ( 2 +  i  + 3 j + 2 k ) v + ( 2 + 3 i + 4 j + 5 k ) w
    3 + 4 i + 2 j +  k  =  ( 3 + 4 i + 5 j + 6 k ) u + ( 4 + 2 i + 3 j + 4 k ) v + ( 1 + 2 i + 3 j + 4 k ) w

-Store these 48 integers in the same order as above into

  R01  R02  R03  R04     R13  R14  R15  R16    R25  R26  R27  R28    R37  R38  R39  R40
  R05  R06  R07  R08     R17  R18  R19  R20    R29  R30  R31  R32    R41  R42  R43  R44
  R09  R10  R11  R12     R21  R22  R23  R24    R33  R34  R35  R36    R45  R46  R47  R48

-This block of registers is moved during the calculation, so execute a SIZE 4n.m+18  or greater ( here, at least SIZE 066 ).

  CF 00   CF 01  CF 02

-There are 3 rows & 4 columns, so:

  3   ENTER^
  4   XEQ "QLS"   >>>>   19782 = R00

-Since R00 # 0 , there is a unique solution ( u , v , w ) in registers R01 thru R12:

  u = -0.398544130 + 0.660903352 i - 0.558993024 j + 0.824992418 k
  v =  1.079162875 - 0.798604792 i + 0.672732788 j - 0.722171672 k              ( rounded to 9D )
  w = 0.462996664 - 0.098422809 i - 0.377767667 j - 0.104185623 k

-Registers R13 to R48 now contain the Identity matrix +/- very small roundoff-errors sometimes...
 

Example2:    If you want to solve the similar system ( same coefficients but different order for the products )

    1 + 2 i + 3 j + 4 k = u ( 2 + 3 i + j + 2 k )  + v ( 1 + 3 i + 4 j + 2 k ) + w ( 3 + 4 i + 2 j + k )
    2 +  i  + 4 j + 3 k  = u ( 1 + 2 i + j + 3 k )  + v ( 2 +  i  + 3 j + 2 k ) + w ( 2 + 3 i + 4 j + 5 k )
    3 + 4 i + 2 j +  k  =  u ( 3 + 4 i + 5 j + 6 k ) + v ( 4 + 2 i + 3 j + 4 k ) + w ( 1 + 2 i + 3 j + 4 k )

-Store the same coefficients into the same registers ( R01 thru R48 )

   SF 02

   3   ENTER^
   4      R/S           >>>>    45542 = R00 # 0

-The solution is also unique and we find in registers  R01 thru R12

  u = -0.204075359 - 0.409029028 i + 0.185499100 j - 0.192569496 k
  v =  0.781696017 + 0.824864960 i - 0.006762988 j - 0.118791445 k              ( rounded to 9D )
  w = 0.170721532 - 0.247441922 i - 0.104189539 j + 0.300623600 k

Note:

-If all the coefficients of the matrix A are integrers, X-output = R00 will be an integer too ( may be with roundoff-errors ).
 

Equation F(q) = 0
 

-"FQ=0" uses the secant method.

-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 "FQ=0" )

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

-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

-This subroutine may be simplified by replacing  STO 09  RDN  STO 10  RDN  STO 11  X<>Y  STO 12  by  QSTO 09  ...

   "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 "FQ=0"  >>>>  yields after many iterations ( the successive x-approximations are displayed ):

          0.808540975   RDN  -1.290564602   RDN  -0.290931376    RDN   0.490521462      LASTX  gives  | f(q) | = 0.000000009

-So, a solution is   q =  0.808540975 - 1.290564602 i  - 0.290931376 j + 0.490521462 k
 

Equation F(q) = q
 

-"FQ=Q"  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 "FQ=Q" )

        R01 = a                        R14 to R17 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 , R14 thru R17 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 R18 or greater can be used ( and altered ) to compute f(q)
-Registers R13 thru R17 can't be modified but can be used too.
 
 
      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 )    or   CLST  QSTO A

    FIX 9
    XEQ "FQ=Q"   The successive x-approximations are displayed,  and eventually:

       it yields                  0.808540975
                        RDN  -1.290564598
                        RDN  -0.290931375
                        RDN   0.490521465

   whence   q =  0.808540975 - 1.290564598 i  - 0.290931375 j + 0.490521465 k  is a root of this polynomial.

Notes:

-The accuracy is controlled by the display settings:  FIX 9 may lead to an infinite loop. Try FIX 8 ...
-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.
 

Systems of n Equations in n unknowns  Fi(q1,....,qn) = qi
 

-"FNQ=Q" uses the same method as above:

-If F is a contraction mapping on a complete metric space, then the equation  F ( X ) = X  has a unique solution
 which is the limit of the sequence:  X(k+1) = F(X(k))  where X(0) is arbitrary.

-This theorem is used to solve a system of n equations in n unknowns which must be re-written under the form:
 

      q1 = f1 ( q1 , ... , qn )
      q2 = f2 ( q1 , ... , qn )
       .............................

      qn = fn ( q1 , ... , qn )
 

Data Registers:          R20 = n = the number of unknowns    ( Registers R20 thru R20+5n are to be initialized before executing "FNQ=Q" )

                                     R21 to R24 = q1                         R21+4n = f1 name            ( q1 , ... , qn )  is an initial guess that you have to choose.
                                     R25 to R28 = q2                         R22+4n = f2 name
                                       ...............................................................................             R00 to R16 are unused

                                     R17+4n to R20+4n = qn             R20+5n   = fn name         R17-R18-R19: temp

Flags:  /
Subroutines:      The n functions   fi  computing   fi ( q1 , ... , qn )  in registers X Y Z T assuming  q1 , ...... ,  qn  are in R21 to R20+4n
 
 
      STACK        INPUTS      OUTPUTS
           T             /             t
           Z             /             z
           Y             /             y
           X             /             x

  with  q = x + y i + z j + t k  = one solution

Example:      Find 3 quaternions  u , v , w  such that

   u = ArcSinh ( u v w )
   v = ( u + v + w2 ) 1/3
   w = ( 1 + 2 i + 3 j + 4 k + u + v + w ) 1/2

-We have to load 3 subroutines - for instance:
 
 
  01  LBL "U"
  02  RCL 21
  03  STO 01         
  04  RCL 22
  05  STO 02 
  06  RCL 23
  07  STO 03 
  08  RCL 24 
  09  STO 04
  10  RCL 25
  11  STO 05 
  12  RCL 26 
  13  STO 06
  14  RCL 27
  15  STO 07
  16  RCL 28
  17  STO 08
  18  XEQ "Q*Q"
  19  STO 01
  20  RDN
  21  STO 02
  22  RDN
  23  STO 03
  24  X<>Y
  25  STO 04
  26  RCL 29
  27  STO 05         
  28  RCL 30 
  29  STO 06
  30  RCL 31
  31  STO 07
  32  RCL 32
  33  STO 08 
  34  XEQ "Q*Q"
  35  XEQ "ASHQ"
  36  RTN
  37  LBL "V"
  38  2
  39  STO 00
  40  RCL 32
  41  RCL 31
  42  RCL 30
  43  RCL 29
  44  XEQ "Q^R"
  45  STO 01
  46  RDN
  47  STO 02
  48  RDN
  49  STO 03         
  50  CLX
  51  RCL 24 
  52  +
  53  RCL 28 
  54  +
  55  RCL 03
  56  RCL 27 
  57  +
  58  RCL 23 
  59  +
  60  3
  61  1/X
  62  STO 00
  63  CLX
  64  RCL 21
  65  RCL 25
  66  +
  67  ST+ 01
  68  CLX
  69  RCL 22
  70  RCL 26
  71  +
  72  ST+ 02
  73  X<> 02
  74  RCL 01         
  75  XEQ "Q^R"
  76  RTN
  77  LBL "W"
  78  RCL 21
  79  RCL 25 
  80  +
  81  RCL 29 
  82  +
  83  1 
  84  +
  85  STO 01
  86  RCL 24
  87  RCL 28
  88  +
  89  RCL 32
  90  +
  91  4
  92  +
  93  RCL 23
  94  RCL 27
  95  +
  96  RCL 31         
  97  +
  98  3
  99  +
100  RCL 22 
101  RCL 26
102  +
103  RCL 30 
104  +
105  .5
106  STO 00
107  1/X
108  +
109  RCL 01
110  XEQ "Q^R"
111  RTN
112  END

 ( These subroutines may be simplified by replacing lines 02 to 17 with  QRCL  21  QSTO  1  QRCL  25  QSTO  5   and similar improvements... )

     n = 3   STO 20   and if we choose  u = v = w = 1 as initial guesses:

     1  STO 21   STO 25   STO 29
     0  STO 22   STO 26   STO 30
     0  STO 23   STO 27   STO 31
     0  STO 24   STO 28   STO 32

( much simpler is  CLST  1  QSTO 21  QSTO 25  QSTO 29 )

     "U"  ASTO 33    "V"   ASTO 34   "W"   ASTO 35

  FIX 9   XEQ "FNQ=Q"   The HP-41 displays the successive approximations and finally, it yields:

                                   >>>>    4.541667578 = R21
                                   RDN     0.195336787 = R22
                                   RDN     0.293005181 = R23
                                   RDN     0.390673575 = R24

-So,  u = 4.541667578 + 0.195336787 i + 0.293005181 j + 0.390673575 k

-Likewise, we find in registers  R25 thru R28 & R29 thru R32:

         v = 2.725387904 + 0.137399822 i + 0.206099734 j + 0.274799645 k     (  QRCL 25 )

        w = 3.590282794 + 0.377430929 i + 0.566146393 j + 0.754861857 k     (  QRCL 29 )
 

Notes:

-The accuracy is controlled by the display settings.
-FIX 9 may lead to infinite loops:  FIX 8 or FIX 7 are often better.
-Press  XEQ 01  if you want to improve the accuracy then.
 

Quaternionic Discrete Fourier Transform
 

-"QDFT" calculates the discrete Fourier Transform of a vector or a matrix of quaternions.

-Due to the non-commutativity of the multiplication, many Fourier Transform may be defined, here we use the right-sided formulae:
 

1-Dimension:     F (p) = ( 1/sqrt(M) ) SUMm=0,1,....,M-1  f(m) exp ( -2PI µ m p / M )

2-Dimensions:    F (p,q) = ( 1/sqrt(M.N) ) SUMm=0,1,....,M-1  SUMn=0,1,....,N-1   f(m,n) exp [ -2PI µ { ( m p / M ) + ( n q / N ) } ]

    where  µ = µ' i + µ" j + µ''' k    is unitary and purely imaginary:  | µ | = 1 & Re(µ) = 0

-The inverse transform is simply obtained by replacing  -2PI  by +2PI in the above formulae.
 

Data Registers:             R00 = M.NNN        ( Registers R00 & R18 thru R20+4M.N are to be initialized before executing "QDFT" )

                                          R01 thru R17: temp    R09 to R12: partial sums

                                        R18 = µ'            R21 to R24 = f(0,0)                          ....................            R21+4M(N-1) to R24+4M(N-1) = f(0,N-1)
                                        R19 = µ"             ..............................                           ....................             ...............................................................
                                        R20 = µ'''          R17+4M to R20+4M = f(M-1,0)      ....................            R17+4M.N to R20+4M.N = f(M-1,N-1)
Flags:   F00 & F01

   CF 00   All the coefficients of the DFT are computed and stored into  R21+4M.N ................... R20+8M.N
   SF 00   Only 1 coefficient  F(p,q)  is computed and returned in the stack.

   CF 01 =  Fourier Transform
   SF 01 =  Inverse Fourier Transform
 
 
 SF 00  STACK        INPUTS      OUTPUTS
           T             /             t
           Z             /             z
           Y             q             y
           X             p             x

      with  F(p,q)  = x + y i + z j + t k       0 <= p < M  ,  0 <= q < N

Example:     µ = ( i + 2 j + 4 k ) / sqrt(21)    and     Matrix A =

   ( 1 + 2 i + 3 j + 4 k )   ( 2 + 3 i + j + 2 k )     ( 1 + 3 i + 4 j + 2 k )     ( 3 + 4 i + 2 j + k )
   ( 2 +  i  + 4 j + 3 k )   ( 1 + 2 i + j + 3 k )      ( 2 +  i  + 3 j + 2 k )      ( 2 + 3 i + 4 j + 5 k )
   ( 3 + 4 i + 2 j +  k )    ( 3 + 4 i + 5 j + 6 k )    ( 4 + 2 i + 3 j + 4 k )     ( 1 + 2 i + 3 j + 4 k )

-M = 3  ,  N = 4   so  3.004  STO 00

-Key in    1   STO 18   2  STO 19    4   STO 20      21  SQRT    ST/ 18  ST/ 19  ST/ 20     ( the program does not check that µ is unitary )

-Store the 12 quaternions ( 48 real numbers ) into R21 thru R68 as shown below:

  R21  R22  R23  R24     R33  R34  R35  R36    R45  R46  R47  R48    R57  R58  R59  R60
  R25  R26  R27  R28     R37  R38  R39  R40    R49  R50  R51  R52    R61  R62  R63  R64
  R29  R30  R31  R32     R41  R42  R43  R44    R53  R54  R55  R56    R65  R66  R67  R68

-To compute, for instance,  F(2;3)

        SF 00   CF 01

    3  ENTER^
    2  XEQ "QDFT"  >>>>   -0.182134053
                                RDN   -0.119263223
                                RDN    1.977063795
                                RDN    1.221466790

  Whence    F(2;3) = -0.182134053 - 0.119263223 i + 1.977063795 j + 1.221466790 k
 

Notes:

-If N = 1 ( vector )  you may also store  M  in R00  instead of  M.001  and  q  is arbitrary  provided it is not an alpha string.
-Register R17 is unused if flag F00 is set.
 
 
 CF 00  STACK        INPUTS      OUTPUTS
           X             /        bbb.eee

      where bbb.eee is the control number of the DFT

Example:     With the same µ and A, which are still in the same registers  ( SIZE 117 or greater )

   CF 00  CF 01  XEQ "QDFT"   >>>>   69.116   the DFT is in registers R69 to R116 at the end

-The first column ( in registers R69 to R80 ) is

  7.217 + 8.949 i + 10.104 j + 10.681 k
 -1.346 + 0.940 i - 1.267 j - 0.080 k                    ( rounded to 3D )
  0.241 + 0.503 i - 0.176 j - 2.807 k

-The last element in registers R113 to R116 - which has been computed earlier -
  is   F(2;3) = -0.182134053 - 0.119263223 i + 1.977063795 j + 1.221466790 k

Notes:

-If N = 1 ( vector ),  you may also store  M  in R00  instead of  M.001

-May be you will want to calculate the Inverse DFT from these results:
-Key in  69.021048  REGMOVE  SF 01  R/S  >>>>  69.116
-Compare the numbers in R69 to R116 to the initial matrix to get an idea of the roundoff-errors...
 

References:

[1]   Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]   http://functions.wolfram.com/
[3]  Georgi Georgiev, Ivan Ivanov, Milena Mihaylova, Tsvetelina Dinkova - "An algorithm for solving a Sylvester quaternion equation"
[4]  Drahoslava Janovsk´a, Gerhard Opfer - "Linear equations in quaternionic variables"
[5]  Patrice Denis - Thèse - "Quaternions et Algèbres Géométriques, de nouveaux outils pour les images numériques couleur" ( in French )