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)