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<Y?
 61  GTO 01
 62  RCL N
 63  REGMOVE
 64  RCL 00
 65   E3
 66  /
 67  ISG X
 68  CLA
 69  CLD
 70  END

 
        ( 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  X<Y?
  79  X<>Y
  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
038   READATA
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
038  READATA
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  X<Y?
  55  X<>Y
  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 - p2 ]  1F2 ( 1 ; (m-p+3)/2 , (m+p+3)/2 ; -a2/4 )

    s(2)m,p(a)  = am+1 / [ (m+1)2 - p2 ]  1F2 ( 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 !