hp41programs

Anionic Special Functions

Anionic Special Functions(I) for the HP-41


Overview
 

 0°)  A Few Subroutines

  0-a)  Focal Programs
  0-b)  M-Code Routines: ST*A  ST/A  DSA  DS*A  AMOVE ASWAP

 1°)  Gamma Function
 2°)  Digamma Function
 3°)  Error Function
 4°)  Exponential Integral, Sine & Cosine Integral

   a) Program#1
   b) Em(a)
   c) An Asymptotic Expansion for Em(a)

 5°)  Kummer's Functions
 6°)  Parabolic Cylinder Function

   a)  Via Kummer's Function
   b)  Asymptotic Expansion

 7°)  Laguerre's Functions
 8°)  Hypergeometric Functions
 9°)  Associated Legrendre Functions

  9-a)  Via Hypergeometric Functions
  9-b)  Functions of the 1st kind - Integer Order and Degree

10°) Struve Functions
11°) Bessel Functions
12°) Regular Coulomb Wave Function
13°) UltraSpherical Functions
14°) Jacobi Functions
15°) Toronto Functions
16°) Weber & Anger Functions
 

-This page generalizes most of the programs listed in "Quaternionic Special Functions for the HP-41"
-They can be used with complexes, quaternions, octonions, sedenions, 32-ons ... and so on ...  provided there are enough registers !

-Always store n = 2 ( complexes ) or 4 ( quaternions ) or 8 ( octonions ) .... in register R00
-The variable itself is to be stored in registers  R01 thru Rnn
-If there are parameters - restricted to real values - they must be stored in the stack.

-In all the examples - except one in §12 - the arguments are quaternions ( n = 4  STO 00 ).
-But the execution times will increase like n increases.
-The required SIZE is about 3n+1 or 4n+1, so these routines may be used for 64-ons but not for 128-ons...
  unless many components of the imaginary part equal zero.

-Most of these programs employ ascending series which return accurate results for "small" arguments only.

-The execution times have been measured with the M-Code routines listed in paragraph 0°)b) below and in "Anions for the HP-41" paragraph 0°)c)
-Otherwise, the programs would run much slower...
 

0°)  A Few Subroutines
 

     a) Focal Programs
 

-Several subroutines are useful in the programs that are listed on this page:

  "ST*A"  multiplies registers  R01 thru Rnn  by  the content of X-register
  "ST/A"    divides   registers  R01 thru Rnn  by  the content of X-register

  "DSA"   adds registers  R01 thru Rnn  to the registers  R2n+1 thru R3n
               and returns in X the difference between the previous and the new contents of registers R2n+1 thru R3n

-Useful to compute an "infinite" sum:   if this difference equals zero, the sum is assumed to be correct...
  except in §12: a flag is also used to avoid a premature stop.
 
 
 

 01  LBL "ST*A"
 02  RCL 00
 03  RDN
 04  LBL 01
 05  ST* IND T 
 06  DSE T
 07  GTO 01
 08  RTN
 09  LBL "ST/A" 
 10  RCL 00
 11  RDN
 12  LBL 02
 13  ST/ IND T 
 14  DSE T
 15  GTO 02
 16  RTN
 17  LBL "DSA" 
 18  RCL 00
 19  STO Q
 20  PI
 21  INT
 22  *
 23  ENTER^
 24  CLX
 25  LBL 03
 26  RCL IND Q
 27  RCL IND Z
 28  +
 29  STO IND Z
 30  LASTX
 31  -
 32  ABS
 33  +
 34  DSE Y
 35  DSE Q
 36  GTO 03      
 37  END

 
     ( 72 bytes )
 

Notes:

-All these subroutines assume that  R00 = n

-Of course, the programs will run faster if these loops are inserted in the programs themselves.
-But using M-Code routines is even better !
 

     b) M-Code Routines
 

 INITIALIZATION:    this subroutine is already listed in "Anions for the HP-41" paragraph 0°)c)1)
 

260  SETHEX         @EFC8  in my ROM
378  C=c
03C  RCR 3
106  A=C S&X
130  LDI S&X
200  200h                           correct value for an HP-41CV/CX or an HP-41C with a QUAD memory module
306  ?A<C S&X
381  ?NCGO
00A  NEXIST
0A6  A<>C S&X
270   RAMSELECT
038   READATA
38D  ?NCXQ
008   N->S&X
2E6  ?C#0 S&X
0B5  ?NCGO
0A2  ERROR
05A  C=0 M
23A  C=C+1 M
106  A=C S&X
1BC  RCR 11
0A6  A<>C S&X
1E6  C=C+C S&X
10E  A=C ALL
378  C=c
03C  RCR 3
0E6  B<>C S&X
378  C=c
0C6  C=B S&X
1BC  RCR 11
0C6  C=B S&X
20E  C=A+C ALL
24C  ?FSET 9
013  JNC+02
03C  RCR 3
070  N=C ALL
10E  A=C ALL
130  LDI S&X
200  200h
306  ?A<C S&X
381  ?NCGO
00A  NEXIST
3E0  RTN         @EFF2  in my ROM

  ( 43 words )
 

-To call this subroutine:    ?NCXQ  EFC8  =  321   3BC    in my ROM
-Change these 2 words if you have loaded this routine at another address...
 

ST*A ,  ST/A  ,  DSA  ,  DS*A
 

081  "A"
02A  "*"
014  "T"
013  "S"
284   CLRF 7
033   JNC+06
081  "A"
02F  "/"
014  "T"
013  "S"
288   SETF 7
248   SETF 9
321   ?NCXQ                               calls the "initialization routine" above
3BC   EFC8                                  change these 2 words according to your own ROM
2A0   SETDEC                             LOOP
0B0   C=N ALL
270    RAMSLCT
038    READATA
10E    A=C ALL
046    C=0 S&X
270    RAMSLCT
0F8    C=X
28C   ?FSET 7
261    ?CXQ
061    A/C
28C   ?FSET 7
135    ?NCXQ
060    A*C
10E   A=C ALL
0B0   C=N ALL
270   RAMSLCT
0AE  A<>C ALL
2F0   WRITDATA
260   SETHEX
0B0   C=N ALL
266   C=C-1 S&X
070   N=C ALL
106   A=C S&X
03C  RCR 3
306   ?A<C S&X
333   JNC -26d                            GOTO LOOP
3E0   RTN

 ( 41 words )
 

-It's often useful to multiply the numbers in registers R01 thru Rnn ( without modifying these registers ) by a real number x
  before adding them to the partial sum in registers R2n+1 thru R3n
-DS*A does this job.
-Place x in X-register, DS*A  saves x in L-register
 

081   "A"
02A  "*"
013   "S"
004   "D"
284   CLRF 7
02B   JNC+05
081   "A"
013   "S"
004   "D"
288   SETF 7
244   CLRF 9
321   ?NCXQ                               calls the "initialization routine" above
3BC   EFC8                                  change these 2 words according to your own ROM
0B0   C=N ALL
106    A=C S&X
03C   RCR 3
146   A=A+C S&X
1BC  RCR 11
11A  A=C M
378   C=c
03C  RCR 3
1C6  A=A-C S&X
130   LDI S&X
200  200h                                     correct value for an HP-41CV/CX or an HP-41C with a QUAD memory module
306  ?A<C S&X
381  ?NCGO                               tests if R3n does exist
00A  NEXIST
0AE  A<>C ALL
070   N=C ALL
0F8   C=X
128   L=C
04E   C=0 ALL
0E8   X=C
2A0   SETDEC                                 LOOP
138   C=L
10E   A=C ALL
0B0   C=N ALL
03C  RCR 3
270   RAMSLCT
038   READATA
28C  ?FSET 7
135   ?NCXQ
060   C=A*C
10E   A=C ALL
0B0   C=N ALL
270   RAMSLCT
038   READATA
01D   C=
060   A+C
10E   A=C ALL
0B0   C=N ALL
270   RAMSLCT
038   READATA
2BE  C=-C
0AE  A<>C ALL
2F0   WRITDATA
01D   C=
060   A+C
05E   C=| C |
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
0F8   C=X
01D  C=
060  A+C
0E8  X=C
260  SETHEX
0B0  C=N ALL
266   C=C-1 S&X
27A  C=C-1 M
070   N=C ALL
03C  RCR 3
106  A=C S&X
03C  RCR 3
306  ?A<C S&X
2B3  JNC-42d                               GOTO LOOP
3E0   RTN

 ( 77 words )
 

-These M-Code routines are to be used in the same way as the corresponding focal programs  "ST*A"  "ST/A"  "DSA"  listed in praragraph 0°)a)
-They run of course much faster...
 

COPYING / SWAPPING 2 ANIONS:
 

-To copy an anion from registers R01 thru Rnn to registers Rn+1 thru R2n, we can use the following instructions, assuming R00 = n:

   RCL 00     +        +
   .1               E3     REGMOVE
   %              /
   ISG X       1

-The following M-Code routines  AMOVE & ASWAP  can be used instead.
-In the above example, simply key in  12  AMOVE
-To swap 2 anions in registers Rn+1 thru R2n & R2n+1 thru R3n,  key in  23   ASWAP  ( or 32  ASWAP )
-It works with an HP-41CX
 

085   "E"
016   "V"
00F   "O"
00D  "M"
001   "A"
044   CLRF 4
03B   JNC+07
090   "P"
001   "A"
017  "W"
013  "S"
001  "A"
048  SETF 4
0F8  C=X
361  ?NCXQ
050   chk alpha & SETDEC
128   L=C
00E   A=0 ALL
35C  PT=12
102  A=C @PT
1A2  A=A-1 @PT
0AE  A<>C ALL
38E   RSHFA
38E   RSHFA
25C  PT=9
1A2  A=A-1 @PT
0A2  A<>C @PT
15C  PT=6
222   C=C+1 @PT
070   N=C ALL
378   C=c
03C  RCR 3
106   A=C S&X
130   LDI S&X
200   200h                                  200h is correct for an HP41-CX
306   ?A<C S&X
381   ?NCGO
00A   NONEXISTENT              Cheks that R00 does exist
0A6   A<>C S&X
270    RAMSLCT
038    READATA
361  ?NCXQ
050   chk alpha & SETDEC
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
0B0  C=N ALL
135   C=
060  A*C
04E  C
35C
050   =
25C
050  1.001
025  C=
060  AB+C
0E8  X=C
311  ?NCGO                                 the last 2 words are correct for an HP-41CX
0FA  REGMOVE/REGSWAP       they should be modified with an X-Functions module.

( 59 words )
 
 

      STACK        INPUTS      OUTPUTS
           X           pq     bbb.eeennn
           L            /            pq

   where the control number  bbb.eeennn  ( needed for REGMOVE & REGSWAP ) is calculated by

  bbb = n ( p-1 ) + 1  ,   eee = n ( q-1 ) + 1     assuming that  R00 = n

Notes:

-In fact, only the first 2 digits of the mantissa are taken into account.
-I've not used AMOVE or ASWAP in the programs below ( except "PMNA" ).
-But many bytes could be saved if you employ these functions...

-Here is an ( almost ) equivalent focal program:
 
 

 01  LBL "AMOVE"
 02  CF 05
 03  GTO 00
 04  LBL "ASWAP"
 05  SF 05
 06  LBL 00
 07  INT
 08  11
 09  -
 10  10
 11  /
 12  INT 
 13  ST- L
 14   E2
 15  ST/ L
 16  X<> L
 17  +
 18   E-6
 19  +
 20  RCL 00
 21  *
 22  1.001
 23  +
 24  FC? 05
 25  REGMOVE
 26  FS?C 05
 27  REGSWAP
 28  END
 
 

1°)  Gamma Function
 

-Here, the Gamma function is computed by a continued fraction:

    Gam(a) =  exp [ (a-1/2) Ln a + Ln (2.PI)1/2  - a + ( 1/12 )/( a + ( 1/30 )/( a + ( 53/210)/( a + (195/371)/( a + ...  )))) ]
 

-The relation   Gam(a+1) = a Gam(a)  is used recursively if Re(a) < 5  until   Re(a+1+.......+1)  > 5
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "GAMA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R4n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flag:  F04
Subroutines:  "A*A1"  "1/A"   "LNA"  "E^A"  "ST*A"
 
 
 

  01  LBL "GAMA"
  02  RCL 00
  03  .1
  04  %
  05  ISG X
  06  +
  07   E3
  08  /
  09  1
  10  +
  11  REGMOVE 
  12  RCL 00
  13   E3
  14  /
  15  ISG X
  16  CLRGX
  17  SIGN
  18  STO 01
  19  LBL 00
  20  XEQ "A*A1"
  21  RCL 00
  22  1
  23  ST+ Y
  24  ST+ IND Y
  25  RCL IND Y
  26  5
  27  X>Y?
  28  GTO 00
  29  RCL 00
  30   E3
  31  /
  32  STO Y
  33  ISG X
  34  RCL 00
  35  3
  36  *
  37  +
  38   E3
  39  /
  40  1
  41  +
  42  REGMOVE 
  43  X<>Y
  44  ENTER^
  45  ST+ Y
  46  ISG X
  47   E3
  48  /
  49  RCL 00
  50  +
  51  1
  52  +
  53  REGMOVE
  54  +
  55  REGMOVE
  56  XEQ "LNA"
  57  SIGN
  58  RCL 00
  59  +
  60  RCL IND X
  61  STO M
  62  .5
  63  ST- IND Z
  64  XEQ "A*A1"
  65  RCL 00
  66  .1
  67  %
  68  ISG X
  69  +
  70  RCL M
  71  STO IND Y
  72  CLX
  73  RCL 00
  74  +
  75   E3
  76  /
  77  1
  78  +
  79  REGSWAP 
  80  371
  81  195
  82  /
  83  CF 04
  84  XEQ 01
  85  210
  86  53
  87  /
  88  XEQ 01
  89  30
  90  XEQ 01
  91  12
  92  SF 04
  93  LBL 01
  94  XEQ "ST*A"
  95  XEQ "1/A"
  96  RCL 00
  97  ENTER^
  98  ST+ Y
  99  LBL 02
100  RCL IND Y 
101  FS? 04
102  CHS
103  ST+ IND Y
104  RDN
105  DSE Y
106  DSE X
107  GTO 02
108  FC?C 04
109  RTN
110  PI
111  ST+ X
112  SQRT
113  LN
114  ST+ 01
115  RCL 00
116  3
117  *
118  RCL 00
119  XEQ 02
120  XEQ "E^A"
121  RCL 00
122  .1
123  %
124  STO Z
125  ST+ Z
126  ISG X
127  +
128   E3
129  /
130  1
131  +
132  REGMOVE 
133  +
134  REGSWAP
135  XEQ "1/A"
136  XEQ "A*A1"
137  CLA
138  END

 
      ( 252 bytes / SIZE 4n+1 )
 
 

      STACK        INPUT      OUTPUT
           X             /         1.nnn

   where   1.nnn   is the control number of the result.

Example:     With the quaternion         q = 4 + 3 i + 2 j + k

-For quaternions,  n = 4  STO 00

      4  STO 01
      3  STO 02
      2  STO 03
      1  STO 04

      XEQ "GAMA"   >>>>    1.004                     ---Execution time = 22s---

            R01 =  0.541968820
            R02 = -0.740334196
            R03 = -0.493556131
            R04 = -0.246778065

     Gam ( 4 + 3 i + 2 j + k ) =  0.541968820 - 0.740334196 i - 0.493556131 j - 0.246778065 k
 

2°)  Digamma Function
 

Formula:             Psi(a) ~ Ln a - 1/(2a) -1/(12a2) + 1/(120a4) - 1/(252a6) + 1/(240a8)    is used if  Re(a) > 8

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

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "PSIA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R4n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:  "A*A1"  "1/A"   "LNA"  "E^A"  "ST/A"

-Line 36 is a synthetic  TEXT0  instruction. It can be replaced by  ABS  ( or another  NOP  ).
 
 

  01  LBL "PSIA"
  02  RCL 00
  03  .1
  04  %
  05  STO Z
  06  ISG X
  07  +
  08   E3
  09   /
  10  1
  11  +
  12  REGMOVE 
  13  STO M
  14  +
  15  REGMOVE 
  16  RCL 00
  17  4
  18  *
  19   E3
  20  /
  21  ISG X
  22  RCL 00
  23  3
  24  *
  25  +
  26  STO N
  27  CLRGX
  28  LBL 00
  29  XEQ "1/A"
  30  RCL N
  31  LBL 01
  32  RCL IND Y
  33  ST- IND Y
  34  RDN
  35  ISG Y
  36  ""
  37  ISG X
  38  GTO 01
  39  SIGN 
  40  ST+ IND Y
  41  X<>Y
  42  INT
  43  RCL 00
  44   E3
  45  /
  46  ISG X
  47  LASTX
  48  / 
  49  +
  50  REGMOVE 
  51  RCL IND X
  52  8
  53  X>Y?
  54  GTO 00
  55  RCL N
  56  RCL 00
  57  -
  58  RCL 01 
  59  STO IND Y
  60  XEQ "A*A1"
  61  XEQ "1/A"
  62  RCL M
  63  REGMOVE
  64  20
  65  XEQ "ST/A"
  66  21
  67  1/X
  68  ST- 01
  69  XEQ "A*A1"
  70  .1
  71  ST+ 01
  72  XEQ "A*A1"
  73  SIGN
  74  ST- 01 
  75  XEQ "A*A1"
  76  6
  77  XEQ "ST/A"
  78  RCL M
  79  REGMOVE
  80  SIGN
  81  RCL 00 
  82   E3
  83  /
  84  ISG X
  85  LASTX
  86  /
  87  +
  88  RCL 00
  89  ST+ X
  90  +
  91  STO N
  92  REGMOVE 
  93  XEQ "1/A"
  94  RCL 00
  95  ENTER^
  96  ST+ Y
  97  LBL 02
  98  RCL IND Y 
  99  RCL IND Y
100  -
101  2
102  /
103  STO IND Y 
104  RDN
105  DSE Y
106  DSE X
107  GTO 02
108  RCL N
109  REGSWAP
110  XEQ "LNA"
111  RCL 00
112  4
113  *
114  STO M
115  RCL 00
116  3
117  *
118  RCL 00
119  LBL 03
120  RCL IND M
121  RCL IND Z
122  +
123  ST+ IND Y
124  RDN
125  DSE M
126  DSE Y
127  DSE X
128  GTO 03
129  RCL 00
130   E3
131  /
132  ISG X
133  CLA
134  END

 
      ( 248 bytes / SIZE 4n+1 )
 
 

      STACK        INPUT      OUTPUT
           X             /         1.nnn

   where   1.nnn   is the control number of the result.

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

     n = 4  STO 00

      1  STO 01
      2  STO 02
      3  STO 03
      4  STO 04

      XEQ "PSIA"   >>>>     1.004                     ---Execution time = 31s---

            R01 =  1.686531557
            R02 =  0.548896352
            R03 =  0.823344528
            R04 =  1.097792703

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

3°)  Error Function
 

Formula:     erf a  = (2/pi1/2)  SUMn=0,1,2,.....  (-1)n a2n+1 / (n! (2n+1))
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "ERFA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R3n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:     "A*A1"                ( cf "Anions for the HP-41" )
                            "ST*A"    "DSA"   ( cf paragraph 0°) above )
 
 

 01  LBL "ERFA"
 02  RCL 00
 03  .1
 04  %
 05  STO Z
 06  ISG X
 07  +
 08   E3
 09  /
 10  1
 11  +
 12  STO M
 13  REGMOVE
 14  +
 15  REGMOVE
 16  XEQ "A*A1"
 17  0
 18  X<> M
 19  REGSWAP
 20  LBL 01
 21  SIGN
 22  RCL M
 23  +
 24  STO M
 25  ENTER^
 26  ST+ X
 27  STO Z
 28  DSE Z
 29  ISG X
 30  CLX
 31  *
 32  /
 33  CHS
 34  XEQ "ST*A"
 35  XEQ "A*A1"
 36  XEQ "DSA"
 37  X#0?
 38  GTO 01 
 39  RCL 00 
 40   E3
 41  /
 42  ISG X
 43  STO Y
 44   E3
 45  /
 46  1
 47  +
 48  RCL 00       
 49  ST+ X
 50  +
 51  REGMOVE
 52  2
 53  PI
 54  SQRT
 55  /
 56  XEQ "ST*A"
 57  X<> Z
 58  CLA
 59  END

 
       ( 115 bytes / SIZE 3n+1 )
 
 

      STACK        INPUT      OUTPUT
           X             /         1.nnn

   where   1.nnn   is the control number of the result.

Example:            a = 1 + 1.1 i + 1.2 j + 1.3 k

    4    STO 00         ( quaternion )

    1    STO 01
   1.1  STO 02
   1.2  STO 03
   1.3  STO 04        XEQ "ERFA"   >>>>    1.004                            ---Execution time = 38s---

           R01 =  -2.355991406
           R02 =  -3.358261142
           R03 =  -3.663557604
           R04 =  -3.968854073

Whence,   Erf ( 1 + 1.1 i + 1.2 j + 1.3 k ) =  -2.355991406 - 3.358261142 i - 3.663557604 j - 3.968854073 k
 

4°)  Exponential Integral, Sine & Cosine Integral
 

     a) Programs#1
 

-The following programs caculates
 

  "EIA"       Ei(a)  = C + Ln(a) + Sumn=1,2,.....  an/(a.a!)                where  C = 0.5772156649... = Euler's constant.

  "SIA"       Si(a)  = Summ=0,1,2,..... (-1)m a2m+1/((2m+1).(2m+1)!)           if  flag  F01  is clear
                 Shi(a) = Summ=0,1,2,.....  a2m+1/((2m+1).(2m+1)!)                   if  flag  F01  is set

  "CIA"      Ci(a)  = C + ln(a) + Summ=1,2,..... (-1)m a2m/(2m.(2m)!)        if  flag  F01  is clear
                 Chi(a)= C + ln(a) + Summ=1,2,.....  a2m/(2m.(2m)!)                 if  flag  F01  is set

  "ENA"    Em(a) = am-1 Gam(1-m) - [1/(1-m)] M( 1-m , 2-m ; -a )      where    M = Kummer's function

                               with   m # 1 , 2 , 3 , ................
 
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing these programs )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R3n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flag:  F01
Subroutines:  "A*A1"   "LNA"                   ( cf "Anions for the HP-41" )
                         "ST*A"   "ST/A"  "DSA"       ( cf paragraph 0 above )
                         "KUMA"                              ( cf paragraph 5 below )
                         "1/G+"                                  ( cf "Gamma Function for the HP-41" )
 
 

  01  LBL "EIA"
  02  XEQ 00
  03  CLX
  04  X<> M
  05  CLRGX
  06  SIGN
  07  STO 01
  08  LBL 01
  09  XEQ "A*A1"
  10  RCL M
  11  ENTER^
  12  ISG X
  13  CLX
  14  STO M
  15  X^2
  16  /
  17  X#0?
  18  XEQ "ST*A"
  19  XEQ "DSA"
  20  X#0?
  21  GTO 01
  22  GTO 04
  23  LBL "SIA"
  24  RCL 00
  25  .1
  26  %
  27  STO N
  28  ISG X
  29  +
  30   E3
  31  /
  32  1
  33  +
  34  REGMOVE
  35  STO M
  36  XEQ "A*A1"
  37  RCL M 
  38  REGSWAP
  39  RCL N
  40  +
  41  REGMOVE
  42  CLA
  43  LBL 02 
  44  RCL M
  45  ISG X
  46  CLX
  47  STO M 
  48  ST+ X
  49  ENTER^
  50  ENTER^
  51  DSE Z
  52  ISG X
  53  CLX
  54  X^2
  55  *
  56  /
  57  FC? 01
  58  CHS
  59  XEQ "ST*A"
  60  XEQ "A*A1"
  61  XEQ "DSA"
  62  X#0?
  63  GTO 02
  64  GTO 04
  65  LBL "CIA"
  66  XEQ 00
  67  RCL M
  68   E3
  69  /
  70  1
  71  +
  72  RCL 00
  73  +
  74  REGMOVE
  75  XEQ "A*A1"
  76  RCL N
  77  REGMOVE
  78  CLX
  79  X<> M
  80  CLRGX
  81  SIGN
  82  STO 01
  83  LBL 03
  84  RCL M
  85  ENTER^
  86  ST+ Y
  87  ISG X
  88  CLX
  89  STO M
  90  ST+ X
  91  X^2
  92  LASTX
  93  DSE X
  94  *
  95  /
  96  X=0?
  97  .25
  98  FC? 01
  99  CHS
100  XEQ "ST*A"
101  XEQ "A*A1"
102  XEQ "DSA"
103  X#0?
104  GTO 03
105  LBL 04
106  RCL 00
107   E3
108  /
109  ISG X
110  STO M
111   E3
112  /
113  1
114  +
115  RCL 00
116  ST+ X
117  +
118  REGMOVE
119  X<> M
120  CLA
121  GTO 06
122  LBL 00
123  RCL 00
124  .1
125  %
126  ISG X
127  STO M
128  +
129   E3
130  /
131  1
132  +
133  STO N
134  REGMOVE
135  XEQ "LNA"
136  FRC
137  RCL N
138  +
139  REGMOVE
140  FRC
141   E3
142  *
143  .5772156649
144  ST+ IND Y
145  RTN
146  LBL "ENA"
147  CHS
148  1
149  +
150  STO Y
151  1
152  ST+ Y
153  CHS
154  XEQ "ST*A"
155  RDN
156  XEQ "KUMA"
157  X<> Z
158  ENTER^
159  CHS
160  STO M
161  XEQ "ST/A"
162  RCL 00
163  .1
164  %
165  ISG X
166  +
167   E3
168  /
169  1
170  +
171  REGSWAP
172  SIGN
173  CHS
174  XEQ "ST*A"
175  RCL M
176  XEQ "A^X"
177  CLX
178  X<> M
179  CHS
180  XEQ "1/G+"
181  XEQ "ST/A"
182  RCL 00
183  ST+ X
184  RCL 00
185  LBL 05
186  RCL IND Y
187  ST+ IND Y
188  RDN
189  DSE Y
190  DSE X
191  GTO 05
192  LBL 06
193  RCL 00
194   E3
195  /
196  ISG X
197  END

 
    ( 406 bytes / SIZE 3n+1 )
 
 

      STACK        INPUT      OUTPUT
           X             /         1.nnn

   where   1.nnn   is the control number of the result.

-For  Em(a), place m in X-register before executing "ENA"

Examples:               a = 1 + 2 i + 3 j + 4 k                      4  STO 00    ( quaternions )
 

a) Exponential Integral Ei(a):

  1   STO 01
  2   STO 02
  3   STO 03
  4   STO 04      XEQ "EIA"   >>>>   1.004

                     R01 =  -0.381065061
                     R02 =   1.052055476
                     R03 =   1.578083216
                     R04 =   2.104110953

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

b) Sine Integral  Si(a):

  1   STO 01
  2   STO 02
  3   STO 03
  4   STO 04      CF 01    XEQ "SIA"   >>>>   1.004

                     R01 =  17.98954626
                     R02 =   7.075096290
                     R03 =  10.61264444
                     R04 =  14.15019258

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

c) Hyperbolic Sine Integral  Shi(a):
 

  1   STO 01
  2   STO 02
  3   STO 03
  4   STO 04      SF 01    XEQ "SIA"   >>>>   1.004

                     R01 =  -0.160779308
                     R02 =   0.522153864
                     R03 =   0.783230796
                     R04 =   1.044307726

  Whence:   Shi ( 1 + 2 i + 3 j + 4 k ) =  -0.160779308 + 0.522153864 i + 0.783230796 j + 1.044307726 k
 

d) Cosine Integral  Ci(a):
 

  1   STO 01
  2   STO 02
  3   STO 03
  4   STO 04      CF 01    XEQ "CIA"   >>>>   1.004

                     R01 =  19.04999179
                     R02 =  -6.098016771
                     R03 =  -9.147025156
                     R04 =  -12.19603355

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

e) Hyperbolic Cosine Integral  Chi(a):
 

  1   STO 01
  2   STO 02
  3   STO 03
  4   STO 04      SF 01    XEQ "CIA"   >>>>   1.004

                     R01 =  -0.220285752
                     R02 =   0.529901612
                     R03 =   0.794852420
                     R04 =   1.059803224

   So,   Chi ( 1 + 2 i + 3 j + 4 k ) =  -0.220285752 + 0.529901612 i + 0.794852420 j + 1.059803224 k
 

f) Exponential Integral Em(a):
 

-Compute:  EPI( 1 + i/2 + j/3 + k/4 )

      m = PI     a = 1 + i/2 + j/3 + k/4

 ( 4   STO 00 )

   1   STO 01
   2   1/X  STO 02
   3   1/X  STO 03
   4   1/X  STO 04

  PI   XEQ "ENA"   >>>>   1.004

            R01 =  0.066444329
            R02 = -0.059916132
            R03 = -0.039944088
            R04 = -0.029958065

   EPI( 1 + i/2 + j/3 + k/4 ) = 0.066444329 - 0.059916132 i - 0.039944088 j - 0.029958065 k
 

     b) Em(a)

-The program above does not work if m is a positive integer.
-"ENA" listed hereunder works in this case too.
 

Formulae:

     Em(a) = (-a)m-1 ( -Ln a - gamma + Sumk=1,...,m-1 1/k ) / (m-1)!  - Sumk#m-1 (-a)k / (k-m+1) / k!      where  gamma = Euler's Constant = 0.5772156649...

     and   E0(a) = (1/a).exp(-a)
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "ENA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R3n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:  "A*A1"   "LNA"  "E^A"  "A^X"                 ( cf "Anions for the HP-41" )
                         "ST*A"   "ST/A"  "DSA"                            ( cf paragraph 0 above )
                         "KUMA"                                                   ( cf paragraph 5 below )
                         "1/G+"                                                       ( cf "Gamma Function for the HP-41" )

-Line 43 may be replaced by  1.001  XEQ "HGFB"    ( cf "Anionic Special Functions(II)" )
-Line 79 is a three-byte  GTO 04
-Lines 124 to 145 may be replaced by the M-Code routine  DS*A
 
 

  01  LBL "ENA"
  02  X#0?
  03  GTO 00
  04  RCL 00
  05  .1
  06  %
  07  ISG X
  08  +
  09   E3
  10  /
  11  1
  12  +
  13  STO M
  14  REGMOVE
  15  SIGN
  16  CHS
  17  XEQ "ST*A"
  18  XEQ "E^A"
  19  CLX
  20  X<> M
  21  REGSWAP
  22  XEQ "1/A"
  23  XEQ "A*A1" 
  24  RTN
  25  LBL 00
  26  STO M
  27  X<0?
  28  GTO 04
  29  FRC
  30  X=0?
  31  GTO 00
  32  LBL 04
  33  RCL M
  34  CHS
  35  1
  36  +
  37  STO Y
  38  1
  39  ST+ Y
  40  CHS
  41  XEQ "ST*A"
  42  RDN
  43  XEQ "KUMA"
  44  X<> Z
  45  ENTER^
  46  CHS
  47  STO M
  48  XEQ "ST/A"
  49  RCL 00
  50  .1
  51  %
  52  ISG X
  53  +
  54   E3
  55  /
  56  1
  57  +
  58  REGSWAP
  59  SIGN
  60  CHS
  61  XEQ "ST*A"
  62  RCL M
  63  XEQ "A^X"
  64  CLX
  65  X<> M
  66  CHS
  67  XEQ "1/G+"
  68  XEQ "ST/A"
  69  RCL 00
  70  ENTER^
  71  ST+ Y
  72  LBL 01
  73  RCL IND Y
  74  ST+ IND Y
  75  RDN
  76  DSE Y
  77  DSE X
  78  GTO 01
  79  GTO 04
  80  LBL 00
  81  SIGN
  82  CHS
  83  XEQ "ST*A"
  84  RCL 00
  85  .1
  86  %
  87  ISG X
  88  STO Z
  89  +
  90   E3
  91  /
  92  1
  93  +
  94  REGMOVE 
  95  X<>Y
  96  CLRGX
  97  FRC
  98  +
  99  REGMOVE
100  SIGN
101  STO 01
102  STO N
103  RCL 00
104  ST+ X
105  +
106  RCL M
107  1
108  -
109  X#0?
110  1/X
111  STO IND Y
112  LBL 02
113  XEQ "A*A1" 
114  RCL M
115  RCL N
116  XEQ "ST/A"
117  -
118  1
119  ST+ N
120  X=Y?
121  GTO 02
122  -
123  1/X
124  STO O
125  3
126  RCL 00
127  *
128  STO P
129  LASTX
130  ENTER^
131  CLX
132  LBL 05
133  RCL IND Y
134  RCL O
135  *
136  RCL IND P
137  +
138  STO IND P
139  LASTX
140  -
141  ABS
142  +
143  DSE P
144  DSE Y
145  GTO 05
146  X#0?
147  GTO 02
148  RCL 00
149   E3
150  /
151  ISG X
152  LASTX
153  /
154  1
155  +
156  RCL 00
157  +
158  STO O
159  REGMOVE
160  SIGN
161  CHS
162  XEQ "ST*A" 
163  XEQ "LNA"
164  RCL M
165  1
166  -
167  STO M
168  SIGN
169  CLST
170  LBL 03
171  CLX
172  LASTX
173 X#0?
174   1/X
175  ST+ Y
176  DSE L
177  GTO 03
178  X<>Y
179  .5772156649
180  -
181  ST- 01
182  RCL M
183  FACT
184  CHS
185  XEQ "ST/A"
186  RCL O
187  REGSWAP
188  RCL M
189  XEQ "A^X"
190  XEQ "A*A1" 
191  XEQ "DSA"
192  RCL O
193  RCL 00
194  +
195  REGMOVE
196  LBL 04
197  RCL 00
198    E3
199  / 
200  ISG X
201  CLA
202  END

 
        ( 387 bytes / SIZE 3n+1 )
 
 

      STACK        INPUT      OUTPUT
           X            m         1.nnn

   where   1.nnn   is the control number of the result.

Example1:           m = PI   ,    a = 1 + i/2 + j/3 + k/4                          ( quaternion ->  4  STO 00 )

   1   STO 01
   2   1/X  STO 02
   3   1/X  STO 03
   4   1/X  STO 04

   PI   XEQ "ENA"  >>>>  1.004                                    ---Execution time = 24s---
 

            R01 =  0.066444329
            R02 = -0.059916132
            R03 = -0.039944088
            R04 = -0.029958065

   EPI( 1 + i/2 + j/3 + k/4 ) = 0.066444329 - 0.059916132 i - 0.039944088 j - 0.029958065 k

Example2:          m = 2   ,    a = 1 + i/2 + j/3 + k/4                          ( quaternion ->  4  STO 00 )

   1   STO 01
   2   1/X  STO 02
   3   1/X  STO 03
   4   1/X  STO 04

   2   XEQ "ENA"  >>>>  1.004                                    ---Execution time = 28s---
 

            R01 =  0.082262051
            R02 = -0.087391298
            R03 = -0.058260865
            R04 = -0.043695649

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

     c) An Asymptotic Expansion for Em(a)
 

-If a is very "large", "ENA" will produce poor accuracy or even meaningless results.
-"AENA" is preferable.
 

Formula:     Em(a) ~  (1/a) exp(-a)  2F0(1,m;;-1/a)
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "AENA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R4n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:  "A*A1"   "E^A"   "1/A"       ( cf "Anions for the HP-41" )
                         "ST*A"                               ( cf paragraph 0 above )
                         "2F0A "                              ( cf "Anionic Special Functions(II)" - paragraph 6°)d) )

-Line 22 may be replaced by  2   XEQ "HGFB"  cf "Anionic Special Functions(II)" - paragraph 6°)b)
 
 

 01  LBL "AENA"
 02  STO M
 03  1
 04  CHS
 05  XEQ "ST*A"
 06  RCL 00
 07   E3
 08  /
 09  ISG X
 10  RCL 00
 11  3
 12  *
 13  +
 14   E3
 15  /
 16  1
 17  +
 18  REGMOVE 
 19  XEQ "1/A"
 20  SIGN
 21  RCL M
 22  XEQ "2F0A"
 23  XEQ "A*A1"
 24  SIGN
 25  CHS
 26  XEQ "ST*A"
 27  RCL 00
 28  .1
 29  %
 30  ISG X
 31  STO Z
 32  +
 33   E3
 34  /
 35  1
 36  +
 37  REGMOVE 
 38  CLX
 39   E3
 40  /
 41  1
 42  +
 43  RCL 00
 44  3
 45  *
 46  +
 47  REGMOVE
 48  XEQ "E^A"
 49  XEQ "A*A1"
 50  END

 
        ( 105 bytes / SIZE 4n+1 )
 
 

      STACK        INPUT      OUTPUT
           X             m         1.nnn

   where   1.nnn   is the control number of the result.

Example:        m = PI   ,    a = 41 + 42 i + 43 j + 44 k                          ( quaternion ->  4  STO 00 )

   41  STO 01
   42  STO 02
   43  STO 03
   44  STO 04

   PI   XEQ "AENA"  >>>>  1.004                                    ---Execution time = 20s---

   R01 =  1.789481142 E-20
   R02 = -1.315072310 E-21
   R03 = -1.346383554 E-21
   R04 = -1.377694801 E-21

-So,   EPI( 41 + 42 i + 43 j + 44 k ) = 1.789481142 10 -20 - 1.315072310 10 -21 i - 1.346383554 10 -21 j - 1.377694801 10 -21 k

Notes:

-If the argument a is too "small" , the series will diverge too soon.
-However, it may also work with small arguments if  m is a negative integer.
-For instance, it returns correctly:

         E -2( 1 + 2 i + 3 j + 4 k ) = 0.047465815 - 0.020205534 i - 0.030308300 j - 0.040411067 k
 

5°)  Kummer's Functions
 
 

Formula:        M(p;q;a) = 1 +  (p)1/(q)1. a1/1! + ............. +  (p)k/(q)k . ak/k! + ..........  where (p)k = p.(p+1)(p+2) ...... (p+k-1)

    and   p & q  are real numbers
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "KUMA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R3n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:     "A*A1"  ( cf "Anions for the HP-41" )   "ST*A"    "DSA"   ( cf paragraph 0°) above )
 
 
 

 01  LBL "KUMA"
 02  CLA
 03  STO O
 04  X<>Y
 05  STO N
 06  RCL 00
 07  .1
 08  %
 09  STO Z
 10  ISG X
 11  STO T
 12  +
 13   E3
 14  /
 15  1
 16  +
 17  REGMOVE 
 18  +
 19  X<>Y
 20  CLRGX
 21  SIGN
 22  STO 01
 23  X<>Y
 24  REGMOVE
 25  LBL 01
 26  XEQ "A*A1" 
 27  RCL N
 28  RCL O
 29  RCL M
 30  ST+ Y
 31  ST+ Z
 32  ISG X
 33  CLX
 34  STO M
 35  *
 36  /
 37  XEQ "ST*A" 
 38  XEQ "DSA"
 39  X#0?
 40  GTO 01
 41  RCL N
 42  RCL O
 43  RCL 00
 44   E3
 45  /
 46  ISG X
 47  STO M
 48   E3
 49  /
 50  1
 51  +
 52  RCL 00
 53  ST+ X
 54  +
 55  REGMOVE 
 56  X<> M
 57  CLA
 58  END

 
    ( 108 bytes / SIZE 3n+1 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             /            p
           Y             p            q
           X             q         1.nnn

   where   1.nnn   is the control number of the result.

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

   4   STO 00        ( Quaternion )

   1   STO 01
   2   STO 02
   3   STO 03
   4   STO 04

   2   SQRT
  PI   XEQ "KUMA"   >>>>   1.004                     ---Execution time =36s---

      R01 =  -0.533892202
      R02 =    0.068805302
      R03 =    0.103207953
      R04 =    0.137610603

Whence,   M ( sqrt(2) , PI ; 1 + 2 i + 3 j + 4 k ) =  -0.533892203 + 0.068805302 i + 0.103207953 j + 0.137610603 k
 

6°)  Parabolic Cylinder Functions
 

     a)  Via Kummer's Functions
 

Formula:         Dm(a) = 2m/2 Pi1/2 exp(-a2/4) [ 1/Gam((1-m)/2) M( -m/2 , 1/2 , a2/2 ) - 21/2 ( a / Gam(-m/2) ) M [ (1-m)/2 , 3/2 , a2/2 ]
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "DNA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R4n:  temp    R4n+1 = m

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:  "A*A1"   "E^A"          ( cf "Anions for the HP-41" )
                         "ST*A"   "ST/A"         ( cf paragraph 0 above )
                         "KUMA"                    ( cf paragraph 5 above )
                         "1/G+"                        ( cf "Gamma Function for the HP-41" )
 
 
 

  01  LBL "DNA"
  02  CHS
  03  2
  04  /
  05  STO M
  06  RCL 00
  07  4
  08  *
  09  1
  10  +
  11  X<>Y
  12  STO IND Y 
  13  RCL 00
  14  .1
  15  %
  16  STO Z
  17  ST+ Z
  18  ISG X
  19  +
  20   E3
  21  /
  22  1
  23  +
  24  REGMOVE
  25  +
  26  REGMOVE
  27  XEQ "A*A1"
  28  .5
  29  XEQ "ST*A"
  30  RCL M
  31  +
  32  1.5
  33  XEQ "KUMA"
  34  RCL 00
  35  4
  36  *
  37  1
  38  +
  39  RCL IND X
  40  STO N
  41  XEQ "1/G+"
  42  2
  43  SQRT
  44  *
  45  XEQ "ST*A"
  46  RCL 00
  47   E3
  48  /
  49  ISG X
  50  RCL 00
  51  3
  52  *
  53  +
  54   E3
  55  /
  56  RCL 00
  57  +
  58  1
  59  +
  60  REGSWAP
  61  XEQ "A*A1" 
  62   E3
  63  /
  64  RCL 00
  65  3
  66  *
  67  +
  68  1
  69  +
  70  REGSWAP
  71  RCL N
  72  .5
  73  XEQ "KUMA"
  74  X<> Z
  75  STO N
  76  +
  77  XEQ "1/G+"
  78  XEQ "ST*A"
  79  RCL 00
  80  4
  81  *
  82  RCL 00
  83  LBL 01
  84  RCL IND Y
  85  ST- IND Y
  86  RDN
  87  DSE Y
  88  DSE X
  89  GTO 01
  90  RCL 00
  91   E3
  92  /
  93  ISG X
  94   E3
  95  /
  96  RCL 00
  97  +
  98  1
  99  +
100  REGSWAP 
101  2
102  CHS
103  XEQ "ST/A"
104  XEQ "E^A"
105  XEQ "A*A1" 
106  PI
107  SQRT
108  2
109  RCL N
110  Y^X
111  /
112  XEQ "ST*A" 
113  X<>Y
114  CLA
115  END

 
    ( 222 bytes / SIZE 4n+2 )
 
 

      STACK        INPUT      OUTPUT
           X             m         1.nnn

   where   1.nnn   is the control number of the result.

Example:       m = 0.4      a = 1 + 1.1 i + 1.2 j + 1.3 k

    4    STO 00

    1    STO 01
   1.1  STO 02
   1.2  STO 03
   1.3  STO 04

   0.4  XEQ "DNA"   >>>>    1.004                           ---Execution time = 68s---

             R01 =  2.606105102
             R02 = -0.969116886
             R03 = -1.057218421
             R04 = -1.145319956

   D0.4( 1 + 1.1 i + 1.2 j + 1.3 k ) = 2.606105102 - 0.969116886 i - 1.057218421 j - 1.145319956 k
 

Note:

-For large arguments, an asymptotic expansion is preferable:
 

     b)  Asymptotic Expansion
 

Formula:    Dm(a)  ~  am exp(-a2/4) [ 1 - m(m-1) / ( 2 a2 ) + m(m-1)(m-2)(m-3) / ( 2 ( 4 a4 ) ) - ....... ]
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "AEDNA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R3n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:  "A*A1"   "E^A"   "1/A"   "A^X"       ( cf "Anions for the HP-41" )     ->  For "A*A1" , use a version that does not alter synthetic register P
                         "ST*A"   "ST/A"                                 ( cf paragraph 0 above )

-Lines 59 to 79 may be replaced by 2 lines:  RCL M   DS*A  where  DS*A  is an M-Code routine listed in paragraph 0°)b)
 
 

 01  LBL "AEDNA"
 02  STO N
 03  RCL 00
 04  .1
 05  %
 06  STO Z
 07  ISG X
 08  +
 09   E3
 10  /
 11  1
 12  +
 13  STO M
 14  REGMOVE
 15  +
 16  STO O
 17  XEQ "A*A1"
 18  RCL M
 19  REGSWAP
 20  RCL N
 21  XEQ "A^X"
 22  RCL M
 23  REGSWAP
 24  RCL O
 25  REGMOVE 
 26  4
 27  CHS
 28  XEQ "ST/A"
 29  XEQ "E^A"
 30  XEQ "A*A1" 
 31  RCL O
 32  REGSWAP
 33  XEQ "1/A"
 34  RCL M
 35  REGMOVE
 36  RCL O
 37  REGSWAP
 38  REGMOVE
 39  SIGN
 40  STO M
 41  CLX
 42  STO P
 43  LBL 01
 44  XEQ "A*A1" 
 45  RCL N
 46  STO Y
 47  RCL P
 48  ST- Z
 49  ISG X
 50  CLX
 51  ST- Y
 52  ISG X
 53  CLX
 54  STO P
 55  /
 56  *
 57  CHS
 58  ST* M
 59  RCL 00
 60  PI
 61  INT
 62  *
 63  STO Q
 64  RCL 00
 65  ENTER^
 66  CLX
 67  LBL 02
 68  RCL IND Y 
 69  RCL M
 70  *
 71  RCL IND Q 
 72  +
 73  STO IND Q 
 74  LASTX
 75  -
 76  ABS
 77  +
 78  DSE Q
 79  DSE Y
 80  GTO 02
 81  X#0?
 82  GTO 01
 83  RCL O
 84  REGSWAP 
 85  RCL 00
 86  E3
 87  /
 88  ISG X
 89  CLA
 90  END

 
         ( 178 bytes / SIZE 3n+1 )
 
 

      STACK        INPUT      OUTPUT
           X             m         1.nnn

   where   1.nnn   is the control number of the result.

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

  4   STO 00     ( quaternion )

  1  STO 01
  2  STO 02
  3  STO 03
  4  STO 04

 PI    XEQ "AEDNA"   >>>> 1.004                                 ---Execution time =21s---        ( with  DS*A )

 R01 = -33138.87292
 R02 =   93318.15134
 R03 =  139977.2269
 R04 =  186636.3028

  DPI( 1 + 2 i + 3 j + 4 k ) =  -33138.87292 + 93318.15134 i + 139977.2269 j + 186636.3028 k

Notes:

-As usual, infinte loops will occur if a is too small.
-Add  RND  after line 80 to get a less accurate result in this case:
-The number of exact digits will be controlled by the display settings.

-This program may also be used if a is relatively small when m is a positive integer.
 

7°)  Laguerre's Functions
 

Formula:      Lm(b)(q) = [ Gam(m+b+1) / Gam(m+1) / Gam(b+1) ]  M ( -n , a+1 , q )           m , b # -1 , -2 , -3 , ....

   where  Gam = Gamma function
     and      M   = Kummer's function
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "LANA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R4n:  temp    R4n+1 = m

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:  "ST*A"          ( cf paragraph 0 above )
                         "KUMA"       ( cf paragraph 4 above )
                         "1/G+"           ( cf "Gamma Function for the HP-41" )
 
 
 

 01  LBL "LANA"
 02  CHS
 03  X<>Y
 04  1
 05  +
 06  XEQ "KUMA"
 07  RDN
 08  STO M           
 09  X<>Y
 10  CHS
 11  STO N 
 12  1
 13  +
 14  XEQ "1/G+" 
 15  X<> M
 16  ST+ N
 17  XEQ "1/G+"   
 18  ST* M
 19  RCL N
 20  XEQ "1/G+" 
 21  RCL M
 22  X<>Y
 23  /
 24  XEQ "ST*A"  
 25  RCL 00          
 26   E3
 27  /
 28  ISG X
 29  CLA
 30  END

 
    ( 86 bytes / SIZE 4n+1 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             b            /
           X             m         1.nnn

   where   1.nnn   is the control number of the result.

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

   4  STO 00

   1   STO 01
   2   STO 02
   3   STO 03
   4   STO 04

   3   SQRT
  PI  XEQ "LANA"   >>>>  1.004                          ---Execution time = 28s---

           R01 = -59.37663947
           R02 =   3.135355698
           R03 =   4.703033534
           R04 =   6.270711379

        Lpisqrt(3) ( 1 + 2 i + 3 j + 4 k ) =  -59.37663947 + 3.135355698 i + 4.703033534 j + 6.270711379 k
 

8°)  Hypergeometric Functions
 

-"HGFA"  computes  1F2(p;q,s;a) = 1 +  (p)1/(q)1/(s)1. a1/1! + ............. +  (p)n/(q)n/(s)n . ak/k! + ..........               if flag F04 is set

                            or  2F1(p,q;s;a) = 1 +  (p)1(q)1/(s)1. a1/1! + ............. +  (p)n(q)n/(s)n . ak/k! + ..........                 if flag F04 is cleared
 

-However, if CF 04 & if  2F1  does not exist, the regularized function  2F1 tilde is returned.  -> Flag F00 is set in that case.

-The parameters must be real numbers.
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "HGFA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R3n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  F00 - F04
Subroutines:  "ST*A"   "DSA"        ( cf paragraph 0 above )
                         "A*A1"  "A^X"     ( cf "Anions for the HP-41" )             ->  For "A*A1" , use a version that does not alter synthetic register P
 
 

  01  LBL "HGFA"
  02  STO M
  03  RDN
  04  STO N
  05  X<>Y
  06  STO O
  07  RCL 00
  08   E3
  09  /
  10  ISG X
  11  RCL 00
  12  +
  13   E3
  14  /
  15  1
  16  +
  17  REGMOVE
  18  CF 00
  19  FS? 04
  20  GTO 02
  21  RCL M
  22  ENTER^
  23  FRC
  24  X=0?
  25  X<Y?
  26  GTO 02
  27  RCL O
  28  ENTER^
  29  FRC
  30  X=0?
  31  X<Y?
  32  GTO 00
  33  RCL M
  34  LASTX
  35  X>Y?
  36  GTO 02
  37  LBL 00
  38  RCL N
  39  ENTER^
  40  FRC
  41  X=0?
  42  X<Y?
  43  GTO 00
  44  RCL M
  45  LASTX
  46  X>Y?
  47  GTO 02 
  48  LBL 00
  49  SF 00
  50  1
  51  RCL M
  52  -
  53  STO P
  54  XEQ "A^X" 
  55  RCL P
  56  X<>Y
  57  SIGN
  58  LBL 01
  59  RCL O
  60  ENTER^
  61  CLX
  62  SIGN
  63  -
  64  R^
  65  +
  66  *
  67  RCL N
  68  R^
  69  +
  70  ENTER^
  71  CLX
  72  SIGN
  73  -
  74  *
  75  R^
  76  /
  77  DSE Y
  78  GTO 01
  79  XEQ "ST*A"
  80  LBL 02
  81  RCL 00
  82  PI
  83  INT
  84  10^X
  85  /
  86  ISG X
  87  FS? 00
  88  GTO 00
  89  CLRGX
  90  SIGN
  91  STO 01
  92  LASTX
  93  LBL 00
  94  RCL 00
  95  ST+ X
  96  +
  97  PI
  98  INT
  99  10^X
100  /
101  ENTER^
102  SIGN
103  +
104  REGMOVE
105  RCL P
106  FC? 00
107  CLX
108  STO P
109  LBL 03
110  XEQ "A*A1"
111  RCL O
112  RCL N
113  RCL M
114  RCL P
115  ST+ Y
116  ST+ Z
117  ST+ T
118  ISG X
119  CLX
120  STO P
121  *
122  RDN
123  FS? 04
124  1/X
125  *
126  R^
127  /
128  XEQ "ST*A"
129  XEQ "DSA"
130  X#0?
131  GTO 03
132  RCL 00
133   E3
134  /
135  ISG X
136  STO Y
137  LASTX
138  /
139  1
140  +
141  RCL 00
142  ST+ X
143  +
144  REGMOVE
145  X<> O
146  RCL N
147  RCL M
148  R^
149  CLA
150  END

 
   ( 238 bytes / SIZE 3n+1 )
 
 

      STACK        INPUTS      OUTPUTS
           T             /            p
           Z             p            q
           Y             q            s
           X             s         1.nnn

   where   1.nnn   is the control number of the result.

Example1:       2F1       a = 0.1 + 0.2 i + 0.3 j + 0.4 k           p = 1.1   q = 1.2   s = 1.3

   ( 4  STO 00 )   CF 04

     0.1  STO 01
     0.2  STO 02
     0.3  STO 03
     0.4  STO 04

     1.1  ENTER^
     1.2  ENTER^
     1.3  XEQ "HGFA"   >>>>   1.004                           ---Execution time = 56s---

            R01 =  0.814236592
            R02 =  0.184421233
            R03 =  0.276631850
            R04 =  0.368842467

-Flag F00 is clear, so:

    2F1 ( 1.1 , 1.2 ; 1.3 ;  0.1 + 0.2 i + 0.3 j + 0.4 k ) = 0.814236592 + 0.184421233 i + 0.276631850 j + 0.368842467 k
 

Example2:      2F1        a = 0.1 + 0.2 i + 0.3 j + 0.4 k           p = 1   q = 2   s = -4

   ( 4  STO 00 )    CF 04

     0.1  STO 01
     0.2  STO 02
     0.3  STO 03
     0.4  STO 04

     1     ENTER^
     2     ENTER^
    -4    XEQ "HGFA"   >>>>   1.004                           ---Execution time = 108s---

            R01 =  -7.152496122
            R02 =  -9.061272199
            R03 =  -13.59190768
            R04 =  -18.12254405

-Flag F00 is set, so F tilde has been calculated:

    2F~1 ( 1 , 2 ; -4 ;  0.1 + 0.2 i + 0.3 j + 0.4 k ) = -7.152496122 - 9.061272199 i - 13.59190768 j - 18.12254405 k
 

Example3:     1F2         a = 0.1 + 0.2 i + 0.3 j + 0.4 k           p = 1.1   q = 1.2   s = 1.3

   ( 4  STO 00 )   SF 04

     0.1  STO 01
     0.2  STO 02
     0.3  STO 03
     0.4  STO 04

     1.1  ENTER^
     1.2  ENTER^
     1.3  XEQ "HGFA"   >>>>   1.004                           ---Execution time = 86s---

            R01 =  1.028367021
            R02 =  0.146116085
            R03 =  0.219174127
            R04 =  0.292232170

-Whence,       1F2 ( 1.1 ; 1.2 , 1.3 ;  0.1 + 0.2 i + 0.3 j + 0.4 k ) = 1.028367021 + 0.146116085 i + 0.219174127 j + 0.292232170 k
 

9°)  Associated Legendre Functions
 

    a)  Via Hypergeometric Functions
 

 "ALFA" computes 4 kinds of associated Legendre functions:

-If  CF 02  &  CF 03   we get the associated Legendre function of the 1st kind, type 2
-If  CF 02  &  SF 03    -------------------------------------------------------, type 3

-If  SF 02   &  CF 03   we get the associated Legendre function of the 2nd kind, type 2
-If  SF 02   &  SF 03    --------------------------------------------------------, type 3
 

Formulae:    1st Kind
 

  Type 2      Prm(a) = [ (a+1)/(1-a) ]m/2  2F1(-r , r+1 ; 1-m ; (1-a)/2 ) / Gamma(1-m)        (  a # 1 )                | a - 1 | < 2

  Type 3      Prm(a) = [ (a+1)/(a-1) ]m/2  2F1(-r , r+1 ; 1-m ; (1-a)/2 ) / Gamma(1-m)        (  a # 1 )                | a - 1 | < 2
 

Formulae:    2nd Kind
 

  Type 2      Qrm(a) = 2m pi1/2  (1-a2)-m/2 [ -Gam((1+m+r)/2)/(2.Gam((2-m+r)/2)) . sin pi(m+r)/2 .  2F1(-r/2-m/2 ; 1/2+r/2-m/2 ; 1/2 ; a2 )          | a | < 1
                                                                + a  Gam((2+r+m)/2) / Gam((1+r-m)/2) . cos pi(m+r)/2 . 2F1((1-m-r)/2 ; (2+r-m)/2 ; 3/2 ; a2 )  ]

  Type 3        Qrm(a) =  exp( I (m.PI) ) 2 -r-1 sqrt(PI) Gam(m+r+1) a -m-r-1 (a2-1)m/2  2F~1( (2+m+r)/2 , (1+m+r)/2 ; r+3/2 ; 1/a2 )                      | a | > 1
 

                     where  I = Im(a) / | Im(a) |  multiplied by the sign of the 1st imaginary component of the anion  (  or  I = i  if  Im(a) = 0 )
 
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "ALFA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R4n+1:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  F00-F02-F03-F04

Subroutines:  "ST*A"                                       ( cf paragraph 0 above )
                         "HGFA"                                      ( cf paragraph 8 above )
                         "A*A1"  "A^X"  "1/A"  "E^A"    ( cf "Anions for the HP-41" )
                         "1/G+"                                        ( cf "Gamma Function for the HP-41" )
 

-Lines 04 and 73 are three-byte  GTO 02   ( not really important... )
 
 

  01  LBL "ALFA"
  02  CF 04
  03  FS? 02
  04  GTO 02
  05  .5
  06  CHS
  07  XEQ "ST*A"
  08  ST- 01
  09  RDN
  10  CHS
  11  1
  12  ST- Z
  13  RCL Y
  14  -
  15  R^
  16  CHS
  17  XEQ "HGFA"
  18  RDN
  19  STO M
  20  FS?C 00
  21  GTO 01
  22  XEQ "1/G+"
  23  XEQ "ST*A"
  24  LBL 01
  25  RCL 00
  26   E3
  27  /
  28  ISG X
  29  STO Y
  30  RCL 00
  31  3
  32  *
  33  +
  34   E3
  35  /
  36  1
  37  +
  38  STO O
  39  REGMOVE
  40  X<>Y
  41   E3
  42  /
  43  RCL 00
  44  +
  45  1
  46  ST- M
  47  +
  48  STO N
  49  REGMOVE
  50  XEQ "1/A"
  51  RCL N
  52  REGSWAP
  53  SIGN
  54  ST- 01
  55  CHS
  56  FC? 03
  57  XEQ "ST*A"
  58  XEQ "A*A1"
  59  RCL M
  60  CHS
  61  2
  62  /
  63  XEQ "A^X"
  64  RCL O
  65  RCL 00
  66  +
  67  REGSWAP
  68  XEQ "A*A1"
  69  CLA
  70  RTN
  71  LBL 02
  72  FC? 03
  73  GTO 02
  74  RCL 00
  75  .1
  76  %
  77  STO M
  78  ISG X
  79  +
  80   E3
  81  /
  82  1
  83  +
  84  REGMOVE 
  85  RCL M
  86  ST+ X
  87  +
  88  REGMOVE 
  89  CLX
  90  RCL 00
  91  4
  92  *
  93  1
  94  +
  95  X<> T
  96  STO IND T
  97  RCL Y
  98  +
  99  1
100  X<>Y
101  ST+ Y
102  2
103  ST/ Z
104  ST+ Y
105  /
106  X<> Z
107  1.5
108  +
109  STO M
110  RDN
111  STO N
112  X<>Y
113  STO O
114  XEQ "A*A1"
115  XEQ "1/A"
116  RCL O
117  RCL N
118  RCL M
119  XEQ "HGFA"
120  RDN
121  STO M
122  RDN
123  STO N
124  X<>Y
125  STO O
126  FS?C 00
127  GTO 01
128  R^
129  XEQ "1/G+"
130  XEQ "ST*A"
131  LBL 01
132  RCL N
133  ST+ X
134  STO N
135  XEQ "1/G+"
136  2
137  RCL M
138  .5
139  -
140  Y^X
141  *
142  PI
143  SQRT
144  X<>Y
145  /
146  XEQ "ST*A"
147  RCL 00
148  .1
149  %
150  ISG X
151  STO Z
152  +
153   E3
154  /
155  1
156  +
157  STO M
158  REGMOVE 
159  X<>Y
160   E3
161  /
162  1
163  +
164  RCL 00
165  3
166  *
167  +
168  STO O
169  REGMOVE
170  RCL N
171  CHS
172  XEQ "A^X"
173  XEQ "A*A1"
174  RCL 00
175  ST+ X
176  +
177   E3
178  /
179  1
180  +
181  REGMOVE
182  RCL O
183  REGMOVE
184  RCL M
185  REGMOVE
186  XEQ "A*A1"
187  SIGN
188  ST- 01
189  RCL 00
190  4
191  *
192  1
193  +
194  RCL IND X
195  STO N
196  2
197  /
198  XEQ "A^X"
199  RCL 00
200  ST+ X
201  RCL M
202  +
203  REGMOVE
204  XEQ "A*A1"
205  RCL M
206  REGMOVE
207  RCL O
208  REGMOVE 
209  CLX
210  STO 01
211  RCL N
212  PI
213  *
214  RCL 00
215  X<>Y
216  0
217  LBL 03
218  RCL IND Z 
219  X^2
220  +
221  DSE Z
222  GTO 03
223  SQRT
224  X#0?
225  GTO 01
226  SIGN
227  STO 02
228  LBL 01
229  RCL 02
230  SIGN
231  *
232  /
233  XEQ "ST*A"
234  XEQ "E^A"
235  XEQ "A*A1"
236  CLA
237  RTN
238  LBL 02
239  STO N
240  X<>Y
241  STO M
242  RCL 00
243  4
244  *
245  1
246  +
247  X<>Y
248  STO IND Y
249  RCL 00
250  .1
251  %
252  STO Z
253  ST+ Z
254  ISG X
255  +
256   E3
257  /
258  1
259  +
260  REGMOVE
261  +
262  REGMOVE
263  XEQ "A*A1"
264  SIGN
265  RCL M
266  -
267  RCL N
268  -
269  2
270  LASTX
271  +
272  RCL M
273  -
274  2
275  ST/ Z
276  /
277  1.5
278  XEQ "HGFA"
279  RDN
280  STO M
281  RDN
282  STO N
283  X<>Y
284  STO O
285  1
286  CHS
287  ACOS
288  *
289  SIN
290  XEQ "ST*A" 
291  RCL N
292  .5
293  -
294  XEQ "1/G+"
295  XEQ "ST*A"
296  RCL M
297  RCL O
298  -
299  XEQ "1/G+"
300  XEQ "ST/A"
301  RCL 00
302  .1
303  %
304  STO Z
305  ST+ Z
306  ISG X
307  +
308   E3
309  /
310  1
311  +
312  +
313  STO P
314  RCL 00
315  +
316  REGSWAP
317  XEQ "A*A1"
318  RCL P
319  REGSWAP
320  .5
321  ST- O
322  ST- N
323  RCL O
324  RCL N
325  RCL M
326  1
327  -
328  XEQ "HGFA"
329  RDN
330  STO M
331  RDN
332  STO N
333  X<>Y
334  STO O
335  CLX
336  .5
337  +
338  XEQ "1/G+"
339  RCL O
340  1
341  CHS
342  ACOS
343  *
344  SIN
345  *
346  XEQ "ST*A"
347  .5
348  RCL O
349  -
350  XEQ "1/G+"
351  ST+ X
352  XEQ "ST/A"
353  4
354  RCL 00
355  ST* Y
356  LBL 04
357  RCL IND X 
358  ST+ IND Z
359  RDN
360  DSE Y
361  DSE X
362  GTO 04
363  RCL 00
364   E3
365  /
366  ISG X
367   E3
368  /
369  1
370  +
371  RCL 00
372  +
373  REGMOVE 
374  SIGN
375  CHS
376  XEQ "ST*A"
377  ST- 01
378  RCL 00
379  4
380  *
381  1
382  +
383  RCL IND X
384  STO M
385  2
386  /
387  CHS
388  XEQ "A^X"
389  PI
390  SQRT
391  2
392  RCL M
393  Y^X
394  *
395  XEQ "ST*A"
396  RCL 00
397  .1
398  %
399  ISG X
400  +
401   E3
402  /
403  1
404  +
405  RCL 00
406  3
407  *
408  +
409  REGMOVE
410  XEQ "A*A1"
411  CLA
412  END

 
      ( 779 bytes / SIZE 4n+2 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             m            /
           X             r         1.nnn

   where   1.nnn   is the control number of the result.

Example1:     Associated Legendre functions of the 1st kind     CF 02
 

  •  CF 02   CF 03   Legendre Function 1st kind - type 2

      m = sqrt(2)    r = sqrt(3)       a = 0.4 + 0.5 i + 0.6 j + 0.7 k

      0.4   STO 01
      0.5   STO 02
      0.6   STO 03
      0.7   STO 04

       2   SQRT
       3   SQRT   XEQ "ALFA"   >>>>   1.004                               ---Execution time = 80s---

       R01 =  -0.866899245
       R02 =  -1.808960474
       R03 =  -2.170752568
       R04 =  -2.532544662

>>>     Psqrt(2)sqrt(3) ( 0.4 + 0.5 i + 0.6 j + 0.7 k ) = -0.866899245 - 1.808960474 i - 2.170752568 j - 2.532544662 k
 

  •  CF 02  SF 03   Legendre Function 1st kind - type 3

      m = sqrt(2)    r = sqrt(3)       a = 0.4 + 0.5 i + 0.6 j + 0.7 k

      0.4   STO 01
      0.5   STO 02
      0.6   STO 03
      0.7   STO 04

       2   SQRT
       3   SQRT   XEQ "ALFA"   >>>>   1.004                              ---Execution time = 80s---

       R01 =  -2.494183067
       R02 =   1.424529612
       R03 =   1.709435534
       R04 =   1.994341458

>>>     Psqrt(2)sqrt(3) ( 0.4 + 0.5 i + 0.6 j + 0.7 k ) = -2.494183067 + 1.424529612 i + 1.709435534 j + 1.994341458 k
 

Example2:        Associated Legendre functions of the 2nd kind    SF 02

  •  SF 02   CF 03   Legendre Function 2nd kind - type 2

    m = 0.4    r = 1.3       a = 0.1 + 0.2 i + 0.3 j + 0.4 k

      0.1   STO 01
      0.2   STO 02
      0.3   STO 03
      0.4   STO 04

      0.4   ENTER^
      1.3   XEQ "ALFA"   >>>>   1.004                               ---Execution time = 71s---

      R01 =  -0.944581516
      R02 =  -0.368022565
      R03 =  -0.552033847
      R04 =  -0.736045129

    Q1.30.4( 0.1 + 0.2 i + 0.3 j + 0.4 k ) = -0.944581516 - 0.368022565 i - 0.552033847 j - 0.736045129 k

  •  SF 02   SF 03   Legendre Function 2nd kind - type 3

   m = sqrt(2)    r = - sqrt(3)       a = 1 + 2 i + 3 j + 4 k

    1   STO 01
    2   STO 02
    3   STO 03
    4   STO 04

    2   SQRT
    3   SQRT  CHS    XEQ "ALFA"  >>>>     1.004                               ---Execution time = 37s---

    R01 = -1.926294207
    R02 =  0.742062001
    R03 =  1.113093001
    R04 =  1.484124002

     Qsqrt(2)-sqrt(3) ( 1 + 2 i + 3 j + 4 k ) = -1.926294207 + 0.742062001 i + 1.113093001 j + 1.484124002 k

Note:

-In the last example, the results are different from those given in "Quaternionic Special Functions"
  because here, I've multiplied by  exp( I (m.PI) ) insead of  exp( i (m.PI) )

                     where  I = Im(a) / | Im(a) |  multiplied by the sign of the 1st imaginary component of the anion  (  or  I = i  if  Im(a) = 0 )

-I think it's a better generalization of the formula for the associated Legendre function of the 2nd kind - type 3.
-It would be probably even better to multiply Im(a) / | Im(a) | by the sign of the 1st imaginary component # 0
 

    b)  Functions of the 1st kind - Integer Order & Degree
 

 "PMNA" compute the functions of type 2 - if CF 03 - or the functions of type 3 - if  SF 03.
 

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

     Type 2                Pmm(a) = (-1)m (2m-1)!!  ( 1-a )m/2 ( 1+a )m/2                    where (2m-1)!! = (2m-1)(2m-3)(2m-5).......5.3.1
     Type 3                Pmm(a) =   (2m-1)!!  ( a-1 )m/2 ( a+1 )m/2
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "PMNA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R4n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flag:  F03
Subroutines:   ST*A   ST/A    AMOVE   ASWAP     ( cf paragraph 0 above )    X/E3 =  E3 /       X/2 =  2  /      X-1 =  1  -
                          A*A1   "A^X"      ( cf "Anions for the HP-41" )
 
 
 

  01  LBL "PMNA"
  02  STO N
  03  X<>Y
  04  STO M
  05  12
  06  AMOVE
  07  13
  08  AMOVE
  09  SIGN
  10  ST+ 01
  11  RCL M
  12  X/2
  13  XEQ "A^X"
  14  12
  15  ASWAP
  16  SIGN
  17  ST- 01
  18  CHS
  19  FC? 03
  20  ST*A
  21  RCL M
  22  X/2
  23  XEQ "A^X"
  24  A*A1
  25  32
  26  AMOVE       
  27  RCL M
  28  ST+ X
  29  1
  30  ST+ Y
  31  X<>Y
  32  LBL 01
  33  2
  34  -
  35  X>0?
  36  ST* Y
  37  X>0?
  38  GTO 01
  39  CLX
  40  SIGN
  41  FC? 03
  42  CHS
  43  RCL M
  44  Y^X
  45  *
  46  ST*A
  47  RCL 00
  48  4
  49  *
  50  .1
  51  %
  52  ISG X
  53  +
  54  RCL 00
  55  -
  56  CLRGX
  57  13
  58  AMOVE       
  59  RCL N
  60  X/E3
  61  RCL M
  62  +
  63  STO O
  64  LBL 02
  65  31
  66  AMOVE
  67  ISG O
  68  X=0?
  69  GTO 00
  70  RCL O
  71  INT
  72  ST+ X
  73  X-1
  74  ST*A
  75  A*A1
  76  RCL 00
  77  ST+ X
  78  ST+ X
  79  RCL 00
  80  ENTER
  81  SIGN
  82  RCL M        
  83  -
  84  RCL O
  85  INT
  86  -
  87  SIGN
  88  LBL 03
  89  CLX
  90  RCL IND Z
  91  LASTX
  92  *
  93  ST+ IND Y
  94  DSE Z
  95  DSE Y
  96  GTO 03
  97  RCL O
  98  INT
  99  RCL M
100  -
101  ST/A
102  34
103  AMOVE       
104  13
105  AMOVE
106  GTO 02
107  LBL 00
108  RCL 00
109  X/E3
110  ISG X
111  CLA
112  END

 
      ( 193 bytes / SIZE 4n+1 )
 
 

      STACK        INPUTS      OUTPUTS
           Y            m            /
           X            n         1.nnn

   where   1.nnn   is the control number of the result.

Example:      m = 3  ,  n = 7  ,  a = 1 + i / 2 + j / 3 + k / 4             ( quaternion ->  4  STO 00 )

  •  CF 03   Legendre Function of type 2

    1   STO 01
    2   1/X   STO 02
    3   1/X   STO 03
    4   1/X   STO 04

    3   ENTER^
    7   XEQ "PMNA"  >>>>  1.004                               ---Execution time = 20s---

   R01 = -12188.53940
   R02 = -8670.127580
   R03 = -5780.085053
   R04 = -4335.063790

     P37 ( 1 + i/2 + j/3 + k/4 ) = -12188.53940 - 8670.127580 i - 5780.085053 j - 4335.063790 k
 

  •  SF 03   Legendre Function of type 3

    1   STO 01
    2   1/X   STO 02
    3   1/X   STO 03
    4   1/X   STO 04

    3   ENTER^
    7   XEQ "PMNA"  >>>>  1.004                               ---Execution time = 20s---

   R01 = 11285.97688
   R02 = -9363.495330
   R03 = -6242.330218
   R04 = -4681.747665

     P37 ( 1 + i/2 + j/3 + k/4 ) = 11285.97688 - 9363.495330 i - 6242.330218 j - 4681.747665 k

Notes:

-This program doesn't check if m & n are integers that satisfy 0 <= m <= n
  Pn-1m(a) is in registers R3n+1 thru R4n.
 

10°)  Struve Functions
 

Formulae:

  CF 01        Hm(q) = (a/2)m+1 1F2( 1 ; 3/2 , m + 3/2 ; - a2/4 ) / Gam( m+3/2 )                with    m+3/2 # 0 , -1 , -2 , ...................
  SF 01         Lm(q) = (a/2)m+1 1F2( 1 ; 3/2 , m + 3/2 ;  a2/4 ) / Gam(m+3/2)
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "STRA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R3n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  F01-F04
Subroutines:  "ST*A"  "ST/A"     ( cf paragraph 0 above )
                         "HGFA"                ( cf paragraph 8 above )
                         "A*A1"  "A^X"     ( cf "Anions for the HP-41" )
                         "1/G+"                   ( cf "Gamma Function for the HP-41" )
 

-Line 33 is  1/Gam(3/2)
 
 

 01  LBL "STRA"
 02  SF 04
 03  STO M
 04  2
 05  XEQ "ST/A" 
 06  RCL 00
 07  .1
 08  %
 09  STO Z
 10  ST+ Z
 11  ISG X
 12  +
 13   E3
 14  /
 15  1
 16  +
 17  REGMOVE 
 18  +
 19  REGMOVE
 20  SIGN
 21  CHS
 22  FC? 01
 23  XEQ "ST*A"
 24  XEQ "A*A1"
 25  SIGN
 26  1.5
 27  ST+ M
 28  RCL M
 29  XEQ "HGFA"
 30  RDN
 31  STO M
 32  XEQ "1/G+"
 33  1.128379167
 34  *
 35  XEQ "ST*A"
 36  RCL 00
 37   E3
 38  /
 39  ISG X
 40   E3
 41  /
 42  1
 43  +
 44  RCL 00 
 45  3
 46  *
 47  +
 48  STO N
 49  REGSWAP
 50  RCL M
 51  .5
 52  -
 53  XEQ "A^X"
 54  FRC
 55  RCL N
 56  +
 57  REGMOVE
 58  XEQ "A*A1"
 59  CLA
 60  END

 
    ( 142 bytes / SIZE 3n+1 )
 
 

      STACK        INPUT      OUTPUT
           X             m         1.nnn

   where   1.nnn   is the control number of the result.

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

   •   CF 01

    1   STO 01
    2   STO 02
    3   STO 03
    4   STO 04

   PI   XEQ "STRA"   >>>>  1.004                                         ---Execution time = 29s---

    R01 =  8.546633232
    R02 = -4.085079933
    R03 = -6.127619898
    R04 =  -8.170159864

   So,   Hpi( 1 + 2 i + 3 j + 4 k ) = 8.546633232 - 4.085079933 i - 6.127619898 j - 8.170159864 k

   •   SF 01

   1   STO 01
    2   STO 02
    3   STO 03
    4   STO 04

   PI   XEQ "STRA"   >>>>  1.004                                         ---Execution time = 30s---

   R01 =  1.864582784
   R02 = -0.116220198
   R03 = -0.174330297
   R04 = -0.232440397

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

Notes:

-Lines 29 to 35 may be replaced by   1.002  CHS  XEQ "HGFB",
 where "HGFB" is listed in "Anionic Special Functions(II)" paragraph 6°b)
-In  this case, "STRA" will also work if  m+3/2 = 0 , -1 , -2 , ...................
 

11°)  Bessel Functions
 

- "BSLA" computes the Bessel functions of the 1st kind if flag  F02 is clear and the Bessel functions of the 2nd kind if flag  F02  is set.

-Set flag  F03  to get the modified  Bessel functions .
 

Formulae:
 

  Bessel Functions of the 1st kind

       Jm(a) = (a/2)m [ 1/Gam(m+1)  +  (-a2/4)1/ (1! Gam(m+2) )  + .... + (-a2/4)k/ (k! Gam(m+k+1) ) + ....  ]      m # -1 ; -2 ; -3 ; ....

  Modified Bessel Functions of the 1st kind

       Im(a) = (a/2)m [ 1/Gam(m+1)  +  (a2/4)1/ (1! Gam(m+2) )  + .... + (a2/4)k/ (k! Gam(m+k+1) ) + ....  ]          m # -1 ; -2 ; -3 ; ....

  Bessel Functions & Modified Bessel Functions of the 2nd kind - non-integer order

        Ym(a) = ( Jm(a) cos(m(pi)) - J-m(a) ) / sin(m(pi))      ;      Km(a) = (pi/2) ( I-m(a) - Im(a) ) / sin(m(pi))           m # .... -3 ; -2 ; -1 ; 0 ; 1 ; 2 ; 3 ....

  Bessel Functions & Modified Bessel Functions of the 2nd kind - non-negative integer order

         Ym(a) = -(1/pi) (a/2)-m SUMk=0,1,....,m-1  (m-k-1)!/(k!) (a2/4)k + (2/pi) ln(a/2) Jm(a)
                       - (1/pi) (a/2)m SUMk=0,1,.....  ( psi(k+1) + psi(m+k+1) ) (-a2/4)k / (k!(m+k)!)                         where   psi = the digamma function

         Km(a) = (1/2) (a/2)-m SUMk=0,1,..,m-1  (m-k-1)!/(k!) (-a2/4)k - (-1)m ln(a/2) Im(a)
                    + (1/2) (-1)m (a/2)m SUMk=0,1,...( psi(k+1) + psi(m+k+1) ) (a2/4)k / (k!(m+k)!)
 
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "BSLA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R3n:  temp    ( R3n+1 to R4n are also used for the functions of the second kind )

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  F02-F03
Subroutines:  "ST*A"  "ST/A"  "DSA"    ( cf paragraph 0 above )
                         "A*A1"  "A^X"  "LNA"   ( cf "Anions for the HP-41" )
                         "1/G+"                             ( cf "Gamma Function for the HP-41" )

-Lines 52 & 130 are three-byte  GTO 09
-Line 207 = 2 x Euler's gamma constant

-Lines 218 to 238 may be replaced by  RCL M   DS*A  where  DS*A  is an M-Code routine listed in paragraph 0°)b)
 
 

  01  LBL "BSLA"
  02  2
  03  XEQ "ST/A"
  04  X<>Y
  05  FS? 02
  06  GTO 02
  07  LBL 10
  08  STO N
  09  RCL 00
  10  .1
  11  %
  12  STO O
  13  ISG X
  14  +
  15   E3
  16  /
  17  1
  18  +
  19  STO M
  20  ST+ O
  21  REGMOVE
  22  XEQ "A*A1"
  23  CLX
  24  X<> M
  25  REGSWAP
  26  RCL N
  27  XEQ "A^X"
  28  RCL N
  29  1
  30  +
  31  XEQ "1/G+"
  32  XEQ "ST*A"
  33  RCL O
  34  REGMOVE
  35  LBL 01
  36  RCL N
  37  RCL M
  38  1
  39  +
  40  STO M
  41  ST+ Y
  42  *
  43  FC? 03
  44  CHS
  45  XEQ "ST/A"
  46  XEQ "A*A1"
  47  XEQ "DSA"
  48  X#0?
  49  GTO 01
  50  RCL O
  51  REGSWAP
  52  GTO 09
  53  LBL 02
  54  RCL 00
  55   E3
  56  /
  57  ISG X
  58  RCL 00
  59  3
  60  *
  61  +
  62   E3
  63  /
  64  1
  65  +
  66  REGMOVE
  67  X<>Y
  68  XEQ 10
  69  X<>Y
  70  STO N
  71  FRC
  72  X=0?
  73  GTO 04
  74  FS? 03
  75  GTO 00
  76  LASTX
  77  1
  78  CHS
  79  ACOS
  80  *
  81  COS
  82  XEQ "ST*A"
  83  LBL 00
  84  RCL 00
  85   E3
  86  /
  87  ISG X
  88  RCL 00
  89  3
  90  *
  91  +
  92   E3
  93  /
  94  1
  95  +
  96  REGSWAP
  97  RCL N
  98  CHS
  99  XEQ 10
100  X<>Y
101  CHS
102  STO N
103  RCL 00
104  4
105  *
106  RCL 00
107  LBL 03
108  RCL IND Y
109  RCL IND Y
110  -
111  STO IND Y
112  RDN
113  DSE Y
114  DSE X
115  GTO 03
116  PI
117  2
118  /
119  CHS
120  FC? 03
121  1
122  RCL N
123  1
124  CHS
125  ACOS
126  *
127  SIN
128  /
129  XEQ "ST*A"
130  GTO 09
131  LBL 04
132  RCL 00
133  .1
134  %
135  ISG X
136  STO Z
137  +
138   E3
139  /
140  1
141  +
142  REGMOVE
143  CLX
144   E3
145  /
146  1
147  +
148  RCL 00
149  3
150  *
151  +
152  REGMOVE
153  XEQ "LNA"
154  1
155  FS? 03
156  CHS
157  RCL N
158  STO M
159  Y^X
160  ST+ X
161  XEQ "ST*A"
162  XEQ "A*A1"
163  STO Y
164  RCL 00
165  ST+ X
166  +
167   E3
168  /
169  1
170  +
171  REGMOVE
172  X<>Y
173  FRC
174  LASTX
175   E3
176  /
177  RCL 00
178  3
179  *
180  +
181  1
182  +
183  REGMOVE
184  +
185  REGMOVE
186  XEQ "A*A1"
187  RCL 00
188  +
189   E3
190  /
191  1
192  +
193  REGSWAP
194  RCL N
195  XEQ "A^X"
196  RCL N
197  FACT
198  1
199  FS? 03
200  CHS
201  LASTX
202  Y^X
203  *
204  XEQ "ST/A"
205  CLX
206  STO O
207  1.15443133
208  X<> M
209  ABS
210  LBL 05
211  X#0?
212  1/X
213  ST- M
214  X<> L
215  DSE X
216  GTO 05
217  LBL 06
218  RCL 00
219  3
220  *
221  STO P
222  RCL 00
223  ENTER^
224  CLX
225  LBL 07
226  RCL IND Y
227  RCL M
228  *
229  RCL IND P
230  +
231  STO IND P
232  LASTX
233  -
234  ABS
235  +
236  DSE P
237  DSE Y
238  GTO 07
239  X=0?
240  GTO 00
241  XEQ "A*A1"
242  SIGN
243  RCL O
244  +
245  STO O
246  STO Y
247  RCL N
248  +
249  *
250  FC? 03
251  CHS
252  XEQ "ST/A"
253  RCL M
254  LASTX
255  1/X
256  -
257  RCL O
258  1/X
259  -
260  STO M
261  GTO 06
262  LBL 00
263  RCL 00
264   E3
265  /
266  ISG X
267   E3
268  /
269  1
270  +
271  RCL 00
272  3
273  *
274  +
275  REGMOVE
276  RCL N
277  X=0?
278  GTO 00
279  CHS
280  XEQ "A^X"
281  RCL N
282  1
283  -
284  FACT
285  CHS
286  XEQ "ST*A"
287  CLX
288  STO O
289  XEQ "DSA"
290  LBL 08
291  RCL N
292  RCL O
293  1
294  +
295  STO O
296  X=Y?
297  GTO 00
298  STO Z
299  -
300  *
301  FS? 03
302  CHS
303  XEQ "ST/A"
304  XEQ "A*A1"
305  XEQ "DSA"
306  GTO 08
307  LBL 00
308  RCL 00
309   E3
310  /
311  ISG X
312   E3
313  /
314  1
315  +
316  RCL 00
317  ST+ X
318  +
319  REGMOVE
320  PI
321  FS? 03
322  -2
323  XEQ "ST/A"
324  LBL 09
325  RCL N
326  RCL 00
327   E3
328  /
329  ISG X
330  CLA
331  END

 
      ( 593 bytes / SIZE 3n+1 or 4n+1 )
 
 

      STACK        INPUT      OUTPUT
           X             m         1.nnn

   where   1.nnn   is the control number of the result.

Example1:         m = PI     a = 1 + 2 i + 3 j + 4 k        ( quaternion ->  4  STO 00 )

   •   CF 02  CF 03  Bessel Functions of the 1st kind    J

   1   STO 01
   2   STO 02
   3   STO 03
   4   STO 04

   PI   XEQ "BSLA"   >>>>   1.004                                            ---Execution time = 20s---

   R01 = -11.22987514
   R02 =  -3.570801501
   R03 =  -5.356202256
   R04 =  -7.141603008

-So,   JPI( 1 + 2 i + 3 j + 4 k ) = -11.22987514 - 3.570801501 i - 5.356202256 j - 7.141603008 k

   •   CF 02  SF 03  Modified  Bessel Functions of the 1st kind    I

   1   STO 01
   2   STO 02
   3   STO 03
   4   STO 04

   PI   XEQ "BSLA"   >>>>   1.004                                            ---Execution time = 22s---

   R01 =  0.315197913
   R02 = -0.129988660
   R03 = -0.194982989
   R04 = -0.259977320

-Whence,   IPI( 1 + 2 i + 3 j + 4 k ) = 0.315197913 - 0.129988660 i - 0.194982989 j - 0.259977320 k
 

Example2:         m = PI     a = 1 + 2 i + 3 j + 4 k

   •   SF 02  CF 03  Bessel Functions of the 2nd kind   Y

   1   STO 01
   2   STO 02
   3   STO 03
   4   STO 04

   PI   XEQ "BSLA"   >>>>   1.004                                            ---Execution time = 50s---

   R01 =   9.617564021
   R02 = -4.171358834
   R03 = -6.257038267
   R04 = -8.342717698

   So,   YPI( 1 + 2 i + 3 j + 4 k ) = 9.617564021 - 4.171358834 i - 6.257038267 j - 8.342717698 k

   •   SF 02  SF 03  Modified Bessel Functions of the 2nd kind   K

   1   STO 01
   2   STO 02
   3   STO 03
   4   STO 04

   PI   XEQ "BSLA"   >>>>   1.004                                            ---Execution time = 51s---

   R01 =  0.203067071
   R02 = -0.055030894
   R03 = -0.082546339
   R04 = -0.110061787

-Whence,   KPI( 1 + 2 i + 3 j + 4 k ) = 0.203067071 - 0.055030894 i - 0.082546339 j - 0.110061787 k
 

Example3:         m = 3     a = 1 + 2 i + 3 j + 4 k

   •   SF 02  CF 03  Bessel Functions of the 2nd kind   Y

   1   STO 01
   2   STO 02
   3   STO 03
   4   STO 04

   3   XEQ "BSLA"   >>>>   1.004                                            ---Execution time = 80s---

   R01 =  7.690027219
   R02 = -5.224812672
   R03 = -7.837219001
   R04 = -10.44962534

   So,   Y3( 1 + 2 i + 3 j + 4 k ) = 7.690027219 - 5.224812672 i - 7.837219001 j - 10.44962534 k

   •   SF 02  SF 03  Modified Bessel Functions of the 2nd kind   K

   1   STO 01
   2   STO 02
   3   STO 03
   4   STO 04

   3   XEQ "BSLA"   >>>>   1.004                                            ---Execution time = 84s---

   R01 =  0.208767801
   R02 = -0.047983196
   R03 = -0.071974795
   R04 = -0.095966394

-Whence,   K3( 1 + 2 i + 3 j + 4 k ) = 0.208767801 - 0.047983196 i - 0.071974795 j - 0.095966394 k

-Likewise, we find:

   Y0( 1 + 2 i + 3 j + 4 k ) = 29.92977302 + 8.767829488 i + 13.15174422 j + 17.53565898 k    ( in 75 seconds )
   K0( 1 + 2 i + 3 j + 4 k ) = 0.190884048 + 0.016289912 i + 0.024434868 j + 0.032579827 k    ( in 85 seconds )

Notes:

-If you have used the M-Code routine  DS*A , the execution times become about 60 seconds instead of about 80s in example(s) 3.
-"BSLA" does not work if m is a negative integer, but we can employ the relations:

      Jm  =  (-1)m  J-m                         Im  =  I-m
      Ym =  (-1)m  Y-m                       Km =  K-m

-SIZE 3n+1 is enough for the functions of the 1st kind.
-The functions of the 2nd kind require SIZE 4n+1.
 

12°)  Regular Coulomb Wave Function
 

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

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

-"RCWFA" also uses   | 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:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "RCWFA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R3n+1:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result

Flag:  F10
Subroutines:  "A*A1"  "A^X"   ( cf "Anions for the HP-41" )             ->  For "A*A1" , use a version that does not alter synthetic register P
 

-Lines 113 to 134 may be replaced by  RCL O   DS*A  where  DS*A  is listed in paragraph 0°)b)
 
 
 

  01  LBL "RCWFA"
  02  STO N
  03  X<>Y
  04  STO M
  05  ENTER^
  06  X=0?
  07  GTO 00
  08  SIGN
  09  LBL 01
  10  RCL N
  11  X^2
  12  RCL Z
  13  X^2
  14  +
  15  *
  16  DSE Y
  17  GTO 01
  18  LBL 00
  19  X=0?
  20  SIGN
  21  RCL N
  22  PI
  23  *
  24  STO O
  25  ENTER^
  26  X=0?
  27  ISG Y
  28  LBL 00
  29  E^X-1
  30  LASTX 
  31  CHS
  32  E^X-1
  33  -
  34  2
  35  /
  36  X=0?
  37  SIGN
  38  /
  39  *
  40  SQRT
  41  2
  42  ST/ O
  43  RCL M           
  44  Y^X
  45  *
  46  RCL O
  47  E^X
  48  /
  49  RCL M
  50  ST+ X
  51  1
  52  +
  53  FACT
  54  /
  55  STO O
  56  RCL 00
  57  .1
  58  %
  59  STO Z
  60  +
  61  STO Z
  62  ST+ Z
  63  1
  64  +
  65   E3
  66  /
  67  1
  68  +
  69  REGMOVE     
  70  RDN
  71  ISG X
  72  CLRGX
  73  +
  74  CLRGX
  75  RCL 00
  76  3
  77  *
  78  1
  79  STO 01
  80  ST+ Y
  81  STO IND Y
  82  RCL O
  83  STO IND T
  84  CLX
  85  STO P
  86  LBL 02
  87  SF 10
  88  LBL 03
  89  RCL 00
  90  PI
  91  INT
  92  *
  93  ISG X
  94  CLX
  95  RCL N
  96  ST+ X
  97  RCL O
  98  ST* Y
  99  X<> P
100  -
101  RCL IND Y    
102  ST/ Y
103  RCL M
104  ST+ X
105  +
106  PI
107  SIGN
108  ST+ IND T
109  +
110  /
111  STO O
112  XEQ "A*A1"
113  RCL 00
114  PI
115  INT
116  *
117  STO Q
118  RCL 00
119  ENTER^
120  CLX
121  LBL 04
122  RCL IND Y
123  RCL O
124  *
125  RCL IND Q 
126  +
127  STO IND Q    
128  LASTX
129  -
130  ABS
131  +
132  DSE Q
133  DSE Y
134  GTO 04
135  X#0?
136  GTO 02
137  FS?C 10
138  GTO 03
139  RCL 00
140  ST+ X
141   E3
142  /
143  ISG X
144   E3
145  /
146  1
147  +
148  RCL 00
149  +
150  REGMOVE    
151  RCL M
152  1
153  +
154  XEQ "A^X"
155  XEQ "A*A1" 
156  CLA
157  END

 
      ( 241 bytes / SIZE 3n+2 )
 
 

      STACK        INPUTS      OUTPUTS
           Y            L            /
           X            h         1.nnn

   where   1.nnn   is the control number of the result.

Example:         L = 2  ,   h = 0.75  ;   a = 1 + 1.2 i + 1.4 j + 1.6 k        ( quaternion ->  4  STO 00 )

   1    STO 01
  1.2  STO 02
  1.4  STO 03
  1.6  STO 04

   2     ENTER^
0.75   XEQ "RCWFA"   >>>>  1.004                                              ---Execution time = 69s---

  R01 =  -0.478767860
  R02 =  -0.179694771
  R03 =  -0.209643900
  R04 =  -0.239593028

-So,   F2( 0.75 ;  1 + 1.2 i + 1.4 j + 1.6 k ) = -0.478767860 - 0.179694771 i - 0.209643900 j - 0.239593028 k

Notes:

  L must be a non-negative integer.
  If you employ the M-Code routine  DS*A , the execution time becomes 39 seconds.

An example with an octonion:

        L = 3  ,  h = 0.75  ;   a = 1 + 1.1 e1 + 1.2 e2 + 1.3 e3  + 1.4 e4 + 1.5 e5 + 1.6 e6  + 1.7 e7

   8   STO 00

   1   STO 01   1.1  STO 02  .........................  1.7  STO 08

   3     ENTER^
0.75   XEQ "RCWFA"   >>>>  1.008                                              ---Execution time = 60s---      ( with DS*A )

  R01 =   1.025562455
  R02 =  -0.294428914
  R03 =  -0.321195179
  R04 =  -0.347961444
  R05 =  -0.374727708
  R06 =  -0.401493973
  R07 =  -0.428260238
  R08 =  -0.455026503

-So,   F3( 0.75 ;  1 + 1.1 e1 + 1.2 e2 + 1.3 e3  + 1.4 e4 + 1.5 e5 + 1.6 e6  + 1.7 e7 ) =

   1.025562455 - 0.294428914 e1 - 0.321195179 e2 - 0.347961444 e3 - 0.374727708 e4 - 0.401493973 e5 - 0.428260238 e6 - 0.455026503 e7
 

13°) UltraSpherical Functions
 
 

 Formula:    Assuming  l # 0

          Cm(l)(q) = [ Gam(m+2l) / Gam(m+1) / Gam(2.l) ]  2F1( -m , m+2l , m+1/2 , (1-a)/2 )

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

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "USFA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R3n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flag:  F04
Subroutines:  "ST*A"      ( cf paragraph 0 above )
                         "HGFA"    ( cf paragraph 8 above )
                         "1/G+"      ( cf "Gamma Function for the HP-41" )
 
 

 01  LBL "USFA"
 02  .5
 03  CHS
 04  XEQ "ST*A" 
 05  ST- 01
 06  RDN
 07  ENTER^
 08  CHS
 09  X<> Z
 10  ST+ Y
 11  ST+ Y
 12  .5
 13  +
 14  CF 04
 15  XEQ "HGFA"
 16  RDN
 17  STO M
 18  RDN
 19  STO N
 20  1
 21  RCL Z
 22  -
 23  XEQ "1/G+" 
 24  X<> M
 25  ST+ X
 26  1
 27  -
 28  XEQ "1/G+" 
 29  ST* M
 30  RCL N
 31  XEQ "1/G+" 
 32  ST/ M
 33  RCL M
 34  XEQ "ST*A"
 35  RCL 00
 36   E3
 37  /
 38  ISG X
 39  CLA
 40  END

 
    ( 97 bytes / SIZE 3n+1 )
 
 

      STACK        INPUTS      OUTPUTS
           Y            l            /
           X            m         1.nnn

   where   1.nnn   is the control number of the result.

Example:      l = sqrt(2)    m = sqrt(3)    a = 1 + i / 2 + j / 3 + k / 4             ( quaternion ->  4  STO 00 )

    1   STO 01
    2   1/X   STO 02
    3   1/X   STO 03
    4   1/X   STO 04

    2   SQRT
    3   SQRT   XEQ "USFA"  >>>>  1.004                               ---Execution time = 34s---

    R01 = 3.241115513
    R02 = 4.848277899
    R03 = 3.232185264
    R04 = 2.424138948

-Whence      Csqrt(3)sqrt(2) ( 1 + i/2 + j/3 + k/4 ) =  3.241115513 + 4.848277899 i + 3.232185264 j + 2.424138948 k
 

14°)  Jacobi Functions
 

Formula:           Pm(a;b) (a) = [ Gam(a+m+1) / Gam(m+1) ]  2F~1 ( -m , a+b+m+1 , a+1 , (1-a)/2 )

                            where  2F1 tilde is the regularized hypergeometric function
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "JCFA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R3n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  F00-F04
Subroutines:  "ST*A"      ( cf paragraph 0 above )
                         "HGFA"    ( cf paragraph 8 above )
                         "1/G+"      ( cf "Gamma Function for the HP-41" )
 
 

 01  LBL "JCFA"
 02  STO M
 03  CLX
 04  .5
 05  CHS
 06  XEQ "ST*A" 
 07  ST- 01
 08  X<> M
 09  CHS
 10  STO T
 11  -
 12  X<>Y
 13  1
 14  +
 15  ST+ Y
 16  CF 04
 17  XEQ "HGFA"
 18  RDN
 19  STO M
 20  RDN
 21  X<>Y
 22  STO O
 23  FS?C 00
 24  GTO 00
 25  R^
 26  XEQ "1/G+"
 27  XEQ "ST*A" 
 28  LBL 00
 29  1
 30  RCL O
 31  ST- M
 32  -
 33  XEQ "1/G+" 
 34  X<> M
 35  XEQ "1/G+"
 36  ST/ M
 37  RCL M
 38  XEQ "ST*A" 
 39  RCL 00
 40   E3
 41  /
 42  ISG X
 43  CLA
 44  END

 
       ( 106 bytes / SIZE 3n+1 )
 
 

      STACK        INPUTS      OUTPUTS
           Z            a            /
           Y            b            /
           X            m         1.nnn

   where   1.nnn   is the control number of the result.

Example:      a = sqrt(2)    b = sqrt(3) m = PI                 a = 1 + i / 2 + j / 3 + k / 4             ( quaternion ->  4  STO 00 )

    1   STO 01
    2   1/X   STO 02
    3   1/X   STO 03
    4   1/X   STO 04

    2   SQRT
    3   SQRT
   PI   XEQ "JCFA"  >>>>  1.004                               ---Execution time = 33s---

    R01 = -10.14397527
    R02 =  11.74139999
    R03 =   7.827599996
    R04 =   5.870699997

-So,        PPIsqrt(2);sqrt(3) ( 1 + i/2 + j/3 + k/4 ) =  -10.14397527 + 11.74139999 i + 7.827599996 j + 5.870699997 k
 

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

Notes:

-Unless the function is a polynomial, the series is convergent if  | 1 - a | < 2
-Lines 23 thru 28 may be deleted if you replace lines 16-17 by  2.001  CHS  XEQ "HGFB"
  where "HGFB" is listed in "Anionic Special Functions(II)" paragraph 6°b)
 

15°) Toronto Functions
 

-The Toronto functions are defined by

   T(m,n;a) = exp(-a2) [ Gam((m+1)/2) / n ! ] a2n-m+1  M( (m+1)/2 , n+1 , a2 )                    where M = Kummer's function
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "TORA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R4n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  /
Subroutines:  "ST*A"                           ( cf paragraph 0 above )
                         "A*A1"  "A^X"  "E^A"   ( cf "Anions for the HP-41" )
                         "KUMA"                        ( cf paragraph 5 above )
                         "1/G+"                            ( cf "Gamma Function for the HP-41" )
 
 

 01  LBL "TORA" 
 02  1
 03  ST+ Z
 04  +
 05  STO M
 06  X<>Y
 07  2
 08  /
 09  STO N
 10  RCL 00
 11  .1
 12  %
 13  STO Z
 14  ST+ Z
 15  ISG X
 16  +
 17   E3
 18  /
 19  1
 20  +
 21  REGMOVE
 22  +
 23  REGMOVE
 24  XEQ "A*A1"
 25  RCL N
 26  RCL M
 27  XEQ "KUMA"
 28  RDN
 29  STO M
 30  X<>Y
 31  STO N
 32  XEQ "1/G+"
 33  STO O
 34  RCL M
 35  XEQ "1/G+" 
 36  RCL O
 37  /
 38  XEQ "ST*A" 
 39  RCL 00
 40  .1
 41  %
 42  ISG X
 43  +
 44   E3
 45  /
 46  1
 47  +
 48  STO O
 49  REGSWAP
 50  SIGN
 51  CHS
 52  XEQ "ST*A" 
 53  XEQ "E^A"
 54  XEQ "A*A1" 
 55  RCL O
 56  REGMOVE
 57  RCL 00
 58   E3
 59  /
 60  ST+ X
 61  +
 62  REGSWAP
 63  RCL M
 64  RCL N
 65  -
 66  ST+ X
 67  XEQ "A^X"
 68  XEQ "A*A1" 
 69  CLA
 70  END

 
    ( 157 bytes / SIZE 4n+1 )
 
 

      STACK        INPUTS      OUTPUTS
           Y            m            /
           X            n         1.nnn

   where   1.nnn   is the control number of the result.

Example:      m = sqrt(2)    n = sqrt(3)                  a = 1 + i / 2 + j / 3 + k / 4             ( quaternion ->  4  STO 00 )

    1   STO 01
    2   1/X   STO 02
    3   1/X   STO 03
    4   1/X   STO 04

    2   SQRT
    3   SQRT   XEQ "USFA"  >>>>  1.004                               ---Execution time = 33s---

    R01 = 0.321293450
    R02 = 0.465223023
    R03 = 0.310148682
    R04 = 0.232611511

-Whence          T( sqrt(2) ,sqrt(3) ; 1 + i/2 + j/3 + k/4 )  = 0.321293450 + 0.465223023 i + 0.310148682 j + 0.232611511 k

Note:

-"TORA" does not work if  n  is a negative integer.
-But it will work if you replace:

  lines 33 to 38  by  XEQ "ST/A"
  line  27  by  1.001  CHS  XEQ "HGFB"     ( cf "Anionic Special Functions(II)" paragraph 6°)b) )
 

16°) Weber & Anger Functions
 

 "WEBA" calculates Weber functions if flag F01 is clear and Anger's functions if flag F01 is set
 

Formulae:
 

   Em(a) = - (a/2) cos ( 90°m )  1F~2( 1 ; (3-m)/2 , (3+m)/2 ; -a2/4 )         Weber's functions

                    + sin ( 90°m ) 1F~2( 1 ; (2-m)/2 , (2+m)/2 ; -a2/4 )

  Jm(a) = + (a/2) sin ( 90°m ) 1F~2( 1 ; (3-m)/2 , (3+m)/2 ; -a2/4 )      Anger's functions

                  + cos ( 90°m ) 1F~2( 1 ; (2-m)/2 , (2+m)/2 ; -a2/4 )
 

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "WEBA" )

                                      •  R01   ......  •  Rnn = the n components of the anion

                                         Rn+1 ..........  R4n:  temp

  >>>  When the program stops:     R01   ......  Rnn = the n components of the result
 

Flags:  F00-F01-F04                 CF 01 for Weber's function
                                                    SF 01 for Anger's function

Subroutines:  "ST*A"  "ST/A"    ( cf paragraph 0 above )
                         "HGFA"               ( cf paragraph 8 above )
                         "A*A1"              ( cf "Anions for the HP-41" )
                         "1/G+"       ( cf "Gamma Function for the HP-41" )
 
 

  01  LBL "WEBA"
  02  STO N
  03  2
  04  XEQ "ST/A"
  05  RCL 00
  06  .1
  07  %
  08  STO Z
  09  ST+ Z
  10  ISG X
  11  +
  12   E3
  13  /
  14  1
  15  +
  16  REGMOVE
  17  +
  18  REGMOVE
  19  SIGN
  20  CHS
  21  XEQ "ST*A"
  22  XEQ "A*A1"
  23  SIGN
  24  RCL N
  25  ENTER^
  26  CHS
  27  3
  28  ST+ Z
  29  +
  30  2
  31  ST/ Z
  32  /
  33  SF 04
  34  XEQ "HGFA"
  35  RDN
  36  STO M
  37  X<>Y
  38  STO N
  39  X<>Y
  40  -
  41  1
  42  ASIN
  43  *
  44  FS? 01
  45  SIN
  46  FC? 01
  47  COS
  48  STO O
  49  RCL M
  50  XEQ "1/G+"
  51  ST* O
  52  RCL N
  53  XEQ "1/G+"
  54  RCL O
  55  *
  56  FC? 01
  57  CHS
  58  XEQ "ST*A" 
  59  RCL 00
  60   E3
  61  /
  62  ISG X
  63  RCL 00
  64  3
  65  *
  66  +
  67   E3
  68  /
  69  1
  70  +
  71  STO O
  72  RCL 00
  73  +
  74  REGSWAP
  75  XEQ "A*A1"
  76  RCL O
  77  REGSWAP
  78  SIGN
  79  RCL N
  80  RCL M
  81  .5
  82  ST- Z
  83  -
  84  XEQ "HGFA"
  85  RDN
  86  STO M
  87  X<>Y
  88  STO N
  89  X<>Y
  90  -
  91  1 
  92  ASIN
  93  *
  94  FS? 01
  95  COS
  96  FC? 01
  97  SIN
  98  STO O
  99  RCL M
100  XEQ "1/G+"
101  ST* O
102  RCL N
103  XEQ "1/G+"
104  RCL O
105  *
106  XEQ "ST*A" 
107  4
108  RCL 00
109  ST* Y
110  LBL 01
111  RCL IND Y
112  ST+ IND Y
113  RDN
114  DSE Y
115  DSE X
116  GTO 01
117  RCL 00
118   E3
119  /
120  ISG X
121  CLA
122  CF 04
123  END

 
       ( 243 bytes / SIZE 4n+1 )
 
 

      STACK        INPUT      OUTPUT
           X             m         1.nnn

   where   1.nnn   is the control number of the result.

Example:         m = PI     a = 1 + 2 i + 3 j + 4 k        ( quaternion ->  4  STO 00 )

   •   CF 01  Weber's function

   1   STO 01
   2   STO 02
   3   STO 03
   4   STO 04

   PI   XEQ "WEBA"   >>>>   1.004                                            ---Execution time = 58s---

   R01 = -9.566817335
   R02 =  4.177204885
   R03 =  6.265807323
   R04 =  8.354409766

-So,     EPi( 1 + 2 i + 3 j + 4 k ) = -9.566817335 + 4.177204885 i + 6.265807323 j + 8.354409766 k

   •   SF 01  Anger's function

   1   STO 01
   2   STO 02
   3   STO 03
   4   STO 04

   PI   XEQ "WEBA"   >>>>   1.004                                            ---Execution time = 58s---

   R01 = -11.24234891
   R02 = -3.564968423
   R03 = -5.347452633
   R04 = -7.129936844

-Whence,   JPi( 1 + 2 i + 3 j + 4 k ) = -11.24234891 - 3.564968423 i - 5.347452633 j - 7.129936844 k

Notes:

-This routine does not work if  m = +/-2 , +/-3 , +/-4 , .............................

-However, it will work if you delete lines 98 to 105, if you replace line 84 by  1.002  CHS  XEQ "HGFB",
 if you delete lines 48 to 55 and if you replace lines 33-34 by    1.002  CHS  XEQ "HGFB"
 where "HGFB" is listed in "Anionic Special Functions(II)" paragraph 6°b)