Bionic Functions2

Bionic Special Functions(II) for the HP-41

Overview

0°) 5 Extra-Registers & M-Code Routines
1°) Associated Legendre Functions
2°) Chebyshev Functions
3°) Incomplete Beta Function
4°) Laguerre's Functions
5°) Lerch Transcendent Function
6°) Lommel Functions
7°) Toronto Functions
8°) Ultraspherical Functions
9°) Whittaker's M- & W-Functions
10°) Riemann Zeta Function

0°) 5 Extra Registers & M-Code Routines

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

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

-Other M-Code routines are employed:

Z*Z 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°) Associated Legendre Functions

"ALFB" 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 P_r^m(b) = [ (b+1)/(1-b) ]^m/2 ₂F₁(-r , r+1 ; 1-m ; (1-b)/2 ) / Gamma(1-m) ( b # 1 )

Type 3 P_r^m(a) = [ (b+1)/(b-1) ]^m/2 ₂F₁(-r , r+1 ; 1-m ; (1-b)/2 ) / Gamma(1-m) ( b # 1 )

Formulae: 2nd Kind

Type 2 Q_r^m(b) = 2^m pi^1/2 (1-b²)^-m/2 [ -Gam((1+m+r)/2)/(2.Gam((2-m+r)/2)) . sin pi(m+r)/2 . ₂F₁(-r/2-m/2 ; 1/2+r/2-m/2 ; 1/2 ; b² )
+ b Gam((2+r+m)/2) / Gam((1+r-m)/2) . cos pi(m+r)/2 . ₂F₁((1-m-r)/2 ; (2+r-m)/2 ; 3/2 ; b² ) ]

Type 3 Q_r^m(b) = exp( i (m.PI) ) 2^-r-1 sqrt(PI) Gam(m+r+1) b^-m-r-1 (b²-1)^m/2 ₂F^~₁( (2+m+r)/2 , (1+m+r)/2 ; r+3/2 ; 1/b² )

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

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

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

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

Flags: F00-F02-F03-F04

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

-Lines 09-81-86-211 are three-byte GTO's

01 LBL "ALFB"
02 STO U
03 RDN
04 STO V
05 RDN
06 STO R
07 X<>Y
08 STO S
09 FS? 02
10 GTO 00
11 .5
12 CHS
13 ST*A
14 ST- 01
15 12
16 AMOVE
17 CLX
18 STO W
19 ST*A
20 SIGN
21 STO 01
22 13
23 AMOVE
24 LBL 01
25 B*B1
26 RCL V
27 CHS
28 RCL U
29 CHS
30 RCL U
31 X+1
32 RCL W
33 ST+ Z
34 +
35 R^
36 CHS
37 X<>Y
38 Z*Z
39 RCL R
40 CHS
41 RCL W
42 X+1
43 STO W
44 +
45 RCL S
46 CHS
47 X<>Y
48 Z/Z
49 RCL W
50 ST/ Z
51 /
52 Z*B
53 DSA
54 X#0?
55 GTO 01

56 21
57 AMOVE
58 SIGN
59 ST- 01
60 CHS
61 FC? 03
62 ST*A
63 12
64 ASWAP
65 XEQ "1/B"
66 B*B1
67 32
68 AMOVE
69 RCL S
70 RCL R
71 2
72 ST/ Z
73 /
74 XEQ "B^Z"
75 RCL S
76 CHS
77 RCL R
78 CHS
79 X+1
80 XEQ "GAMZ"
81 1/Z
82 GTO 04
83 LBL 00
84 12
85 AMOVE
86 FC? 03
87 GTO 00
88 13
89 AMOVE
90 B*B1
91 12
92 AMOVE
93 XEQ "1/B"
94 12
95 ASWAP
96 SIGN
97 ST- 01
98 RCL S
99 RCL R
100 2
101 ST/ Z
102 /
103 XEQ "B^Z"
104 23
105 ASWAP
106 12
107 ASWAP
108 RCL V
109 CHS
110 RCL S

111 -
112 RCL R
113 CHS
114 RCL U
115 -
116 X-1
117 XEQ "B^Z"
118 B*B1
119 32
120 AMOVE
121 13
122 AMOVE
123 CLX
124 STO W
125 LBL 02
126 B*B1
127 RCL V
128 RCL S
129 +
130 RCL U
131 RCL R
132 +
133 X+1
134 ENTER^
135 X+1
136 2
137 ST/ Z
138 ST/ T
139 /
140 RCL W
141 ST+ Z
142 +
143 Y<>Z
144 Z*Z
145 RCL U
146 .5
147 +
148 RCL W
149 X+1
150 STO W
151 +
152 RCL V
153 X<>Y
154 Z/Z
155 RCL W
156 ST/ Z
157 /
158 Z*B
159 DSA
160 X#0?
161 GTO 02
162 31
163 AMOVE
164 RCL V
165 RCL S

166 +
167 RCL U
168 RCL R
169 +
170 X+1
171 XEQ "GAMZ"
172 Z*B
173 RCL V
174 RCL U
175 1.5
176 +
177 XEQ "GAMZ"
178 1/Z
179 Z*B
180 1
181 CHS
182 ACOS
183 RCL R
184 *
185 RCL S
186 PI
187 *
188 CHS
189 E^X
190 PI
191 SQRT
192 *
193 X/2
194 P-R
195 Z*B
196 RCL V
197 2
198 LN
199 CHS
200 *
201 RCL U
202 LASTX
203 *
204 E^X
205 RAD
206 P-R
207 Z*B
208 RCL 00
209 X/E3
210 ISG X
211 GTO 05
212 LBL 00
213 B*B1
214 12
215 ASWAP
216 SF 04
217 XEQ 00
218 ST*T
219 RDN
220 *

221 CHS
222 X<>Y
223 Z*B
224 14
225 AMOVE
226 CLX
227 ST*A
228 SIGN
229 STO 01
230 LBL 00
231 CLX
232 STO W
233 13
234 AMOVE
235 LBL 03
236 B*B1
237 RCL V
238 RCL S
239 +
240 CHS
241 X/2
242 RCL U
243 RCL R
244 +
245 CHS
246 FS? 04
247 X+1
248 X/2
249 RCL U
250 RCL R
251 -
252 2
253 FC? 04
254 SIGN
255 +
256 X/2
257 RCL W
258 ST+ Z
259 +
260 RCL V
261 SIGN
262 RDN
263 RCL S
264 ST- L
265 X<> L
266 X/2
267 X<>Y
268 Z*Z
269 .5
270 FS? 04
271 X+1
272 RCL W
273 ST+ Y
274 X+1
275 STO W

276 *
277 ST/ Z
278 /
279 Z*B
280 DSA
281 X#0?
282 GTO 03
283 31
284 AMOVE
285 RCL S
286 RCL V
287 +
288 X/2
289 RCL R
290 RCL U
291 +
292 2
293 FC? 04
294 SIGN
295 +
296 X/2
297 XEQ "GAMZ"
298 Z*B
299 RCL V
300 RCL S
301 -
302 X/2
303 RCL U
304 RCL R
305 -
306 2
307 FS? 04
308 SIGN
309 +
310 X/2
311 XEQ "GAMZ"
312 2
313 CHS
314 FS? 04
315 ST/ X
316 ST* Z
317 *
318 1/Z
319 Z*B
320 RCL S
321 RCL V
322 +
323 PI
324 X/2
325 *
326 COSH
327 LASTX
328 SINH
329 RCL R
330 RCL U

331 +
332 PI
333 X/2
334 *
335 RAD
336 1
337 P-R
338 FS?C 04
339 RTN
340 ST* Z
341 X<> T
342 *
343 Z*B
344 24
345 ASWAP
346 A+A
347 12
348 AMOVE
349 41
350 AMOVE
351 SIGN
352 CHS
353 ST*A
354 ST- 01
355 RCL S
356 RCL R
357 2
358 ST/ Z
359 /
360 XEQ "B^Z"
361 XEQ "1/B"
362 RCL S
363 2
364 LN
365 *
366 RAD
367 2
368 RCL R
369 Y^X
370 PI
371 SQRT
372 *
373 P-R
374 LBL 04
375 Z*B
376 B*B1
377 LBL 05
378 DEG
379 END

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

STACK	INPUTS	OUTPUTS
T	m₂	/
Z	m₁	/
Y	r₂	/
X	r₁	1.2n

where m = m₁ + i m₂ , r = r₁ + i r₂
and 1.2n 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 = 1 + 2 i , r = 3 + 4 i , b = ( 1.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( -0.1 - 0.2 i ) e₂ + ( -0.3 - 0.4 i ) e₃ ( biquaternion -> 8 STO 00 )

1.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 -0.1 STO 05 -0.2 STO6 -0.3 STO 07 -0.4 STO 08

       4 ENTER^
       3 ENTER^
       2 ENTER^
       1 XEQ "ALFB"   >>>>   1.008                               ---Execution time = 3m20s---

And we get in R01 thru R08:

>>> P_r^m(b) = ( 93.10375888 - 105.0475343 i ) + ( 75.56897000 + 91.98871986 i ) e₁
+ ( -25.89125050 - 46.52055456 i ) e₂ + ( -75.56897000 - 91.98871986 i ) e₃

m # 1 , 2 , 3 , .......... & b # 1

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

m = 1 + 2 i , r = 3 + 4 i , b = ( 1.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e₁ + ( -0.1 - 0.2 i ) e₂ + ( -0.3 - 0.4 i ) e₃ ( biquaternion -> 8 STO 00 )

1.1 STO 01 0.2 STO 02 0.3 STO 03 0.4 STO 04 -0.1 STO 05 -0.2 STO 06 -0.3 STO 07 -0.4 STO 08

       4 ENTER^
       3 ENTER^
       2 ENTER^
       1 XEQ "ALFB"   >>>>   1.008                               ---Execution time = 3m20s---

And we get in R01 thru R08:

>>> P_r^m(b) = ( 57.14529246 + 522.8336891 i ) + ( 346.1203818 - 37.93305254 i ) e₁
+ ( -155.3276125 - 10.99908746 i ) e₂ + ( -346.1203818 + 37.93305254 i ) e₃

m # 1 , 2 , 3 , .......... & b # 1

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

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

m = 0.2 + 0.3 i , r = 1.4 + 1.6 i , b = ( 0.05 + 0.1 i ) + ( 0.15 + 0.2 i ) e₁ + ( 0.25 + 0.3 i ) e₂ + ( 0.35 + 0.4 i ) e₃ ( biquaternion -> 8 STO 00 )

0.05 STO 01 0.1 STO 02 0.15 STO 03 0.2 STO 04 0.25 STO 05 0.3 STO6 0.35 STO 07 0.4 STO 08

      0.3   ENTER^
      0.2   ENTER^
      1.6   ENTER^
      1.4   XEQ "ALFB"   >>>>   1.008                               ---Execution time = 8m06s---

And we get in R01 thru R08:

Q_r^m(b) = ( -0.537350197 - 1.483228607 i ) + ( 1.318605388 - 1.639974390 i ) e₁
+ ( 1.925826457 - 2.663848476 i ) e₂ + ( 2.533047525 - 3.687722566 i ) e₃

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

m = 1.2 + 1.3 i , r = 1.4 + 1.5 , b = ( 1.1 + 1.2 i ) + ( 1.3 + 1.4 i ) e₁ + ( 1.5 + 1.6 i ) e₂ + ( 1.7 + 1.8 i ) e₃ ( biquaternion -> 8 STO 00 )

1.1 STO 01 1.2 STO 02 1.3 STO 03 1.4 STO 04 1.5 STO 05 1.6 STO 06 1.7 STO 07 1.8 STO 08

    1.5 ENTER^
    1.4 ENTER^
    1.3 ENTER^
    1.2 XEQ "ALFB" >>>>     1.008                               ---Execution time = 2m03s---

And we get in R01 thru R08:

Q_r^m(b) = ( 0.094724134 + 0.042366477 i ) + ( 0.021457752 - 0.047321593 i ) e₁
+ ( 0.024373027 - 0.054440173 i ) e₂ + ( 0.027288302 - 0.061558753 i ) e₃

>>> We must have r + 3/2 # 0 , 1 , 2 , ....................

Note:

-See also "Bionic Orthogonal Polynomials for the HP-41" , paragraph 7

2°) Chebyshev Functions

Formulae:

CF 02 T_m(cos B) = cos m.B ( first kind ) with cos B = b
SF 02 U_m(cos B) = [ sin (m+1).B ] / sin B ( second kind )

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

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

R2n+1 .......... R4n or R6n temp

>>> When the program stops: R01 ...... R2n = the n components of f(b)

Flags: F01-F02 CF 02 for the Chebyshev Functions of the 1st kind
SF 02 ----------------------------------- 2nd kind

Subroutines: AMOVE ASWAP Z*B STO U RCL U STO V RCL V X+1
"ACOSB" "COSB" "SINB" "1/B" B*B1 ( cf "Bions for the HP-41" )

01 LBL "CHBB"
02 STO U
03 X<>Y
04 STO V
05 XEQ "ACOSB"
06 13
07 FS? 02
08 AMOVE
09 RCL V
10 RCL U
11 FS? 02
12 X+1
13 Z*B
14 FS? 02
15 GTO 00
16 XEQ "COSB"
17 RTN
18 LBL 00
19 XEQ "SINB"
20 13
21 ASWAP
22 XEQ "SINB"
23 XEQ "1/B"
24 32
25 AMOVE
26 B*B1
27 END

( 78 bytes / SIZE 4n+1 or 6n+1 )

STACK	INPUTS	OUTPUTS
Y	m₂	/
X	m₁	1.2n

where m = m₁ + i m₂
and 1.2n is the control number of the result.

Example: m = 1.2 + 1.3 i

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 )

• CF 02 Chebyshev Function 1st kind

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

1.3 ENTER^
1.2 XEQ "CHBB" >>>> 1.008 ---Execution time = 51s---

And in R01 thru R08

T_m(b) = ( 1.070193727 + 1.688865678 i ) + ( -0.040122708 + 0.878077603 i ) e₁
+ ( 0.007654783 + 1.373010877 i ) e₂ + ( 0.055432275 + 1.867944152 i ) e₃

• SF 02 Chebyshev Function 2nd kind

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

1.3 ENTER^
1.2 XEQ "CHBB" >>>> 1.008 ---Execution time = 84s---

And in R01 thru R08

U_m(b) = ( 3.332558956 + 2.316362967 i ) + ( -0.806753418 + 1.096755418 i ) e₁
+ ( -1.170794897 + 1.775478725 i ) e₂ + ( -1.534836377 + 2.454202031 i ) e₃

3°) Incomplete Beta Function

Formula:

B_b(p,q) = ( b^p / p ) F(p,1-q;p+1;b) where F = the hypergeometric function

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

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

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

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

Flags: /

Subroutines: B*B1 "B^Z" Z*B DSA 1/Z Z*Z Z/Z AMOVE X+1 X/E3
STO R RCL R STO S RCL S STO U RCL U STO V RCL V STO W RCL W

01 LBL "IBFB"
02 STO R
03 RDN
04 STO S
05 RDN
06 STO U
07 X<>Y
08 STO V
09 12
10 AMOVE
11 X<> Z

12 XEQ "B^Z"
13 RCL V
14 RCL U
15 1/Z
16 Z*B
17 13
18 AMOVE
19 CLX
20 STO W
21 LBL 01
22 B*B1

23 RCL V
24 RCL U
25 RCL R
26 CHS
27 RCL W
28 ST+ Z
29 X+1
30 STO W
31 +
32 RCL S
33 CHS

34 X<>Y
35 Z*Z
36 RCL U
37 RCL W
38 ST+ Y
39 *
40 RCL V
41 ST* L
42 X<> L
43 X<>Y
44 Z/Z

45 Z*B
46 DSA
47 X#0?
48 GTO 01
49 31
50 AMOVE
51 RCL 00
52 X/E3
53 ISG X
54 END

( 105 bytes / SIZE 6n+1 )

STACK	INPUTS	OUTPUTS
T	p₂	/
Z	p₁	/
Y	q₂	/
X	q₁	1.2n

where p = p₁ + i p₂ , q = q₁ + i q₂
and 1.2n is the control number of the result.

Example: p = 0.2 + 0.3 i , q = 1.4 + 1.6 i ,

b = ( 0.05 + 0.1 i ) + ( 0.15 + 0.2 i ) e₁ + ( 0.25 + 0.3 i ) e₂ + ( 0.35 + 0.4 i ) e₃ ( biquaternion -> 8 STO 00 )

0.05 STO 01 0.1 STO 02 0.15 STO 03 0.2 STO 04 0.25 STO 05 0.3 STO6 0.35 STO 07 0.4 STO 08

      0.3   ENTER^
      0.2   ENTER^
      1.6   ENTER^
      1.4   XEQ "IBFB"   >>>>   1.008                               ---Execution time = 4m10s---

And we get in R01 thru R08:

B_b(p,q) = ( 0.781441740 - 1.586229595 i ) + ( 0.445532319 - 0.184816869 i ) e₁
+ ( 0.680245068 - 0.323956902 i ) e₂ + ( 0.914957818 - 0.463096935 i ) e₃

Note:

-Another formula may also be used: B_b(p,q) = ( b^p / p ) ( 1 - b )^q F(1,p+q;p+1;b)
-It gives better accuracy if q is a large real posive number, but with complexes, it's not sure...

01 LBL "IBFB"
02 STO R
03 RDN
04 STO S
05 RDN
06 STO U
07 X<>Y
08 STO V
09 12
10 AMOVE
11 13
12 AMOVE

13 RDN
14 X<> Z
15 XEQ "B^Z"
16 RCL V
17 RCL U
18 1/Z
19 Z*B
20 12
21 ASWAP
22 SIGN
23 CHS
24 ST*A

25 ST- 01
26 RCL S
27 RCL R
28 XEQ "B^Z"
29 B*B1
30 32
31 AMOVE
32 13
33 AMOVE
34 CLX
35 STO W
36 LBL 01

37 B*B1
38 RCL S
39 RCL R
40 RCL U
41 ST+ Y
42 RCL W
43 ST+ Z
44 X+1
45 STO W
46 +
47 RCL V
48 ST+ T

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

( 121 bytes / SIZE 6n+1 )

STACK	INPUTS	OUTPUTS
T	p₂	/
Z	p₁	/
Y	q₂	/
X	q₁	1.2n

where p = p₁ + i p₂ , q = q₁ + i q₂
and 1.2n is the control number of the result.

Example: With the same example p = 0.2 + 0.3 i , q = 1.4 + 1.6 i ,

b = ( 0.05 + 0.1 i ) + ( 0.15 + 0.2 i ) e₁ + ( 0.25 + 0.3 i ) e₂ + ( 0.35 + 0.4 i ) e₃ ( biquaternion -> 8 STO 00 )

0.05 STO 01 0.1 STO 02 0.15 STO 03 0.2 STO 04 0.25 STO 05 0.3 STO6 0.35 STO 07 0.4 STO 08

      0.3   ENTER^
      0.2   ENTER^
      1.6   ENTER^
      1.4   XEQ "IBFB"   >>>>   1.008                               ---Execution time = 6m48s---

And we get in R01 thru R08:

B_b(p,q) = ( 0.781441737 - 1.586229596 i ) + ( 0.445532322 - 0.184816870 i ) e₁
+ ( 0.680245069 - 0.323956905 i ) e₂ + ( 0.914957821 - 0.463096937 i ) e₃

Note:

-With this example, the 2nd version is slower and the precision is not better !

4°) Laguerre's Functions

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

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

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

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

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

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

Flags: /
Subroutines: B*B1 Z*B   ST*A DSA 1/Z Z/Z   AMOVE X+1 X/E3
                        STO R RCL R STO S RCL S 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 "LANB"
02 STO R
03 RDN
04 STO S
05 RDN
06 STO U
07 X<>Y
08 STO V
09 12
10 AMOVE
11 CLX

12 STO W
13 ST*A
14 SIGN
15 STO 01
16 13
17 AMOVE
18 LBL 01
19 B*B1
20 RCL S
21 CHS
22 RCL R

23 CHS
24 RCL U
25 RCL W
26 ST+ Z
27 X+1
28 STO W
29 +
30 RCL V
31 X<>Y
32 Z/Z
33 RCL W

34 ST/ Z
35 /
36 Z*B
37 DSA
38 X#0?
39 GTO 01
40 31
41 AMOVE
42 RCL S
43 RCL V
44 +

45 RCL R
46 RCL U
47 +
48 X+1
49 XEQ "GAMZ"
50 Z*B
51 RCL S
52 RCL R
53 X+1
54 XEQ "GAMZ"
55 1/Z

56 Z*B
57 RCL V
58 RCL U
59 X+1
60 XEQ "GAMZ"
61 1/Z
62 Z*B
63 RCL 00
64 X/E3
65 ISG X
66 END

( 135 bytes / SIZE 6n+1 )

STACK	INPUTS	OUTPUTS
T	a₂	/
Z	a₁	/
Y	m₂	/
X	m₁	1.2n

where a = a₁ + i a₂ , m = m₁ + i m₂
and 1.2n is the control number of the result.

Example: a = 0.2 + 0.3 i , m = 1.4 + 1.6 i ,

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

      0.3   ENTER^
      0.2   ENTER^
      1.6   ENTER^
      1.4   XEQ "LANB"   >>>>   1.008                               ---Execution time = 1m31s---

And we get in R01 thru R08:

L_m^(a)(b) = ( 2.195071146 + 1.392323359 i ) + ( 0.503589446 - 0.948842409 i ) e₁
+ ( 0.709692142 - 1.520481314 i ) e₂ + ( 0.915794840 - 2.092120220 i ) e₃

5°) Lerch Transcendent Function

Formula:

F( b , s , a ) = SUM_k=0,1,2,.... b^k / ( a + k )^s ( F is the greek letter "PHI" )

where s & a are complex numbers

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

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

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

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

Flags: /
Subroutines: B*B1 Z*B ST*A DSA Z*Z AMOVE X+1 X/E3
STO R RCL R STO S RCL S STO U RCL U STO V RCL V STO W RCL W

01 LBL "LERCHB"
02 RAD
03 STO R
04 RDN
05 STO S
06 RDN
07 CHS
08 STO U
09 X<>Y
10 CHS
11 STO V

12 12
13 AMOVE
14 CLX
15 STO W
16 ST*A
17 13
18 AMOVE
19 SIGN
20 STO 01
21 GTO 02
22 LBL 01

23 B*B1
24 RCL W
25 X+1
26 STO W
27 LBL 02
28 14
29 AMOVE
30 RCL S
31 RCL R
32 RCM W
33 +

34 R-P
35 LN
36 RCL V
37 RCL U
38 Z*Z
39 E^X
40 P-R
41 Z*B
42 DSA
43 41
44 AMOVE

45 X<>Y
46 X#0?
47 GTO 01
48 31
49 AMOVE
50 RCL 00
51 X/E3
52 ISG X
53 DEG
54 END

( 97 bytes / SIZE 8n+1 )

STACK	INPUTS	OUTPUTS
T	s₂	/
Z	s₁	/
Y	a₂	/
X	a₁	1.2n

where a = a₁ + i a₂ , s = s₁ + i s₂
and 1.2n is the control number of the result.

Example: s = 3.1 + 3.2 i , a = 1.2 + 1.4 i ,

b = ( 0.05 + 0.1 i ) + ( 0.15 + 0.2 i ) e₁ + ( 0.25 + 0.3 i ) e₂ + ( 0.31 + 0.4 i ) e₃ ( biquaternion -> 8 STO 00 )

0.05 STO 01 0.1 STO 02 0.15 STO 03 0.2 STO 04 0.25 STO 05 0.3 STO6 0.35 STO 07 0.4 STO 08

      3.1   ENTER^
      3.2   ENTER^
      1.2   ENTER^
      1.4   XEQ "LERCHB"   >>>>   1.008                               ---Execution time = 4m24s---

And we get in R01 thru R08:

F( b , s , a ) = ( -0.186440789 + 2.364004295 i ) + ( -0.058004148 + 0.054445111 i ) e₁
+ ( -0.086130862 + 0.089574705 i ) e₂ + ( -0.114257576 + 0.124704299 i ) e₃

6°) Lommel Functions

-"LOMB" calculates the Lommel functions of the 1st kind if flag F02 is clear and the Lommel functions of the 2nd kind if flag F02 is set.
-It uses the formulae:

s⁽¹⁾_m,p(b) = b^m+1 / [ (m+1)² - p² ] ₁F₂ ( 1 ; (m-p+3)/2 , (m+p+3)/2 ; -b²/4 )

    s⁽²⁾_m,p(b) = b^m+1 / [ (m+1)² - p² ] ₁F₂ ( 1 ; (m-p+3)/2 , (m+p+3)/2 ; b²/4 )
                      + 2^m+p-1 Gam(p) Gam((m+p+1)/2) b^-p / Gam((-m+p+1)/2) ₀F₁ ( ; 1-p ; -b²/4 )
                      + 2^m-p-1 Gam(-p) Gam((m-p+1)/2) b^p / Gam((-m-p+1)/2) ₀F₁ ( ; 1+p ; -b²/4 )

where _pF_q is the generalized hypergeometric function and Gam is Euler Gamma function.

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

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

R2n+1 .......... R6n or R10n temp

>>> When the program stops: R01 ...... R2n = the n components of f(b)

Flags: F01-F02 CF 02 for the Lommel Functions of the 1st kind
SF 02 ------------------------------- 2nd kind

Subroutines:   B*B1 Z*B "B^Z" ST*A ST/A A+A DSA 1/Z Z*Z Z^2 AMOVE ASWAP X+1 X-1 X/2 X/E3
                         STO R RCL R STO S RCL S 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 78 is a three-byte GTO 04

01 LBL "LOMB"
02 STO R
03 RDN
04 STO S
05 RDN
06 STO U
07 X<>Y
08 STO V
09 14
10 FS? 02
11 AMOVE
12 12
13 AMOVE
14 B*B1
15 4
16 ST/A
17 SIGN
18 CHS
19 FC? 02
20 ST*A
21 12
22 ASWAP
23 RCL V
24 RCL U
25 X+1
26 XEQ "B^Z"
27 RCL V
28 RCL U
29 X+1
30 Z^2
31 RCL S

32 RCL R
33 Z^2
34 Z-Z
35 1/Z
36 Z*B
37 13
38 AMOVE
39 CLX
40 STO W
41 LBL 01
42 B*B1
43 RCL V
44 RCL S
45 -
46 RCL U
47 X+1
48 RCL X
49 RCL R
50 ST- Z
51 +
52 2
53 ST/ T
54 ST/ Z
55 /
56 RCL W
57 X+1
58 STO W
59 ST+ Z
60 +
61 RCL V
62 SIGN

63 RDN
64 RCL S
65 ST+ L
66 X<> L
67 X/2
68 X<>Y
69 Z*Z
70 1/Z
71 Z*B
72 DSA
73 X#0?
74 GTO 01
75 31
76 AMOVE
77 FC? 02
78 GTO 04
79 15
80 AMOVE
81 21
82 AMOVE
83 SIGN
84 CHS
85 ST*A
86 12
87 AMOVE
88 SF 01
89 XEQ 02
90 RCL S
91 CHS
92 STO S
93 RCL R

94 CHS
95 STO R
96 LBL 02
97 41
98 AMOVE
99 RCL S
100 RCL R
101 XEQ "B^Z"
102 RCL V
103 RCL S
104 -
105 RCL U
106 RCL R
107 -
108 X+1
109 2
110 ST/ Z
111 /
112 XEQ "GAMZ"
113 Z*B
114 RCL S
115 CHS
116 RCL R
117 CHS
118 XEQ "GAMZ"
119 Z*B
120 RCL V
121 RCL S
122 -
123 2
124 LN

125 *
126 2
127 RCL U
128 RCL R
129 -
130 X-1
131 Y^X
132 RAD
133 P-R
134 Z*B
135 RCL V
136 RCL S
137 +
138 CHS
139 RCL U
140 RCL R
141 +
142 CHS
143 X+1
144 2
145 ST/ Z
146 /
147 XEQ "GAMZ"
148 1/Z
149 Z*B
150 13
151 AMOVE
152 CLX
153 STO W
154 LBL 03
155 B*B1

156 RCL S
157 RCL R
158 RCL W
159 X+1
160 STO W
161 ST+ Y
162 ST* Z
163 *
164 1/Z
165 Z*B
166 DSA
167 X#0?
168 GTO 03
169 31
170 AMOVE
171 25
172 ASWAP
173 A+A
174 52
175 AMOVE
176 15
177 AMOVE
178 FS?C 01
179 RTN
180 LBL 04
181 RCL 00
182 X/E3
183 ISG X
184 END

( 345 bytes / SIZE 6n+1 or 10n+1 )

STACK	INPUTS	OUTPUTS
T	m₂	/
Z	m₁	/
Y	p₂	/
X	p₁	1.2n

where m = m₁ + i m₂ , p = p₁ + i p₂
and 1.2n is the control number of the result.

Example: m = 0.2 + 0.3 i , p = 1.4 + 1.6 i

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 )

• CF 02 Lommel Function 1st kind

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

      0.3   ENTER^
      0.2   ENTER^
      1.6   ENTER^
      1.4   XEQ "TORB"   >>>>   1.008                               ---Execution time = 56s---

And we get in R01 thru R08:

s⁽¹⁾_m,p(b) = ( 0.093585318 + 0.033014799 i ) + ( -0.071915563 + 0.063783295 i ) e₁
+ ( -0.107085614 + 0.105255185 i ) e₂ + ( -0.142255666 + 0.146727076 i ) e₃

• SF 02 Lommel Function 2nd kind

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

      0.3   ENTER^
      0.2   ENTER^
      1.6   ENTER^
      1.4   XEQ "TORB"   >>>>   1.008                               ---Execution time = 3m39s---

And in R01 thru R08:

s⁽²⁾_m,p(b) = ( -1.431572986 - 2.691843964 i ) + ( -1.074362126 + 0.548149519 i ) e₁
+ ( -1.632152955 + 0.941062218 i ) e₂ + ( -2.189943785 + 1.333974919 i ) e₃

7°) Toronto Functions

Formula:

T(p,q;b) = exp(-b²) [ Gam((p+1)/2) / Gam(q+1) ] b^2q-p+1 M( (p+1)/2 , q+1 , b² ) where M = Kummer's function

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

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

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

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

Flags: /
Subroutines: B*B1 Z*B "B^Z" "E^B" ST*A DSA 1/Z Z/Z   AMOVE ASWAP X+1 X/2 X/E3
                        STO R RCL R STO S RCL S 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"TORB"
02 X+1
03 STO R
04 ST+ X
05 RDN
06 STO S
07 ST+ X
08 RDN
09 X+1
10 STO U
11 ST- Z
12 RDN
13 STO V
14 ST- Z
15 CLX
16 12

17 AMOVE
18 RDN
19 XEQ "B^Z"
20 13
21 AMOVE
22 21
23 AMOVE
24 B*B1
25 12
26 AMOVE
27 SIGN
28 CHS
29 ST*A
30 XEQ "E^B"
31 23
32 ASWAP

33 B*B1
34 32
35 AMOVE
36 13
37 AMOVE
38 CLX
39 STO W
40 LBL 01
41 B*B1
42 RCL V
43 X/2
44 RCL U
45 X/2
46 RCL R
47 RCL W
48 ST+ Z

49 +
50 RCL S
51 X<>Y
52 Z/Z
53 RCL W
54 X+1
55 STO W
56 ST/ Z
57 /
58 Z*B
59 DSA
60 X#0?
61 GTO 01
62 31
63 AMOVE
64 RCL V

65 RCL U
66 2
67 ST/ Z
68 /
69 XEQ "GAMZ"
70 Z*B
71 RCL S
72 RCL R
73 XEQ "GAMZ"
74 1/Z
75 Z*B
76 RCL 00
77 X/E3
78 ISG X
79 END

( 163 bytes / SIZE 6n+1 )

STACK	INPUTS	OUTPUTS
T	p₂	/
Z	p₁	/
Y	q₂	/
X	q₁	1.2n

where p = p₁ + i p₂ , q = q₁ + i q₂
and 1.2n is the control number of the result.

Example: p = 0.2 + 0.3 i , q = 1.4 + 1.6 i ,

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

      0.3   ENTER^
      0.2   ENTER^
      1.6   ENTER^
      1.4   XEQ "TORB"   >>>>   1.008                               ---Execution time = 2m02s---

And we get in R01 thru R08:

T(p,q;b) = ( 2.358850025 - 5.982167499 i ) + ( 2.067337951 + 0.962846884 i ) e₁
+ ( 3.302074949 + 1.336654106 i ) e₂ + ( 4.536811954 + 1.710461323 i ) e₃

8°) Ultraspherical Functions

Formula: Assuming p # 0

C_q^(p)(b) = [ Gam(q+2p) / Gam(q+1) / Gam(2.p) ] ₂F₁( -q , q+2p , q+1/2 , (1-b)/2 )

where ₂F₁ is the hypergeometric function and Gam = Gamma function

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

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

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

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

Flags: /
Subroutines: B*B1 Z*B   ST*A DSA 1/Z Z*Z Z/Z   AMOVE X+1 X/E3
                        STO R RCL R STO S RCL S 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 "USFB"
02 STO R
03 RDN
04 STO S
05 RDN
06 STO U
07 X<>Y
08 STO V
09 .5
10 CHS
11 ST*A
12 ST- 01
13 12
14 AMOVE
15 CLX
16 STO W
17 ST*A

18 SIGN
19 STO 01
20 13
21 AMOVE
22 LBL 01
23 B*B1
24 RCL V
25 ST+ X
26 RCL U
27 ST+ X
28 RCL R
29 ST+ Y
30 CHS
31 RCL W
32 ST+ Z
33 +
34 RCL S

35 ST+ T
36 CHS
37 X<>Y
38 Z*Z
39 RCL U
40 .5
41 +
42 RCL W
43 ST+ Y
44 X+1
45 STO W
46 *
47 RCL V
48 ST* L
49 X<> L
50 X<>Y
51 Z/Z

52 Z*B
53 DSA
54 X#0?
55 GTO 01
56 31
57 AMOVE
58 RCL V
59 ST+ X
60 RCL S
61 +
62 RCL U
63 ST+ X
64 RCL R
65 +
66 XEQ "GAMZ"
67 Z*B
68 RCL S

69 RCL R
70 X+1
71 XEQ "GAMZ"
72 1/Z
73 Z*B
74 RCL V
75 ST+ X
76 RCL U
77 ST+ X
78 XEQ "GAMZ"
79 1/Z
80 Z*B
81 RCL 00
82 X/E3
83 ISG X
84 END

( 168 bytes / SIZE 6n+1 )

STACK	INPUTS	OUTPUTS
T	p₂	/
Z	p₁	/
Y	q₂	/
X	q₁	1.2n

where p = p₁ + i p₂ , q = q₁ + i q₂
and 1.2n is the control number of the result.

Example: p = 0.2 + 0.3 i , q = 1.4 + 1.6 i ,

b = ( 1.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 )

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

      0.3   ENTER^
      0.2   ENTER^
      1.6   ENTER^
      1.4   XEQ "USFB"   >>>>   1.008                               ---Execution time = 6m11s---

And we get in R01 thru R08:

C_q^(p)(b) = ( -1.082467436 + 0.708001467 i ) + ( -0.196556331 - 0.294814043 i ) e₁
+ ( -0.330212998 - 0.444185398 i ) e₂ + ( -0.463869666 - 0.593556755 i ) e₃

9°) Whittaker's M- & W-Functions

Formulae:

M_q,p(b) = exp(-b/2) b^p+1/2 M( p-q+1/2 , 1+2p , b ) where M = Kummer's functions

W_q,p(b) = [ Gam(2p) / Gam(p-q+1/2) ] M_q,-p(b) + [ Gam(-2p) / Gam(-p-q+1/2) ] M_q,p(b) assuming 2p is not an integer.

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

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

R2n+1 .......... R6n+1 or R8n+1: temp

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

Flags: /
Subroutines: B*B1 Z*B "B^Z" "E^B" A+A ST/A DSA 1/Z Z/Z   AMOVE ASWAP X+1 X/E3
                        STO R RCL R STO S RCL S 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 "WHIWB"
02 XEQ 01
03 14
04 AMOVE
05 21
06 AMOVE
07 RCL V
08 RCL U
09 RCL S
10 CHS
11 RCL R
12 CHS
13 XEQ 01
14 42
15 AMOVE
16 A+A
17 RTN
18 LBL 01
19 XEQ 02
20 RCL S
21 RCL R

22 2
23 CHS
24 ST* Z
25 *
26 XEQ "GAMZ"
27 Z*B
28 RCL V
29 RCL S
30 +
31 CHS
32 .5
33 RCL U
34 -
35 RCL R
36 -
37 XEQ "GAMZ"
38 1/Z
39 Z*B
40 RTN
41 LBL 02
42 LBL "WHIMB"

43 STO R
44 RDN
45 STO S
46 RDN
47 STO U
48 X<>Y
49 STO V
50 12
51 AMOVE
52 13
53 AMOVE
54 RCL S
55 RCL R
56 .5
57 +
58 XEQ "B^Z"
59 12
60 ASWAP
61 2
62 CHS
63 ST/A

64 XEQ "E^B"
65 B*B1
66 32
67 AMOVE
68 13
69 AMOVE
70 CLX
71 STO W
72 LBL 03
73 B*B1
74 RCL V
75 CHS
76 .5
77 RCL U
78 -
79 RCL R
80 ST+ Y
81 ST+ X
82 RCL W
83 ST+ Z
84 X+1

85 STO W
86 +
87 RCL S
88 ST+ T
89 ST+ X
90 X<>Y
91 Z/Z
92 RCL W
93 ST/ Z
94 /
95 Z*B
96 DSA
97 X#0?
98 GTO 03
99 31
100 AMOVE
101 RCL 00
102 X/E3
103 ISG X
104 END

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

STACK	INPUTS	OUTPUTS
T	q₂	/
Z	q₁	/
Y	p₂	/
X	p₁	1.2n

where q = q₁ + i q₂ , p = p₁ + i p₂
and 1.2n is the control number of the result.

Example: q = 0.2 + 0.3 i , p = 1.4 + 1.6 i ,

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 )

• Whittaker M-function

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

      0.3   ENTER^
      0.2   ENTER^
      1.6   ENTER^
      1.4   XEQ "WHIMB"   >>>>   1.008                               ---Execution time = 1m35s---

And in R01 thru R08:

M_q,p(b) = ( 1.697417279 - 1.038735336 i ) + ( 0.312886027 + 0.617812026 i ) e₁
+ ( 0.537527164 + 0.938755878 i ) e₂ + ( 0.762168301 + 1.259699731 i ) e₃

• Whittaker W-function

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

      0.3   ENTER^
      0.2   ENTER^
      1.6   ENTER^
      1.4   XEQ "WHIWB"   >>>>   1.008                               ---Execution time = 4m00s---

And we get in R01 thru R08:

W_q,p(b) = ( 7.064436874 + 2.782608838 i ) + ( 1.108472812 - 2.517738290 i ) e₁
+ ( 1.527798530 - 4.016349561 i ) e₂ + ( 1.947124243 - 5.514960825 i ) e₃

10°) Riemann Zeta Function

-The following program employs the method given by P. Borwein in "An Efficient Algorithm for the Riemann Zeta Function" if ReRe(b) >= 1/2
-If ReRe(b) < 1/2, it uses: Zeta(b) = Zeta(1-b) Pi^b-1/2 Gamma((1-b)/2) / Gamma(b/2)

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

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

R2n+1 .......... R6n+1 or R12n+1: temp

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

Flags: /
Subroutines: B*B1 "1/B" "X^B" ST*A   ST/A DSA 1/Z Z/Z   AMOVE ASWAP X+1 X-1 X/E3 3X
                        STO R RCL R STO S RCL S STO U RCL U STO V RCL V
                        "GAMB" ( cf "Bionic Special Functions(I) for the HP-41" )

01 LBL "ZETAB"
02 XEQ 13
03 RCL 01
04 X^2
05 +
06 X#0?
07 GTO 00
08 .5
09 CHS
10 STO 01
11 RCL 00
12 X/E3
13 ISG X
14 RTN
15 LBL 00
16 RCL 01
17 .5
18 X<=Y?
19 GTO 00
20 16
21 AMOVE
22 X<>Y
23 ST- 01
24 PI
25 XEQ "X^B"
26 14
27 AMOVE
28 61
29 AMOVE
30 SIGN
31 CHS
32 ST*A
33 ST- 01

34 XEQ 00
35 42
36 AMOVE
37 B*B1
38 15
39 AMOVE
40 61
41 AMOVE
42 2
43 CHS
44 ST/A
45 1/X
46 ST- 01
47 XEQ "GAMB"
48 52
49 AMOVE
50 B*B1
51 15
52 AMOVE
53 61
54 AMOVE
55 2
56 ST/A
57 XEQ "GAMB"
58 XEQ "1/B"
59 52
60 AMOVE
61 B*B
62 RTN
63 LBL 13
64 RCL 00
65 E-3
66 +

67 0
68 LBL 14
69 RCL IND Y
70 X^2
71 +
72 DSE Y
73 GTO 14
74 SQRT
75 RTN
76 LBL 00
77 12
78 AMOVE
79 RCL 00
80 3X
81 .1
82 %
83 X+1
84 +
85 RCL 00
86 -
87 CLRGX
88 XEQ 13
89 PI
90 X/2
91 X<>Y
92 ST* Y
93 ST+ X
94 LN1+X
95 +
96 3 E10
97 LN
98 +
99 8

100 SQRT
101 3
102 +
103 LN
104 /
105 INT
106 1
107 STO U
108 STO V
109 +
110 STO R
111 STO S
112 LBL 01
113 21
114 AMOVE
115 RCL R
116 1/X
117 XEQ "X^B"
118 SIGN
119 CHS
120 RCL R
121 Y^X
122 RCL U
123 *
124 ST*A
125 DSA
126 RCL R
127 ENTER^
128 ST+ Y
129 ST* Y
130 -
131 RCL V
132 *

133 RCL S
134 X^2
135 RCL R
136 X-1
137 STO R
138 X^2
139 -
140 ST+ X
141 /
142 STO V
143 RCL U
144 +
145 STO U
146 RCL R
147 X#0?
148 GTO 01
149 21
150 AMOVE
151 SIGN
152 CHS
153 ST*A
154 ST- 01
155 2
156 XEQ "X^B"
157 SIGN
158 ST- 01
159 XEQ "1/B"
160 32
161 AMOVE
162 RCL U
163 ST/A
164 B*B1
165 END

( 306 bytes / SIZE 6n+1 or 12n+1 )

STACK	INPUT	OUTPUT
X	/	1.2n

where 1.2n is the control number of the result.

Example1: b = ( 1.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 )

1.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 "ZETAB" >>>> 1.008 ---Execution time = 3m31s---

And we get in R01 thru R08

Zeta(b) = ( 0.670436770 - 0.041923691 i ) + ( -0.145445652 + 0.208534940 i ) e₁
+ ( -0.210212423 + 0.336950158 i ) e₂ + ( -0.274979192 + 0.465365376 i ) e₃

Example2: 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 "ZETAB" >>>> 1.008 ---Execution time = 6m06s---

And in R01 thru R08

Zeta(b) = ( 0.478160113 + 0.577429237 i ) + ( -0.256274162 + 0.142614011 i ) e₁
+ ( -0.388378571 + 0.242979789 i ) e₂ + ( -0.520482981 + 0.343345568 i ) e₃

hp41programs

Bionic Special Functions(II) for the HP-41