Bionic Functions

Bionic Special Functions for the HP-41

Overview

0°) Extra Registers & M-Code Routines
1°) Gamma, Digamma, Polygamma Functions
2°) Catalan Numbers
3°) Error Function
4°) Exponential, Sine and Cosine Integrals
5°) Exponential Integral Em

a) Ascending Series
b) Asymptotic Expansion

6°) Bessel Functions
7°) Struve Functions
8°) Jacobian Elliptic Functions
9°) Weber & Anger Functions
10°) Incomplete Gamma Function
11°) Airy Functions
12°) Parabolic Cylinder Functions

a) Ascending Series
b) Asymptotic Expansion

13°) Hermite Functions
14°) Harmonic Numbers
15°) Generalized Error Functions
16°) Lambert-W Function

-These programs generalize similar routines listed in "Anionic Functions for the HP-41" to "bi-onic" arguments i-e complexified anions.
-As usual, the result are not very accurate for large arguments when ascending series are used.

0°) Extra Registers & M-Code Routines

-We need 3 extra-registers that I've called: U V W
-STO W & RCL W are listed in "A few M-Code Routines for the HP-41"
-Registers U and V are coded in the same way.

-Of course, if n is not too large, you can replace these registers by R97 R98 R99

-Other M-Code routines are employed:

Z*Z Z^2 SQRTZ 1/Z Z/Z cf "A few M-Code Routines for the HP-41"

and X+1 X-1 SINH COSH -------------------------------------------

-And M-Code routines that are useful for anions may also be used for bi-ons without any modification:

A+A A-A ST*A ST/A AMOVE ASWAP DSA DS*A cf "Anions" and "Anionic Functions for the HP-41"

1°) Gamma, Digamma, Polygamma Functions

Formula1:

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

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

Formula2: Psi(b) ~ Ln b - 1/(2b) -1/(12b²) + 1/(120b⁴) - 1/(252b⁶) + 1/(240b⁸) is used if Re(b) > 8

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

Formula3:

• If m > 0 , y^(m) (b) ~ (-1)^m-1 [ (m-1)! / b^m + m! / (2.b^m+1) + SUM_k=1,2,.... B_2k (2k+m-1)! / (2k)! / b^2k+m where B_2k are Bernoulli numbers

• If m = 0 , y (b) ~ Ln b - 1/(2b) - SUM_k=1,2,.... B_2k / (2k) / b^2k ( digamma function )

-So, the digamma function may be computed with "PSINB" too.

and the recurrence relation: y^(m) (a+p) = y^(m) (a) + (-1)^m m! [ 1/a^m+1 + ....... + 1/(a+p-1)^m+1 ] where p is a positive integer

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "GAMB" , "PSIB" or "PSINB" )

• R01 ...... • R2n = the n components of the bion

R2n+1 .......... R8n: temp

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

Flag: F10
Subroutines: B*B1 "1/B" "LNB" "E^B" ST*A ST/A A+A A-A DSA AMOVE ASWAP STO W RCL W X+1 X-1 X/E3

-Line 109 = TEXT0
-Line 171 is a three-byte GTO 09

01 LBL "GAMB"
02 12
03 AMOVE
04 CLX
05 ST*A
06 SIGN
07 STO 01
08 LBL 00
09 B*B1
10 RCL 00
11 1
12 ST+ Y
13 ST+ IND Y
14 RCL IND Y
15 5
16 X>Y?
17 GTO 00
18 14
19 AMOVE
20 21
21 AMOVE
22 13
23 AMOVE
24 XEQ "LNB"
25 RCL 00
26 X+1
27 RCL IND X
28 STO W
29 .5
30 ST- IND Z
31 B*B1
32 RCL 00
33 X+1
34 RCL W
35 STO IND Y
36 13
37 ASWAP
38 371
39 195
40 /
41 CF 10
42 XEQ 01
43 210
44 53

45 /
46 XEQ 01
47 30
48 XEQ 01
49 12
50 SF 10
51 LBL 01
52 ST*A
53 XEQ "1/B"
54 RCL 00
55 ENTER^
56 ST+ Y
57 LBL 02
58 RCL IND Y
59 FS? 10
60 CHS
61 ST+ IND Y
62 RDN
63 DSE Y
64 DSE X
65 GTO 02
66 FC?C 10
67 RTN
68 PI
69 ST+ X
70 SQRT
71 LN
72 ST+ 01
73 3
74 RCL 00
75 ST* Y
76 XEQ 02
77 XEQ "E^B"
78 12
79 AMOVE
80 41
81 AMOVE
82 XEQ "1/B"
83 B*B1
84 RTN
85 LBL "PSIB"
86 12
87 AMOVE
88 23

89 AMOVE
90 RCL 00
91 4
92 *
93 .1
94 %
95 ISG X
96 +
97 RCL 00
98 -
99 STO W
100 CLRGX
101 LBL 03
102 XEQ "1/B"
103 RCL W
104 LBL 04
105 RCL IND Y
106 ST- IND Y
107 RDN
108 ISG Y
109 ""
110 ISG X
111 GTO 04
112 SIGN
113 ST+ IND Y
114 21
115 AMOVE
116 RCL 01
117 8
118 X>Y?
119 GTO 03
120 RCL W
121 RCL 00
122 -
123 RCL 01
124 STO IND Y
125 B*B1
126 XEQ "1/B"
127 12
128 AMOVE
129 20
130 ST/A
131 21
132 1/X

133 ST- 01
134 B*B1
135 .1
136 ST+ 01
137 B*B1
138 SIGN
139 ST- 01
140 B*B1
141 6
142 ST/A
143 12
144 AMOVE
145 31
146 AMOVE
147 XEQ "1/B"
148 A-A
149 2
150 CHS
151 ST/A
152 13
153 ASWAP
154 XEQ "LNB"
155 RCL 00
156 4
157 *
158 STO M
159 RCL 00
160 ST- Y
161 LBL 05
162 RCL IND M
163 RCL IND Z
164 +
165 ST+ IND Y
166 RDN
167 DSE M
168 DSE Y
169 DSE X
170 GTO 05
171 GTO 09
172 LBL "PSINB"
173 STO W
174 12
175 AMOVE
176 CLX

177 ST*A
178 13
179 AMOVE
180 LBL 06
181 21
182 AMOVE
183 8
184 RCL W
185 +
186 RCL 01
187 X>Y?
188 GTO 07
189 1
190 LASTX
191 +
192 CHS
193 XEQ "B^X"
194 DSA
195 RCL 00
196 1
197 ST+ Y
198 ST+ IND Y
199 GTO 06
200 LBL 07
201 14
202 AMOVE
203 B*B1
204 XEQ "1/B"
205 12
206 AMOVE
207 RCL W
208 9
209 +
210 FACT
211 39.6
212 /
213 ST*A
214 B*B1
215 RCL W
216 7
217 FACT
218 +
219 ST- 01
220 B*B1

221 40
222 ST/A
223 RCL W
224 5
225 +
226 FACT
227 ST+ 01
228 B*B1
229 42
230 ST/A
231 RCL W
232 3
233 +
234 FACT
235 ST- 01
236 B*B1
237 60
238 ST/A
239 RCL W
240 X+1
241 FACT
242 ST+ 01
243 B*B1
244 6
245 ST/A
246 13
247 ASWAP
248 RCL W
249 FACT
250 ST*A
251 13
252 ASWAP
253 12
254 AMOVE
255 41
256 AMOVE
257 XEQ "1/B"
258 RCL W
259 FACT
260 ST*A
261 A+A
262 2
263 ST/A
264 RCL W

265 X=0?
266 GTO 07
267 X-1
268 FACT
269 ST+ 01
270 12
271 AMOVE
272 41
273 AMOVE
274 RCL W
275 CHS
276 XEQ "B^X"
277 B*B1
278 GTO 08
279 LBL 07
280 12
281 AMOVE
282 41
283 AMOVE
284 XEQ "LNB"
285 A-A
286 SIGN
287 CHS
288 ST*A
289 LBL 08
290 DSA
291 31
292 AMOVE
293 SIGN
294 CHS
295 RCL W
296 Y^X
297 CHS
298 ST*A
299 LBL 09
300 RCL 00
301 X/E3
302 ISG X
303 CLA
304 END

( 579 bytes / SIZE 8n+1 )

STACK	INPUT	OUTPUT
X	/	1.2n

where 1.2n is the control number of the result.

Exception: For Psi^(m) (b) , place m in X-register ( m = non-negative integer )

Example: b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( 0.5 + 0.6 i ) e₂ + ( 0.7 + 0.8 i ) e₃ ( biquaternion -> 8 STO 00 )

0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08

• XEQ "GAMB" >>>> 1.008 ---Execution time = 70s---

And we find in registers R01 thru R08

Gam(b) = ( 0.258452256 + 0.166123162 i ) + ( 0.014416808 + 0.150189863 i ) e₁
+ ( 0.034505408 + 0.233142842 i ) e₂ + ( 0.054594009 + 0.316095822 ) e₃

• Likewise, store again b into R01 to R08 XEQ "PSIB" >>>> 1.008 ---Execution time = 82s---

And the result is in registers R01 thru R08

Psi(b) = ( 0.300919280 + 0.872730291 i ) + ( 0.587618165 - 0.077794679 i ) e₁
+ ( 0.910460763 - 0.168369153 i ) e₂ + ( 1.233303362 - 0.258943626 ) e₃

• To calculate, say Psi⁽³⁾ (b) , store b into R01 to R08 , 3 XEQ "PSINB" >>>> 1.008 ---Execution time = 4m15s---

R01 thru R08 give:

Psi⁽³⁾ (b) = ( 0.349772972 + 0.790426125 i ) + ( -0.281836208 - 0.432414491 i ) e₁
+ ( -0.474257644 - 0.652019709 i ) e₂ + ( -0.666679081 - 0.871624929 ) e₃

Note:

0 XEQ "PSINB" is another way to calculate the digamma function, though it will be slower ( 3m04s instead of 82 seconds in the example above )

2°) Catalan Numbers

Formula: C(b) = 4^b Gam(b+1/2) / Gam(b+2) / sqrt(PI)

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "GAMB" , "PSIB" or "PSINB" )

• R01 ...... • R2n = the n components of the bion

R2n+1 .......... R12n: temp

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

Flag: /
Subroutines: B*B1 "1/B" "X^B" ST*A AMOVE ASWAP "GAMB"

01 LBL "CATB"
02 15
03 AMOVE
04 16
05 AMOVE
06 .5
07 ST+ 01
08 XEQ "GAMB"
09 15
10 ASWAP
11 4
12 XEQ "X^B"
13 52
14 AMOVE
15 B*B1
16 16
17 ASWAP
18 2
19 ST+ 01
20 XEQ "GAMB"
21 PI
22 SQRT
23 ST*A
24 XEQ "1/B"
25 62
26 AMOVE
27 B*B1
28 END

( 72 bytes / SIZE 12n+1)

STACK	INPUT	OUTPUT
X	/	1.2n

where 1.2n is the control number of the result.

Example: b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( 0.5 + 0.6 i ) e₂ + ( 0.7 + 0.8 i ) e₃ ( biquaternion -> 8 STO 00 )

0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08

XEQ "CATB" >>>> 1.008 ---Execution time = 2m23s---

Cat(b) is in registers R01 thru R08:

Cat(b) = ( 0.474992330 - 0.263534256 i ) + ( 0.048669671 + 0.153014311 i ) e₁
+ ( 0.088165831 + 0.234808752 i ) e₂ + ( 0.127661988 + 0.316603194 ) e₃

3°) Error Function

Formula: erf a = (2/pi^1/2) SUM_{n=0,1,2,.....} (-1)ⁿ a²ⁿ⁺¹ / (n! (2n+1))

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "ERFB" )

• R01 ...... • R2n = the n components of the bion

Rn+1 .......... R6n: temp

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

Flags: /
Subroutines: B*B1 ST*A DSA STO W RCL W X+1

01 LBL "ERFB"
02 12
03 AMOVE
04 13
05 AMOVE
06 B*B1
07 12
08 ASWAP

09 CLX
10 STO W
11 LBL 01
12 B*B1
13 RCL W
14 ST+ X
15 X+1
16 RCL W

17 X+1
18 STO W
19 ENTER^
20 ST+ X
21 X+1
22 *
23 /
24 CHS

25 ST*A
26 DSA
27 X#0?
28 GTO 01
29 31
30 AMOVE
31 2
32 PI

33 SQRT
34 /
35 ST*A
36 RCL 00
37 E3
38 /
39 ISG X
40 END

( 74 bytes / SIZE 6n+1 )

STACK	INPUT	OUTPUT
X	/	1.2n

where 1.2n is the control number of the result.

Example: b = ( 1 + 0.9 i ) + ( 0.8 + 0.7 i ) e₁ + ( 0.6 + 0.5 i ) e₂ + ( 0.4 + 0.3 i ) e₃ ( biquaternion -> 8 STO 00 )

1 STO 01 0.9 STO 02 0.8 STO 03 0.7 STO 04 0.6 STO 05 0.5 STO 06 0.4 STO 07 0.3 STO 08

XEQ "ERFB" >>>> 1.008 ---Execution time = 79s---

R01 thru R08 yields:

Erf(b) = ( 2.952728868 + 8.149265562 i ) + ( 6.172244022 - 1.372589195 i ) e₁
+ ( 4.509301553 - 1.117428240 i ) e₂ + ( 2.846359083 - 0.862267287 ) e₃

4°) Exponential, Sine and Cosine Integrals

-The following programs caculates

"EIB" Ei(b) = C + Ln(b) + Sum_n=1,2,..... bⁿ/(n.n!) where C = 0.5772156649... = Euler's constant.

"SIB" Si(b) = Sum_{m=0,1,2,.....} (-1)^m b^2m+1/((2m+1).(2m+1)!) if flag F01 is clear
Shi(b) = Sum_{m=0,1,2,.....} b^2m+1/((2m+1).(2m+1)!) if flag F01 is set

"CIB" Ci(b) = C + ln(b) + Sum_m=1,2,..... (-1)^m b^2m/(2m.(2m)!) if flag F01 is clear
Chi(b)= C + ln(b) + Sum_m=1,2,..... b^2m/(2m.(2m)!) if flag F01 is set

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing these programs )

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

Rn+1 .......... R6n: temp

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

Flag: F01
Subroutines: B*B1 "LNB" ST*A ST/A DSA A+A AMOVE ASWAP STO W RCL W GEU X+1

-Lines 07-45 may be replaced by .5772156649 ( Gamma-Euler Constant )

01 LBL "EIB"
02 12
03 AMOVE
04 13
05 AMOVE
06 XEQ "LNB"
07 GEU
08 ST+ 01
09 A+A
10 13
11 ASWAP
12 SIGN
13 STO W
14 LBL 01
15 B*B1
16 RCL W

17 ENTER^
18 X+1
19 STO W
20 X^2
21 /
22 ST*A
23 DSA
24 X#0?
25 GTO 01
26 GTO 00
27 LBL "SIB"
28 12
29 AMOVE
30 13
31 AMOVE
32 B*B1

33 12
34 ASWAP
35 SIGN
36 STO W
37 GTO 02
38 LBL "CIB"
39 12
40 AMOVE
41 B*B1
42 12
43 ASWAP
44 XEQ "LNB"
45 GEU
46 ST+ 01
47 13
48 AMOVE

49 CLX
50 STO W
51 ST*A
52 SIGN
53 STO 01
54 LBL 02
55 B*B1
56 RCL W
57 ENTER^
58 X+1
59 ENTER^
60 X+1
61 STO W
62 X^2
63 *
64 /

65 X=0?
66 .25
67 FC? 01
68 CHS
69 ST*A
70 DSA
71 X#0?
72 GTO 02
73 LBL 00
74 31
75 AMOVE
76 RCL 00
77 E3
78 /
79 ISG X
80 END

( 162 bytes / SIZE 6n+1 )

STACK	INPUT	OUTPUT
X	/	1.2n

where 1.2n is the control number of the result.

Example: b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( 0.5 + 0.6 i ) e₂ + ( 0.7 + 0.8 i ) e₃ ( biquaternion -> 8 STO 00 )

0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08

• XEQ "EIB" >>>> 1.008 ---Execution time = 52s---

Ei(b) = ( 1.140687091 + 0.557945078 i ) + ( 0.829134747 + 0.443634358 i ) e₁
+ ( 1.328940953 + 0.625738819 i ) e₂ + ( 1.828747160 + 0.808843279 ) e₃

• CF 01 XEQ "SIB" >>>> 1.008 ---Execution time = 25s---

Si(b) = ( 0.029007647 + 0.214825821 i ) + ( 0.254033031 + 0.422973668 i ) e₁
+ ( 0.430129422 + 0.639516280 i ) e₂ + ( 0.606225812 + 0.856058892 ) e₃

• SF 01 XEQ "SIB" >>>> 1.008 ---Execution time = 25s---

Shi(b) = ( 0.168831684 + 0.171286637 i ) + ( 0.343698222 + 0.372373658 i ) e₁
+ ( 0.565959119 + 0.553407049 i ) e₂ + ( 0.788220016 + 0.734440440 ) e₃

• CF 01 XEQ "CIB" >>>> 1.008 ---Execution time = 36s---

Ci(b) = ( 0.822497127 + 1.344232914 i ) + ( 0.534656969 - 0.027912750 i ) e₁
+ ( 0.831831852 - 0.086316448 i ) e₂ + ( 1.129006736 - 0.144720145 ) e₃

• SF 01 XEQ "CIB" >>>> 1.008 ---Execution time = 36s---

Chi(b) = ( 0.971855407 + 0.386658442 i ) + ( 0.485436524 + 0.071260700 i ) e₁
+ ( 0.762981834 + 0.072331770 i ) e₂ + ( 1.040527144 + 0.073402840 ) e₃

5°) Exponential Integral Em

a) Ascending Series

Formulae:

• If m # 1 , 2 , 3 , ..... E_m(b) = b^m-1 Gam(1-m) - [1/(1-m)] M( 1-m , 2-m ; -b ) where M = Kummer's function

• If m = 1 , 2 , 3 , ..... E_m(b) = (-b)^m-1 ( -Ln b - gamma + Sum_k=1,...,m-1 1/k ) / (m-1)! - Sum_k#m-1 (-b)^k / (k-m+1) / k!

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

• E₀(b) = (1/b).exp(-b)

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "ENB" )

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

Rn+1 .......... R6n: temp

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

Flags: /
Subroutines: B*B1   "LNB" "E^B" "B^Z" "B^X" Z*B ST*A   ST/A   DSA DS*A A+A 1/Z Z/Z   AMOVE ASWAP X+1 X-1
                        STO U RCL U STO V RCL V STO W RCL W
                        "GAMZ"    ( a version that does not use any data register - cf "Gamma Function for the HP-41" paragraph 1°) h-3) )

-Line 147 ( GEU ) may be replaced with .5772156649 ( Euler's constant )
-Line 163 may be replaced by E3 /

01 LBL "ENB"
02 X=Y?
03 X#0?
04 GTO 00
05 12
06 AMOVE
07 SIGN
08 CHS
09 ST*A
10 XEQ "E^B"
11 12
12 ASWAP
13 XEQ "1/B"
14 B*B1
15 RTN
16 LBL 00
17 STO U
18 X<>Y
19 STO V
20 X#0?
21 GTO 02
22 X<>Y
23 X<0?
24 GTO 02
25 FRC
26 X=0?
27 GTO 00
28 LBL 02
29 1
30 CHS
31 ST*A
32 12
33 AMOVE

34 CLX
35 STO W
36 ST*A
37 SIGN
38 STO 01
39 13
40 AMOVE
41 LBL 01
42 B*B1
43 RCL V
44 CHS
45 RCL U
46 CHS
47 STO Z
48 X+1
49 RCL W
50 X+1
51 STO W
52 ST+ Y
53 ST+ T
54 ST* Z
55 *
56 RCL V
57 CHS
58 RDN
59 Z/Z
60 Z*B
61 DSA
62 X#0?
63 GTO 01
64 31
65 AMOVE
66 RCL V

67 RCL U
68 X-1
69 STO U
70 1/Z
71 Z*B
72 12
73 ASWAP
74 SIGN
75 CHS
76 ST*A
77 RCL V
78 RCL U
79 XEQ "B^Z"
80 RCL V
81 CHS
82 RCL U
83 CHS
84 XEQ "GAMZ"
85 Z*B
86 A+A
87 RTN
88 LBL 00
89 SIGN
90 CHS
91 ST*A
92 12
93 AMOVE
94 CLX
95 ST*A
96 13
97 AMOVE
98 SIGN
99 STO V

100 STO 01
101 RCL 00
102 ST+ X
103 +
104 RCL U
105 X-1
106 X#0?
107 1/X
108 STO IND Y
109 LBL 03
110 B*B1
111 RCL U
112 RCL V
113 ST/A
114 -
115 LASTX
116 X+1
117 STO V
118 CLX
119 SIGN
120 X=Y?
121 GTO 03
122 -
123 1/X
124 DS*A
125 X#0?
126 GTO 03
127 21
128 AMOVE
129 SIGN
130 CHS
131 ST*A
132 XEQ "LNB"

133 RCL U
134 X-1
135 STO U
136 SIGN
137 CLST
138 LBL 04
139 CLX
140 LASTX
141 X#0?
142 1/X
143 ST+ Y
144 DSE L
145 GTO 04
146 X<>Y
147 GEU
148 -
149 ST- 01
150 RCL U
151 FACT
152 CHS
153 ST/A
154 12
155 ASWAP
156 RCL U
157 XEQ "B^X"
158 B*B1
159 DSA
160 31
161 AMOVE
162 RCL 00
163 X/E3
164 ISG X
165 END

( 292 bytes / ZIZE 6n+1 )

STACK	INPUT	OUTPUT
Y	q	/
X	p	1.2n

where 1.2n is the control number of E_p+i.q (b)

Example: b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( 0.5 + 0.6 i ) e₂ + ( 0.7 + 0.8 i ) e₃ ( biquaternion -> 8 STO 00 )

0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08

-With m = 2 + 3.i

3 ENTER^
2 XEQ "ENB" >>>> 1.008 ---Execution time = 95s---

-And we get in R01 thru R08:

E_2+3.i (b) = ( -0.214869070 - 0.150734076 i ) + ( -0.062956984 + 0.108669920 i ) e₁
+ ( -0.089519302 + 0.174561633 i ) e₂ + ( -0.116081620 + 0.240453347 ) e₃

-Likewise, in 82 seconds: E₃ (b) = ( -0.2411685505 - 0.522115501 i ) + ( -0.214772638 + 0.090964781 i ) e₁
+ ( -0.327768133 + 0.159086869 i ) e₂ + ( -0.440763628 + 0.227208957 ) e₃

Note:

-Don't forget to place 0 in Y-register when m is a real number.

b) Asymptotic Expansion

Formula: E_m(b) ~ (1/b) exp(-b) ₂F₀(1,m;;-1/b)

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "AENB" )

• R01 ...... • R2n = the n components of the bion

Rn+1 .......... R6n: temp

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

Flags: /
Subroutines: B*B1 "1/B" "E^B" Z*B ST*A DSA AMOVE ASWAP X+1
STO U RCL U STO V RCL V STO W RCL W

01 LBL "AENB"
02 STO U
03 X<>Y
04 STO V
05 12
06 AMOVE
07 XEQ "1/B"
08 12

09 ASWAP
10 SIGN
11 CHS
12 ST*A
13 XEQ "E^B"
14 B*B1
15 13
16 AMOVE

17 CLX
18 STO W
19 LBL 01
20 B*B1
21 RCL V
22 CHS
23 RCL U
24 RCL W

25 ST+ Y
26 X+1
27 STO W
28 RDN
29 CHS
30 Z*B
31 DSA
32 X#0?

33 GTO 01
34 31
35 AMOVE
36 RCL 00
37 E3
38 /
39 ISGX
40 END

( 82 bytes / SIZE 6n+1 )

STACK	INPUT	OUTPUT
Y	q	/
X	p	1.2n

where 1.2n is the control number of E_p+i.q (b)

Example: b = ( 30 + 31 i ) + ( 32 + 33 i ) e₁ + ( 34 + 35 i ) e₂ + ( 36 + 37 i ) e₃ ( biquaternion -> 8 STO 00 )

30 STO 01 31 STO 02 32 STO 03 33 STO 04 34 STO 05 35 STO 06 36 STO 07 37 STO 08

-With m = 3.14 + 2.718 i

2.718 ENTER^
3.14 XEQ "ENB" >>>> 1.008 ---Execution time = 56s---

-And we get in R01 thru R08:

E_m (b) = ( -8.3244596 E10 + 7.1723931 E10 i ) + ( 3.8924375 E10 + 4.5261116 E10 i ) e₁
+ ( 4.1361995 E10 + 4.8008914 E10 i ) e₂ + ( 4.3799615 E10 + 5.0756712 E10 ) e₃

6°) Bessel Functions

- "BSLB" 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 F01 to get the modified Bessel functions .

Formulae:

Bessel Functions of the 1st kind

J_m(b) = (b/2)^m [ 1/Gam(m+1) + (-b²/4)¹/ (1! Gam(m+2) ) + .... + (-b²/4)^k/ (k! Gam(m+k+1) ) + .... ] m # -1 ; -2 ; -3 ; ....

Modified Bessel Functions of the 1st kind

I_m(b) = (b/2)^m [ 1/Gam(m+1) + (b²/4)¹/ (1! Gam(m+2) ) + .... + (b²/4)^k/ (k! Gam(m+k+1) ) + .... ] m # -1 ; -2 ; -3 ; ....

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

Y_m(b) = ( J_m(b) cos(m(pi)) - J_-m(b) ) / sin(m(pi)) ; K_m(b) = (pi/2) ( I_-m(b) - I_m(b) ) / sin(m(pi)) m # .... -3 ; -2 ; -1 ; 0 ; 1 ; 2 ; 3 ....

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

Y_m(b) = -(1/pi) (b/2)^-m SUM_{k=0,1,....,m-1} (m-k-1)!/(k!) (b²/4)^k + (2/pi) Ln(b/2) J_m(b)
- (1/pi) (b/2)^m SUM_k=0,1,..... ( psi(k+1) + psi(m+k+1) ) (-b²/4)^k / (k!(m+k)!) where psi = the digamma function

K_m(b) = (1/2) (b/2)^-m SUM_k=0,1,..,m-1 (m-k-1)!/(k!) (-b²/4)^k - (-1)m Ln(b/2) I_m(b)
+ (1/2) (-1)^m (b/2)^m SUM_k=0,1,...( psi(k+1) + psi(m+k+1) ) (b²/4)^k / (k!(m+k)!)

>>> Here, m may be a complex number p + i.q.

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "BSLB" )

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

R2n+1 .......... R6n: temp ( R6n+1 to R8n are also used for the functions of the second kind )

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

Flags: F01-F02

Subroutines: B*B1   "LNB" "E^B" "B^Z" "B^X" Z*B ST*A   ST/A   DSA DS*A A-A 1/Z Z/Z   AMOVE ASWAP X+1 X-1 X/2 X/E3
                        STO U RCL U STO V RCL V STO W RCL W SINH COSH
                        "GAMZ"    ( a version that does not use any data register - cf "Gamma Function for the HP-41" paragraph 1°) h-3) )

-Lines 48-127 are three-byte GTO 09
-Line 164 ( GEU ) may be replaced with .5772156649 ( Euler's constant )

01 LBL "BSLB"
02 2
03 ST/ A
04 RDN
05 FS? 02
06 GTO 02
07 LBL 10
08 STO U
09 X<>Y
10 STO V
11 12
12 AMOVE
13 B*B1
14 12
15 ASWAP
16 RCL V
17 RCL U
18 XEQ "B^Z"
19 RCL V
20 RCL U
21 X+1
22 XEQ "GAMZ"
23 1/Z
24 Z*B
25 13
26 AMOVE
27 CLX
28 STO W
29 LBL 01
30 B*B1
31 RCL V
32 RCL U
33 RCL W
34 X+1
35 STO W

36 ST+ Y
37 FC? 01
38 CHS
39 ST* Z
40 *
41 1/Z
42 Z*B
43 DSA
44 X#0?
45 GTO 01
46 13
47 ASWAP
48 GTO 09
49 LBL 02
50 14
51 AMOVE
52 RDN
53 XEQ 10
54 RCL U
55 RCL V
56 X#0?
57 GTO 03
58 X<>Y
59 FRC
60 X=0?
61 GTO 04
62 X<> L
63 X<>Y
64 LBL 03
65 FS? 01
66 GTO 00
67 X<>Y
68 PI
69 ST* Z
70 *

71 RAD
72 1
73 P-R
74 X<>Y
75 RCL Z
76 SINH
77 *
78 CHS
79 X<>Y
80 R^
81 COSH
82 *
83 Z*B
84 LBL 00
85 14
86 ASWAP
87 RCL V
88 CHS
89 RCL U
90 CHS
91 XEQ 10
92 RCL V
93 CHS
94 STO V
95 RCL U
96 CHS
97 STO U
98 42
99 AMOVE
100 A-A
101 R^
102 R^
103 PI
104 ST* Z
105 *

106 RAD
107 1
108 P-R
109 RCL Z
110 SINH
111 *
112 X<>Y
113 R^
114 COSH
115 *
116 1/Z
117 PI
118 X/2
119 CHS
120 FC? 01
121 ST/ X
122 CHS
123 ST* Z
124 *
125 Z*B
126 DEG
127 GTO 09
128 LBL 04
129 12
130 AMOVE
131 41
132 AMOVE
133 XEQ "LNB"
134 SIGN
135 FS? 01
136 CHS
137 RCL U
138 Y^X
139 ST+ X
140 ST*A

141 B*B1
142 13
143 AMOVE
144 41
145 AMOVE
146 42
147 AMOVE
148 B*B1
149 12
150 ASWAP
151 RCL U
152 XEQ "B^X"
153 RCL U
154 FACT
155 1
156 FS? 01
157 CHS
158 LASTX
159 Y^X
160 *
161 ST/A
162 CLX
163 STO W
164 GEU
165 ST+ X
166 STO M
167 RCL U
168 ABS
169 LBL 05
170 X#0?
171 1/X
172 ST- M
173 X<> L
174 DSE X
175 GTO 05

176 RCL M
177 GTO 00
178 LBL 06
179 B*B1
180 RCL W
181 X+1
182 ENTER^
183 STO W
184 RCL U
185 +
186 *
187 FC? 01
188 CHS
189 ST/A
190 RCL V
191 LASTX
192 1/X
193 -
194 RCL W
195 1/X
196 -
197 LBL 00
198 STO V
199 DS*A
200 X#0?
201 GTO 06
202 41
203 AMOVE
204 RCL U
205 X=0?
206 GTO 00
207 CHS
208 XEQ "B^X"
209 RCL U
210 X-1

211 FACT
212 CHS
213 ST*A
214 CLX
215 STO W
216 DSA
217 LBL 08
218 RCL U
219 RCL W
220 X+1
221 STO W
222 X=Y?
223 GTO 00
224 STO Z
225 -
226 *
227 FS? 01
228 CHS
229 ST/A
230 B*B1
231 DSA
232 GTO 08
233 LBL 00
234 31
235 AMOVE
236 PI
237 FS? 01
238 -2
239 ST/A
240 LBL 09
241 RCL 00
242 X/E3
243 ISG X
244 CLA
245 END

( 427 bytes / SIZE 6n+1 or 8n+1 )

STACK	INPUT	OUTPUT
Y	q	/
X	p	1.2n

where 1.2n is the control number of B_p+i.q (b)

B = J if CF 01 CF 02 B = Y if CF 01 SF 02
B = I if SF 01 CF 02 B = K if SF 01 SF 02

Example: b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( 0.5 + 0.6 i ) e₂ + ( 0.7 + 0.8 i ) e₃ ( biquaternion -> 8 STO 00 )

0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08

• With m = 2 + 3.i CF 01 CF 02 JNB

3 ENTER^
2 XEQ "BSLB" >>>> 1.008 ---Execution time = 56s---

-And we get in R01 thru R08:

J_2+3.i (b) = ( 1.761978801 + 2.178620406 i ) + ( -0.807601648 + 0.577907527 i ) e₁
+ ( -1.213625969 + 0.966143875 i ) e₂ + ( -1.619650291 + 1.354380222 ) e₃

• With m = 2 + 3.i SF 01 SF 02 KNB

3 ENTER^
2 XEQ "BSLB" >>>> 1.008 ---Execution time = 121s---

-And we get in R01 thru R08:

K_2+3.i (b) = ( 14.38068939 - 63.26852081 i ) + ( -22.09224564 - 6.403916002 i ) e₁
+ ( -34.97621648 - 8.222729323 i ) e₂ + ( -47.86018735 - 10.04154263 ) e₃

• With m = 3 CF 01 SF 02 YNB

0 ENTER^
3 XEQ "BSLB" >>>> 1.008 ---Execution time = 138s---

-And in R01 thru R08:

Y₃ (b) = ( -0.712755221 - 0.475908557 i ) + ( 0.575617774 + 0.146363744 i ) e₁
+ ( 0.909672826 + 0.182278018 i ) e₂ + ( 1.243727879 + 0.218192293 ) e₃

-Likewise, SF 01 CF 02 to calculate the modified Bessel functions of the 1st kind INB

Notes:

-"BSLB" does not work if m is a negative integer, but we can employ the relations:

J_m = (-1)^m J_-m I_m = I_-m
Y_m = (-1)^m Y_-m K_m = K_-m

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

7°) Struve Functions

Formulae:

CF 01 H_m(b) = (b/2)^m+1 ₁F₂( 1 ; 3/2 , m + 3/2 ; - b²/4 ) / Gam( m+3/2 ) with m+3/2 # 0 , -1 , -2 , ...................
SF 01 L_m(b) = (b/2)^m+1 ₁F₂( 1 ; 3/2 , m + 3/2 ; b²/4 ) / Gam(m+3/2)

>>> m may be a complex number.

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "STRB" )

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

R2n+1 .......... R8n: temp

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

Flags: F01-F04
Subroutines: B*B1   "B^Z" Z*B ST*A   ST/A   DSA   1/Z   AMOVE   X+1
                        STO U RCL U STO V RCL V STO W RCL W
                        "GAMZ"    ( a version that does not use any data register - cf "Gamma Function for the HP-41" paragraph 1°) h-3) )

-Line 56 is Gam(3/2)

01 LBL "STRB"
02 1.5
03 +
04 STO U
05 X<>Y
06 STO V
07 2
08 ST/A
09 12
10 AMOVE
11 14
12 AMOVE
13 SIGN

14 CHS
15 FC? 01
16 ST*A
17 B*B1
18 12
19 AMOVE
20 CLX
21 STO W
22 ST*A
23 SIGN
24 STO 01
25 13
26 AMOVE

27 LBL 01
28 B*B1
29 RCL V
30 RCL U
31 1.5
32 RCL W
33 ST+ Y
34 ST+ Z
35 X+1
36 STO W
37 RDN
38 ST* Z
39 *

40 1/Z
41 Z*B
42 DSA
43 X#0?
44 GTO 01
45 41
46 AMOVE
47 RCL V
48 RCL U
49 2
50 1/X
51 -
52 XEQ "B^Z"

53 RCL V
54 RCL U
55 XEQ "GAMZ"
56 .8862269255
57 ST* Z
58 *
59 1/Z
60 Z*B
61 32
62 AMOVE
63 B*B1
64 END

( 137 bytes / SIZE 8n+1 )

STACK	INPUT	OUTPUT
Y	q
X	p	1.2n

where 1.2n is the control number of H_p+i.q (b) if CF 01 or L_p+i.q (b) if SF 01

Example: b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( 0.5 + 0.6 i ) e₂ + ( 0.7 + 0.8 i ) e₃ ( biquaternion -> 8 STO 00 )

0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08

• With m = 2 + 3.i CF 01

3 ENTER^
2 XEQ "STRB" >>>> 1.008 ---Execution time = 61s---

And in R01 thru R08:

H_2+3.i (b) = ( 1.051041908 - 0.109956748 i ) + ( 0.017273466 + 0.374543488 i ) e₁
+ ( 0.056910085 + 0.582905964 i ) e₂ + ( 0.096546705 + 0.791268440 ) e₃

• With m = 2 + 3.i SF 01

3 ENTER^
2 XEQ "STRB" >>>> 1.008 ---Execution time = 61s---

And in R01 thru R08:

L_2+3.i (b) = ( 0.996410287 - 0.241120848 i ) + ( 0.064842486 + 0.357891440 i ) e₁
+ ( 0.129785593 + 0.553123248 i ) e₂ + ( 0.194728701 + 0.748355056 ) e₃

Note:

-"STRB" does not work if m = -3/2 , -5/2 , ....

8°) Jacobian Elliptic Functions

"JEFB" employs Gauss' transformation to calculate sn ( | m ) , cn ( b | m ) & dn ( b | m )

-If m # 1 , let m' = 1-m , µ = [ ( 1-sqrt(m') / ( 1+sqrt(m') ]² and v = b / ( 1+sqrt(µ) ] , then:

-These formulas are applied recursively until µ is small enough to use

   sn ( v | 0 ) = Sin v
   cn ( v | 0 ) = Cos v
   dn ( v | 0 ) = 1

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

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "JEFB" )

• R01 ...... • R2n = the n components of the bion

R2n+1 .......... R8n: temp R8n+1 , R8n+2 ..... temp

>>> When the program stops:     R01   ...... R2n = the n components of sn ( b | m )
                                                      R2n+1 .... R4n = -------------------- cn ( b | m )
                                                      R4n+1 .... R6n = -------------------- dn ( b | m )

Flag: F01
Subroutines: ST*A ST/A B*B1 Z*B "1/B" "CHB" "THB" "COSB" "SINB" AMOVE ASWAP Z*Z Z/Z Z^2 SQRTZ X+1 X-1 X/2
STO W RCL W STO U RCL U STO V RCL V

-Line 165 is a three-byte GTO 03

01 LBL "JEFB"
02 STO U
03 X<>Y
04 STO V
05 X#0?
06 GTO 01
07 SIGN
08 X#Y?
09 GTO 01
10 13
11 AMOVE
12 XEQ "CHB"
13 XEQ "1/B"
14 13
15 ASWAP
16 XEQ "THB"
17 32
18 AMOVE
19 X<>Y
20 RTN
21 LBL 01
22 RCL 00
23 4
24 *
25 .1

26 %
27 +
28 CLA
29 STO O
30 STO P
31 SIGN
32 STO M
33 RCL V
34 RCL U
35 ENTER^
36 LBL 02
37 RDN
38 RCL Y
39 CHS
40 RCL Y
41 CHS
42 X+1
43 SQRTZ
44 X+1
45 X<>Y
46 STO V
47 X<>Y
48 STO U
49 Z^2
50 Z/Z

51 ISG O
52 CLX
53 STO IND O
54 X<>Y
55 ISG O
56 CLX
57 STO IND O
58 X<>Y
59 Z^2
60 X<> M
61 X<>Y
62 X<> N
63 X<>Y
64 RCL V
65 X/2
66 RCL U
67 X/2
68 Z*Z
69 X<> M
70 X<>Y
71 X<> N
72 X<>Y
73 RCL Y
74 ABS
75 RCL Y

76 ABS
77 +
78 X#0?
79 GTO 02
80 RCL P
81 X+1
82 RCL 00
83 -
84 CLRGX
85 SIGN
86 STO IND L
87 RCL O
88 STO W
89 RCL N
90 RCL M
91 Z*B
92 13
93 AMOVE
94 XEQ "COSB"
95 13
96 ASWAP
97 XEQ "SINB"
98 12
99 AMOVE
100 LBL 03

101 41
102 AMOVE
103 23
104 ASWAP
105 B*B1
106 13
107 ASWAP
108 12
109 AMOVE
110 B*B1
111 RCL W
112 RCL IND X
113 DSE Y
114 RCL IND Y
115 Z*B
116 23
117 ASWAP
118 14
119 AMOVE
120 SIGN
121 ST+ 01
122 XEQ "1/B"
123 B*B1
124 13
125 ASWAP

126 RCL W
127 ENTER^
128 X-1
129 STO W
130 RCL IND Y
131 RCL IND Y
132 X+1
133 Z*B
134 12
135 AMOVE
136 41
137 AMOVE
138 SIGN
139 ST+ 01
140 XEQ "1/B"
141 B*B1
142 14
143 ASWAP
144 12
145 AMOVE
146 SIGN
147 CHS
148 ST- 01
149 ST*A
150 12

151 ASWAP
152 SIGN
153 ST+ 01
154 XEQ "1/B"
155 B*B1
156 14
157 ASWAP
158 12
159 AMOVE
160 RCL W
161 ENTER^
162 X-1
163 STO W
164 DSE Y
165 GTO 03
166 21
167 AMOVE
168 32
169 AMOVE
170 43
171 AMOVE
172 RCL 00
173 X/E3
174 ISG X
175 END

( 339 bytes / ZIZE 8n+1+ ??? )

STACK	INPUT	OUTPUT
Y	q
X	p	1.2n

where 1.2n is the control number of sn ( b | p+q.i )

2n+1.4n ------------------------ cn ( b | p+q.i )
4n+1.6n ------------------------ dn ( b | p+q.i )

Example: b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( 0.5 + 0.6 i ) e₂ + ( 0.7 + 0.8 i ) e₃ ( biquaternion -> 8 STO 00 )

0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08

• With m = 2 + 3.i

3 ENTER^
2 XEQ "JEFB" >>>> 1.008 ---Execution time = 4m28s---

In R01 thru R08:

sn ( b | 2+3i ) = ( -0.119297321 - 0.145921679 i ) + ( 0.001430948 + 0.135563787 i ) e₁
+ ( 0.013077382 + 0.211365032 i ) e₂ + ( 0.024723817 + 0.287166277 i ) e₃

In R09 thru R16:

cn ( b | 2+3i ) = ( 0.928961221 - 0.011137456 i ) + ( -0.021319014 + 0.017378296 i ) e₁
+ ( -0.031867398 + 0.028815663 i ) e₂ + ( -0.042415782 + 0.040253031 i ) e₃

And in R17 thru R24:

dn ( b | 2+3i ) = ( -0.926505041 + 0.207185359 i ) + ( 0.084840287 + 0.047110613 i ) e₁
+ ( 0.136119697 + 0.066705334 i ) e₂ + ( 0.187399107 + 0.086300054 i ) e₃

[ dn ( b | 2+3i ) is also stored in R25 thru R32 ]

9°) Weber & Anger Functions

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

Formulae:

E_m(b) = - (b/2) cos ( PI.m/2 ) ₁F^~₂( 1 ; (3-m)/2 , (3+m)/2 ; -b²/4 ) Weber's functions

+ sin ( PI.m/2 ) ₁F^~₂( 1 ; (2-m)/2 , (2+m)/2 ; -b²/4 )

J_m(b) = + (b/2) sin ( PI.m/2 ) ₁F^~₂( 1 ; (3-m)/2 , (3+m)/2 ; -b²/4 ) Anger's functions

+ cos ( PI.m/2 ) ₁F^~₂( 1 ; (2-m)/2 , (2+m)/2 ; -b²/4 )

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "WEBB" )

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

R2n+1 .......... R8n: temp

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

Flags: F01-F02-F03 CF 01 for Weber's function
SF 01 for Anger's function

Subroutines: B*B1   "LNB" "E^B" "B^Z" "B^X" Z*B ST*A   ST/A   DSA DS*A A+A 1/Z Z*Z   AMOVE ASWAP X+1 X/2
                        STO U RCL U STO V RCL V STO W RCL W SINH COSH
                        "GAMZ"    ( a version that does not use any data register - cf "Gamma Function for the HP-41" paragraph 1°) h-3) )

01 LBL "WEBB"
02 STO U
03 X<>Y
04 STO V
05 2
06 ST/A
07 12
08 AMOVE
09 14
10 AMOVE
11 SIGN
12 CHS
13 ST*A
14 B*B1
15 12
16 AMOVE
17 SF 02
18 XEQ 01
19 14
20 AMOVE
21 LBL 01
22 CLX
23 STO W
24 ST*A
25 SIGN
26 STO 01
27 13

28 AMOVE
29 LBL 02
30 B*B1
31 RCL U
32 ENTER^
33 CHS
34 2
35 FS? 02
36 X+1
37 ST+ Z
38 +
39 2
40 ST/ Z
41 /
42 RCL W
43 ST+ Y
44 ST+ Z
45 X+1
46 STO W
47 RDN
48 RCL V
49 X/2
50 STO T
51 CHS
52 X<>Y
53 Z*Z
54 1/Z

55 Z*B
56 DSA
57 X#0?
58 GTO 02
59 31
60 AMOVE
61 FC? 02
62 GTO 01
63 24
64 SWAP
65 SIGN
66 CHS
67 FC? 01
68 ST*A
69 B*B1
70 24
71 SWAP
72 LBL 01
73 RCL V
74 RCL U
75 PI
76 X/2
77 ST* Z
78 *
79 RAD
80 1
81 P-R

82 X<> Z
83 SINH
84 LASTX
85 COSH
86 CF 03
87 FS? 01
88 FS? 02
89 SF 03
90 FS? 02
91 FS? 01
92 FS? 30
93 CF 03
94 FS? 03
95 X<>Y
96 ST* T
97 RDN
98 *
99 FC? 03
100 CHS
101 FC?C 03
102 X<>Y
103 Z*B
104 RCL V
105 CHS
106 RCL U
107 CHS
108 2

109 FS? 02
110 X+1
111 +
112 2
113 ST/ Z
114 /
115 XEQ "GAMZ"
116 1/Z
117 Z*B
118 RCL V
119 RCL U
120 2
121 FS? 02
122 X+1
123 +
124 2
125 ST/ Z
126 /
127 XEQ "GAMZ"
128 1/Z
129 Z*B
130 FS?C 02
131 RTN
132 42
133 AMOVE
134 A+A
135 END

( 242 bytes / SIZE 8n+1 )

STACK	INPUT	OUTPUT
Y	q	/
X	p	1.2n

where 1.2n is the control number of E_p+i.q (b) if CF 01 ( Weber ) or J_p+i.q (b) if SF 01 ( Anger )

Example: b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( 0.5 + 0.6 i ) e₂ + ( 0.7 + 0.8 i ) e₃ ( biquaternion -> 8 STO 00 )

0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08

• With m = 0.4 + 0.7 i CF 01

0.7 ENTER^
0.4 XEQ "WEBB" >>>> 1.008 ---Execution time = 125s---

And we get in registers R01 thru R08:

E_m (b) = ( 0.638513905 + 1.155323367 i ) + ( -0.381725526 - 0.246679413 i ) e₁
+ ( -0.615226174 - 0.354281842 i ) e₂ + ( -0.848726822 - 0.461884271 i ) e₃

• With m = 0.4 + 0.7 i SF 01

0.7 ENTER^
0.4 XEQ "WEBB" >>>> 1.008 ---Execution time = 125s---

And we get in registers R01 thru R08:

J_m (b) = ( 1.425286898 - 0.337000377 i ) + ( -0.117808154 + 0.381684479 i ) e₁
+ ( -0.153245961 + 0.604852439 i ) e₂ + ( -0.188683769 + 0.828020400 i ) e₃

Note:

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

10°) Incomplete Gamma Function

Formula: g(m,b) = ( b^m / m ) exp(-b) M(1,m+1;b) where M = Kummer's function

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "IGFB" )

• R01 ........ • R2n = the n components of the bion

R2n+1 .......... R6n: temp

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

Flags: /
Subroutines: B*B1 "E^B" "B^Z" Z*B ST*A DSA 1/Z AMOVE ASWAP X+1
STO U RCL U STO V RCL V STO W RCL W

01 LBL "IGFB"
02 STO U
03 X<>Y
04 STO V
05 12
06 AMOVE
07 CLX
08 STO W

09 ST*A
10 SIGN
11 STO 01
12 13
13 AMOVE
14 LBL 01
15 B*B1
16 RCL V

17 RCL U
18 RCL W
19 X+1
20 STO W
21 +
22 1/Z
23 Z*B
24 DSA

25 X#0?
26 GTO 01
27 21
28 AMOVE
29 RCL V
30 RCL U
31 XEQ "B^Z"
32 12

33 ASWAP
34 SIGN
35 CHS
36 ST*A
37 XEQ "E^B"
38 B*B1
39 RCL V
40 RCL U

41 1/Z
42 Z*B
43 32
44 MOVE
45 B*B1
46 END

( 96 bytes / SIZE 6n+1 )

STACK	INPUT	OUTPUT
Y	q	/
X	p	1.2n

where 1.2n is the control number of the result

Example: b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( 0.5 + 0.6 i ) e₂ + ( 0.7 + 0.8 i ) e₃ ( biquaternion -> 8 STO 00 )

0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08

• With m = 2 + 3 i CF 01

3 ENTER^
2 XEQ "IGFB" >>>> 1.008 ---Execution time = 86s---

And we get in R01 thru R08:

g(m,b) = ( 0.316021540 - 0.294439242 i ) + ( 0.097233572 + 0.118344541 i ) e₁
+ ( 0.161151935 + 0.176838799 i ) e₂ + ( 0.225070299 + 0.235333057 i ) e₃

11°) Airy Functions

Formulae:

Ai(b) = f(b) - g(b) with f(b) = [ 3^-2/3 / Gamma(2/3) ] ₀F₁( 2/3 ; b³/9 )
Bi(b) = [ f(b) + g(b) ] sqrt(3) and g(b) = [ 3^-1/3 / Gamma(1/3) ] ₀F₁( 4/3 ; b³/9 ) b

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "AIRYB" )

• R01 ...... • R2n = the n components of the bion b

Rn+1 .......... R8n: temp

>>> When the program stops: R01 ...... R2n = the n components of Ai(b)
and R2n+1..... R4n = -------------------- Bi(b)

Flag: F01 is cleared
Subroutines: AMOVE ASWAP B*B1 ST*A ST/A DSA A+A A-A X+1
STO U RCL U STO W RCL W

01 LBL "AIRYB"
02 12
03 AMOVE
04 14
05 AMOVE
06 4
07 3
08 /
09 STO U
10 B*B1
11 B*B1
12 9

13 ST/A
14 12
15 AMOVE
16 SF 01
17 XEQ 01
18 24
19 ASWAP
20 B*B1
21 24
22 ASWAP
23 .2588194038
24 ST*A

25 14
26 AMOVE
27 2
28 3
29 /
30 STO U
31 LBL 01
32 CLX
33 STO W
34 ST*A
35 SIGN
36 STO 01

37 13
38 AMOVE
39 LBL 02
40 B*B1
41 RCL U
42 RCL W
43 ST+ Y
44 X+1
45 STO W
46 *
47 ST/A
48 DSA

49 X#0?
50 GTO 02
51 31
52 AMOVE
53 FS?C 01
54 RTN
55 .3550280539
56 ST*A
57 13
58 AMOVE
59 42
60 AMOVE

61 A+A
62 3
63 SQRT
64 ST*A
65 13
66 ASWAP
67 A-A
68 32
69 AMOVE
70 X<>Y
71 END

( 153 bytes / SIZE 8n+1 )

STACK	INPUT	OUTPUT
X	/	1.2n

where 1.2n is the control number of the result

Example: b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( 0.5 + 0.6 i ) e₂ + ( 0.7 + 0.8 i ) e₃ ( biquaternion -> 8 STO 00 )

0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08

XEQ "AIRYB" >>>> 1.008 ---Execution time = 67s---

We get in R01 thru R08:

Ai(b) = ( 0.478468884 - 0.045373965 i ) + ( -0.040623036 - 0.145492499 i ) e₁
+ ( -0.075011336 - 0.223718455 i ) e₂ + ( -0.109399637 - 0.301944411 i ) e₃

And in R09 thru R16:

Bi(b) = ( 0.650545601 + 0.036453982 i ) + ( 0.247761441 + 0.146528307 i ) e₁
+ ( 0.398230113 + 0.208763244 i ) e₂ + ( 0.548698784 + 0.270998181 i ) e₃

12°) Parabolic Cylinder Functions

a) Ascending Series

Formula:

D_m(b) = 2^m/2 Pi^1/2 exp(-b²/4) [ 1/Gam((1-m)/2) M( -m/2 , 1/2 , b²/2 ) - 2^1/2 ( b / Gam(-m/2) ) M [ (1-m)/2 , 3/2 , b²/2 ]

where M = Kummer's functions.

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "DNB" )

• R01 ...... • R2n = the n components of the bion b

Rn+1 .......... R8n: temp

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

Flag: F01 is cleared
Subroutines: AMOVE ASWAP B*B1 Z*B 1/Z "E^B" ST*A ST/A DSA A-A X+1
STO U RCL U STO V RCL V STO W RCL W
"GAMZ" ( a version that does not use any data register - cf "Gamma Function for the HP-41" paragraph 1°) h-3) )

01 LBL "DNB"
02 2
03 CHS
04 ST/ Z
05 /
06 STO U
07 X<>Y
08 STO V
09 12
10 AMOVE
11 14
12 AMOVE
13 B*B1
14 2
15 ST/A
16 12
17 AMOVE
18 SF 01
19 XEQ 01
20 24

21 ASWAP
22 B*B1
23 24
24 ASWAP
25 2
26 SQRT
27 ST*A
28 14
29 AMOVE
30 LBL 01
31 CLX
32 STO W
33 ST*A
34 SIGN
35 STO 01
36 13
37 AMOVE
38 LBL 02
39 B*B1
40 RCL V

41 RCL U
42 .5
43 FS? 01
44 ST+ Y
45 FS? 01
46 X+1
47 RCL W
48 ST+ Y
49 ST+ Z
50 X+1
51 STO W
52 *
53 ST/ Z
54 /
55 Z*B
56 DSA
57 X#0?
58 GTO 02
59 31
60 AMOVE

61 RCL V
62 RCL U
63 .5
64 FS? 01
65 CLX
66 +
67 XEQ "GAMZ"
68 1/Z
69 Z*B
70 FS?C 01
71 RTN
72 24
73 ASWAP
74 A-A
75 12
76 AMOVE
77 41
78 AMOVE
79 2
80 CHS

81 ST/A
82 XEQ "E^B"
83 RCL V
84 CHS
85 2
86 LN
87 *
88 RAD
89 2
90 RCL U
91 CHS
92 Y^X
93 PI
94 SQRT
95 *
96 P-R
97 Z*B
98 B*B1
99 DEG
100 END

( 181 bytes / SIZE 8n+1 )

STACK	INPUT	OUTPUT
Y	q	/
X	p	1.2n

where 1.2n is the control number of the result , m = p + q i

Example: b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( 0.5 + 0.6 i ) e₂ + ( 0.7 + 0.8 i ) e₃ ( biquaternion -> 8 STO 00 )

0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08

• With m = 0.4 + 0.7 i

0.7 ENTER^
0.4 XEQ "DNB" >>>> 1.008 ---Execution time = 149s---

And we get in registers R01 thru R08:

D_m (b) = ( 0.409729789 + 0.317378267 i ) + ( -0.162419934 + 0.271785448 i ) e₁
+ ( -0.231632221 + 0.436979673 i ) e₂ + ( -0.300844508 + 0.602173399 i ) e₃

b) Asymptotic Expansion

Formula: D_m(b) ~ b^m exp(-b²/4) [ 1 - m(m-1) / ( 2 b² ) + m(m-1)(m-2)(m-3) / ( 2 ( 4 b⁴ ) ) - ....... ]

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "AEDNB" )

• R01 ...... • R2n = the n components of the bion b

Rn+1 .......... R6n: temp

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

Flags: /
Subroutines: AMOVE ASWAP B*B1 Z*B Z*Z "E^B" "B^Z" ST/A DSA X-1
STO U RCL U STO V RCL V STO W RCL W

01 LBL "AEDNB"
02 STO U
03 X<>Y
04 STO V
05 12
06 AMOVE
07 X<> Z
08 XEQ "B^Z"
09 13
10 AMOVE
11 21
12 AMOVE

13 B*B1
14 12
15 AMOVE
16 4
17 CHS
18 ST/A
19 XEQ "E^B"
20 23
21 ASWAP
22 B*B1
23 13
24 ASWAP

25 XEQ "1/B"
26 12
27 AMOVE
28 CLX
29 STO W
30 31
31 AMOVE
32 LBL 01
33 B*B1
34 RCL V
35 RCL U
36 RCL W

37 -
38 ENTER^
39 X-1
40 RCL V
41 X<>Y
42 Z*Z
43 RCL W
44 2
45 +
46 STO W
47 CHS
48 ST/ Z

49 /
50 Z*B
51 DSA
52 X#0?
53 GTO 01
54 31
55 AMOVE
56 RCL 00
57 E3
58 /
59 ISG X
60 END

( 122 bytes / SIZE 6n+1 )

STACK	INPUT	OUTPUT
Y	q	/
X	p	1.2n

where 1.2n is the control number of the result , m = p + q i

Example: b = ( 5 + 6 i ) + ( 7 + 8 i ) e₁ + ( 9 + 10 i ) e₂ + ( 11 + 12 i ) e₃ ( biquaternion -> 8 STO 00 )

5 STO 01 6 STO 02 7 STO 03 8 STO 04 9 STO 05 10 STO 06 11 STO 07 12 STO 08

• With m = 3.14 + 2.718 i

2.718 ENTER^
3.14 XEQ "DNB" >>>> 1.008 ---Execution time = 60s---

And we get in registers R01 thru R08:

D_m (b) = ( -4.3259184 E34 + 2.1466036 E36 i ) + ( 9.6478579 E35 + 3.4397076 E34 i ) e₁
+ ( 1.2215323 E36 + 2.6453388 E34 i ) e₂ + ( 1.4782789 E36 + 1.8509600 E34 i ) e₃

Note:

-The asymptotic series diverge too soon if b is "small" ... unless m is a positive real integer !

13°) Hermite Functions

Formula:

H_m(b) = 2^m sqrt(PI) [ (1/Gam((1-m)/2)) M(-m/2,1/2,b²) - ( 2.b / Gam(-m/2) ) M((1-m)/2,3/2,b²) ]

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

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "HMTB" )

• R01 ...... • R2n = the n components of the bion

R2n+1 .......... R8n: temp

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

Flag: F01
Subroutines: AMOVE ASWAP B*B1 Z*B 1/Z ST*A ST/A DSA A-A X+1 X/E3
STO U RCL U STO V RCL V STO W RCL W
"GAMZ" ( a version that does not use any data register - cf "Gamma Function for the HP-41" paragraph 1°) h-3) )

01 LBL "HMTB"
02 2
03 CHS
04 ST/ Z
05 /
06 STO U
07 X<>Y
08 STO V
09 12
10 AMOVE
11 14
12 AMOVE
13 B*B1
14 12
15 AMOVE
16 SF 01
17 XEQ 01
18 24
19 ASWAP

20 B*B1
21 24
22 ASWAP
23 2
24 ST*A
25 14
26 AMOVE
27 LBL 01
28 CLX
29 STO W
30 ST*A
31 SIGN
32 STO 01
33 13
34 AMOVE
35 LBL 02
36 B*B1
37 RCL V
38 RCL U

39 .5
40 FS? 01
41 ST+ Y
42 FS? 01
43 X+1
44 RCL W
45 ST+ Y
46 ST+ Z
47 X+1
48 STO W
49 *
50 ST/ Z
51 /
52 Z*B
53 DSA
54 X#0?
55 GTO 02
56 31
57 AMOVE

58 RCL V
59 RCL U
60 .5
61 FS? 01
62 CLX
63 +
64 XEQ "GAMZ"
65 1/Z
66 Z*B
67 FS?C 01
68 RTN
69 24
70 ASWAP
71 A-A
72 RCL V
73 CHS
74 ST+ X
75 2
76 LN

77 *
78 RAD
79 4
80 RCL U
81 CHS
82 Y^X
83 PI
84 SQRT
85 *
86 P-R
87 Z*B
88 RCL 00
89 X/E3
90 ISG X
91 DEG
92 END

( 96 bytes / SIZE 8n+1 )

STACK	INPUT	OUTPUT
Y	q	/
X	p	1.2n

where 1.2n is the control number of the result , m = p + q i

Example: b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( 0.5 + 0.6 i ) e₂ + ( 0.7 + 0.8 i ) e₃ ( biquaternion -> 8 STO 00 )

0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08

• With m = 0.4 + 0.7 i

0.7 ENTER^
0.4 XEQ "HMTB" >>>> 1.008 ---Execution time = 3mn00s---

And we get in registers R01 thru R08:

H_m (b) = ( 0.817434194 + 0.870134917 i ) + ( -0.140724908 + 0.379062438 i ) e₁
+ ( -0.189205861 + 0.602595396 i ) e₂ + ( -0.237686815 + 0.826128354 i ) e₃

14°) Harmonic Numbers

"HRMB" calculates Hr_N(b) = SUM_k=1,...,N k^(-b)

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "HRMB" )

• R01 ...... • R2n = the n components of the bion

R2n+1 .......... R6n: temp

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

Flags: /
Subroutines: ST*A DSA AMOVE "X^B" X-1 STO W RCL W

01 LBL "HRMB"
02 STO W
03 SIGN
04 CHS
05 ST*A
06 12

07 AMOVE
08 CLX
09 ST*A
10 13
11 AMOVE
12 LBL 01

13 21
14 AMOVE
15 RCL W
16 XEQ "X^B"
17 DSA
18 RCL W

19 ENTER^
20 X-1
21 STO W
22 DSE Y
23 GTO 01
24 31

25 AMOVE
26 RCL 00
27 E3
28 /
29 ISG X
30 END

( 63 bytes / SIZE 6n+1 )

STACK	INPUT	OUTPUT
X	N	1.2n

where 1.2n is the control number of the result , N = a positive integer

Example: b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( 0.5 + 0.6 i ) e₂ + ( 0.7 + 0.8 i ) e₃ ( biquaternion -> 8 STO 00 )

0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08

• With N = 12

12 XEQ "HRMB" >>>> 1.008 ---Execution time = 127s---

And in registers R01 thru R08:

Hr_N (b) = ( -20.51713285 - 23.81908134 i ) + ( -9.341156837 + 7.380203802 i ) e₁
+ ( -13.98178837 + 12.26041048 i ) e₂ + ( -18.62241988 + 17.14061716 i ) e₃

15°) Generalized Error Functions

Formula:

Erf_m(b) = b exp(-b^m) M( 1 ; 1+1/m ; b^m ) where M = Kummer's function

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "GERFB" )

• R01 ...... • R2n = the n components of the bion

R2n+1 .......... R6n: temp

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

Flag: /
Subroutines: AMOVE ASWAP B*B1 Z*B "B^Z" "E^B" 1/Z ST*A DSA X+1
STO U RCL U STO V RCL V STO W RCL W

01 LBL "GERFB"
02 STO U
03 X<>Y
04 STO V
05 13
06 AMOVE
07 X<> Z
08 XEQ "B^Z"
09 RCL V
10 RCL U

11 1/Z
12 STO U
13 X<>Y
14 STO V
15 12
16 AMOVE
17 SIGN
18 CHS
19 ST*A
20 XEQ "E^B"

21 23
22 ASWAP
23 B*B1
24 32
25 AMOVE
26 13
27 AMOVE
28 CLX
29 STO W
30 LBL 01

31 B*B1
32 RCL V
33 RCL U
34 RCL W
35 X+1
36 ST+ Y
37 STO W
38 RDN
39 1/Z
40 Z*B

41 DSA
42 X#0?
43 GTO 01
44 31
45 AMOVE
46 RCL 00
47 E3
48 /
49 ISG X
50 END

( 104 bytes SIZE 6n+1 )

STACK	INPUT	OUTPUT
Y	q	/
X	p	1.2n

where 1.2n is the control number of the result , m = p + q i

Example: b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( 0.5 + 0.6 i ) e₂ + ( 0.7 + 0.8 i ) e₃ ( biquaternion -> 8 STO 00 )

0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08

• With m = 0.4 + 0.7 i

0.7 ENTER^
0.4 XEQ "GERFB" >>>> 1.008 ---Execution time = 101s---

And we get in registers R01 thru R08:

Erf_m (b) = ( -0.073280236 + 0.530416043 i ) + ( 0.176820829 + 0.257161065 i ) e₁
+ ( 0.296413379 + 0.387025596 i ) e₂ + ( 0.416005929 + 0.516890126 i ) e₃

16°) Lambert-W Function

-We solve b = w(b) exp [ w(b) ]
-"LBWB" employs Newton's method.

Data Registers: • R00 = 2n > 3 ( Registers R00 thru R2n are to be initialized before executing "LBWB" )

• R01 ...... • R2n = the n components of the bion

R2n+1 .......... R8n: temp

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

Flag: /
Subroutines: AMOVE ASWAP B*B1 "1/B" "LNB" "E^B" A-A ST*A DSA

01 LBL "LBWB"
02 14
03 AMOVE
04 SIGN
05 ST+ 01
06 XEQ "LNB"
07 13
08 AMOVE

09 LBL 01
10 31
11 AMOVE
12 VIEW 01
13 SIGN
14 CHS
15 ST*A
16 XEQ "E^B"

17 42
18 AMOVE
19 B*B1
20 32
21 AMOVE
22 A-A
23 12
24 ASWAP

25 SIGN
26 ST+ 01
27 XEQ "1/B"
28 B*B1
29 DSA
30 E-8
31 X<Y?
32 GTO 01

33 31
34 AMOVE
35 RCL 00
36 E3
37 /
38 ISG X
39 CLD
40 END

( 88 bytes / SIZE 8n+1 )

STACK	INPUT	OUTPUT
X	/	1.2n

where 1.2n is the control number of the result

Example: b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( 0.5 + 0.6 i ) e₂ + ( 0.7 + 0.8 i ) e₃ ( biquaternion -> 8 STO 00 )

0.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 0.5 STO 05 0.6 STO 06 0.7 STO 07 0.8 STO 08

XEQ "LBWB" >>>> 1.008 ---Execution time = 121s---

And in registers R01 thru R08:

W (b) = ( 0.494380927 + 0.401088979 i ) + ( 0.217079229 + 0.074534310 i ) e₁
+ ( 0.344606343 + 0.098907186 i ) e₂ + ( 0.472133456 + 0.123280061 i ) e₃

Notes:

-This routine does not work for b = -1, we have W(-1) = -0.3181315052 + 1.337235701 i
-Even in the neighborhood of -1 , the 1st approximation ( Ln(1+b) , cf line 06 ) is not a very good one !
-Change lines 04-05-06 if you prefer a better guess.

hp41programs

Bionic Special Functions for the HP-41