hp41programs

Bionic Functions3

Bionic Special Functions(III) for the HP-41


Overview
 

 0°)  2 Non-Merged Functions & Other M-Code Routines
 1°)  Spheroidal Wave Functions
 2°)  Jacobi Functions
 3°)  Coulomb Wave Functions

    a)  Regular Coulomb Wave Function
    b)  Irregular Coulomb Wave Function
    c)  Asymptotic Expansions

 4°)  Weierstrass Elliptic Functions

    a)  Laurent Series
    b)  Duplication Formula

 5°)  Hypergeometric Functions
 6°)  Generalized Hypergeometric Functions
 
 
 

0°)  2 Non-Merged Functions & Other M-Code Routines
 

-We could create many extra-registers like R S U V W ....
-Though I've still used some of them, it's also useful to create functions that use data registers
  but the addresses of these registers may depend on the content of R00

-For example, if registers R01 thru R6n are used, we may have to store parameters into R6n+1 , R6n+2 , ... and so on
-Assuming R00 = 2n , two M-Code functions ( SSTO - - - & RRCL - - - ) are listed below
-The address - for instance 6n+1 - is extracted from the following line ( 301 ) in a program:

-The most important lines below were written by William Wilder ( and thanks to Ángel Martin who sent me these instructions ! )
 

1A0  A=B=C=0         @DFA0   in my ROM
158   M=C ALL
141   ?NCXQ                  get
0A4   2950                       PC
01D   ?NCXQ                 get
0B4   2D07                      nxtbyte
056    C=0 XS
2E6    ?C#0 S&X
3E3    JNC -04
39C   PT=0
06E   A<>B ALL
106   A=C S&X
130   LDI S&X
01A  26d
306   ?A<C S&X
06B  JNC+13d
042   C=0 @PT
306   ?A<C S&X
057   JC+10d
1D8  C<>M ALL
2FC  RCR 13
0A2  A<>C @PT
1D8  C<>M ALL
019   ?NCXQ                 nbytab
0B4   2D06
17E   A=A+1 S&X
056   C=0 XS
373   JNC -18d
06E   A<>B ALL
31D  .NCXQ
0A4   29C7
346   ?A#0 S&X
0BD  ?CXQ                   put
08D   232F                    PC X
198   C=M ALL
05A  C=0 M
05E   C=0 MS
07C  RCR 4
070   N=C ALL
2A0  SETDEC
046  C=0 S&X
270  RAMSLCT
378  C=c
03C  RCR 3
106   A=C S&X
130  LDI S&X
200  512d
306  ?A<C S&X
381   ?NCGO
00A   NEXIST
0AE  A<>C ALL
270   RAMSLCT
038   READATA
10E  A=C ALL
0B0  C=N ALL
19C  PT=11
04A  C=0 PT<-
135   C=
060   A*C
0B0  C=N ALL
35C  PT=12
042  C=0 @PT
2FC  RCR 13
2E2  ?C#0 @PT
01F  JC+03
2FC  RCR 13
013  JNC+02
226  C=C+1 S&X
025  C=
060  AB+C
260  SETHEX
38D  C->
008  S&X
106  A=C S&X
378  C=c
03C  RCR 3
146  A=A+C S&X         the address of the register is now in CPU register  A  S&X
130  LDI S&X
200  512d
306  ?A<C S&X
381   ?NCGO
00A   NEXIST
3E0   RTN                @DFF2   in my ROM
 

   ?NCXQ   DFA0   =     281   37C         Change these instructions if the routine is coded elsewhere in your own ROM
 

-Now, we just need 2 short routines to store or recall a number in X-register:
 

08F  "O"
014  "T"
113  "S"
113  "S"
281  ?NCXQ
37C  DFA0
046  C=0 S&X
270  RAMSLCT
0F8  C=X
0AE  A<>C ALL
270  RAMSLCT
0AE  A<>C ALL
2F0  WRITDATA
3E0  RTN
08C  "L"
003  "C"
112  "R"
112  "R"
281  ?NCXQ
37C  DFA0
0AE  A<>C ALL
270   RAMSLCT
038   READATA
10E  A=C ALL
046  C=0 S&X
270  RAMSLCT
0AE  A<>C ALL
028  T=C
3B5  ?NCGO
052   R^

-If, for example, in a program, you have the lines   SSTO   301    whereas R00 = 2n = 8  ( biquaternions )
 the content of X-register will be stored into R( 3x8 + 01 ) = R25
-Likewise, SSTO   329  will be equivalent to  STO 53  because  3 x 8 + 29 = 53

-RRCL works similarly.
-However,  CLX  RRCL - - -  will place  0  in Y-register  ( use  RDN  RRCL - - -  if you want to overwrite register X )

Notes:

-These M-Code functions are of course slower than the standard  STO & RCL and also slower than  STO W   RCL W ...
- SSTO - - -  requires about  115 ms
-Each instruction actually employs 2 lines and the second line is a three-digit number so, it occupies 5 bytes.

-The existence of the registers is checked but there is no check for alpha data.

-Provided n is not too large - you can use  R99  R98  R97 .....

-Other M-Code routines are also employed:

  Z*Z  1/Z  Z/Z   SINH  COSH                cf "A few M-Code Routines for the HP-41"

 and   X+1   X-1   X/2  3X  X/E3             -------------------------------------------

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

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

1°)  Spheroidal Wave Functions
 

-"SWFB" computes the angular spheroidal wave function of the first kind.
-Given  m , n  and c2 and the corresponding eigenvalue  l

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

     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 = 2n > 3         ( Registers R00 thru R2n & R8n+1 thru R8n+8 are to be initialized before executing "SWFB" )

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

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

                                       •  R8n+1-R8n+2 = m   •  R8n+3-R8n+4 = n   •  R8n+5-R8n+6 = c2   •  R8n+7-R8n+8 = l

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

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

-Lines 78-162-212  are three-byte GTO's
 
 

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

 
          ( 495 bytes / SIZE 8n+13 )
 
 

      STACK        INPUT      OUTPUT
           X             /         1.2n

   where       1.2n   is the control number of the result.

Example:        m = 0.2 + 0.3 i  ,  n = 0.6 + 0.7 i  ,  c2 = 1.7 + 1.8 i   ,   l = 1.462455887 + 2.197449718 i

                        b = ( 0.05 + 0.1 i ) + ( 0.15 + 0.2 i ) e1 + ( 0.25 + 0.3 i ) e2 + ( 0.35 + 0.4 i ) e3          ( biquaternion ->  8  STO 00 )

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

      0.2  STO 33   0.3 STO 34         1.7  STO 37   1.8  STO 38
      0.6  STO 35   0.7  STO 36         1.462455887  STO 39   2.197449718   STO 40
 

      XEQ "SWFB"    >>>>    1.008                             ---Execution time = 10m21s---

And we find in R01 thru R08

  Smn(b) = ( 0.628023288 - 0.491662142 i ) + ( -0.021846567 + 0.346286109 i ) e1
                                                                      + ( -0.006377756 + 0.541954055 i ) e2 + ( 0.009091056 + 0.737622001 i ) e3
 

2°)  Jacobi Functions
 

Formula:

     Pm(a;b) (b) = [ Gam(a+m+1) / Gam(m+1) ]  2F~1 ( -m , a+b+m+1 , a+1 , (1-b)/2 )

     where  2F1 tilde is the regularized hypergeometric function
 

Data Registers:           •  R00 = 2n > 3         ( Registers R00 thru R2n & R6n+1 thru R6n+6 are to be initialized before executing "JCFB" )

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

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

                                       •  R6n+1-R6n+2 = a   •  R6n+3-R6n+4 = b   •  R6n+5-R6n+6 = m

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

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

 01  LBL "JCFB"
 02  .5
 03  CHS
 04  ST*A
 05  ST- 01
 06  12
 07  AMOVE        
 08  CLX
 09  STO W
 10  ST*A
 11  SIGN
 12  STO 01 
 13  13
 14  AMOVE 
 15  LBL 01
 16  B*B1
 17  RRCL 
 18  306
 19  CHS
 20  RRCL
 21  305
 22  CHS
 23  RRCL
 24  303
 25  RRCL 
 26  301
 27  +
 28  RCL Y
 29  -
 30  RCL W         
 31  ST+ Z
 32  X+1
 33  STO W 
 34  +
 35  STO M
 36  RDN
 37  RRCL
 38  302
 39  RRCL
 40  304
 41  +
 42  R^
 43  -
 44  RCL M
 45  Z*Z
 46  RRCL
 47  301
 48  RCL W         
 49  +
 50  RRCL
 51  302
 52  X<>Y
 53  Z/Z
 54  RCL W
 55  ST/ Z
 56  /
 57  Z*B
 58  DSA
 59  X#0?
 60  GTO 01 
 61  31
 62  AMOVE
 63  RRCL
 64  302
 65  RRCL
 66  306
 67  +
 68  RRCL
 69  301
 70  RRCL
 71  305
 72  +
 73  X+1
 74  XEQ "GAMZ"
 75  Z*B
 76  RRCL
 77  302
 78  RRCL
 79  301
 80  X+1
 81  XEQ "GAMZ"
 82  1/Z
 83  Z*B
 84  RRCL
 85  306
 86  RRCL
 87  305
 88  X+1
 89  XEQ "GAMZ"
 90  1/Z
 91  Z*B
 92  RCL 00
 93  X/E3
 94  ISG X
 95  END

 
     ( 204 Bytes / SIZE 6n+7 )
 
 

      STACK        INPUT      OUTPUT
           X             /         1.2n

   where       1.2n   is the control number of the result.

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

                        b = ( 1.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3          ( biquaternion ->  8  STO 00 )

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

      0.2  STO 25    0.3  STO 26
      0.6  STO 27    0.7  STO 28
      1.4  STO 29    1.6  STO 30

      XEQ "JCFB"    >>>>   1.008                                 ---Execution time = 9m24s---

And we find in R01 thru R08

  Pm(a;b) (b) = ( -1.594720719 + 3.464906169 i ) + ( -1.081150690 - 0.614094392 i ) e1
                                                                              + ( -1.735722626 - 0.871495197 i ) e2 + ( -2.390294571 - 1.128896005 i ) e3
 

3°)  Coulomb Wave Functions
 

     a)  Regular Coulomb Wave Function
 

Formulae:         FL(h,b) = CL(h) b L+1 exp(-i.b)  M ( L+1-i.h ; 2L+2 ; 2i.b )       where   M = Kummer's function

    with   CL(h) = 2L exp [ (1/2) { -PI.h + Lngamma(L+1+i.h) + Lngamma(L+1-i.h) } - Lngamma(2L+2) ]
 

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

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

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

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

Flag:  /
Subroutines:  B*B1  Z*B  "B^Z"  ST*A   DSA  Z*Z  Z/Z  Z+Z  Z-Z  AMOVE  ASWAP  X+1  X/E3  X/2
                        STO R  RCL R  STO S  RCL S  STO U  RCL U  STO V  RCL V  STO W  RCL W  SSTO ---  RRCL ---
                        "LNGZ"    ( cf "Gamma Function for the HP-41" paragraph 1°) h-4) )
 
 

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

 
          ( 294 bytes / SIZE 8n+1 )
 
 

      STACK        INPUTS      OUTPUTS
           T          Im(L)            /
           Z          Re(L)            /
           Y          Im(h)            /
           X          Re(h)          1.2n

   where       1.2n   is the control number of the result.

Example:        L = 2 + 3 i  ,  h = -0.6 - 0.7 i

         b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3               ( biquaternion ->  2n = 8  STO 00 )

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

          3  ENTER^
          2  ENTER^
         0.7  CHS  ENTER^
         0.6  CHS  XEQ "RCWFB"  >>>>   1.008                                     ---Execution time = 2m30s---

And we find in R01 thru R08

        FL(h,b) = ( 1.432663948 + 2.724686924 i ) + ( -0.994443197 + 0.451696564 i ) e1
                                                                              + ( -1.515195662 + 0.784202097 i ) e2 + ( -2.035948129 + 1.116707629 i ) e3

Notes:

 "RCWFB"  does not work if  2.L = -1 , -2 , -3 , ...............
-Lines 24 to 39 multiply b by i  in R01 thru R2n
 

     b)  Irregular Coulomb Wave Function
 

Formulae:         GL(h,b) = [ FL(h,b) Cos c - F-L-1(h,b) ] / sin c

    with   c = sL(h) - s-L-1(h) - ( 2.L + 1 ) PI / 2
    and   sL(h) = [ Lngamma ( 1 + L + i.h )  -  Lngamma ( 1 + L - i.h ) ] / ( 2.i )
 
 

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

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

                                         R2n+1 ..........  R10n:  temp
 

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

Flag:  /
Subroutines:  Z*B  ST*A   DSA  A-A  Z*Z  1/Z  Z+Z  Z-Z  AMOVE  X+1  X-1  X/E3  X/2  SINH  COSH
                        STO R  RCL R  STO S  RCL S  STO U  RCL U  STO V  RCL V  STO W  RCL W  SSTO ---  RRCL ---
                        "LNGZ"    ( cf "Gamma Function for the HP-41" paragraph 1°) h-4) )
                        "RCWFB"  ( cf paragraph 3°) a) above )
 
 

  01  LBL "ICWF"
  02  R^
  03  CHS
  04  R^
  05  CHS
  06  X-1
  07  R^
  08  R^
  09  XEQ "RCWFB"
  10  15
  11  AMOVE
  12  41
  13  AMOVE
  14  RCL V 
  15  CHS
  16  RCL U
  17  CHS
  18  X-1
  19  RCL S
  20  RCL R
  21  XEQ "RCWFB"
  22  RCL V
  23  RCL R
  24  +
  25  RCL U
  26  RCL S
  27  -
  28  X+1
  29  XEQ "LNGZ"  
  30  SSTO
  31  301
  32  X<>Y
  33  STO W
  34  RCL V 
  35  RCL R
  36  +
  37  CHS
  38  RCL S
  39  RCL U
  40  -
  41  XEQ "LNGZ" 
  42  RCL W
  43  RRCL
  44 301
  45  Z+Z
  46  SSTO
  47  301
  48  X<>Y
  49  STO W
  50  RCL R
  51  RCL V 
  52  -
  53  RCL U
  54  RCL S
  55  +
  56  CHS
  57  XEQ "LNGZ"  
  58  RCL W
  59  RRCL
  60  301
  61  Z-Z
  62  SSTO
  63  301
  64  X<>Y
  65  STO W
  66  RCL V 
  67  RCL R
  68  -
  69  RCL U
  70  RCL S
  71  +
  72  X+1
  73  XEQ "LNGZ"  
  74  RCL W
  75  RRCL
  76  301
  77  Z+Z
  78  RCL V
  79  ST+ X
  80  PI
  81  *
  82  -
  83  X/2
  84  STO N
  85  SINH
  86  X<>Y
  87  CHS
  88  RCL U 
  89  ST+ X
  90  X+1
  91  PI
  92  *
  93  -
  94  X/2
  95  STO M            
  96  RAD
  97  SIN
  98  *
  99  CHS
100  RCL M
101  COS
102  RCL N
103  COSH
104  *
105  Z*B
106  52
107  AMOVE
108  A-A
109  RCL N
110  SINH
111  RCL M
112  COS
113  *
114  RCL N
115  COSH
116  RCL M            
117  SIN
118  *
119  1/Z
120  Z*B
121  RCL 00
122  X/E3
123  ISG X
124  DEG
125  CLA
126  END

 
      ( 250 bytes / SIZE 10n+1 )
 
 

      STACK        INPUTS      OUTPUTS
           T          Im(L)            /
           Z          Re(L)            /
           Y          Im(h)            /
           X          Re(h)          1.2n

   where       1.2n   is the control number of the result.

Example:        L = 2 + 3 i  ,  h = -0.6 - 0.7 i

         b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3               ( biquaternion ->  2n = 8  STO 00 )

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

          3  ENTER^
          2  ENTER^
         0.7  CHS  ENTER^
         0.6  CHS  XEQ "ICWFB"  >>>>   1.008                                     ---Execution time = 6m09s---

And we find in R01 thru R08

        GL(h,b) = ( -4.805817330 - 31.48658762 i ) + ( -13.14440341 - 0.948059260 i ) e1
                                                                               + ( -19.02111405 - 0.507420170 i ) e2 + ( -25.89782473 - 0.066781108 i ) e3

Note:

 "ICWFB"  does not work if  2.L is an integer
 

     c)  Asymptotic Expansion
 

Formulae:

    with   q = b - h Ln 2.b - L PI / 2 + [ Lngamma ( 1 + L + i.h )  -  Lngamma ( 1 + L - i.h ) ] / ( 2.i )

      HL+(h,b)  ~  exp ( i.q2F0 ( - L + i.h , 1 + L + i.h , 1 / 2.i.b )
      HL-(h,b)  ~  exp ( -i.q2F0 ( - L - i.h , 1 + L - i.h , -1 / 2.i.b )

   then,   FL(h,b) = [ HL+(h,b) -  HL-(h,b) ] / (2.i)               GL(h,b) = [ HL+(h,b) +  HL-(h,b) ] / 2
 
 

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

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

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

  >>>  When the program stops:     R01   ......  R2n = the n components of FL(h,b)
                                                      R2n+1 ....  R4n = the n components of GL(h,b)

Flag:  F04
Subroutines:  B*B1  Z*B  "LNB"  "E^B"  "1/B"  ST*A   DSA  A+A  A-A  Z*Z  AMOVE  ASWAP  X+1  X/E3  X/2
                        STO R  RCL R  STO S  RCL S  STO U  RCL U  STO V  RCL V  STO W  RCL W
                        "LNGZ"    ( cf "Gamma Function for the HP-41" paragraph 1°) h-4) )
 
 

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

 
            ( 320 bytes / SIZE 8n+1 )
 
 

      STACK        INPUTS      OUTPUTS
           T          Im(L)            /
           Z          Re(L)            /
           Y          Im(h)            /
           X          Re(h)          1.2n

   where       1.2n   is the control number of  FL(h,b)    [ the control number of  GL(h,b)  is  2n+1.4n

Example:        L = 1.2 + 1.3 i  ,  h = 0.6 + 0.7 i

         b = ( 5.1 + 5.2 i ) + ( 5.3 + 5.4 i ) e1 + ( 5.5 + 5.6 i ) e2 + ( 5.7 + 5.8 i ) e3               ( biquaternion ->  2n = 8  STO 00 )

       5.1  STO 01   5.2  STO 02   5.3  STO 03   5.4  STO 04   5.5  STO 05   5.6  STO 06   5.7  STO 07   5.8  STO 08

         1.3  ENTER^
         1.2  ENTER^
         0.7  ENTER^
         0.6  XEQ "AECWB"  >>>>   1.008                                     ---Execution time = 3m57s---

We find in R01 thru R08

        FL(h,b) = ( -4338.130616 + 4323.462574 i ) + ( 2250.237631 + 1939.122345 i ) e1
                                                                               + ( 2335.028740 + 2010.820458 i ) e2 + ( 2419.819849 + 2082.518573 i ) e3

and in R09 thru R16

        GL(h,b) = ( 4046.165344 + 3484.318834 i ) + ( 2414.155964 - 2404.336065 i ) e1
                                                                               + ( 2503.556959 - 2495.053311 i ) e2 + ( 2592.957955 - 2585.770556 i ) e3

Notes:

-As usual with asymptotic expansions, the series will diverge too soon if b is too "small"
-Add  RND  after line 149 to control the precision by the display settings in this case.
-Lines 94 to 109 multiply b by i  in registers R01 thru R2n
 

4°)  Weierstrass Elliptic Functions
 

     a)  Laurent Series
 

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

      P(b;g2;g3) = 1/ b2 + c2.b2 + c3.b4 + ...... + ck.b2k-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  R8n+1  R8n+2  ......
 

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

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

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

                                       •  R8n+1-R8n+2 = c2  •  R8n+3-R8n+4 = c3  •  R8n+5-R8n+6 = c4             ......................

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

Flags:  /
Subroutines:  B*B1  Z*B  "1/B"  ST*A   DSA  Z*Z  AMOVE  X+1  X-1  X/E3  X/3
                        STO U   RCL U   STO V   RCL V   STO W   RCL W
                        SSTO - - -   RRCL - - -  ( cf paragraph 0 above )

-Lines 121 & 126 are three-byte GTO's
 
 

  01  LBL "WEFB"
  02  STO M
  03  RDN
  04  STO N 
  05  CLX
  06  20
  07  ST/ Z
  08  /
  09  SSTO
  10  401
  11  X<>Y
  12  SSTO
  13  402
  14  28
  15  ST/ M
  16  ST/ N
  17  RCL M
  18  SSTO
  19  403
  20  RCL N
  21  SSTO
  22  404
  23  12
  24  AMOVE
  25  B*B1
  26  12
  27  AMOVE
  28  XEQ "1/B"
  29  13
  30  AMOVE 
  31  21
  32  AMOVE
  33  RRCL
  34  402
  35  RRCL
  36  401
  37  Z*B
  38  DSA
  39  21 
  40  AMOVE      
  41  B*B1
  42  14
  43  AMOVE
  44  RRCL
  45  404
  46  RRCL
  47  403
  48  Z*B
  49  DSA
  50  41
  51  AMOVE 
  52  CLX
  53  STO U 
  54  LBL 01
  55  3
  56  STO V
  57  LBL 02
  58  B*B1
  59  14
  60  AMOVE 
  61  RCL U 
  62  X+1
  63  STO U
  64  STO M
  65  ST+ X
  66  RCL 00       
  67  4
  68  *
  69  ST+ Y
  70  X+1
  71  STO N 
  72  X<>Y
  73  STO O
  74  LASTX
  75  +
  76  STO W
  77  CLX
  78  STO P
  79  STO Q
  80  LBL 03
  81  RCL IND O 
  82  DSE O
  83  RCL IND O
  84  RCL IND N
  85  ISG N
  86  CLX
  87  RCL IND N 
  88  X<>Y
  89  Z*Z
  90  ST+ P
  91  X<>Y
  92  ST+ Q
  93  ISG N
  94  CLX
  95  DSE O
  96  DSE M
  97  GTO 03
  98  RCL P
  99  RCL Q
100  RCL U
101  ST+ X
102  7
103  +
104  RCL U
105  *
106  X/3
107  ST/ Z
108  /
109  RCL W
110  RDN
111  STO IND T
112  DSE T
113  X<>Y
114  STO IND T 
115  Z*B
116  DSA
117  41
118  AMOVE
119  X<>Y
120  X#0?
121  GTO 01
122  RCL V 
123  X-1
124  STO V
125  X#0?
126  GTO 02
127  31
128  AMOVE
129  RCL 00
130  X/E3
131  ISG X
132  CLA
133  END

 
      ( 257 bytes / SIZE 8n+ ??? )
 
 

      STACK        INPUTS      OUTPUTS
           T          Im(g2)            /
           Z          Re(g2)            /
           Y          Im(g3)            /
           X          Re(g3)          1.2n

   where       1.2n   is the control number of the result.

Example:        g2 = 1.2 + 1.3 i  ,  g3 = 1.6 + 1.7 i

         b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3               ( biquaternion ->  2n = 8  STO 00 )

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

        1.3  ENTER^
        1.2  ENTER^
        1.7  ENTER^
        1.6  XEQ "WEFB"  >>>>   1.008                                     ---Execution time = 8m56s---

And we find in R01 thru R08

        P(b;g2;g3) = ( 0.089457020 + 0.121607346 i ) + ( -0.005461551 + 0.122834455 i ) e1
                                                                                 + (  0.001306738 + 0.192058673 i ) e2 + ( 0.008075026 + 0.261282892 i ) e3

Note:

 SIZE 095 ( at least ) in this example !
 

     b)  Duplication Formula
 

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

-"WF2B" takes P(b) in registers R01 .......... R2n  and returns P(2b) in the same registers
 

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

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

                                      •  R01   ......  •  R2n = the n components of P(b)

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

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

Flags:  /
Subroutines:  B*B1  "1/B"  ST*A  A+A  AMOVE  ASWAP  X/2
                        STO R   RCL R   STO S   RCL S   STO U   RCL U   STO V   RCL V
 
 

 01  LBL "WF2B"
 02  STO R 
 03  RDN
 04  STO S
 05  RDN
 06  STO U 
 07  X<>Y
 08  STO V
 09  12
 10  AMOVE
 11  13
 12  AMOVE
 13  B*B1
 14  14
 15  AMOVE      
 16  6
 17  ST*A
 18  RCL U 
 19  X/2
 20  ST- 01
 21  RCL V
 22  X/2
 23  ST- 02
 24  12
 25  AMOVE      
 26  B*B1
 27  14
 28  ASWAP
 29  4
 30  ST*A
 31  RCL U 
 32  ST- 01
 33  RCL V
 34  ST- 02
 35  32
 36  AMOVE      
 37  B*B1
 38  RCL R 
 39  ST- 01
 40  RCL S
 41  ST- 02
 42  4
 43  ST*A
 44  XEQ "1/B"
 45  42
 46  AMOVE      
 47  B*B1
 48  12
 49  AMOVE
 50  31
 51  AMOVE
 52  -2
 53  ST*A 
 54  A+A
 55  END

 
       ( 114 bytes / SIZE 8n+1 )
 
 

      STACK        INPUTS      OUTPUTS
           T          Im(g2)            /
           Z          Re(g2)            /
           Y          Im(g3)            /
           X          Re(g3)          1.2n

   where       1.2n   is the control number of the result.

Example:       g2 = 1.2 + 1.3 i  ,  g3 = 1.6 + 1.7 i   ,  first calculate   P(b;g2;g3)

    where  b = ( 0.05 + 0.1 i ) + ( 0.15 + 0.2 i ) e1 + ( 0.25 + 0.3 i ) e2 + ( 0.35 + 0.4 i ) e3     then deduce  P(2b;g2;g3)
 

 >>>  "WEFB"  gives in  2m06s

           P(b;g2;g3) = ( 0.403054401 + 1.804707647 i ) + ( -0.027262872 + 0.221241952 i ) e1
                                                                                    + ( -0.024830724 + 0.347318475 i ) e2 + ( -0.022398576 + 0.473394997 i ) e3

>>>  Then    1.3  ENTER^
                    1.2  ENTER^
                    1.7  ENTER^
                    1.6  XEQ "WEFB"  >>>>   1.008                                     ---Execution time = 18s---

And we find in R01 thru R08

        P(2b;g2;g3) = ( 0.089457022 + 0.121607344 i ) + ( -0.005461551 + 0.122834455 i ) e1
                                                                                   + (  0.001306737 + 0.192058672 i ) e2 + ( 0.008075025 + 0.261282893 i ) e3

Notes:

-We have saved more than 6 minutes
-SIZE 059 is enough instead of SIZE 095
 

5°)  Hypergeometric Functions
 

-"HGFB"  computes  2F1(p;q,s;b) = 1 +  (p)1/(q)1/(s)1. b1/1! + ............. +  (p)n/(q)n/(s)n . bk/k! + ..........

-However, if  2F1  does not exist ( it may happen if c is a negative integer ),
 the regularized function  2F1 tilde is returned.                 -> Flag F00 is set in that case.
 
 

Data Registers:           •  R00 = 2n > 3         ( Registers R00 thru R2n & R6n+1 thru R6n+6 are to be initialized before executing "HGFB" )

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

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

                                       •  R6n+1-R6n+2 = p   •  R6n+3-R6n+4 = q   •  R6n+5-R6n+6 = s

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

Flag:   F00
Subroutines:  B*B1  Z*B  ST*A   DSA  Z*Z  Z/Z  AMOVE  ASWAP  X+1  X-1  X/E3
                        STO V   RCL V   STO W   RCL W
                        SSTO - - -   RRCL - - -  ( cf paragraph 0 above )

-Lines 10-17 are three-byte GTO 02
 
 

  01  LBL "HGFB"
  02  12
  03  AMOVE
  04  CLX
  05  STO W 
  06  CF 00
  07  RRCL
  08  306
  09  X#0?
  10  GTO 02
  11  RRCL
  12  305
  13  ENTER^
  14  FRC
  15  X=0?
  16  X<Y?
  17  GTO 02
  18  RRCL 
  19  302
  20  X#0?
  21  GTO 00
  22  RRCL
  23  301
  24  ENTER^
  25  FRC
  26  X=0?
  27  X<Y?
  28  GTO 00
  29  RRCL
  30  305
  31  LASTX
  32  X>Y?
  33  GTO 02
  34  LBL 00 
  35  RRCL
  36  304
  37  X#0?
  38  GTO 00       
  39  RRCL
  40  303
  41  ENTER^
  42  FRC
  43  X=0?
  44  X<Y?
  45  GTO 00 
  46  RRCL
  47  305
  48  LASTX
  49  X>Y?
  50  GTO 02
  51  LBL 00
  52  SF 00
  53  1
  54  RRCL
  55  305
  56  -
  57  STO V
  58  STO W
  59  XEQ "B^X"
  60  CLST
  61  SIGN
  62  ENTER^
  63  LBL 03        
  64  RDN
  65  RRCL
  66  301
  67  X-1
  68  RCL V 
  69  +
  70  RRCL
  71  302
  72  X<>Y
  73  Z*Z
  74  RRCL 
  75  303
  76  X-1
  77  RCL V
  78  +
  79  RRCL
  80  304
  81  X<>Y
  82  Z*Z
  83  RCL V
  84  ST/ Y
  85  ST/ Z
  86  X-1
  87  STO V
  88  X#0?
  89  GTO 03       
  90  RDN
  91  Z*B
  92  LBL 02 
  93  FS? 00
  94  GTO 00 
  95  CLX
  96  ST*A
  97  SIGN
  98  STO 01
  99  LBL 00
100  13
101  AMOVE 
102  LBL 01
103  B*B1
104  RRCL
105  302
106  RRCL
107  301
108  RRCL
109  303
110  RCL W
111  ST+ Z
112  + 
113  RRCL
114  304
115  X<>Y
116  Z*Z
117  RRCL
118  305
119  RCL W        
120  + 
121  RRCL
122  306
123  X<>Y
124  Z/Z
125  RCL W
126  X+1
127  STO W 
128  ST/ Z
129  /
130  Z*B
131  DSA
132  X#0?
133  GTO 01 
134  31
135  AMOVE
136  RCL 00
137  X/E3
138  ISG X
139  END

 
     ( 261 bytes / SIZE 6n+7 )
 
 

      STACK        INPUT      OUTPUT
           X             /         1.2n

   where       1.2n   is the control number of the result.

Example1:        p = 0.2 + 0.3 i  ,  q = 0.6 + 0.7 i  ,  s = 1.4 + 1.6 i

                         b = ( 0.05 + 0.1 i ) + ( 0.15 + 0.2 i ) e1 + ( 0.25 + 0.3 i ) e2 + ( 0.35 + 0.4 i ) e3          ( biquaternion ->  8  STO 00 )

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

      0.2  STO 25    0.3  STO 26
      0.6  STO 27    0.7  STO 28
      1.4  STO 29    1.6  STO 30

      XEQ "HGFB"    >>>>   1.008                                 ---Execution time = 5m01s---

>>>   F00  is clear and we find in R01 thru R08

      F (b) = ( 1.020242672 + 0.003637406 i ) + ( -0.013019430 + 0.038951549 i ) e1
                                                                        + ( -0.017194186 + 0.061805972 i ) e2 + ( -0.021368943 + 0.084660394 i ) e3

Example2:        p = 0.2 + 0.3 i  ,  q = 0.6 + 0.7 i  ,  s = -4

                         b = ( 0.05 + 0.1 i ) + ( 0.15 + 0.2 i ) e1 + ( 0.25 + 0.3 i ) e2 + ( 0.35 + 0.4 i ) e3          ( biquaternion ->  8  STO 00 )

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

      0.2  STO 25    0.3  STO 26
      0.6  STO 27    0.7  STO 28
       -4  STO 29      0   STO 30

      XEQ "HGFB"    >>>>   1.008                                 ---Execution time = 10m17s---

>>>   F00 is set and we find in R01 thru R08

      Ftilde (b) = ( -13.42936413 - 18.84144453 i ) + ( 6.922905074 - 4.319886604 i ) e1
                                                                              + ( 10.45414095 - 7.292855522 i ) e2 + ( 13.98537690 - 10.26582445 i ) e3

Example3:        p = 0.2 + 0.3 i  ,  q = -3  ,  s = -4

                         b = ( 0.05 + 0.1 i ) + ( 0.15 + 0.2 i ) e1 + ( 0.25 + 0.3 i ) e2 + ( 0.35 + 0.4 i ) e3          ( biquaternion ->  8  STO 00 )

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

      0.2  STO 25    0.3  STO 26
       -3  STO 27      0   STO 28
       -4  STO 29      0   STO 30

      XEQ "HGFB"    >>>>   1.008                                 ---Execution time = 20s---

>>>   F00  is clear and we get in R01 thru R08

      F (b) = ( 1.041533875 + 0.022087750 i ) + ( -0.025527625 + 0.066593375 i ) e1
                                                                        + ( -0.034495625 + 0.105927875 i ) e2 + ( -0.043463625 + 0.145262375 i ) e3

Note:

-In the 3rd example, the function reduces to a polynomial and there is no roundoff-error !
 

6°)  Generalized Hypergeometric Functions
 

-This program computes   pFq( a1,a2,....,ap ; b1,b2,....,bq ; b ) =  SUMk=0,1,2,.....    [(a1)k(a2)k.....(ap)k] / [(b1)k(b2)k.....(bq)k] . bk/k!         if   X > 0

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

           ai & bj  are complexes    and    b is a "bi-on"

>>>   or the regularized function F tilde:

           pF~q( a1,a2,....,ap ; b1,b2,....,bq ; b ) =  SUMk=0,1,2,.....    [ (a1)k(a2)k.....(ap)k ] / [Gam(k+b1) Gam(k+b2).....Gam(k+bq)] . bk/k!      if  X < 0

    where Gam = Euler's Gamma function.
 
 
 

Data Registers:           •  R00 = 2n > 3         ( Registers R00 thru R2n & R6n+1 ......  are to be initialized before executing "GHGFB" )

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

                                         R2n+1 ........  R6n:  temp         •  R6n+1-R6n+2 = a1 ......  •  R6n+2p+2q-1-R6n+2p+2q = bq

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

Flags:  /
Subroutines:  B*B1   "B^X"  Z*B  ST*A   DSA  1/Z  Z*Z  Z/Z   AMOVE  X+1  X-1  3X  X*E3  X/E3
                        STO R  RCL R  STO S  RCL S  STO U  RCL U  STO V  RCL V  STO W  RCL W
                        SSTO - - -   RRCL - - -  ( cf paragraph 0 above )
                        "GAMZ"    ( a version that does not use any data register - cf "Gamma Function for the HP-41" paragraph 1°) h-3) )
                         "1/G"    M-Code routine  ( cf "Gamma Function for the HP-41" paragraph 1-f) )

-Line 42 is a three-byte  GTO 07
 
 

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

 
           ( 416 bytes / SIZE 6n+1+2p+2q )
 
 

      STACK        INPUT      OUTPUT
           X       +/- p.qqq         1.2n

   where       1.2n   is the control number of the result.

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

Example1:        a1 = 0.2 + 0.3 i  ,  b1 = 0.6 + 0.7 i  ,  b3 = 1.4 + 1.6 i              -> Calculate  2F3( a1,a2 ; b1,b2,b3 ; b )     with
                          a2 = 0.4 + 0.5 i  ,  b2 = 1.2 + 1.3 i

         b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3               ( biquaternion ->  2n = 8  STO 00 )

       0.2  STO 25   0.3  STO 26      0.6  STO 29   0.7  STO 30     1.4   STO 33   1.6  STO 34       ( 6n+1 = 25 )
       0.4  STO 27   0.5  STO 28      1.2  STO 31   1.3  STO 32
 

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

       2.003   XEQ "GHGFB"  >>>>   1.008                                     ---Execution time = 68s---

And we find in R01 thru R08

      2F3(b) = ( 1.003991856 + 0.005192296 i ) + ( 0.032853131 + 0.009144529 i ) e1
                                                                        + ( 0.051982446 + 0.011637214 i ) e2 + ( 0.071111761 + 0.014129900 i ) e3

Example2:      With the same arguments, calculate the regularized function F tilde:   2F~3( a1,a2 ; b1,b2,b3 ; b )

-Registers R25 thru R34 are unchanged, so:
 

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

       2.003  CHS   XEQ "GHGFB"  >>>>   1.008                                     ---Execution time = 1m53s---

And we get in R01 thru R08

      2F~3(b) = ( 9.112735981 + 3.989347215 i ) + ( 0.262953245 + 0.212181889 i ) e1
                                                                         + ( 0.427181613 + 0.309967488 i ) e2 + ( 0.591409981 + 0.407753086 i ) e3

-The result of the 1st example has been simply divided by  [ Gam(b1) Gam(b2) Gam(b3) ]
 

Example3:      Calculate again the regularized function F tilde:    2F~3( a1,a2 ; b1,b2,b3 ; b )    but with

                          a1 = 0.2 + 0.3 i  ,  b1 = 0.6 + 0.7 i  ,  b3 = -4
                          a2 = 0.4 + 0.5 i  ,  b2 = -3.14

         b = ( 0.1 + 0.2 i ) + ( 0.3 + 0.4 i ) e1 + ( 0.5 + 0.6 i ) e2 + ( 0.7 + 0.8 i ) e3

       0.2  STO 25   0.3  STO 26      0.6  STO 29   0.7  STO 30     -4   STO 33   0  STO 34
       0.4  STO 27   0.5  STO 28    -3.14  STO 31   0   STO 32
 

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

       2.003  CHS   XEQ "GHGFB"  >>>>   1.008                                     ---Execution time = 2m06s---

And we get in R01 thru R08

      2F~3(b) = ( 0.185476239 - 0.060667558 i ) + ( 0.133911078 + 0.071283991 i ) e1
                                                                          + ( 0.214604002 + 0.100490141 i ) e2 + ( 0.295296925 + 0.129696289 i ) e3
 

Notes:

-The M-Code routines  SSTO - - - & RRCL - - - may be avoided here:

-Replace lines 94 to 97 by         RCL 00  X+1  ST+ X  SIGN  CLX  RCL IND L  DSE L  RCL IND L
-And replace lines 80 to 85 by   RCL 00  X+1  ST+ X  X<> T  STO IND T  DSE T  X<> Z  STO IND T