# hp41programs

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      Prm(b) = [ (b+1)/(1-b) ]m/2  2F1(-r , r+1 ; 1-m ; (1-b)/2 ) / Gamma(1-m)        (  b # 1 )

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

Formulae:    2nd Kind

Type 2      Qrm(b) = 2m pi1/2  (1-b2)-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 ; b2 )
+ b  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 ; b2 )  ]

Type 3        Qrm(b) =  exp( i (m.PI) ) 2 -r-1 sqrt(PI) Gam(m+r+1) b -m-r-1 (b2-1)m/2  2F~1( (2+m+r)/2 , (1+m+r)/2 ; r+3/2 ; 1/b2 )

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 m2 / Z m1 / Y r2 / X r1 1.2n

where        m = m1 + i m2  ,   r = r1 + i r2
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 ) e1 + ( -0.1 - 0.2 i ) e2 + ( -0.3 - 0.4 i ) e3                 ( 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:

>>>     Prm(b) = ( 93.10375888 - 105.0475343 i )  + ( 75.56897000 + 91.98871986 i ) e1
+ ( -25.89125050 - 46.52055456 i ) e2 + ( -75.56897000 - 91.98871986 i ) e3

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 ) e1 + ( -0.1 - 0.2 i ) e2 + ( -0.3 - 0.4 i ) e3                 ( 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:

>>>     Prm(b) = ( 57.14529246 + 522.8336891 i )  + ( 346.1203818 - 37.93305254 i ) e1
+ ( -155.3276125 - 10.99908746 i ) e2 + ( -346.1203818 + 37.93305254 i ) e3

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 ) e1 + ( 0.25 + 0.3 i ) e2 + ( 0.35 + 0.4 i ) e3          ( 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:

Qrm(b) = ( -0.537350197 - 1.483228607 i )  + ( 1.318605388 - 1.639974390 i ) e1
+ ( 1.925826457 - 2.663848476 i ) e2 + ( 2.533047525 - 3.687722566 i ) e3

•  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 ) e1 + ( 1.5 + 1.6 i ) e2 + ( 1.7 + 1.8 i ) e3               ( 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:

Qrm(b) = ( 0.094724134 + 0.042366477 i )  + ( 0.021457752 - 0.047321593 i ) e1
+ ( 0.024373027 - 0.054440173 i ) e2 + ( 0.027288302 - 0.061558753 i ) e3

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

Note:

2°)  Chebyshev Functions

Formulae:

CF 02          Tm(cos B) = cos m.B                          ( first kind )                           with       cos B = b
SF 02          Um(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 m2 / X m1 1.2n

where       m = m1 + i m2
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 ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3               ( 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

Tm(b) =  ( 1.070193727 + 1.688865678 i )  + ( -0.040122708 + 0.878077603 i ) e1
+ (  0.007654783 + 1.373010877 i ) e2 + ( 0.055432275 + 1.867944152 i ) e3

•  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

Um(b) =  ( 3.332558956 + 2.316362967 i )  + ( -0.806753418 + 1.096755418 i ) e1
+ ( -1.170794897 + 1.775478725 i ) e2 + ( -1.534836377 + 2.454202031 i ) e3

3°)  Incomplete Beta Function

Formula:

Bb(p,q) = ( bp / 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 p2 / Z p1 / Y q2 / X q1 1.2n

where        p = p1 + i p2  ,   q = q1 + i q2
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 ) e1 + ( 0.25 + 0.3 i ) e2 + ( 0.35 + 0.4 i ) e3          ( 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:

Bb(p,q) =  ( 0.781441740 - 1.586229595 i )  + ( 0.445532319 - 0.184816869 i ) e1
+ ( 0.680245068 - 0.323956902 i ) e2 + ( 0.914957818 - 0.463096935 i ) e3

Note:

-Another formula may also be used:   Bb(p,q) = ( bp / 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 p2 / Z p1 / Y q2 / X q1 1.2n

where        p = p1 + i p2  ,   q = q1 + i q2
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 ) e1 + ( 0.25 + 0.3 i ) e2 + ( 0.35 + 0.4 i ) e3          ( 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:

Bb(p,q) =  ( 0.781441737 - 1.586229596 i )  + ( 0.445532322 - 0.184816870 i ) e1
+ ( 0.680245069 - 0.323956905 i ) e2 + ( 0.914957821 - 0.463096937 i ) e3

Note:

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

4°)  Laguerre's Functions

Formula:      Lm(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 a2 / Z a1 / Y m2 / X m1 1.2n

where        a = a1 + i a2  ,   m = m1 + i m2
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 ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3               ( 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:

Lm(a)(b) =  ( 2.195071146 + 1.392323359 i )  + ( 0.503589446 - 0.948842409 i ) e1
+ ( 0.709692142 - 1.520481314 i ) e2 + ( 0.915794840 - 2.092120220 i ) e3

5°)  Lerch Transcendent Function

Formula:

F( b , s , a ) = SUM k=0,1,2,....   bk / ( 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 s2 / Z s1 / Y a2 / X a1 1.2n

where        a = a1 + i a2  ,   s = s1 + i s2
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 ) e1 + ( 0.25 + 0.3 i ) e2 + ( 0.31 + 0.4 i ) e3               ( 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 ) e1
+ ( -0.086130862 + 0.089574705 i ) e2 + ( -0.114257576 + 0.124704299 i ) e3

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(1)m,p(b)  = bm+1 / [ (m+1)2 - p21F2 ( 1 ; (m-p+3)/2 , (m+p+3)/2 ; -b2/4 )

s(2)m,p(b)  = bm+1 / [ (m+1)2 - p21F2 ( 1 ; (m-p+3)/2 , (m+p+3)/2 ; b2/4 )
+ 2m+p-1 Gam(p) Gam((m+p+1)/2) b -p / Gam((-m+p+1)/2) 0F1 ( ; 1-p ; -b2/4 )
+ 2m-p-1 Gam(-p) Gam((m-p+1)/2) bp / Gam((-m-p+1)/2) 0F1 ( ; 1+p ; -b2/4 )

where  pFq 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 m2 / Z m1 / Y p2 / X p1 1.2n

where        m = m1 + i m2  ,   p = p1 + i p2
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 ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3               ( 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(1)m,p(b)  =  ( 0.093585318 + 0.033014799 i )  + ( -0.071915563 + 0.063783295 i ) e1
+ ( -0.107085614 + 0.105255185 i ) e2 + ( -0.142255666 + 0.146727076 i ) e3

•  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(2)m,p(b)  =  ( -1.431572986 - 2.691843964 i )  + ( -1.074362126 + 0.548149519 i ) e1
+ ( -1.632152955 + 0.941062218 i ) e2 + ( -2.189943785 + 1.333974919 i ) e3

7°)  Toronto Functions

Formula:

T(p,q;b) = exp(-b2) [ Gam((p+1)/2) / Gam(q+1) ] b2q-p+1  M( (p+1)/2 , q+1 , b2 )                    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 p2 / Z p1 / Y q2 / X q1 1.2n

where        p = p1 + i p2  ,   q = q1 + i q2
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 ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3               ( 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 ) e1
+ ( 3.302074949 + 1.336654106 i ) e2 + ( 4.536811954 + 1.710461323 i ) e3

8°)  Ultraspherical Functions

Formula:          Assuming  p # 0

Cq(p)(b) = [ Gam(q+2p) / Gam(q+1) / Gam(2.p) ]  2F1( -q , q+2p , q+1/2 , (1-b)/2 )

where  2F1 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 p2 / Z p1 / Y q2 / X q1 1.2n

where        p = p1 + i p2  ,   q = q1 + i q2
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 ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3               ( 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:

Cq(p)(b) =  ( -1.082467436 + 0.708001467 i )  + ( -0.196556331 - 0.294814043 i ) e1
+ ( -0.330212998 - 0.444185398 i ) e2 + ( -0.463869666 - 0.593556755 i ) e3

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

Formulae:

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

Wq,p(b) = [ Gam(2p) / Gam(p-q+1/2) ]  Mq,-p(b) + [ Gam(-2p) / Gam(-p-q+1/2) ]  Mq,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 q2 / Z q1 / Y p2 / X p1 1.2n

where        q = q1 + i q2   ,   p = p1 + i p2
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 ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3               ( 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:

Mq,p(b) = ( 1.697417279 - 1.038735336 i )  + ( 0.312886027 + 0.617812026 i ) e1
+ ( 0.537527164 + 0.938755878 i ) e2 + ( 0.762168301 + 1.259699731 i ) e3

•  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:

Wq,p(b) = ( 7.064436874 + 2.782608838 i )  + ( 1.108472812 - 2.517738290 i ) e1
+ ( 1.527798530 - 4.016349561 i ) e2 + ( 1.947124243 - 5.514960825 i ) e3

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 ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3               ( 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 ) e1
+ ( -0.210212423 + 0.336950158 i ) e2 + ( -0.274979192 + 0.465365376 i ) e3

Example2:          b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3               ( 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 ) e1
+ ( -0.388378571 + 0.242979789 i ) e2 + ( -0.520482981 + 0.343345568 i ) e3