# hp41programs

Anionic Functions2

# Anionic Special Functions(II) for the HP-41

Overview

1°)  Spheroidal Wave Functions
2°)  Jacobian Elliptic Functions

a)  Sn( a | m ) , Cn( a | m ) , Dn( a | m )
b)  Sn( a | m ) only

3°)  Catalan Numbers
4°)  Lambert W Function
5°)  Weierstrass Elliptic Functions

a)  Laurent Series
b)  Duplication Formula

6°)  Generalized Hypergeometric Functions

a)  General Case
b)  p+q < 4
c)  0F1 ( 2 versions )
d)  2F0

7°)  Airy Functions
8°)  Hermite Functions
9°)  Chebyshev Functions
10°) Riemann Zeta Function

a)  If  Re(a) > 1
b)  Borwein Algorithm

11°) Lommel Functions

1°)  Spheroidal Wave Functions

-"SWFA" computes the angular spheroidal wave function of the first kind.
-Given  m , n  and c2 , the corresponding eigenvalue  l  may be calculated by a program listed in "Spheroidal Harmonics for the HP-41"

-We assume that  | a | <= 1  and  Smn(a) is computed by      Smn(a) = ( 1 - a2 ) m/2   Sumk=0,1,.... dk ak

with     (k+1)(k+2) dk+2 - [ k ( k + 2m + 1 ) - l + m ( m + 1 ) ] dk  - c2  dk-2 = 0

Flammer's Scheme:           the coefficients are normalized as follows:

d0 = Pnm(0)  =  2m sqrt(PI) / [ Gam((1-m-n)/2) Gam((2-m+n)/2 ]
d1 = P'nm(0) = ( m + n ) 2m sqrt(PI) / [ Gam((2-m-n)/2) Gam((1-m+n)/2 ]

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "SWFA" )

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

Rn+1 ..........  R3n:  temp   R3n+1 = m ,  R3n+2 = n ,  R3n+3 = c2 ,  R3n+4 = l

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

Flag:  F01
Subroutines:  "ST*A"    ( cf  "Anionic Special Functions(I) "paragraph 0 )
"A*A1"    ( cf "Anions for the HP-41" )
"1/G+"      ( cf "Gamma Function for the HP-41" )

-Lines 118 to 138 may be replaced by a single instruction:   DS*A  ( an M-Code routine listed in "Anionic Special Functions" §0°)c)

 01  LBL "SWFA"   02  R^   03  STO M   04  X<> T   05  STO N   06  CLX   07  RCL 00   08  3   09  *   10  4    11  +    12  RDN   13  STO IND T   14  RDN   15  DSE Z   16  STO IND Z   17  CLX   18  RCL N   19  DSE Z   20  STO IND Z   21  RCL M   22  DSE T   23  STO IND T    24  +   25  CHS   26  2   27  ST+ Y   28  /   29  XEQ "1/G+"   30  STO O   31  1   32  RCL M    33   -   34  RCL N   35  ST+ M 36  +   37  2   38  /    39  XEQ "1/G+"   40  RCL O   41  *   42  RCL M   43  *   44  STO N   45  RCL 00   46  .1   47  %   48  ISG X   49  +   50   E3   51  /   52  1   53  +    54  STO M   55  REGMOVE   56  XEQ "A*A1"   57  CLX    58  X<> M   59  REGSWAP   60  SIGN    61  STO O   62  RCL 00   63  3   64  *   65  .1   66  %   67  ISG X   68  +   69  RCL 00   70  - 71  CLRGX   72  SF 01   73  RCL N   74  GTO 02   75  LBL 01   76  XEQ "A*A1"   77  RCL 00   78  3   79  *   80  1   81  +   82  RCL IND X   83  ENTER^   84  ST+ X   85  1   86  +   87  RCL O   88  ST+ Y   89  *   90  3   91  ST+ T   92  CLX   93  RCL IND T   94  -   95  X<>Y   96  X^2   97  LASTX   98  +   99  + 100  RCL N 101  ST* Y 102  X<> M 103  DSE Z 104  RCL IND Z 105  * 106  + 107  RCL O 108  1 109  + 110  RCL X 111  LASTX 112  + 113  STO O 114  * 115  / 116  LBL 02 117  STO N 118  RCL 00 119  3 120  * 121  STO P 122  RCL 00 123  ENTER^ 124  CLX 125  LBL 03 126  RCL IND Y 127  RCL N 128  * 129  RCL IND P 130  + 131  STO IND P 132  LASTX 133  - 134  ABS 135  + 136  DSE P 137  DSE Y 138  GTO 03 139  X#0? 140  GTO 01 141  CLA 142  RCL 00 143  3 144  * 145  2 146  + 147  RCL IND X 148  STO N 149  DSE Y 150  RCL IND Y 151  FC?C 01 152  GTO 00 153  ST+ N 154  - 155  2 156  ST+ Y 157  / 158  XEQ "1/G+" 159  X<> N 160  1 161  X<>Y 162  - 163  2 164  / 165  XEQ "1/G+" 166  RCL N 167  * 168  0 169  XEQ "ST*A" 170  SIGN 171  STO 01 172  X<>Y 173  GTO 02 174  LBL 00 175  STO M 176  RCL 00 177  ST+ X 178   E3 179   / 180  ISG X 181   E3 182  / 183  1 184  + 185  RCL 00 186  + 187  REGMOVE 188  SIGN 189  CHS 190  XEQ "ST*A" 191  ST- 01 192  RCL M 193  2 194  / 195  XEQ "A^X" 196  XEQ "A*A1" 197  2 198  RCL M 199  Y^X 200  PI 201  SQRT 202  * 203  XEQ "ST*A" 204  X<>Y 205  CLA 206  END

( 347 bytes / SIZE 3n+5 )

 STACK INPUTS OUTPUTS T m / Z n / Y c2 / X l 1.nnn

where   1.nnn   is the control number of the result.

Example:    m = 0.2  ,  n = 0.6 ,  c2 = 1.7 ,  l = 2.246866651     a = 0.1 + 0.2 i + 0.3 j + 0.4 k                 ( Quaternion ->   4   STO 00 )

0.1  STO 01
0.2  STO 02
0.3  STO 03
0.4  STO 04

0.2   ENTER^
0.6   ENTER^
1.7   ENTER^
2.246866651   XEQ "SWFA"  >>>>   1.004                             ---Execution time = 86s---            ( with  DS*A )

R01 = 0.391049684
R02 = 0.161100168
R03 = 0.241650252
R04 = 0.322200335

-Whence     S( 0.1 + 0.2 i + 0.3 j + 0.4 k ) = 0.391049684 + 0.161100168 i + 0.241650252 j + 0.322200335 k

2°)  Jacobian Elliptic Functions

a) Sn( a | m ) , Cn( a | m ) , Dn( a | m )

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

-If   m < 1 , we have:

-With  m' = 1-m ,  let be  µ = [ ( 1-sqrt(m') / ( 1+sqrt(m') ]2   and   v = a / ( 1+sqrt(µ) ]  , then:

sn ( a | m ) = [ ( 1 + sqrt(µ) ) sn ( v | µ ) ] / [ 1 + sqrt(µ) sn2 ( v | µ ) ]
cn ( a | m ) = [ cn ( v | µ ) dn ( v | µ ) ] / [ 1 + sqrt(µ) sn2 ( v | µ ) ]
dn ( a | m ) = [ 1 - sqrt(µ) sn2 ( v | µ ) ] / [ 1 + sqrt(µ) sn2 ( v | µ ) ]

-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, the program uses the relations:

sn ( a | m ) = sn ( a m1/2 | 1/m ) / m1/2
cn ( a | m ) = dn ( a m1/2 | 1/m )
dn ( a | m ) = cn ( a m1/2 | 1/m )

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

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "JEFA" )

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

Rn+1 ..........  R4n:  temp   R4n+1 = m ,  R4n+2 ..... temp

>>>  When the program stops:     R01   ......  Rnn = the n components of sn ( a | m )
Rn+1 ......  R2n = --------------------  cn ( a | m )
R2n+1 ....  R3n = --------------------  dn ( a | m )

Flags:  F01-F10
Subroutines:  "ST*A"  "ST/A"                                                        ( cf  "Anionic Special Functions(I) "paragraph 0 )
"A*A1"  "1/A"  "CHA"  "THA"  "SINA"  "COSA"    ( cf "Anions for the HP-41" )

-Lines 35 & 189 are three-byte GTOs

 01  LBL "JEFA"   02  RCL 00   03  4   04  *   05  1   06  +   07  X<>Y   08  STO IND Y   09  RCL 00   10  .1   11  %   12  ISG X   13  +   14   E3   15  /   16  1   17  +   18  STO N   19  REGMOVE   20  CF 10   21  SIGN   22  X>Y?   23  GTO 01   24  X#Y?   25  GTO 00   26  XEQ "CHA"   27  XEQ "1/A"   28  FRC   29  RCL N   30  +   31  REGMOVE   32  LASTX   33  REGSWAP   34  XEQ "THA"   35  GTO 04   36  LBL 00   37  X<>Y   38  SQRT 39  XEQ "ST*A"   40  LASTX   41  1/X   42  ENTER^   43  SF 10   44  LBL 01   45  X<>Y   46  STO M   47  RCL 00   48  4   49  *   50  1   51  +   52  .1   53  %   54  +   55  STO O   56  SIGN   57  X<> M   58  LBL 02   59  ENTER^   60  CHS   61  1   62  +   63  SQRT   64  1   65  ST+ O   66  +   67  RCL X   68  2   69  /   70  ST* M   71  RDN   72  X^2   73  /   74  STO IND O   75  X^2   76  X#0? 77  GTO 02   78  RCL 00   79  4   80  *   81  .1   82  %   83  ISG X   84  +   85  RCL 00   86  -   87  CLRGX   88  SIGN   89  STO IND L   90  RCL M   91  XEQ "ST*A"   92  RCL N   93  REGMOVE   94  XEQ "SINA"   95  RCL N   96  REGSWAP   97  XEQ "COSA"   98  RCL N   99  RCL 00 100   E3 101  / 102  STO M 103  + 104  REGMOVE 105  LBL 03 106  RCL 00 107  3 108  * 109  RCL M 110  - 111  RCL N 112  + 113  REGMOVE 114  LASTX 115  RCL 00 116  ST+ X 117  + 118  REGSWAP 119  XEQ "A*A1" 120  FRC 121  RCL N 122  + 123  REGSWAP 124  LASTX 125  REGMOVE 126  XEQ "A*A1" 127  RCL IND O 128  XEQ "ST*A" 129  RCL 00 130  ST+ X 131  RCL N 132  + 133  REGSWAP 134  LASTX 135  RCL M 136  ST+ X 137  + 138  REGMOVE 139  SIGN 140  ST+ 01 141  XEQ "1/A" 142  XEQ "A*A1" 143  FRC 144  RCL N 145  + 146  REGSWAP 147  SIGN 148  RCL IND O 149  + 150  XEQ "ST*A" 151  RCL N 152  REGMOVE 153  RCL M 154  - 155  RCL 00 156  3 157  * 158  + 159  REGMOVE 160  SIGN 161  ST+ 01 162  XEQ "1/A" 163  XEQ "A*A1" 164  FRC 165  ST+ X 166  RCL N 167  + 168  REGSWAP 169  LASTX 170  REGMOVE 171  SIGN 172  CHS 173  XEQ "ST*A" 174  ST- 01 175  RCL N 176  REGSWAP 177  SIGN 178  ST+ 01 179  XEQ "1/A" 180  XEQ "A*A1" 181  FRC 182  ST+ X 183  RCL N 184  + 185  REGSWAP 186  LASTX 187  REGMOVE 188  DSE O 189  GTO 03 190  RCL M 191  3 192  * 193  ISG X 194   E3 195  / 196  1 197  + 198  RCL 00 199  + 200  REGMOVE 201  FC?C 10 202  GTO 04 203  RCL 00 204  4 205  * 206  1 207  + 208  RCL IND X 209  SQRT 210  XEQ "ST/A" 211  RCL N 212  RCL 00 213  ST+ X 214  + 215  REGSWAP 216  LBL 04 217  CF 01 218  RCL 00 219   E3 220  / 221  ISG X 222  CLA 223  END

( 410 bytes / SIZE 4n+10 )

 STACK INPUTS OUTPUTS X m 1.nnn

where   1.nnn   is the control number of  sn ( a | m )

Example1:      m = 1/PI   q = 1 + 2 i + 3 j + 4 k      ( Quaternion ->   4   STO 00 )

1  STO 01
2  STO 02
3  STO 03
4  STO 04

PI  1/X   XEQ "JEFA"   >>>> 1.004                    ---Execution time = 70s---

R01 =  1.498262248
R02 =  0.205407205
R03 =  0.308110807
R04 =  0.410814410

sn ( 1 + 2 i + 3 j + 4 k | 1/PI ) = 1.498262248 + 0.205407205 i + 0.308110807 j + 0.410814410 k    and

cn ( 1 + 2 i + 3 j + 4 k | 1/PI ) = -0.694940081 + 0.442849491 i + 0.664274236 j + 0.885698980 k   in  R05 to 508

dn ( 1 + 2 i + 3 j + 4 k | 1/PI ) = -0.719248901 + 0.136199160 i + 0.204298741 j + 0.272398321 k    in R09 to R12

Example2:      m = PI   a = 1 + 2 i + 3 j + 4 k      ( Quaternion ->   4   STO 00 )

1  STO 01
2  STO 02
3  STO 03
4  STO 04

PI   XEQ "JEFA"   >>>>  1.004                    ---Execution time = 72s---

R01 =  0.860242322
R02 = -0.005873618
R03 = -0.008810427
R04 = -0.011747237

sn ( 1 + 2 i + 3 j + 4 k | PI ) = 0.860242322 - 0.005873618 i - 0.008810427 j - 0.011747237 k    and

cn ( 1 + 2 i + 3 j + 4 k | PI ) = 0.510825404 + 0.009891315 i + 0.014836973 j + 0.019782630 k   in  R05 to 508

dn ( 1 + 2 i + 3 j + 4 k | PI ) = -0.037125130 - 0.427571169 i - 0.645356753 j - 0.855142337 k    in R09 to R12

Example3:      m = -PI   q = 1 + 2 i + 3 j + 4 k      ( Quaternion ->   4   STO 00 )

1  STO 01
2  STO 02
3  STO 03
4  STO 04

PI  CHS   XEQ "JEFA"   >>>>  1.004                    ---Execution time = 81s---

R01 = -1.473447237
R02 = -0.079007894
R03 = -0.118511841
R04 = -0.158015787

sn ( 1 + 2 i + 3 j + 4 k | -PI ) = -1.473447237 - 0.079007894 i - 0.118511841 j - 0.158015787 k    and

cn ( 1 + 2 i + 3 j + 4 k | -PI ) = 0.285290836 - 0.408053637 i - 0.612080455 j - 0.816107274 k   in  R05 to 508

dn ( 1 + 2 i + 3 j + 4 k | -PI ) = -2.793322221 - 0.130928415 i - 0.196392622 j - 0.261856829 k    in R09 to R12

Notes:

-The subroutine "1/A" is called 3 times instead of once inside the loop LBL 03 ( lines 105 thru 189 )
-This avoids a SIZE = 5n+10 but it also increases the execution time...

b) Sn( a | m ) only

-If you just want to compute sn ( a | m ), the routine may of course be simplified:

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "SNA" )

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

Rn+1 ..........  R2n:  temp   R2n+1 = m ,  R2n+2 ..... temp

>>>  When the program stops:     R01   ......  Rnn = the n components of sn ( a | m )

Flags:  F01-F10
Subroutines:  "ST*A"  "ST/A"                            ( cf  "Anionic Special Functions(I) "paragraph 0 )
"A*A1"  "1/A"  "THA"  "SINA"    ( cf "Anions for the HP-41" )

-Line 26 is three-byte GTO 04

 01  LBL "SNA"   02  RCL 00   03  ST+ X   04  1   05  +   06  X<>Y   07  STO IND Y   08  RCL 00   09  .1   10  %   11  ISG X   12  +   13   E3   14  /   15  1   16  +   17  STO N   18  REGMOVE   19  CF 10   20  SIGN   21  X>Y? 22  GTO 01   23  X#Y?   24  GTO 00   25  XEQ "THA"   26  GTO 04   27  LBL 00   28  X<>Y   29  SQRT   30  XEQ "ST*A"   31  LASTX   32  1/X   33  ENTER^   34  SF 10   35  LBL 01   36  X<>Y   37  STO M   38  RCL 00   39  ST+ X   40  1   41  +   42  .1 43  %   44  +   45  STO O   46  SIGN   47  X<> M    48  LBL 02    49  ENTER^   50  CHS   51  1   52  +   53  SQRT   54  1   55  ST+ O   56  +   57  RCL X           58  2   59  /   60  ST* M   61  RDN   62  X^2   63  / 64  STO IND O   65  X^2   66  X#0?   67  GTO 02   68  RCL M   69  XEQ "ST*A"   70  XEQ "SINA"   71  LBL 03   72  RCL N   73  REGMOVE   74  XEQ "A*A1"   75  RCL IND O   76  XEQ "ST*A"   77  1   78  ST+ 01   79  XEQ "1/A"   80  XEQ "A*A1"   81  SIGN   82  RCL IND O   83  +   84  XEQ "ST*A" 85  DSE O   86  GTO 03   87  FC?C 10   88  GTO 04   89  RCL 00   90  ST+ X   91  1   92  +   93  RCL IND X   94  SQRT   95  XEQ "ST/A"   96  LBL 04   97  CF 01   98  RCL 00   99   E3 100  / 101  ISG X 102  CLA 103  END

( 197 bytes / SIZE 2n+10 )

 STACK INPUTS OUTPUTS X m 1.nnn

where   1.nnn   is the control number of  sn ( a | m )

Example:      m = 1/PI     a = 1 + 2 i + 3 j + 4 k      ( Quaternion ->   4   STO 00 )

1  STO 01
2  STO 02
3  STO 03
4  STO 04

PI  1/X   XEQ "JEFA"   >>>>  1.004                    ---Execution time = 26s---

R01 =  1.498262249
R02 =  0.205407205
R03 =  0.308110807
R04 =  0.410814409

sn ( 1 + 2 i + 3 j + 4 k | 1/PI ) = 1.498262249 + 0.205407205 i + 0.308110807 j + 0.410814409 k

-Similar results with m = PI & m = -PI

Notes:

-"SNA" runs of course faster than "JEFA"
-Moreover, since SIZE 2n+10 is sufficient, the argument may be a 128-on !

3°)  Catalan Numbers

-The Catalan numbers may be defined by  C(n) = 4n Gam(n+1/2) / [ sqrt(PI) Gam(n+2) ]
-This formula is used hereunder after replacing n by the anion a

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "CATA" )

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

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

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

Flag:  F04
Subroutines:  "ST*A"   "GAMA"         ( cf  "Anionic Special Functions(I)" )
"A*A1"  "X^A"  "1/A"    ( cf "Anions for the HP-41" )

 01  LBL "CATA"   02  RCL 00  03   E3  04  /  05  ISG X  06  RCL 00  07  4  08  *  09  +  10   E3  11  /  12  1  13  +  14  STO M 15  REGMOVE  16  4  17  XEQ "X^A"  18  FRC  19  RCL M  20  +  21  REGMOVE  22  LASTX  23  REGSWAP  24  REGMOVE  25  .5  26  ST+ 01  27  XEQ "GAMA"  28  RCL 00 29  .1  30  %  31  ISG X  32  +  33   E3  34  /  35  1  36  +  37  RCL 00  38  5  39  *  40  +  41  STO M  42  REGMOVE 43  XEQ "A*A1"  44  FRC  45  RCL 00  46  +  47  RCL M  48  X<>Y  49  -  50  REGSWAP  51  2  52  ST+ 01  53  XEQ "GAMA"  54  PI  55  SQRT  56  XEQ "ST*A" 57  XEQ "1/A"  58  RCL 00  59  +  60   E3  61  /  62  1  63  +  64  RCL 00  65  4  66  *  67  +  68  REGMOVE     69  XEQ "A*A1"  70  END

( 133 bytes / SIZE 6n+1 )

 STACK INPUTS OUTPUTS X m 1.nnn

where   1.nnn   is the control number of the result

Example:               a = 1 + i/2 + j/3 + k/4                            ( Quaternion ->   4   STO 00 )

1   STO 01
2   STO 02
3   STO 03
4   STO 04   XEQ "CATA"     >>>>    1.004                       ---Execution time = 54s---

R01 = 0.844036224
R02 = 0.239299403
R03 = 0.159532935
R04 = 0.119649702

-Thus,     C ( 1 + i/2 + j/3 + k/4 ) =  0.844036224 + 0.239299403 i + 0.159532935 j + 0.119649702 k

Note:

-Since SIZE = 6n+1 , the argument cannot be a 64-on.

4°)  Lambert W Function

-We have to solve  x = w(x) exp [ w(x) ]
-"LBWA" employs Newton's method.

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "LBWA" )

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

Rn+1 ..........  R4n:  temp

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

Flags:  /
Subroutines:  "ST*A"   "DSA"                          ( cf  "Anionic Special Functions(I)" )
"A*A1"  "1/A"  "LNA"  "E^A"    ( cf "Anions for the HP-41" )

 01  LBL "LBWA"  02  RCL 00  03   E3  04  /  05  ISG X  06   E3  07  /  08  RCL 00  09  3  10  *  11  +  12  1  13  +  14  STO M 15  REGSWAP  16  REGMOVE  17  SIGN  18  ST+ 01  19  XEQ "LNA"  20  RCL M  21  RCL 00  22  -  23  STO N  24  REGSWAP  25  LBL 01  26  RCL N  27  REGMOVE  28  VIEW 01 29  SIGN  30  CHS  31  XEQ "ST*A"  32  XEQ "E^A"  33  FRC  34  RCL M  35  +  36  REGMOVE  37  XEQ "A*A1"  38  3  39  RCL 00  40  ST* Y  41  LBL 02  42  RCL IND Y 43  ST- IND Y  44  RDN  45  DSE Y  46  DSE X  47  GTO 02  48  RCL N  49  RCL 00  50  -  51  REGSWAP  52  RCL N  53  REGMOVE  54  SIGN  55  ST+ 01  56  XEQ "1/A" 57  XEQ "A*A1"  58  XEQ "DSA"  59   E-8  60  X

( 143 bytes / SIZE 4n+1 )

 STACK INPUTS OUTPUTS X / 1.nnn

where   1.nnn   is the control number of the result

Example:                a = 1 + 2 i + 3 j + 4 k                            ( Quaternion ->   4   STO 00 )

1  STO 01
2  STO 02
3  STO 03
4  STO 04  XEQ "LBWA"   >>>>      1.004                       ---Execution time = 52s---

R01 =  1.281459811
R02 =  0.304051227
R03 =  0.456076840
R04 =  0.608102453

Whence,    W(1 + 2 i + 3 j + 4 k) =  1.281459811 + 0.304051227 i + 0.456076840 j + 0.608102453 k

Note:

-This routine doesn't work if a = -1  but   W(-1) =   -0.318131505 + 1.337235701 i

5°)  Weierstrass Elliptic Functions

a) Laurent Series

-The Weierstrass Elliptic Function  P(a;g2,g3) may be calculated by a Laurent series:

P(a;g2;g3) = a -2 + c2.a2 + c3.a4 + ...... + ck.a2k-2 + ....

where   c2 = g2/20  ;   c3 = g3/28  and   ck =  3 ( c2. ck-2 + c3. ck-3 + ....... + ck-2. c2 ) / (( 2k+1 )( k-3 ))      ( k > 3 )

-The successive  ck  are computed and stored into registers  R3n+1  R3n+2  ......

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "WEFA" )

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

Rn+1 ..........  R4n:  temp

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

Flags: /
Subroutines:    "A*A"  "1/A"   ( cf "Anions for the HP-41" )

-If you can use the M-Code routine  DS*A  listed in "Anionic Functions (I) paragrah  0°)b)

Replace lines 111 to 133  by  DS*A

-Likewise, replace lines 41 to 65 by   DS*A   XEQ "A*A1"   RCL N   1   +   RCL IND X   DS*A

 01  LBL "WEFA"   02  RCL 00   03  3   04  *   05  2   06  +   07  X<>Y   08  28   09  /   10  STO IND Y   11  X<> Z   12  20   13  /   14  DSE Y   15  STO IND Y   16  X<>Y   17  STO N   18  RCL 00   19  .1   20  %   21  ISG X   22  +   23   E3   24  /   25  1   26  + 27  STO O   28  REGMOVE   29  XEQ "A*A1"    30  RCL O   31  REGMOVE   32  XEQ "1/A"   33  FRC   34  RCL O   35  +   36  REGMOVE   37  LASTX   38  REGSWAP   39  REGMOVE   40  RCL IND N   41  STO M   42  XEQ 00   43  XEQ "A*A1"   44  RCL N   45  1   46  +   47  RCL IND X   48  STO M   49  XEQ 00   50  GTO 02   51  LBL 00   52  3 53  RCL 00   54  ST* Y   55  LBL 01   56  RCL IND X   57  RCL M   58  *   59  ST+ IND Z   60  RDN   61  DSE Y   62  DSE X   63  GTO 01   64  RTN   65  LBL 02   66  3   67  STO M   68  LBL 03   69  3   70  STO O    71  LBL 04   72  XEQ "A*A1"   73  SIGN   74  ST+ M   75  RCL 00   76  RCL 00   77   E3   78  / 79  +    80  X<>Y   81  -   82  3   83  *   84  RCL M   85  +   86  STO P   87  RCL N   88  ENTER^   89  CLX   90  LBL 05   91  RCL IND Y    92  RCL IND P    93  *   94  +   95  ISG Y   96  CLX   97  DSE P   98  GTO 05   99  RCL M 100  ST+ X 101  1 102  ST+ T 103  + 104  / 105  RCL M 106  3 107  ST* Z 108  - 109  / 110  STO IND Y 111  STO P 112  RCL 00 113  PI 114  INT 115  * 116  STO Q 117  RCL 00 118  ENTER^ 119  CLX 120  LBL 06 121  RCL IND Y  122  RCL P 123  * 124  RCL IND Q 125  + 126  STO IND Q 127  LASTX 128  - 129  ABS 130  + 131  DSE Q 132  DSE Y 133  GTO 06 134  X#0? 135  GTO 03 136  DSE O 137  GTO 04 138  RCL 00 139   E3 140  / 141  ISG X 142  STO Y 143  LASTX 144  / 145  1 146  + 147  RCL 00 148  ST+ X 149  + 150  REGMOVE  151  X<>Y 152  CLA 153  END

( 251 bytes / SIZE 3n+3+??? )

 STACK INPUTS OUTPUTS Y g2 / X g3 1.nnn

where   1.nnn   is the control number of the result

Example:         a = 0.2 + 0.3 i + 0.4 j + 0.5 k  ;   g2 = 0.9  ,   g3 = 1.4                            ( Quaternion ->   4   STO 00 )

0.2   STO 01
0.3   STO 02
0.4   STO 03
0.5   STO 05

0.9   ENTER^
1.4   XEQ "WEFA"   >>>>      1.004                       ---Execution time = 41s---    ( With  "DS*A" )

R01 =  -1.591637275
R02 =  -0.411614082
R03 =  -0.548818776
R04 =  -0.686023470

P( 0.2 + 0.3 i + 0.4 j + 0.5 k  ;  0.9 ,  1.4 ) = -1.591637275 - 0.411614082 i - 0.548818776 j - 0.686023470 k

b) Duplication Formula

-The argument q may be "large" and the Laurent series doesn't converge.
-The duplication formula may be used one or several times in this case.

-"WF2A" takes P(a) in registers R01 .......... Rnn  and returns P(2a) in the same registers

P(2a) = -2 P(a) + ( 6 P2(a) - g2/2 )2 / ( 4 ( 4 P3(a) - g2 P(a) - g3 ) )

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "WF2A" )

•  R01   ......  •  Rnn = the n components of P(a)

Rn+1 ..........  R3n:  temp

>>>  When the program stops:     R01   ......  Rnn = the n components of  P(2a)

Flags:  /
Subroutines:  "ST*A"                ( cf  "Anionic Special Functions(I)" )
"A*A1"  "1/A"     ( cf "Anions for the HP-41" )

 01  LBL "WF2A"  02  STO N  03  X<>Y  04  STO M  05  RCL 00  06  .1  07  %  08  STO Z  09  ISG X  10  +  11   E3  12  /  13  1  14  +  15  STO O  16  REGMOVE  17  + 18  REGMOVE  19  XEQ "A*A1"  20  6  21  XEQ "ST*A"  22  RCL M  23  2  24  /  25  ST- 01  26  RCL O  27  REGMOVE  28  XEQ "A*A1"  29  FRC  30  RCL O  31  +  32  REGSWAP  33  LASTX  34  REGMOVE 35  XEQ "A*A1"  36  XEQ "A*A1"  37  RCL 00  38  ENTER^  39  ST+ Y  40  LBL 01  41  RCL IND Y  42  RCL M  43  *  44  RCL IND Y  45  ST+ X  46  ST+ X  47  X<>Y  48  -  49  STO IND Y  50  RDN  51  DSE Y 52  DSE X  53  GTO 01  54  RCL N  55  ST- 01  56  4  57  XEQ "ST*A"  58  XEQ "1/A"  59  FRC  60  RCL O  61  +  62  RCL 00  63  +  64  REGSWAP  65  XEQ "A*A1"  66  3  67  RCL 00  68  ST* Y 69  LBL 02  70  RCL IND Y   71  ST+ X  72  ST- IND Y  73  RDN  74  DSE Y  75  DSE X  76  GTO 02  77  RCL M  78  RCL N  79  RCL 00  80   E3  81  /  82  ISG X  83  CLA  84  END

( 173 bytes / SIZE 3n+1 )

 STACK INPUTS OUTPUTS Y g2 / X g3 1.nnn

where   1.nnn   is the control number of the result

Example:     We have found  P( 0.2 + 0.3 i + 0.4 j + 0.5 k  ;  0.9 ,  1.4 ) = -1.591637275 - 0.411614082 i - 0.548818776 j - 0.686023470 k

-If the result is still in registers  R01 thru Rnn

0.9   ENTER^
1.4   XEQ "WF2A"   >>>>      1.004                       ---Execution time = 9s---

R01 = -0.370892103
R02 = -0.169841922
R03 = -0.226455893
R04 = -0.283069868

-So,   P( 0.4 + 0.6 i + 0.8 j + k  ;  0.9 ,  1.4 ) = -0.370892103 - 0.169841922 i - 0.226455893 j - 0.283069868 k

Note:

-In fact, the Laurent series still converges for 2a   and a direct evaluation of P(2a) with "WEFA"  yields ( in 2m01s )

P( 0.4 + 0.6 i + 0.8 j + k  ;  0.9 ,  1.4 ) = -0.370892103 - 0.169841920 i - 0.226455894 j - 0.283069867 k

6°)  Generalized Hypergeometric Functions

a)  General Case

-This program computes   pFq( b1,b2,....,bp ; c1,c2,....,cq ; a ) =  SUMk=0,1,2,.....    [(b1)k(b2)k.....(bp)k] / [(c1)k(c2)k.....(cq)k] . ak/k!         if   X > 0

where (bi)k = bi(bi+1)(bi+2) ...... (bi+k-1)   &  (bi)0 = 1    ,    likewise for  (cj)k    ( Pochhammer's symbol )

or the regularized function F tilde:

pF~q( b1,b2,....,bp ; c1,c2,....,cq ; a ) =  SUMk=0,1,2,.....    [ (b1)k(b2)k.....(bp)k ] / [Gam(k+c1) Gam(k+c2).....Gam(k+cq)] . ak/k!      if  X < 0

where Gam = Euler's Gamma function.

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn & R3n+1 ......  are to be initialized before executing "HGFA+" )

•  R01   ......  •  Rnn = the n components of the anion a

Rn+1 ..........  R3n:  temp         •  R3n+1 = b1 ......  •  R3n+p = bp   •  R3n+p+1 = c1  ...........  •  R3n+p+q = cq

>>>  When the program stops:     R01   ......  Rnn = the n components of  f(a)

Flags:  /
Subroutines:  "ST*A"  "DSA"                ( cf  "Anionic Special Functions(I)" )
"A*A1"  "A^X"                ( cf "Anions for the HP-41" )
"1/G"    M-Code routine  ( cf "Gamma Function for the HP-41" paragraph 1-f) )

 01  LBL "HGFA+"   02  CLA   03  STO M   04  RCL 00   05  PI   06  INT   07  10^X   08  /   09  ISG X   10  RCL X   11  RCL 00   12  +   13  PI   14  INT   15  10^X   16  /   17  ENTER^   18  SIGN   19  +   20  REGMOVE   21  RDN   22  CLRGX   23  SIGN   24  STO 01   25  X<>Y   26  ABS   27  INT   28  RCL X   29  PI   30  INT   31  10^X   32  /   33  STO O   34  CLX   35  RCL M   36  FRC   37  PI   38  INT   39  10^X   40  *   41  ABS 42  +   43  RCL X   44  PI   45  INT   46  10^X   47  /   48  RCL 00   49  ENTER^   50  SIGN   51  CHS   52  10^X   53  %   54  +   55  PI   56  INT   57  *   58  ST+ O   59  +   60  STO N   61  CLX   62  RCL O   63  +   64  RCL M   65  X>0?   66  GTO 05   67  CLX   68  SIGN   69  ENTER^   70  LBL 01   71  RCL IND Z     72  FRC   73  X#0?   74  GTO 01   75  X<> L   76  X>0?   77  GTO 01   78  XY   80  GTO 00   81  LBL 01   82  CLX 83  RCL IND T   84  1/G   85  ST* Z   86  LBL 00   87  RDN   88  DSE Z   89  GTO 01   90  X>0?   91  GTO 00   92  1   93  X<>Y   94  -   95  STO P   96  LBL 02   97  ST/ Y   98  ENTER^   99  SIGN 100  - 101  X<>Y 102  ISG N 103  LBL 03 104  RCL Y 105  RCL IND N   106  + 107  ISG O 108  GTO 04 109  STO L 110  X>0? 111  GTO 03 112  FRC 113  X#0? 114  GTO 03 115  CLX 116  SIGN 117  SIGN 118  LBL 03 119  X<> L 120  1/X 121  LBL 04 122  * 123  ISG N 124  GTO 03 125  RCL 00 126  PI 127  INT 128  * 129  RCL N 130  FRC 131  X<>Y 132  + 133  STO N 134  X<> O 135  LASTX 136  X<>Y 137  FRC 138  + 139  STO O 140  X<> Z 141  X#0? 142  GTO 02 143  LBL 00 144  X<>Y 145  RCL 00 146  ST+ X 147  ENTER^ 148  SIGN 149  + 150  X<>Y 151  STO IND Y   152  RCL 00 153  PI 154  INT 155  10^X 156  / 157  ISG X 158  LASTX 159  / 160  ENTER^ 161  SIGN 162  + 163  RCL 00 164  + 165  REGMOVE 166  RCL P 167  XEQ "A^X" 168  RCL 00 169  ST+ X 170  ENTER^ 171  SIGN 172  + 173  RCL IND X  174  XEQ "ST*A"  175  LBL 05 176  RCL 00 177  PI 178  INT 179  10^X 180  / 181  ISG X 182  RCL 00 183  ST+ X 184  + 185  PI 186  INT 187  10^X 188  / 189  ENTER^ 190  SIGN 191  + 192  REGMOVE   193  LBL 06 194  XEQ "A*A1" 195  SIGN 196  ST+ N 197  LBL 07 198  RCL IND N 199  RCL P 200  + 201  ISG O 202  FS? 30 203  1/X 204  * 205  ISG N 206  GTO 07 207  RCL N 208  FRC 209  RCL 00 210  PI 211  INT 212  * 213  + 214  STO N 215  X<> O 216  LASTX 217  X<>Y 218  FRC 219  + 220  STO O 221  SIGN 222  RCL P 223  + 224  STO P 225  / 226  XEQ "ST*A"  227  XEQ "DSA" 228  X#0? 229  GTO 06 230  RCL M 231  RCL 00 232   E3 233  / 234  ISG X 235  RCL X 236  LASTX 237  / 238  1 239  + 240  RCL 00 241  ST+ X 242  + 243  REGMOVE 244  RDN 245  CLA 246  END

( 354 bytes / SIZE 3n+1+p+q )

 STACK INPUTS OUTPUTS Y / +/- p.qqq X +/- p.qqq 1.nnn

where   1.nnn   is the control number of the result.

+ p.qqq   for the hypergeometric function
- p.qqq   for the regularized hypergeometric function

Example1:             b1 = PI  ,   b2 = e  ;  c1 = 1 ,  c2 = 2 ,   c3 = 4  ;   a = 1 + 2 i + 3 j + 4 k              ( Quaternion ->   4   STO 00 )

-Since n = 4 , the parameters must be stored into  R13  R14  R15  ................

PI   STO 13    1  E^X  STO 14  ,   1  STO 15   2  STO 16   4  STO 17

•   1  STO 01   2  STO 02   3  STO 03   4  STO 04

2.003   XEQ "HGFA+"  >>>>      1.004                       ---Execution Time= 45s---

R01 = -6.691126901
R02 =  1.302530584
R03 =  1.953795876
R04 =  2.605061166

2F3( PI , e ; 1 , 2 , 4 ;  1 + 2 i + 3 j + 4 k ) = -6.691126901 + 1.302530584 i + 1.953795876 j + 2.605061166 k

•   1  STO 01   2  STO 02   3  STO 03   4  STO 04

2.003  CHS   XEQ "HGFA+"  >>>>      1.004                       ---Execution Time= 53s---

R01 = -1.115187818
R02 =  0.217088431
R03 =  0.325632646
R04 =  0.434176861

2F~3( PI , e ; 1 , 2 , 4 ;  1 + 2 i + 3 j + 4 k ) = -1.115187818 + 0.217088431 i + 0.325632646 j + 0.434176861 k

-The first result has been simply divided by  Gam(1) Gam(2) Gam(4) = 6

Example2:             b1 = PI  ,   b2 = e  ;  c1 = 1 ,  c2 = -2 ,   c3 = -4  ;   a = 1 + 2 i + 3 j + 4 k

PI   STO 13    1  E^X  STO 14  ,   1  STO 15   -2  STO 16   -4  STO 17

( the non-regularized hypergeometric function does not exist for these parameters )

•    1  STO 01   2  STO 02   3  STO 03   4  STO 04

2.003  CHS   XEQ "HGFA+"  >>>>      1.004                       ---Execution Time= 67s---

R01 = -2910223.711
R02 =  192140.2259
R03 =  288210.3384
R04 =  384280.4514

2F~3( PI , e ; 1 , -2 , -4 ;  1 + 2 i + 3 j + 4 k ) = -2910223.711 + 192140.2259 i + 288210.3384 j + 384280.4514 k

b)  p + q < 4

-If p+q < 4 , the parameters bi & cj may be entered in the stack without using registers R3n+1 ... and so on ...
-However, it requires several M-Code routines + the 2 following ones:  RCLIX  and   @HG

RCLIX  takes an integer m in X-register and replaces it by the content of

M-register if m = 1
N-register if m = 2
O-register if m = 3

@HG  corresponds to the loop  LBL 06  in the previous listing.

-It takes  q in register Y and p+q in register X  and also the required data in registers R00 thru R3n.

-But first  RCLIX:

098  "X"
009  "I"
00C  "L"
003  "C"
012  "R"
0F8  C=X
38D  C->
008  S&X
106  A=C S&X
130  LDI S&X
004  004
146  A=A+C S&X
130   LDI S&X
009   009
306   ?A<C S&X
0B5   ?NCGO
0A2   ERROR
0A6   A<>C S&X
270   RAMSLCT
0E8   X=C
3E0   RTN

Note:

"RCLIX" gives a DATA ERROR message if  X > 4

-Then  @HG

087  "G"
008  "H"
000  "@"
0B8  C=Y
38D  C->
008   S&X
070   N=C ALL
0F8   C=X
38D  C->
008   S&X
106  A=C S&X
0B0  C=N ALL
1BC  RCR 11
0A6  A<>C S&X
10E  A=C ALL
130  LDI S&X
004  004
306  ?A<C S&X
0B5   ?NCGO                                 gives a DATA ERROR message if  p+q > 3
0A2   ERROR
0E6  B<>C S&X
0C6  C=B S&X
1BC  RCR 11
0C6  C=B S&X
20E  C=A+C ALL
1BC  RCR 11
0C6  C=B S&X
226   C=C+1 S&X
03C  RCR 3
268  Q=C  =  addresses  q , p+q , 005
3C8  CLRKEY
119  ?NCXQ                                            LOOP1
384   E146                 E146 = the address of the 1st executable word of "A*A1" in my ROM -> Change these 2 words according to your own ROM
04E  C
35C =                         I've modified the M-Code routine A*A1 so that it does not use synthetic register Q
050  1
0E8  X=C
278  C=Q
070  N=C ALL
0B0  C=N ALL                                        LOOP2
03C  RCR 3
106  A=C S&X
0B0  C=N ALL
244   CLRF 9
306  ?A<C S&X
013  JNC+02
248  SETF 9
270  RAMSLCT
10E  A=C ALL
238  C=P
01D  C=
060  A+C
24C  ?FSET 9
239   ?NCXQ
060   1/AB
0F8   C=X
13D   C=
060   AB*C
0E8  X=C
0B0  C=N ALL
266   C=C-1 S&X
070  N=C ALL
106  A=C S&X
1BC RCR 11
306  ?A<C S&X
32B  JNC -27d                                      GOTO  LOOP2
238  C=P
02E  C
0FA ->
0AE  AB
001   C=
060   AB+1
228  P=C
10E  A=C ALL
0F8  C=X
0AE  A<>C ALL
261   C=
060   A/C
0E8  X=C
3D5  ?NCXQ
388   E2F5 = the 1st executable word of "ST*A" in my ROM.    -> Change these 2 words according to your own ROM
046  C=0 S&X
270  RAMSLCT
2F1  ?NCXQ
384  E1BC = the 1st executable word of "DSA" in my ROM.    -> Change these 2 words according to your own ROM
0F8  C=X
3CC  ?KEY
360   ?C RTN             Press a key to stop an infinite loop or if the routine is too slow.
2EE  ?C#0 ALL
227   JC-60d                                        GOTO  LOOP1
3C1
002   ( END )

( 93 words )

-Only use @HG  with the program below !

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "HGFB" )

•  R01   ......  •  Rnn = the n components of the anion a  which are saved in  Rn+1 .... R2n

Rn+1 ..........  R3n:  temp

>>>  When the program stops:     R01   ......  Rnn = the n components of  f(a)

Flags:  /
Subroutines:  "ST*A"                                                       ( cf  "Anionic Special Functions(I)" )
"A^X"                                                       ( cf "Anions for the HP-41" )
"1/G"                                                         ( cf "Gamma Function for the HP-41" paragraph 1-f) )
"ST<>A"  "X+1"  "X-1"  "X*E3"  "X/E3"  ( cf "A few M-Code routines for the HP-41" )
and  RCLIX & @HG which are listed above

-Line 39 is a three-byte  GTO 04

 01  LBL "HGFB"   02  CLA   03  RDN   04  ST<>A   05  R^   06  RCL 00   07  X/E3   08  ISG X   09  RCL X   10  RCL 00   11  +   12  X/E3   13  X+1   14  REGMOVE   15  RDN   16  CLRGX   17  SIGN   18  STO 01   19  X<>Y   20  CF 00   21  X<0?   22  SF 00   23  ABS   24  FRC   25  STO Y   26  LASTX   27  INT   28  X/E3   29  +   30  RCL 00   31  X+1 32  ST+ X   33  RDN   34  STO IND T   35  DSE T   36  X<>Y   37  STO IND T    38  FC?C 00   39  GTO 04   40  X*E3   41  PI   42  SIGN   43  ENTER^   44  LBL 01   45  RCL Z   46  RCLIX   47  FRC   48  X#0?   49  GTO 01   50  CLX   51  LASTX   52  X>0?   53  GTO 01   54  XY   56  GTO 00   57  LBL 01   58  X<> L   59  1/G   60  ST* Z   61  LBL 00   62  RDN 63  DSE Z   64  GTO 01   65  X>0?   66  GTO 00   67  CHS   68  X+1   69  STO P   70  RCL 00   71  X+1   72  ST+ X   73  STO Q   74  RDN   75  LBL 02   76  ST/ Y   77  X-1   78  X<>Y   79  ISG IND Q   80  LBL 03   81  RCL Y   82  RCL IND Q    83  RCLIX   84  +   85  DSE Q   86  ISG IND Q   87  GTO 02   88  *   89  GTO 03   90  LBL 02   91  STO L   92  X>0?   93  GTO 02 94  FRC   95  X#0?   96  GTO 02   97  CLX   98  SIGN   99  SIGN 100  LBL 02 101  X<> L 102  / 103  LBL 03 104  ISG Q 105  CLX 106  ISG IND Q 107  GTO 03 108  RCL IND Q 109  FRC 110  STO IND Q 111  DSE Q 112  X<> IND Q 113  FRC 114  STO IND Q  115  SIGN 116  ST+ Q 117  X<> Z 118  X#0? 119  GTO 02 120  LBL 00 121  X<>Y 122  RCL 00 123  ST+ X 124  X+1 125  X<>Y 126  X<> IND Y 127  X*E3 128  ISG Y 129  CLX 130  ST+ IND Y 131  RCL 00 132  X/E3 133  ISG X 134  X/E3 135  X+1 136  RCL 00 137  + 138  REGMOVE 139  RCL P 140  XEQ "A^X" 141  RCL 00 142  ST+ X 143  X+1 144  RCL IND X  145  ST*A 146  CLX 147  SIGN 148  + 149  RCL IND X 150  FRC 151  LASTX 152  INT 153  X/E3 154  LBL 04 155  RCL 00 156  X/E3 157  ISG X 158  RCL 00 159  ST+ X 160  + 161  X/E3 162  X+1 163  REGMOVE 164  RDN 165  X*E3 166  X<>Y 167  X*E3 168  @HG 169  RCL 00 170  X/E3 171  ISG X 172  ENTER^ 173  X/E3 174  X+1 175  RCL 00 176  ST+ X 177  + 178  REGMOVE  179  R^ 180  R^ 181  ST<>A 182  R^ 183  CLA 184  END

( 293 bytes /SIZE 3n+1 )

 STACK INPUTS OUTPUTS T bi or cj bi or cj Z bi or cj bi or cj Y bi or cj bi or cj X +/- p.qqq 1.nnn

where   1.nnn   is the control number of the result.

+ p.qqq   for the hypergeometric function
- p.qqq   for the regularized hypergeometric function

>>> Enter the bi first, then the cj and finally  +/- p.qqq

Example1:       2F1       a = 0.1 + 0.2 i + 0.3 j + 0.4 k           p = 1.1   q = 1.2   s = 1.3              ( Quaternion ->   4   STO 00 )

0.1  STO 01
0.2  STO 02
0.3  STO 03
0.4  STO 04

1.1  ENTER^
1.2  ENTER^
1.3  ENTER^
2.001  XEQ "HGFB"   >>>>   1.004                           ---Execution time = 42s---

R01 =  0.814236592
R02 =  0.184421233
R03 =  0.276631850
R04 =  0.368842467

2F1 ( 1.1 , 1.2 ; 1.3 ;  0.1 + 0.2 i + 0.3 j + 0.4 k ) = 0.814236592 + 0.184421233 i + 0.276631850 j + 0.368842467 k

Example2:      2F~1        a = 0.1 + 0.2 i + 0.3 j + 0.4 k           p = 1   q = 2   s = -4                        ( Quaternion ->   4   STO 00 )

0.1  STO 01
0.2  STO 02
0.3  STO 03
0.4  STO 04

1     ENTER^
2     ENTER^
-4    ENTER^
2.001  CHS   XEQ "HGFB"   >>>>   1.004                          ---Execution time = 92s---

R01 =  -7.152496297
R02 =  -9.061272756
R03 =  -13.59190873
R04 =  -18.12254510

2F~1 ( 1 , 2 ; -4 ;  0.1 + 0.2 i + 0.3 j + 0.4 k ) = -7.152496297 - 9.061272756 i - 13.59190873 j - 18.12254510 k

Example3:     1F2         a = 0.1 + 0.2 i + 0.3 j + 0.4 k           p = 1.1   q = 1.2   s = 1.3              ( Quaternion ->   4   STO 00 )

0.1  STO 01
0.2  STO 02
0.3  STO 03
0.4  STO 04

1.1  ENTER^
1.2  ENTER^
1.3  ENTER^
1.002  XEQ "HGFB"   >>>>   1.004                           ---Execution time = 10s---

R01 =  1.028367021
R02 =  0.146116085
R03 =  0.219174127
R04 =  0.292232170

-Whence,       1F2 ( 1.1 ; 1.2 , 1.3 ;  0.1 + 0.2 i + 0.3 j + 0.4 k ) = 1.028367021 + 0.146116085 i + 0.219174127 j + 0.292232170 k

Notes:

-Likewise, for the Kummer's Function 1F1  ,  place  1.001  in X-register
and  -1.001  in X-register for the regularized Kummer's Function  1F~1  and so on ....

c)  0F1

-It's also possible to write programs that don't use M-Code routines to compute all kinds of hypergeometric functions:
-For example, the following one, 0F1 , is useful to calculate Airy's functions and Lommel functions:

Data Registers:           •  R00 = n > 1                          ( Registers R00 thru Rnn are to be initialized before executing "0F1A" )

•  R01   ......  •  Rnn = the n components of the anion a  which are saved in  Rn+1 .... R2n

Rn+1 ..........  R3n:  temp

>>>  When the program stops:     R01   ......  Rnn = the n components of  f(a)

Flags:  /
Subroutines:   "ST/A"  "DSA"           ( cf  "Anionic Special Functions(I)" )
"A*A1"                        ( cf "Anions for the HP-41" )

 01  LBL "0F1A"   02  CLA  03  STO M  04  RCL 00  05   E3  06  /  07  ISG X  08  STO O  09  RCL 00  10  +  11   E3 12  /  13  1  14  +  15  REGMOVE   16  RCL O  17  CLRGX  18  SIGN  19  STO 01  20  LASTX  21  RCL 00  22  ST+ X 23  +  24   E3  25  /  26  +  27  REGMOVE  28  LBL 01  29  XEQ "A*A1"  30  RCL M  31  RCL N  32  ST+ Y  33  ISG X 34  CLX  35  STO N  36  *  37  XEQ "ST/A"  38  XEQ "DSA"  39  X#0?  40  GTO 01  41  SIGN  42  RCL O  43   E3  44  / 45  +  46  RCL 00  47  ST+ X  48  +  49  REGMOVE   50  RCL O  51  CLA  52  END

( 97 bytes / SIZE 3n+1 )

 STACK INPUTS OUTPUTS X m 1.nnn

where   1.nnn   is the control number of  the result

Example:    Calculate  0F1( m ; a )   with   m = PI   a = 1 + 2 i + 3 j + 4 k      ( Quaternion ->   4   STO 00 )

1  STO 01
2  STO 02
3  STO 03
4  STO 04

PI  XEQ "0F1A"   >>>>  1.004                    ---Execution time = 15s---              ( with M-Code routines... )

R01 = 0.106086113
R02 = 0.641727862
R03 = 0.962591792
R04 = 1.283455724

-Whence,  0F1( PI ; 1 + 2 i + 3 j + 4 k ) = 0.106086113 + 0.641727862 i + 0.962591792 j + 1.283455724 k

Notes:

-To get the regularized functions, simply divide the result by the product of the Gamma(cj)
....  but it will not work if a cj is a non-positive integer !
-Lommel functions also call "0F1A" but it is preferable to save registers  Y & Z  in this case ( cf paragraph 11 below )

-For "A*A1" - line 34 - use a version that does not alter synthetic register P

 01  LBL "0F1A"  02  STO M  03  RDN  04  STO N  05  X<>Y  06  STO O  07  RCL 00  08  RCL 00  09   E3  10  /  11  ISG X  12  STO Z  13  + 14   E3  15  /  16  1  17  +  18  REGMOVE   19  X<>Y  20  CLRGX  21  SIGN  22  STO 01  23  LASTX  24  RCL 00  25  ST+ X  26  + 27   E3  28  /  29  +  30  REGMOVE  31  CLX  32  STO P  33  LBL 01  34  XEQ "A*A1"  35  RCL M  36  RCL P  37  ST+ Y  38  ISG X  39  CLX 40  STO P  41  *  42  XEQ "ST/A"  43  XEQ "DSA"  44  X#0?  45  GTO 01  46  RCL 00  47   E3  48  /  49  ISG X  50  STO Y  51  LASTX  52  / 53  1  54  +  55  RCL 00  56  ST+ X  57  +  58  REGMOVE   59  X<> O  60  RCL N  61  RCL M  62  R^  63  CLA  64  END

( 114 bytes / SIZE 3n+1 )

 STACK INPUTS OUTPUTS T / Z Z Z Y Y Y m X m 1.nnn

where   1.nnn   is the control number of  the result

-The same example gives the same result.
-Use this variant to compute Lommel functions of the second kind in paragraph 11.

d)  2F0

-"2F0A" computes the asymptotic function   2F0 (p,q;;a)
-It is used to calculate an asymptotic expansion of the exponential integral  Em(a)

Data Registers:           •  R00 = n > 1                          ( Registers R00 thru Rnn are to be initialized before executing "2F0A" )

•  R01   ......  •  Rnn = the n components of the anion a  which are saved in  Rn+1 .... R2n

Rn+1 ..........  R3n:  temp

>>>  When the program stops:     R01   ......  Rnn = the n components of  f(a)

Flags:  /
Subroutines:   "ST*A"  "DSA"           ( cf  "Anionic Special Functions(I)" )
"A*A1"                        ( cf "Anions for the HP-41" )

 01  LBL "2F0A"   02  STO M  03  X<>Y  04  STO N  05  RCL 00  06  .1  07  %  08  STO Z  09  ISG X  10  +  11   E3  12  / 13  1  14  +  15  REGMOVE  16  +  17  STO O  18  CLX  19  XEQ "ST*A"  20  SIGN  21  STO 01  22  CLX  23  X<> O  24  REGMOVE 25  LBL 01  26  XEQ "A*A1"  27  RCL M  28  RCL N  29  RCL O  30  ST+ Y  31  ST+ Z  32  ISG X  33  CLX  34  STO O  35  /  36  * 37  XEQ "ST*A"  38  XEQ "DSA"  39  X#0?  40  GTO 01  41  RCL N  42  RCL M  43  RCL 00  44   E3  45  /  46  ISG X  47  STO O  48  LASTX 49  /  50  1  51  +  52  RCL 00  53  ST+ X  54  +  55  REGMOVE   56  X<> O  57  CLA  58  END

( 112 bytes / SIZE 3n+1 )

 STACK INPUTS OUTPUTS Z / p Y p q X q 1.nnn

where   1.nnn   is the control number of  the result

Notes:

-See "Anionic Special Functions (I)" paragraph 4°)c) for an application of this routine.
- XEQ "2F0A"  may be replaced by  2  XEQ "HGFB"

7°)  Airy Functions

Formulae:

Ai(a) =   f(a) - g(a)                             with            f(a) = [ 3 -2/3 / Gamma(2/3) ]  0F1( 2/3 ; a3/9 )
Bi(a) = [ f(a) + g(a) ] sqrt(3)               and            g(a) = [ 3 -1/3 / Gamma(1/3) ]  0F1( 4/3 ; a3/9 )   a

Data Registers:           •  R00 = n > 1                          ( Registers R00 thru Rnn are to be initialized before executing "AIRYA" )

•  R01   ......  •  Rnn = the n components of the anion a

Rn+1 ..........  R4n:  temp

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

Flags:  /
Subroutines:  "ST/A"  "ST*A"           ( cf  "Anionic Special Functions(I)" )
"A*A1"                       ( cf "Anions for the HP-41" )
"0F1A"                       ( cf paragraph  6°c) above )

 01  LBL "AIRYA"  02  RCL 00  03  .1  04  %  05  STO Z  06  ST+ Z  07  ISG X  08  +  09   E3  10  /  11  1  12  +  13  REGMOVE  14  + 15  REGMOVE  16  XEQ "A*A1"  17  XEQ "A*A1"  18  9  19  XEQ "ST/A"  20  4  21  3  22  /  23  XEQ "0F1A"   24  RCL 00  25  3  26  *  27  +  28   E3 29  /  30  1  31  +  32  STO M  33  RCL 00  34  +  35  REGSWAP  36  XEQ "A*A1"  37  .2588194038  38  XEQ "ST*A"  39  RCL M  40  REGSWAP  41  2  42  3 43  /  44  XEQ "0F1A"  45 .3550280539  46  XEQ "ST*A"  47  RCL 00  48  STO M  49  ST+ X  50  ENTER^  51  ST+ Y  52  LBL 01  53  RCL IND M  54  RCL IND Z  55  ST- IND M  56  + 57  3  58  SQRT  59  *  60  STO IND Y   61  RDN  62  DSE Y  63  DSE X  64  DSE M  65  GTO 01  66  RCL 00  67   E3  68  /  69  ISG X  70  END

( 167 bytes / SIZE 4n+1 )

 STACK INPUTS OUTPUTS X m 1.nnn

where   1.nnn   is the control number of  Ai(a)     [  Bi(a)  is in registers  Rn+1 ..... R2n ]

Example:             a = 1 + i/2 + j/3 + k/4                                                  ( Quaternion ->   4   STO 00 )

1    STO 01
2   1/X   STO 02
3   1/X   STO 03
4   1/X   STO 04    XEQ "AIRYA"    >>>>  1.004                                      ---Execution time = 27s---

R01 =  0.105299445
R02 = -0.078442074
R03 = -0.052294716
R04 = -0.039221037

Thus,    Ai ( 1 + i/2 + j/3 + k/4 ) = 0.105299445 - 0.078442074 i - 0.052294716 j - 0.039221037 k

and     Bi ( 1 + i/2 + j/3 + k/4 ) = 0.973463977 + 0.394832415 i + 0.263221610 j + 0.197416207 k    in   R05 to R08

Notes:

-This program returns accurate results for "small" arguments only.
-Lines 23 & 44 may be replaced by  E-3  XEQ "HGFB"

8°)  Hermite Functions

Formula:       Hm(a) = 2m sqrt(PI) [ (1/Gam((1-m)/2))  M(-m/2,1/2,a2) - ( 2.a / Gam(-m/2) ) M((1-m)/2,3/2,a2) ]

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

Data Registers:           •  R00 = n > 1         ( Registers R00 thru Rnn are to be initialized before executing "HMTA" )

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

Rn+1 ..........  R4n:  temp   R4n+1 = -m/2

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

Flags:  /
Subroutines:  "ST*A"   "KUMA"  ( cf  "Anionic Special Functions(I) )
"A*A1"                      ( cf "Anions for the HP-41" )
"1/G+"                 ( cf "Gamma Function for the HP-41" )

 01  LBL "HMTA"   02  CHS   03  2   04  /   05  STO M   06  RCL 00   07  4   08  *   09  1   10  +   11  X<>Y   12  STO IND Y    13  RCL 00   14  .1   15  %   16  STO Z   17  ST+ Z   18  ISG X   19  +   20   E3   21  / 22  1   23  +   24  REGMOVE   25  +   26  REGMOVE   27  XEQ "A*A1"   28  RCL M   29  .5   30  ST+ Y   31  3   32  *   33  XEQ "KUMA"   34  SIGN   35  RCL 00   36  4   37  *   38  +   39  RCL IND X   40  STO N   41  XEQ "1/G+"   42  ST+ X 43  XEQ "ST*A"   44  RCL 00   45   E3   46  /   47  ISG X   48  RCL 00   49  3   50  *   51  +   52   E3   53  /   54  RCL 00   55  +   56  1   57  +   58  REGSWAP    59  XEQ "A*A1"    60   E3   61  /   62  RCL 00   63  3 64  *   65  +   66  1   67  +   68  REGSWAP   69  RCL N   70  .5   71  XEQ "KUMA"   72  X<> Z   73  STO N   74  +   75  XEQ "1/G+"   76  XEQ "ST*A"   77  4   78  RCL 00   79  ST* Y   80  LBL 01   81  RCL IND Y   82  ST- IND Y   83  RDN   84  DSE Y 85  DSE X   86   GTO 01   87  2    88  RCL N   89  ST+ X   90  CHS   91  Y^X   92  PI   93  SQRT   94  *   95  XEQ "ST*A"    96  RCL 00   97   E3   98  /   99  ISG X 100  CLA 101  END

( 190 bytes / SIZE 4n+2 )

 STACK INPUTS OUTPUTS X m 1.nnn

where   1.nnn   is the control number of the result

Example:             a = 1 + i/2 + j/3 + k/4       m = PI                                           ( Quaternion ->   4   STO 00 )

1    STO 01
2   1/X   STO 02
3   1/X   STO 03
4   1/X   STO 04

PI   XEQ "HMTA"    >>>>  1.004                                      ---Execution time = 45s---

R01 = -17.81760490
R02 =   2.845496135
R03 =   1.896997426
R04 =   1.422748067

-Whence,     Hpi( 1 + i/2 + j/3 + k/4 ) =  -17.81760490 + 2.845496135 i + 1.896997426 j + 1.422748067 k

Note:

-Lines 71 & 33 may be replaced by   1.001  XEQ "HGFB"

9°)  Chebyshev Functions

Formulae:

CF 02          Tm(cos A) = cos m.A                          ( first kind )                           with       cos A = a
SF 02          Um(cos A) = [ sin (m+1).A ] / sin A   ( second kind )

Data Registers:           •  R00 = n > 1                          ( Registers R00 thru Rnn are to be initialized before executing "CHBA" )

•  R01   ......  •  Rnn = the n components of the anion a

Rn+1 ..........  R2n:  temp   R2n+1 = m

>>>  When the program stops:     R01   ......  Rnn = the n components of  f(a)

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

Subroutines:   "ST*A"                                                              ( cf  "Anionic Special Functions(I)" )
"A*A1"  "ACOSA"  "COSA"  "SINA"  "1/A"           ( cf "Anions for the HP-41" )

 01  LBL "CHBA"  02  RCL 00  03  ST+ X  04  1  05  +  06  X<>Y  07  STO IND Y  08  XEQ "ACOSA"  09  RCL 00 10  STO Z  11  ST+ Z  12  +  13   E3  14  /  15  1  16  +  17  FS? 02  18  REGMOVE 19  SIGN  20  +  21  RCL IND X  22  LASTX  23  FC? 02  24  CLX  25  +  26  XEQ "ST*A"    27  FS? 02 28  GTO 00  29  XEQ "COSA"    30  RTN  31  LBL 00  32  XEQ "SINA"   33  RCL 00  34  +  35   E3  36  / 37  1  38  +  39  REGSWAP  40  XEQ "SINA"    41  XEQ "1/A"  42  XEQ "A*A1"   43  END

( 100 bytes / SIZE 2n+2 )

 STACK INPUTS OUTPUTS X m 1.nnn

where   1.nnn   is the control number of the result

Example:             a = 0.1 + 0.2 i + 0.3 j + 0.4 k   m = PI                                           ( Quaternion ->   4   STO 00 )

•   CF 02      Chebyshev Function 1st kind

0.1   STO 01
0.2   STO 02
0.3   STO 03
0.4   STO 04

PI  XEQ "CHBA"   >>>>    1.004                                      ---Execution time = 15s---

R01 = -0.143159414
R02 = -0.905090619
R03 = -1.357635928
R04 = -1.810181237

TPI( 0.1 + 0.2 i + 0.3 j + 0.4 k ) =  -0.143159414 - 0.905090619 i - 1.357635928 j - 1.810181237 k

•   SF 02      Chebyshev Function 2nd kind

0.1   STO 01
0.2   STO 02
0.3   STO 03
0.4   STO 04

PI  XEQ "CHBA"   >>>>    1.004                                      ---Execution time = 20s---

R01 = -0.386199512
R02 = -1.369695948
R03 = -2.054543921
R04 = -2.739391894

And    UPI( 0.1 + 0.2 i + 0.3 j + 0.4 k ) =  -0.386199512 - 1.369695948 i - 2.054543921 j - 2.739391894 k

Note:

-Here is a variant that employs the M-Code routines  AMOVE & ASWAP  listed in "Anionic Special Function(I)"
and the M-Code routines  X+1   STO W  &  RCL W  listed in "A few M-Code routines for the HP-41"

 01  LBL "CHBA"  02  STO W  03  ACOSA  04  12  05  FS? 02  06  AMOVE  07  RCL W  08  FS? 02  09  X+1  10  ST*A  11  FS? 02  12  GTO 00  13  COSA  14  RTN  15  LBL 00  16  SINA  17  12  18  ASWAP  19  SINA  20  1/A  21  A*A1  22  END

( 49 bytes / SIZE 2n+1 )

-Assuming the focal programs "ACOSA" "COSA" "SINA" "1/A"  are in a HEPAX or a NOVRAM module,  we've saved 51 bytes !
-The results are of course identical.

10°) Riemann Zeta Function

a) If  Re(a) > 1

-"ZETAA" uses  Zeta(a) = 1 + 1/2a + 1/3a + ........ + 1/na + ......

Data Registers:           •  R00 = n > 1                          ( Registers R00 thru Rnn are to be initialized before executing "ZETAA" )

•  R01   ......  •  Rnn = the n components of the anion a

Rn+1 ..........  R3n:  temp

>>>  When the program stops:     R01   ......  Rnn = the n components of  f(a)

Flags:  /
Subroutines:  "DSA"      ( cf  "Anionic Special Functions(I)" )
"X^A"          ( cf "Anions for the HP-41" )

 01  LBL "ZETAA"  02  CLA  03  RCL 00  04  .1  05  %  06  ISG X  07  +  08   E3  09  /  10  1  11  +  12  REGMOVE 13  RCL 00  14  3  15  *  16  .1  17  %  18  ISG X  19  +  20  RCL 00           21  -  22  CLRGX  23  LBL 01  24  RCL 00 25   E3  26  /  27  ISG X  28  LASTX  29  /  30  1  31  +  32  RCL 00  33  +  34  REGMOVE     35  SIGN  36  RCL M 37  +  38  STO M  39  1/X  40  XEQ "X^A"  41  XEQ "DSA"     42  X#0?  43  GTO 01  44  RCL 00  45   E3  46  /  47  ISG X  48  LASTX 49  /  50  1  51  +  52  RCL 00  53  ST+ X  54  +  55  REGMOVE     56  FRC  57   E3  58  *  59  CLA  60  END

( 96 bytes / SIZE 3n+1 )

 STACK INPUT OUTPUT X / 1.nnn

where   1.nnn   is the control number of the result

Example:                a = 6 + 7 i + 8 j + 9 k                          ( Quaternion ->   4   STO 00 )

6  STO 01
7  STO 02
8  STO 03
9  STO 04

XEQ "ZETAA"      >>>>    1.004                                      ---Execution time = 2m53s---

R01 = 0.983702003
R02 = 0.001473512
R03 = 0.001684014
R04 = 0.001894515

Zeta ( 6 + 7 i + 8 j + 9 k ) = 0.983702003 + 0.001473512 i + 0.001684014 j + 0.001894515 k

Notes:

-Add  RND  after line 41 and the accuracy will be controlled by the display settings.
-If Re(a) is close to 1, the execution time becomes prohibitive and a better method must be used:

b) Borwein Algorithm

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

Data Registers:           •  R00 = n > 1                          ( Registers R00 thru Rnn are to be initialized before executing "ZETAA" )

•  R01   ......  •  Rnn = the n components of the anion a

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

>>>  When the program stops:     R01   ......  Rnn = the n components of  f(a)

Flags:  /
Subroutines:   "ST*A"  "ST/A"  "GAMA"  "DSA"                          ( cf  "Anionic Special Functions(I)" )
"A*A1"  "X^A"   "1/A"                                                 ( cf "Anions for the HP-41" )

-Lines 20 & 139 are three-byte GTO's

 01  LBL "ZETAA"   02  XEQ 14   03  RCL 01   04  X^2   05  +   06  X#0?   07  GTO 00   08  .5   09  CHS   10  STO 01   11  RCL 00   12   E3   13  /   14  ISG X   15  RTN   16  LBL 00   17  RCL 01   18  .5   19  X<=Y?   20  GTO 01   21  RCL 00   22   E3   23  /   24  ISG X   25  RCL 00   26  5   27  *   28  +   29   E3   30  /   31  1   32  +   33  STO M   34  REGMOVE   35  X<>Y   36  ST- 01   37  PI   38  XEQ "X^A"   39  RCL 00   40   E3   41  /   42  ISG X   43  RCL 00   44  3   45  *   46  +   47   E3   48  / 49  1   50  +   51  REGMOVE   52  RCL M   53  REGSWAP   54  REGMOVE   55  SIGN   56  CHS   57  XEQ "ST*A"    58  ST- 01   59  XEQ 01   60  RCL 00   61  .1   62  %   63  ISG X   64  +   65   E3   66  /   67  1   68  +   69  RCL 00   70  3   71  *   72  +   73  REGMOVE   74  XEQ "A*A1"    75  FRC   76  LASTX   77  RCL 00   78  4   79  *   80  +   81   E3   82  /   83  1   84  +   85  REGMOVE   86  +   87  REGSWAP   88  REGMOVE   89  2   90  CHS   91  XEQ "ST/A"   92  1/X   93  ST- 01   94  XEQ "GAMA"   95  RCL 00   96  + 97   E3   98  /   99  1 100  + 101  RCL 00 102  4 103  * 104  + 105  REGMOVE 106  XEQ "A*A1"  107  FRC 108  LASTX 109  RCL 00 110  4 111  * 112  + 113   E3 114  / 115  1 116  + 117  REGMOVE 118  + 119  REGSWAP 120  REGMOVE 121  2 122  XEQ "ST/A" 123  XEQ "GAMA" 124  XEQ "1/A" 125  RCL 00 126  .1 127  % 128  ISG X 129  + 130   E3 131  / 132  1 133  + 134  RCL 00 135  4 136  * 137  + 138  REGMOVE 139  GTO 03 140  LBL 14 141  RCL 00 142   E-3 143  + 144  0 145  LBL 13 146  RCL IND Y 147  X^2 148  + 149  DSE Y 150  GTO 13 151  SQRT 152  RTN 153  LBL 01 154  RCL 00 155  .1 156  % 157  ISG X 158  + 159   E3 160  / 161  1 162  + 163  REGMOVE   164  RCL 00 165  3 166  * 167  .1 168  % 169  ISG X 170  + 171  RCL 00 172  - 173  CLRGX 174  XEQ 14 175  PI 176  2 177  / 178  X<>Y 179  ST* Y 180  ST+ X 181  LN1+X 182  + 183  3 E10 184  LN 185  + 186  8 187  SQRT 188  3 189  + 190  LN 191  / 192  INT 193  1 194  STO O 195  STO P 196  + 197  STO M 198  STO N 199  LBL 02 200  RCL 00 201  PI 202  INT 203  10^X 204  / 205  ISG X 206  LASTX 207  / 208  ENTER^ 209  SIGN 210  + 211  RCL 00 212  + 213  REGMOVE   214  RCL M 215  1/X 216  XEQ "X^A" 217  SIGN 218  CHS 219  RCL M 220  Y^X 221  RCL O 222  * 223  XEQ "ST*A" 224  XEQ "DSA" 225  RCL M 226  ENTER^ 227  ST+ Y 228  ST* Y 229  - 230  RCL P 231  * 232  RCL N 233   X^2 234  RCL M 235  ENTER^ 236  SIGN 237  - 238  X^2 239  - 240  ST+ X 241  / 242  ST+ O 243  STO P 244  DSE M 245  GTO 02 246  RCL 00 247   E3 248  / 249  ISG X 250  LASTX 251  / 252  1 253  + 254  RCL 00 255  + 256  REGMOVE 257  SIGN 258  CHS 259  XEQ "ST*A"  260  ST- 01 261  2 262  XEQ "X^A"  263  SIGN 264  ST- 01 265  XEQ "1/A" 266  RCL 00 267  .1 268  % 269  ISG X 270  + 271   E3 272  / 273  1 274  + 275  RCL 00 276  ST+ X 277  + 278  REGMOVE 279  RCL O 280  XEQ "ST/A" 281  LBL 03 282  XEQ "A*A1" 283  CLA 284  END

( 469 bytes / SIZE 6n+1 )

 STACK INPUT OUTPUT X / 1.nnn

where   1.nnn   is the control number of the result

Examples:                               a = -7.6 + 7.7 i + 7.8 j + 7.9 k                          ( Quaternion ->   4   STO 00 )

7.6   CHS  STO 01
7.7   STO 02
7.8   STO 03
7.9   STO 04

XEQ "ZETAA"    >>>>    1.004                                      ---Execution time = 3m00s---

R01 =  762.9852912
R02 = -14.26823510
R03 = -14.45353680
R04 = -14.63883870

-Whence,      Zeta ( -7.6 + 7.7 i + 7.8 j + 7.9 k ) = 762.9852912 - 14.26823510 i - 14.45353680 j - 14.63883870 k

-Likewise,  Zeta( 1.1 + 7 i + 8 j + 9 k ) = 0.389296028 - 0.027737071 i - 0.031699509 j - 0.035661948 k          ( in 2m08s )

Notes:

-If a is very close to 1, there could be a loss of accuracy when 2^(1-a)-1 is executed:
-Replace lines 261 to 264 by

 261  2 262  LN 263  RCL 01 264  * 265  STO 01        266  E^X-1 267  XEQ 14 268  STO M  269  2 270  LN 271  * 272  STO N        273  X<>Y 274  RDN 275  RAD 276  COS 277  ST* Y 278  1 279  - 280  + 281  RCL 01 282  E^X 283  RCL N 284  SIN 285  * 286  RCL M        287  X=0? 288  SIGN 289  / 290  XEQ "ST*A" 291  X<>Y 292  STO 01

-The last decimals may depend on the version of "E^A" , "A*A1" , ....  that you are using.
-If Re(a) is large, it could be preferable to employ the 1st version of "ZETAA".
-In all other cases, Borwein algorithm is the best !

-Using the M-Code routines AMOVE & ASWAP would save many bytes...

11°) Lommel Functions

-"LOMA" 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(a)  = am+1 / [ (m+1)2 - p21F2 ( 1 ; (m-p+3)/2 , (m+p+3)/2 ; -a2/4 )

s(2)m,p(a)  = am+1 / [ (m+1)2 - p21F2 ( 1 ; (m-p+3)/2 , (m+p+3)/2 ; a2/4 )
+ 2m+p-1 Gam(p) Gam((m+p+1)/2) a -p / Gam((-m+p+1)/2) 0F1 ( ; 1-p ; -a2/4 )
+ 2m-p-1 Gam(-p) Gam((m-p+1)/2) ap / Gam((-m-p+1)/2) 0F1 ( ; 1+p ; -a2/4 )

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

Data Registers:           •  R00 = n > 1                          ( Registers R00 thru Rnn are to be initialized before executing "LOMA" )

•  R01   ......  •  Rnn = the n components of the anion a

Rn+1 ..........  R5n:  temp

>>>  When the program stops:     R01   ......  Rnn = the n components of  f(a)

Flags:  F02 & F04               CF 02  to compute  s(1)m,p(a)
SF 02  to compute  s(2)m,p(a)

Subroutines:   "0F1A"                                            ( cf paragraph 6°)c) above - the second version )
"DSA"  "ST*A"  "ST/A"  "HGFA"    ( cf  "Anionic Special Functions(I)" )
"A*A1"   "A^X"                               ( cf "Anions for the HP-41" )
"1/G+"                                              ( cf "Gamma Function for the HP-41" )

-Line 139 may be replaced by         E-3    XEQ "HGFB"    ( cf paragraph 6°)b) above )
-Lines 59-60 may be replaced by  1.002   XEQ "HGFB"

-Lines 102 & 103  are synthetic three-byte GTOs

 01  LBL "LOMA"   02  RCL 00   03  .1   04  %   05  ISG X   06  +   07  RCL 00   08  ST+ X   09  +   10   E3   11  /   12  1   13  +   14  REGMOVE   15  CLX   16  RCL 00   17  4   18  *   19  2   20  +   21  STO T   22  RDN   23  STO IND Z   24  DSE Z   25  STO N   26  X<>Y   27  STO IND Z   28  STO O   29  2   30  XEQ "ST/A"   31  RCL 00   32  .1   33  %   34  ISG X   35  +   36   E3   37  /   38  1   39  +   40  REGMOVE   41  XEQ "A*A1"   42  SIGN   43  CHS   44  FC? 02   45  XEQ "ST*A" 46  CHS   47  RCL O   48  RCL N   49  -   50  RCL O   51  LASTX   52  +   53  3   54  ST+ Z   55  +   56  2   57  ST/ Z   58  /   59  SF 04   60  XEQ "HGFA"   61  STO Y   62  RCL 00   63  +   64   E3   65  /   66  1   67  +   68  REGMOVE   69  CLX   70   E3   71  /   72  1   73  +   74  RCL 00   75  3   76  *   77  +   78  REGMOVE   79  LASTX   80  RCL 00   81  +   82  2   83  +   84  RCL IND X     85  STO N   86  DSE Y   87  RCL IND Y   88  STO O   89  1   90  + 91  XEQ "A^X"   92  XEQ "A*A1"   93  RCL O   94  1   95  +   96  X^2   97  RCL N   98  X^2   99  - 100  XEQ "ST/A" 101  FC? 02 102  GTO 03 103  GTO 02 104  LBL 01 105  RCL 00 106   E3 107  / 108  ISG X 109  LASTX 110  / 111  RCL 00 112  3 113  * 114  + 115  1 116  + 117  REGMOVE 118  2 119  XEQ "ST/A" 120  RCL 00 121  .1 122  % 123  ISG X 124  + 125  E3 126  / 127  1 128  + 129  REGMOVE  130  XEQ "A*A1" 131  SIGN 132  CHS 133  XEQ "ST*A" 134  RCL O 135  RCL N 136  1 137  RCL Y 138  - 139  XEQ "0F1A" 140  RDN 141  STO M 142  RDN 143  STO N 144  X<>Y 145  STO O 146  RCL 00 147  .1 148  % 149  ISG X 150  STO Z 151  + 152   E3 153  ST/ Z 154  / 155  1 156  ST+ Z 157  + 158  REGMOVE 159  CLX 160  RCL 00 161  3 162  * 163  + 164  REGMOVE 165  RCL N 166  CHS 167  XEQ "A^X" 168  XEQ "A*A1" 169  2 170  RCL O 171  RCL M 172  - 173  Y^X 174  STO M 175  RCL N 176  1 177  + 178  RCL O 179  - 180  2 181  / 182  XEQ "1/G+" 183  ST* M 184  RCL N 185  XEQ "1/G+" 186  ST/ M 187  RCL N 188  RCL O 189  + 190  1 191  + 192  2 193  / 194  XEQ "1/G+" 195  ST/ M 196  RCL M 197  XEQ "ST*A"  198  RCL 00 199  .1 200  % 201  ISG X 202  + 203  RCL 00 204  + 205   E3 206  / 207  1 208  + 209  RCL 00 210  4 211  * 212  + 213  REGMOVE 214  XEQ "DSA" 215  RTN 216  LBL 02 217  RCL 00 218   E3 219  / 220  ISG X 221  RCL 00 222  4 223  * 224  + 225   E3 226  / 227  1 228  + 229  REGMOVE 230  XEQ 01 231  RCL 00 232   E3 233  / 234  ISG X 235  RCL 00 236  4 237  * 238  + 239   E3 240  / 241  1 242  + 243  RCL 00 244  2 245  * 246  + 247  REGMOVE   248  SIGN 249  CHS 250  ST* N 251  XEQ 01 252  RCL 00 253   E3 254  / 255  ISG X 256  LASTX 257  / 258  1 259  + 260  RCL 00 261  ST+ X 262  + 263  REGMOVE 264  LBL 03 265  RCL 00 266  E3 267   / 268  ISG X 269  CLA 270  END

( 462 bytes / SIZE 4n+3 or 5n+1 )

 STACK INPUTS OUTPUTS Y m / X p 1.nnn

where   1.nnn   is the control number of the result

Examples:      m = sqrt(2)    ,    p = sqrt(3)    ,         a = 1 + 1.1 i + 1.2 j + 1.3 k                ( Quaternion ->   4   STO 00 )

CF 02  Lommel functions of the 1st kind:

CF 02

1     STO 01
1.1   STO 02
1.2   STO 03
1.3   STO 04

2     SQRT
3     SQRT   XEQ "LOMA"    >>>>    1.004                                      ---Execution time = 19s---

R01 = -2.555663121
R02 =  1.081642020
R03 =  1.179973113
R04 =  1.278304206

-Whence   s(1)sqrt(2),sqrt(3)( 1 + 1.1 i + 1.2 j + 1.3 k ) = -2.555663121 + 1.081642020 i + 1.179973113 j + 1.278304206 k

SF 02  Lommel functions of the 2nd kind:

SF 02

1     STO 01
1.1   STO 02
1.2   STO 03
1.3   STO 04

2     SQRT
3     SQRT   XEQ "LOMA"    >>>>    1.004                                      ---Execution time = 68s---

R01 =  1.461743799
R02 = -0.935738231
R03 = -1.020805344
R04 = -1.105872457

-Whence   s(2)sqrt(2),sqrt(3)( 1 + 1.1 i + 1.2 j + 1.3 k ) = 1.461743799 - 0.935738231 i - 1.020805344 j - 1.105872457 k

Note:

-Using M-Code routines +  AMOVE & ASWAP would reduce the number of bytes to 247 instead of 462 !