# 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
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
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
28C  ?FSET 7
135   ?NCXQ
060   C=A*C
10E   A=C ALL
0B0   C=N ALL
270   RAMSLCT
01D   C=
060   A+C
10E   A=C ALL
0B0   C=N ALL
270   RAMSLCT
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
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  XY?   36  GTO 02   37  LBL 00   38  RCL N   39  ENTER^   40  FRC   41  X=0?   42  XY?   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)