hp41programs

Hypergeo

Hypergeometric Functions for the HP-41


Overview
 

 1°) Gauss' Hypergeometric Function

    a) Program#1
    b) Program#2
    c) Program#3

 2°) Generalized Hypergeometric Functions

    a) Program#1
   b) M-Code Routine
   c) Regularized Generalized Hypergeometric Functions
   d) More Complete M-Code Routine

 3°) 3 Special cases

    a)  0F1(;a;x)
    b)  0F3(;a,b,c;x)
    c)  1F2(a;b,c;x)
 

1°) Gauss' Hypergeometric Funtion
 

-The first program calculates  F(a,b,c,x)  for | x | < 1
-The second routine also computes   Limt tends to c ( F(a,b,t,x) )/Gam(t)   when c is a negative integer.
 

     a) Program#1
 

Formula:      F(a,b;c;x) = 1 +  (a)1(b)1/(c)1. x1/1! + ............. +  (a)n(b)n/(c)n . xn/n! + ..........  where (a)n = a(a+1)(a+2) ...... (a+n-1)    and     | x | < 1
 

Data Registers:            R00 =  x

                                   •   R01 = a
                                   •   R02 = b              registers R01 R02 R03 are to be initialized before executing "HGF"
                                   •   R03 = c
Flags: /
Subroutines:  /
 
 

01  LBL "HGF"
02  STO 00     
03  CLST
04  SIGN
05  ENTER^
06  ENTER^
07  LBL 01
08  RDN
09  X<>Y
10  RCL 01     
11  R^
12  ST+ Y
13  RDN
14  *
15  RCL 02
16  R^
17  ST+ Y
18  RDN
19  *
20  RCL 03     
21  R^
22  ST+ Y
23  ISG X
24  CLX
25  ST* Y
26  RDN
27  /
28  RCL 00     
29  *
30  STO Z
31  X<>Y
32  ST+ Y
33  X#Y?
34  GTO 01     
35  END

 
  ( 51 bytes / SIZE 004 )
 
 

      STACK        INPUTS      OUTPUTS
           X             x       F(a,b;c;x)
           L             /             x

 
Example:   Calculate  F(0.3 , -0.7 ; 0.4 ; 0.49)

  0.3   STO 01    -0.7   STO 02    0.4    STO 03

  0.49   XEQ "HGF"  >>>>  F(0.3 , -0.7 ; 0.4 ; 0.49) =  0.720164356     ( in 18 seconds )
 

     b) Program#2
 

-This second program uses the following properties:

- If    c-a-b > 0  ,  F(a;b;c;1) = Gam(c).Gam(c-a-b)/(Gam(c-a).Gam(c-b))
- F(a;b;b;x) = F(b;a;b;x) = 1/(1-x)a
- If  c is a negative integer and neither a nor b is a negative integer such that  -a<-c , -b<-c
  then F is not defined but "HGF2" sets flag F00 and calculates:

         Limt tends to c ( F(a,b,t,x) )/Gam(t)  = [ (a)1-c(b)1-c/(1-c)! ] x1-c F(a-c+1;b-c+1;2-c;x)

-This may be useful to compute Associated Legendre Functions.
 

Data Registers:             R00:  x

                                   •   R01 = a
                                   •   R02 = b              registers R01 R02 R03 are to be initialized before executing "HGF2"
                                   •   R03 = c

Flags:  F00 & F05
Subroutine:  "GAM"  ( only if  x = 1 )  cf "Gamma Function for the HP-41"

-Lines 100 to 131 may be replaced by  RCL 00   XEQ "HGF"
 
 

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

 
  ( 240 bytes / SIZE 004 )
 
 

      STACK        INPUTS      OUTPUTS
           X             x            f(x)

       f(x) = F(a;b;c;x)                                            if F00 is clear
       f(x) =  Limt -> c ( F(a,b,t,x) )/Gam(t)              if F00 is set

Examples:        1  STO 01   2  STO 02   3  STO 03         0.1  XEQ "HGF2"  >>>>  F(1;2;3;0.1) = 1.072103131  ( in 7 seconds )

and similarly                                F(1;2;7;1)   = 1.5                    ( 14 s )           ( flag F00 is clear )
                                                  F(2;3;3;0.7) = 11.11111111   ( 1.4 s )           ( flag F00 is clear )

          Lim t tends to -7  ( F(2;2;t;0.1)/Gam(t) ) = 0.105229334   ( 18 s )            ( flag F00 is set )
          Lim t tends to -1  ( F(-4;-5;t;1)/Gam(t) ) = 420                  ( 17 s )            ( flag F00 is set )    .....  etc  .....
 

     c) Program#3
 

-"HGF" & "HGF2" assume that | x | doesn't exceed 1
-"HGF3" overcomes this limitation.
 

Formulae:

   If  x < 0     F(a,b,c,x) =  ( 1 - x ) -a  F(a,c-b,c,-x/(1-x))

   If  x > 1    F(a,b,c,x) = [ Gam(c) Gam(b-a) / Gam(b) / Gam(c-a) ] (-x) -a  F(a,1+a-c,1+a-b,1/x)
                                   + [ Gam(c) Gam(a-b) / Gam(a) / Gam(c-b) ] (-x) -b  F(b,1+b-c,1+b-a,1/x)       which may be a complex number ( F04 is set )
 

Data Registers:             R00 and  R04 to R07: temp

                                   •   R01 = a
                                   •   R02 = b                          ( Registers R01 R02 R03 are to be initialized before executing "HGF3" )
                                   •   R03 = c

Flags:  F00 , F04 , F07    If F04 is set at the end, the result is a complex number.
                                         If F04 is clear at the end, the result is real.

                                         If F00 is set at the end, X-output = Limt->c ( F(a,b,t,x) )/Gam(t)

Subroutine:  "GAM"  or  "GAM+" ....   ( cf "Gamma Function for the HP-41" )

-Line 57 is a three-byte GTO 00
 
 

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

 
  ( 335 bytes / SIZE 005 or 008 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             / Im (f(x)) if SF 04
           X             x      Re( f(x) )

      f(x) = F(a;b;c;x)                                            if F00 is clear             If CF 04  f(x) is real
      f(x) =  Limt -> c ( F(a,b,t,x) )/Gam(t)              if F00 is set                If SF 04   f(x) is complex.

Examples:

  •   1.2  STO 01    2.3  STO 02    3.7  STO 03

     0.4  XEQ "HGF3"  >>>>   1.435242953                                  ( in 18 seconds )      F00 is cleared in the first 4 examples.
      1    XEQ "HGF3"  >>>>  16.23332443                                     ( 7 seconds )
     -3   XEQ "HGF3"  >>>>   0.309850661                                    ( 45 seconds )
      4    XEQ "HGF3" >>>>  -0.804024119 - 0.105379849 i          ( 41 seconds )        F04 is set in the 4th example

  •   2  STO 01  STO 02   -7  STO 03

     0.1   XEQ "HGF3"  >>>>    0.105229334          ( 26 seconds )     but  F00 is set.

Whence:   Lim t tends to -7  ( F(2;2;t;0.1)/Gam(t) ) = 0.105229334

Notes:

-I've used a M-Code routine GAM to compute these values, so the last decimal may vary according to the "GAM" version that you are using.
-Execution times may also be different.
-If F04 is set but if Y-output is very small, f(x) is probably real.
-If F00 & F04 are both set, the result is meaningless and line 104 ( SF 99 ) produces an error message: NONEXISTENT.

-If  x > 1 , the result will be often - but not always - complex.
-So, if you only deal with real numbers, you could use a simplified version that works for x <= 1 ( negative arguments too )
-The following one uses the M-Code function  X=N? ( cf "A few M-Code Routines for the HP-41" ) to simplify the tests at the beginning of "HGF2"
 

  •   The initial definition is used if  x > 0
  •   Otherwise, F(a,b,c,x) =  ( 1 - x ) -a  F(a,c-b,c,-x/(1-x))  is employed.

  •   F(a;b;b;x) = F(b;a;b;x) = 1/(1-x)a  is also applied. If you don't want to use it, delete lines 33 to 47 and 30-31
 

Data Registers:             R00 = x

                                   •   R01 = a
                                   •   R02 = b         ( Registers R01 R02 R03 are to be initialized before executing "HGF3" )
                                   •   R03 = c

Flag:     F00       If F00 is set at the end,  X-output = Limt->c ( F(a,b,t,x) )/Gam(t)
Subroutines: /

-Line 06, this M-Code function may be replaced by  1  -
-Line 35, this M-Code function may be replaced by  RCL 01   X<>Y
-Line 116, this M-Code function may be replaced by   ISG X   CLX
 
 

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

 
   ( 167 bytes / SIZE 004 )
 
 

      STACK        INPUTS      OUTPUTS
           X             x           f(x) 

     f(x) = F(a,b,c;x)                                            if F00 is clear
     f(x) =  Limt -> c ( F(a,b,t;x) )/Gam(t)              if F00 is set

Examples:

  •   1.2  STO 01    2.3  STO 02    3.7  STO 03

     0.4  XEQ "HGF3"  >>>>   1.435242954            ( in 16 seconds )      F00 is cleared
     -3   XEQ "HGF3"  >>>>   0.309850661              ( 40 seconds )        F00 is cleared

  •   2  STO 01  STO 02   -7  STO 03

     0.1   XEQ "HGF3"  >>>>    0.105229334          ( 13 seconds )          F00 is set

Whence:   Lim t tends to -7  [ F( 2 , 2 ; t ; 0.1 )/Gam(t) ] = 0.105229334
 

2°) Generalized Hypergeometric Funtions
 

     a) Program#1
 

-The definition of  F(a,b,c;x) may be generalized as follows:    given 2 non-negative integers  m & p  ( at least one of them positive )

        mFp( a1,a2,....,am ; b1,b2,....,bp ; x ) =  SUMk=0,1,2,.....    [(a1)k(a2)k.....(am)k] / [(b1)k(b2)k.....(bp)k] . xk/k!

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

Data Registers:             R00 =  x                            ( Registers R01 thru Rm+p are to be initialized before executing "GHGF"  )

                                    •  R01 = a1    •  R02 = a2   .......................   •  Rmm = am
                                    •  Rm+1 = b1 •  Rm+2 = b2  ....................   •  Rm+p = bp
Flags: /
Subroutines:  /

-Synthetic registers may be replaced by any standard registers,
 but if, for instance, you replace register O by R99,
 delete line 44 and replace line 02 by  0  STO 99  RDN
 
 

01  LBL "GHGF"
02  CLA
03  STO 00
04  X<> Z
05  ST+ Y
06   E3
07  ST/ Z
08  /
09  STO N
10  X<>Y
11  STO M        
12  SIGN
13  STO T
14  LBL 01
15  R^
16  RCL 00
17  *
18  ISG M
19  LBL 02
20  RCL IND M
21  RCL O
22  +
23  ISG N
24  FS? 30
25  1/X
26  *
27  ISG M
28  GTO 02
29  RCL M        
30  FRC
31  STO M
32  X<> N
33  FRC
34  STO N
35  SIGN
36  RCL O
37  +
38  STO O        
39  /
40  STO T
41  +
42  X#Y?
43  GTO 01
44  CLA
45  END

 
  ( 76 bytes / SIZE m+p+1 )
 
 

      STACK      INPUTS                       OUTPUTS
           Z           m                              /
           Y           p        mFp( a1,a2,....,am ; b1,b2,....,bp ; x )
           X           x        mFp( a1,a2,....,am ; b1,b2,....,bp ; x )

Examples:    Calculate   3F4( 1 , 4 , 7 ; 2 , 3 , 6 , 5 ; Pi )   and 4F3( 1 , 4 , 7 , 2 ; 3 , 6 , 5 ; 1/Pi )

    1  STO 01   4  STO 02   7  STO 03   2  STO 04   3  STO 05   6  STO 06   5  STO 07

    3  ENTER^
    4  ENTER^
   PI  XEQ "GHGF"  >>>>   3F4( 1 , 4 , 7 ; 2 , 3 , 6 , 5 ; Pi )  =  1.631019643  ( in 27 seconds )

   4  ENTER^
   3  ENTER^
  PI  1/X  R/S    >>>>   4F3( 1 , 4 , 7 , 2 ; 3 , 6 , 5 ; 1/Pi ) =  1.258050204  ( in 50 seconds )

Remarks:

  2F1  =  "The" hypergeometric Function
  1F1  =  Kummer's Function
 

     b)  M-Code Routine
 

-This M-code routine calculates the sum    SUMk=k0,k0+1,k0+2,.....    [(a1)k(a2)k.....(am)k] / [(b1)k(b2)k.....(bp)k] . xk/k!
-Here, k0 is not necessarily 0
-@HGF checks that register  Rm+p  does exist  ( NONEXISTENT is displayed otherwise )

-It takes   k0   in synthetic register O       and     uk0  in register  Y     ( here, the first term of the sum may be different from 1. This is used in the next paragraph )
                m    in synthetic register N                    x   in register   X
             m+p  in synthetic register M

Warning:    @HGF  do not check for alpha data. So, be sure that X , Y , M , N , O and R01 thru Rm+p  do not contain alpha strings
                  Check also that m+p is not smaller than p
 

086   "F"
007   "G"
008   "H"
000   "@"
1B8   READ 6(N)              first executable word
38D  ?NCXQ
008   02E3
00E   A=0 ALL
106   A=C S&X
0AE  A<>C ALL
1BC  RCR 11
10E  A=C ALL
178  READ 5(M)
38D  ?NCXQ
008   02E3
106   A=C S&X
378   READ 13(c)
03C  RCR 3
0EE  C<>B ALL
0C6  C=B S&X
1BC  RCR 11
0C6  C=B S&X
20E  C=C+A ALL
070  N=C ALL
10E  A=C ALL
130  LDI S&X
200  200h                         200h  is the correct value if you have an HP-41 CX, CV or C with a Quad memory module or 4 memory modules.
306  ?A<C S&X              If you have an HP-41C without any memory module, replace 200h by 100h
381  ?NCGO
00A  02E0                        if register Rm+p  does not exist, the routine stops after displaying "NONEXISTENT"
0F8  READ 3(X)
128  WRIT 4(L)
0B8  READ 2(Y)
0E8  WRIT 3(X)
2A0  SETDEC
378   READ 13(c)
03C  RCR 3
268   WRIT 9(Q)
0B0  C=N ALL
03C  RCR 3
10E  A=C ALL
278   READ 9(Q)
260   SETHEX
226   C=C+1 S&X
306   ?A<C S&X
097   JC +18d=12h
268   WRIT 9(Q)
270   RAMSLCT
038   READDATA
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
2A0  SETDEC
1F8   READ 7(O)
01D  ?NCXQ             C=
060   1807                A+C
10E   A=C ALL
0B8   READ 2(Y)
135   ?NCXQ            C=
060   184D                A*C
0A8   WRIT 2(Y)
34B   JNC -23d=-17h
32B   JNC -27d=-1Bh
0B0   C=N ALL
10E   A=C ALL
278   READ 9(Q)
260   SETHEX
226   C=C+1 S&X
268   WRIT 9(Q)
306   ?A<C S&X
08F   JC +17d=11h
270   RAMSLCT
038   READDATA
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
2A0  SETDEC
1F8   READ 7(O)
01D  ?NCXQ             C=
060   1807                 A+C
10E   A=C ALL
0B8   READ 2(Y)
0AE  A<>C ALL
261   ?NCXQ             C=
060   1898                 A/C
0A8   WRIT 2(Y)
34B   JNC -23d=-17h
2A0   SETDEC
00E   A=0 ALL           A
35C  PT=12                =
162   A=A+1 @PT     1
1F8   READ 7(O)
01D  ?NCXQ             C=
060   1807                 A+C
1E8   WRIT 7(O)
10E   A=C ALL
0B8   READ 2(Y)
0AE  A<>C ALL
261   ?NCXQ             C=
060   1898                 A/C
10E   A=C ALL
138   READ 4(L)
135   ?NCXQ            C=
060   184D                A*C
0A8  WRIT 2(Y)
10E   A=C ALL
0F8   READ 3(X)
01D  ?NCXQ             C=
060   1807                 A+C
10E   A=C ALL
0F8   READ 3(X)
0AE  A<>C ALL
0E8   WRIT 3(X)
3CC  ?KEY
360   ?C RTN
36E   ?A#C ALL
257   JC-54d=-36h                          replace this word by     267   JC -52d=-34h    if you do not key in the 2 words:   3CC    360   written in red above.
3E0   RTN
 

-The routine stops when 2 consecutive sums are equal ( the result is in X-register )
-The 2 words  3CC  360  allows you to stop the routine if there is an infinite loop:  simply press any key for a few seconds to stop the loop.
-Even if the series is convergent, execution time may be prohibitive.
-But if you have a "newest" HP-41, these lines are unuseful: simply press  ENTER^   ON  to stop the loop

-We can now write a much shorter and faster "GHGF"
 
 

 01  LBL "GHGF"
 02  CLA
 03  STO 00
 04  RDN
 05  ABS
 06  X<>Y
 07  ABS
 08  STO N
 09  +
 10  STO M
 11  SIGN
 12  RCL 00
 13  @HGF
 14  CLA
 15  END

 
( 27 bytes / SIZE m+p )
 
 

      STACK      INPUTS                       OUTPUTS
           Z           m                              /
           Y           p                  the first neglected uk
           X           x         mFp( a1,a2,....,am ; b1,b2,....,bp ; x )

 
Examples:

  •   1  STO 01   4  STO 02   7  STO 03   2  STO 04   3  STO 05   6  STO 06   5  STO 07

      3  ENTER^
      4  ENTER^
     PI  XEQ "GHGF"  >>>>  3F4( 1 , 4 , 7 ; 2 , 3 , 6 , 5 ; Pi )  =  1.631019643  ( in 5 seconds instead of 27 )

     4  ENTER^
     3  ENTER^
    PI  1/X  R/S    >>>> 4F3( 1 , 4 , 7 , 2 ; 3 , 6 , 5 ; 1/Pi ) =  1.258050204  ( in 10 seconds instead of 50 )

Notes:

-If m = p = 0, "GHGF" returns exp(x)
-The series is always convergent if  m < p + 1
-If m = p + 1, the series converges if  | x | < 1
-If a  bj is a negative integer, the function is not defined and you'll get DATA ERROR after some time...
  unless 2 consecutive sums were equal before, but the result is meaningless in this case.
 

     c) Regularized Generalized Hypergeometric Functions
 

-The regularized generalized hypergeometric function F tilde is defined by

      mF~p( a1,a2,....,am ; b1,b2,....,bp ; x ) =  SUMk=0,1,2,.....    [ (a1)k(a2)k.....(am)k ] / [Gam(k+b1) Gam(k+b2).....Gam(k+bp)] . xk/k!

-If no bj is a negative integer, we simply divide F by the product    Gam(b1) Gam(b2).....Gam(bp)
-Otherwise,  the calculation is more complex.
-"HGF+" calculates  F or F tilde  according to the sign of Y-input.
 
 

Data Registers:             R00 =  x                            ( Registers R01 thru Rm+p are to be initialized before executing "HGF+"  )

                                    •  R01 = a1    •  R02 = a2   .......................   •  Rmm = am
                                    •  Rm+1 = b1 •  Rm+2 = b2  ....................   •  Rm+p = bp
Flags: /
Subroutines:  /

-This program uses several M-code routines:

     GAM  listed in  "Gamma Function for the HP-41"
     X+1   X-1   X=N? listed in  "A few M-Code Routines for the HP-41"
     and  @HGF listed in paragraph 2-b) above
 
 

01  LBL "HGF+"
02  CLA
03  STO 00
04  RDN
05  STO T
06  ABS
07  X<>Y
08  ABS
09  STO N
10  +
11  STO M
12  R^
13  X<0?
14  GTO 00
15  SIGN
16  GTO 05
17  LBL 00
18  X<> L
19   E3
20  ST/ M
21  ST/ N
22  /
23  +
24  1
25  ENTER^
26  LBL 01
27  RCL IND Z
28  X=N?
29  X>0?
30  GTO 01
31  X<Y?
32  X<>Y
33  GTO 00
34  LBL 01
35  GAM
36  ST/ Z
37  LBL 00
38  RDN
39  DSE Z
40  GTO 01
41  X>0?
42  GTO 00
43  CHS
44  X+1
45  STO O
46  RCL 00       
47  RCL Y
48  Y^X
49  R^
50  *
51  X<>Y
52  LBL 02
53  ST/ Y
54  X-1
55  X<>Y
56  ISG M
57  LBL 03
58  RCL Y
59  RCL IND M
60  +
61  ISG N
62  GTO 04
63  X=N?
64  X>0?
65  GTO 03
66  CLX
67  SIGN
68  LBL 03
69  1/X
70  LBL 04
71  *
72  ISG M
73  GTO 03
74  RCL M
75  FRC
76  STO M
77  X<> N
78  FRC
79  STO N
80  X<> Z
81  X#0?
82  GTO 02       
83  LBL 00
84   E3
85  ST* M
86  ST* N
87  X<> Z
88  LBL 05
89  RCL 00
90  @HGF
91  CLA
92  END

 
  ( 144 bytes / SIZE m+p+1 )
 
 

      STACK      INPUTS        OUTPUTS
           Z           m                /
           Y         +/- p                /
           X           x             f(x)
           L           /               x

  where  f(x) =   mFp( a1,a2,....,am ; b1,b2,....,bp ; x )    if   Y = +p      hypergeometric function
 
   and     f(x) =   mF~p( a1,a2,....,am ; b1,b2,....,bp ; x )    if   Y = -p     regularized hypergeometric function

Examples:

    •    1  STO 01   4  STO 02   7  STO 03   2  STO 04   3  STO 05   6  STO 06   5  STO 07

    3   ENTER^
    4   ENTER^
   PI  XEQ "HGF+"  >>>>   3F4( 1 , 4 , 7 ; 2 , 3 , 6 , 5 ; Pi )  =  1.631019643  ( in 6 seconds )

    3   ENTER^
  -4   ENTER^               ~
   PI     R/S      >>>>   3F4( 1 , 4 , 7 ; 2 , 3 , 6 , 5 ; Pi )  =  0.0002831631328  ( in 8 seconds )

   •    1  STO 01   4  STO 02   -2  STO 03   -5  STO 04

   2    ENTER^
  -2   ENTER^            ~
  0.1    R/S     >>>>  2F2( 1 , 4 ; -2 , -5 ; 0.1 ) = 0.01289656888  ( in 16 seconds )

Notes:

-If m = p = 0  "HGF+"  returns exp(x)
-The program doesn't check if the series are convergent or not.
-Even when they are convergent, execution time may be prohibitive: press any key to stop @HGF
-Remember that @HGF does not check for alpha data...
-The following variant doesn't use @HGF and is of course much slower.

-It seems impossible to compute F tilde without using extra data-registers if you don't have a M-Code function GAM
-Assuming p+q < 16, you can replace lines 25 thru 43 by

    STO 16        LBL 01              FRC                    RCL 17         GTO 00             ST/ 18              RCL 18
    1                  RCL IND 16      X#0?                   X>Y?            LBL 01              LBL 00            RCL 17
    STO 17        X>0?                 GTO 01               X<>Y           RCL IND 16      DSE 16            X>0?
    STO 18        GTO 01             RCL IND 16       STO 17         XEQ "GAM"      GTO 01           GTO 00

-Lines 64 to 74 may be replaced by

          ISG N             GTO 03          CLX           X<> L
          GTO 04           FRC               SIGN          1/X
          STO L             X#0?              SIGN           LBL 04
          X>0?               GTO 03          LBL 03        *
 
 

  01  LBL "HGF+"
  02  CLA
  03  STO 00
  04  RDN
  05  STO T
  06  ABS
  07  X<>Y
  08  ABS
  09  STO N
  10  +
  11  STO M
  12   E3
  13  ST/ M
  14  ST/ N
  15  R^
  16  X<0?
  17  GTO 00
  18  SIGN
  19  STO T
  20  GTO 05
  21  LBL 00 
  22  X<> Z
  23  RCL N
  24  +
  25  1           
  26  ENTER^  
  27  LBL 01 
  28  RCL IND Z
  29  X=N?
  30  X>0?
  31  GTO 01
  32  X<Y?
  33  X<>Y
  34  GTO 00
  35  LBL 01
  36  GAM
  37  ST/ Z
  38  LBL 00
  39  RDN
  40  DSE Z
  41  GTO 01
  42  X>0?
  43  GTO 00
  44  1
  45  X<>Y
  46  -
  47  STO O
  48  RCL 00
  49  RCL Y
  50  Y^X
  51  R^
  52  *
  53  X<>Y
  54  LBL 02
  55  ST/ Y
  56  1
  57  -
  58  X<>Y
  59  ISG M
  60  LBL 03
  61  RCL Y   
  62  RCL IND M
  63  +
  64  ISG N   
  65  GTO 04
  66  X=N?  
  67  X>0?   
  68  GTO 03
  69  CLX
  70  SIGN
  71  LBL 03
  72  1/X
  73  LBL 04
  74  *
  75  ISG M
  76  GTO 03
  77  RCL M
  78  FRC
  79  STO M
  80  X<> N
  81  FRC
  82  STO N
  83  X<> Z
  84  X#0?
  85  GTO 02
  86  LBL 00
  87  X<>Y
  88  STO T
  89  LBL 05
  90  R^
  91  RCL 00
  92  *
  93  ISG M
  94  LBL 06
  95  RCL IND M
  96  RCL O
  97  +
  98  ISG N
  99  FS? 30
100  1/X
101  *
102  ISG M
103  GTO 06
104  RCL M        
105  FRC
106  STO M
107  X<> N
108  FRC
109  STO N
110  SIGN
111  RCL O
112  +
113  STO O
114  /
115  STO T
116  +
117  X#Y?
118  GTO 05
119  CLA
120  END

 
   ( 183 bytes / SIZE m+p+1 )
 

-The instructions are identical, but x is not saved in L-register.
-If  m = p = 0, the result is meaningless.
 

   d) More Complete M-Code Routine
 

-The M-Code program listed hereunder takes  m , +/-p , x  in registers  Z , Y , X  and returns:

     mFp( a1,a2,....,am ; b1,b2,....,bp ; x )          if  Y-input > 0
     mF~p( a1,a2,....,am ; b1,b2,....,bp ; x )         if  Y-input < 0

-It checks that register Rm+p exists but it does not check for alpha data.

-An M-Code program "GAM" is also called as a subroutine ( cf "Gamma Function for the HP-41" §1-c-2 or §1-d-2 )
 
 

Data Registers:             R00 =  unused                            ( Registers R01 thru Rm+p are to be initialized before executing "HGF+"  )

                                    •  R01 = a1    •  R02 = a2   .......................   •  Rmm = am
                                    •  Rm+1 = b1 •  Rm+2 = b2  ....................   •  Rm+p = bp
Flags: /
Subroutine:   An M-Code function that computes Gamma(x)
 

0AB  "+"
006   "F"
007   "G"
008   "H"
2A0   SETDEC
284   CLRF7
078   READ 1(Z)
05E   C=0 MS          C = | C |           m may positive or negative, only the magnitude is taken into account
070   N=C ALL
10E   A=C ALL
0B8   READ 2(Y)
2FE   ?C#0 MS
013   JNC +02
288   SETF 7
05E   C=0 MS          C = | C |
01D   ?NCXQ           C=
060   1807               A+C
0F0   C<>N ALL
260   SETHEX
38D  ?NCXQ
008   02E3
00E   A=0 ALL
106   A=C S&X
0AE  A<>C ALL
1BC  RCR 11
10E   A=C ALL
0B0   C=N ALL
38D  ?NCXQ
008   02E3
106   A=C S&X
378   READ 13(c)
03C  RCR 3
0EE  C<>B ALL
0C6  C=B S&X
1BC  RCR 11
0C6  C=B S&X
20E  C=A+C ALL
228  WRIT 8(P)
10E  A=C ALL
130  LDI S&X
200  200h                         200h  is the correct value if you have an HP-41 CX, CV or C with a Quad memory module or 4 memory modules.
306  ?A<C S&X              If you have an HP-41C without any memory module, replace 200h by 100h
381  ?NCGO
00A  02E0                        If register Rm+p  does not exist, the routine stops after displaying "NONEXISTENT"
0F8   READ 3(X)
0A8  WRIT 2(Y)
2A0  SETDEC
04E  C=0 ALL
068   WRIT 3(Z)
35C   PT=12                   C=
050   LD@PT- 1              1
168   WRIT 5(M)
1A8   WRIT 6(N)
28C   ?FSET 7
249    ?NCGO                 If Y-input > 0 , we jump to the address  FD92  in my ROM ( see below ).
3F6    FD92                     Change these 2 words according to your own ROM
238    READ 8(P)
03C   RCR 3
1E8   WRIT 7(O)
260    SETHEX
238    READ 8(P)
10E   A=C ALL
1F8   READ 7(O)
226   C=C+1 S&X
1E8   WRIT 7(O)
306   ?A<C S&X
147   JC+28h=+40d
270   RAMSLCT
038   READDATA
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
0AE  A<>C ALL
070   N=C ALL
2A0  SETDEC
2EE  ?C#0 ALL
043   JNC+08
2FE  ?C#0 MS
08B  JNC +11h=+17d
084   CLRF 5               C
0ED  ?NCXQ              =
064   193B                frc(C)
2EE  ?C#0 ALL
067   JC +0Ch=12d
178   READ 5(M)
2BE  C=-C-1 MS      C=-C
10E   A=C ALL
0B0   C=N ALL
01D   ?NCXQ           C=
060   1807               A+C
2FE  ?C#0 MS
303   JNC -20h=-32d
0B0  C=N ALL
168   WRIT 5(M)
2EB   JNC-23h=-35d
0B0   C=N ALL
0E8   WRIT 3(X)
059   ?NCXQ                 F216  is the address of the first executable word of a M-Code routine GAM ( Gamma function ) in my ROM.
3C8   F216                     Change these 2 words to call your Gamma M-Code routine.
1B8   READ 6(N)
10E   A=C ALL
0F8   READ 3(X)
261   ?NCXQ                 C=
060   1898                     A/C
1A8   WRIT 6(N)
293    JNC -2Eh=-46d
2A0   SETDEC
00E   A=0 ALL               A
35C   PT=12                   =
162   A=A+1 @PT         1
178   READ 5(M)
36E   ?A#C ALL
249    ?NCGO                                  We skip lines to arrive at the address  FD92  in my ROM ( see below ).
3F6    FD92                                      Change these 2 words according to your own ROM
2BE  C=-C-1 MS      C=-C
01D   ?CXQ             C=
061   1807              A+C
068   WRIT 1(Z)
168   WRIT 5(M)
10E   A=C ALL
0B8  READ 2(Y)
0AE  A<>C ALL
0E8   WRIT 3(X)
3C4  ST=0                 C
045   ?NCXQ            =
06C  1B11               A^C
10E   A=C ALL
1B8   READ 6(N)
135   ?NCXQ           C=
060   184D              A*C
1A8  WRIT 6(N)
1B8   READ 6(N)
10E   A=C ALL
178   READ 5(M)
261   ?NCXQ                 C=
060   1898                     A/C
1A8  WRIT 6(N)
178   READ 5(M)
00E   A=0 ALL               A
1BE  A=A-1 MS             =
35C  PT=12
162   A=A+1@PT         -1
01D   ?NCXQ             C=
060   1807                  A+C
168   WRIT 5(M)
378   READ 13(c)
03C  RCR 3
268   WRIT 9(Q)
238   READ 8(P)
03C   RCR 3
10E   A=C ALL
278   READ 9(Q)
260   SETHEX
226   C=C+1 S&X
306   ?A<C S&X
097   JC+12h=18d
268   WRIT 9(Q)
270   RAMSLCT
038   READDATA
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
2A0  SETDEC
178   READ 5(M)
01D   ?NCXQ             C=
060   1807                  A+C
10E   A=C ALL
1B8   READ 6(N)
135   ?NCXQ           C=
060   184D              A*C
1A8  WRIT 6(N)
34B   JNC-17h=-23d
2BB  JNC-29h=-41d
238   READ 8(P)
10E   A=C ALL
278   READ 9(Q)
260   SETHEX
226   C=C+1 S&X
268   WRIT 9(Q)
306   ?A<C S&X
0DF  JC+1Bh=+27d
270   RAMSLCT
038   READDATA
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
2A0  SETDEC
178   READ 5(M)
01D   ?NCXQ             C=
060   1807                  A+C
2EE   ?C#0 ALL
373   JNC-12h=-18d
0E8   WRIT 3(X)
2FE   ?C#0 MS
033   JNC+06
084   CLRF 5              C
0ED  ?NCXQ             =
064   193B               frc(C)
2EE   ?C#0 ALL
333    JNC-1Ah=-26d
1B8   READ 6(N)
10E   A=C ALL
0F8   READ 3(X)
261   ?NCXQ                 C=
060   1898                     A/C
1A8  WRIT 6(N)
2FB  JNC-21h=-33d
2A0  SETDEC
178   READ 5(M)
2EE  ?C#0 ALL
2D7  JC-26h=-38d
0B8  READ 2(Y)                     this line is at the address  FD92  in my ROM ( cf the ?NCGO written in red above )
128   WRIT 4(L)
1B8   READ 6(N)
0A8  WRIT 2(Y)
0E8  WRIT 3(X)                     The rest of the routine is similar to @HGF
378   READ 13(c)
03C  RCR 3
268   WRIT 9(Q)
238   READ 8(P)
03C  RCR 3
10E  A=C ALL
278   READ 9(Q)
260   SETHEX
226   C=C+1 S&X
306   ?A<C S&X
097   JC +18d=12h
268   WRIT 9(Q)
270   RAMSLCT
038   READDATA
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
2A0  SETDEC
078   READ 1(Z)
01D  ?NCXQ             C=
060   1807                 A+C
10E   A=C ALL
0B8   READ 2(Y)
135   ?NCXQ            C=
060   184D                A*C
0A8   WRIT 2(Y)
34B   JNC -17h=-23d
32B   JNC -1Bh=-27d
238   READ 8(P)
10E   A=C ALL
278   READ 9(Q)
260   SETHEX
226   C=C+1 S&X
268   WRIT 9(Q)
306   ?A<C S&X
08F   JC +17d=11h
270   RAMSLCT
038   READDATA
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
2A0  SETDEC
078   READ 1(Z)
01D  ?NCXQ             C=
060   1807                 A+C
10E   A=C ALL
0B8   READ 2(Y)
0AE  A<>C ALL
261   ?NCXQ             C=
060   1898                 A/C
0A8   WRIT 2(Y)
34B   JNC -17h=-23d
2A0   SETDEC
00E   A=0 ALL           A
35C  PT=12                =
162   A=A+1 @PT     1
078   READ 1(Z)
01D  ?NCXQ             C=
060   1807                 A+C
068   WRIT 1(Z)
10E   A=C ALL
0B8   READ 2(Y)
0AE  A<>C ALL
261   ?NCXQ             C=
060   1898                 A/C
10E   A=C ALL
138   READ 4(L)
135   ?NCXQ            C=
060   184D                A*C
0A8  WRIT 2(Y)
10E   A=C ALL
0F8   READ 3(X)
01D  ?NCXQ             C=
060   1807                 A+C
10E   A=C ALL
0F8   READ 3(X)
0AE  A<>C ALL
0E8   WRIT 3(X)
3CC  ?KEY
360   ?C RTN
36E   ?A#C ALL
257   JC-54d=-36h                          replace this word by     267   JC -52d=-34h    if you do not key in the 2 words:   3CC    360   written in red above.
345   ?NCXQ            xeq
040   10D1               CLA
3C1   ?NCGO                                 4 subroutine levels are used, so  3E0    RTN
002    00F0                                                      must be replaced by   3C1   002     i-e   ?NCGO   00F0
 

-The routine stops when 2 consecutive sums are equal ( the result is in X-register )
-The 2 words  3CC  360  allows you to stop the routine if there is an infinite loop:  simply press any key for a few seconds to stop the loop.
-Even if the series is convergent, execution time may be prohibitive.
-But if you have a "newest" HP-41, these lines are unuseful: simply press  ENTER^   ON  to stop the loop.

-The words  345  040  clear alpha "register"
 
 

      STACK      INPUTS        OUTPUTS
           T            T               T
           Z           m       last k-value
           Y         +/- p   1st neglected term
           X           x             f(x)
           L           /               x

  where  f(x) =   mFp( a1,a2,....,am ; b1,b2,....,bp ; x )    if   Y = + p      hypergeometric function
 
   and     f(x) =   mF~p( a1,a2,....,am ; b1,b2,....,bp ; x )    if   Y = - p     regularized hypergeometric function

Examples:

    •    1  STO 01   4  STO 02   7  STO 03   2  STO 04   3  STO 05   6  STO 06   5  STO 07

    3   ENTER^
    4   ENTER^
   PI  XEQ "HGF+"  >>>>   3F4( 1 , 4 , 7 ; 2 , 3 , 6 , 5 ; Pi )  =  1.631019643  ( in 5 seconds )

    3   ENTER^
  -4   ENTER^                         ~
   PI  XEQ "HGF+"   >>>>   3F4( 1 , 4 , 7 ; 2 , 3 , 6 , 5 ; Pi )  =  0.0002831631326  ( in 6 seconds )

   •    1  STO 01   4  STO 02   -2  STO 03   -5  STO 04

   2    ENTER^
  -2   ENTER^                         ~
  0.1  XEQ "HGF+"   >>>>   2F2( 1 , 4 ; -2 , -5 ; 0.1 ) = 0.01289656888  ( in 4 seconds )

Notes:

-If m = p = 0 , HGF+  returns exp(x)
-The program doesn't check if the series are convergent or not.
-Even when they are convergent, execution time may be prohibitive: press any key to stop HGF+
-Remember that HGF+ does not check for alpha data...
-Stack register T is saved.
 
 

3°)  3 Special Cases
 

   0F1 is useful to compute Bessel and Airy functions
   0F3 facilitates the calculation of Kelvin functions.
   1F2 is used to compute Anger and Weber functions.

-Of course, one can use "GHGF" but the following routines are faster! ( but slower than "HGF+" with the M-code routine @HGF )
 

     a)   0F1(;a;x)
 

Data Registers:        R00 = x     •  R01 = a      ( Register R01 is to be initialized before executing "0F1" )
Flags: /
Subroutines: /
 
 

01  LBL "0F1"
02  STO 00    
03  CLST
04  SIGN
05  ENTER^
06  ENTER^
07  LBL 01
08  RDN
09  X<>Y
10  RCL 00    
11  *
12  RCL 01    
13  R^
14  ST+ Y
15  ISG X
16  CLX
17  ST* Y
18  RDN
19  /
20  STO Z    
21  X<>Y
22  ST+ Y
23  X#Y?
24  GTO 01    
25  END

 
  ( 39 bytes / SIZE 002 )
 
 

      STACK        INPUTS      OUTPUTS
           X             x      0F1(;a;x)

Example:

   2  SQRT   STO 01

  PI  XEQ "0F1"   >>>>   0F1(;21/2;PI) = 5.198957195
 

     b)    0F3(;a,b,c;x)
 

Data Registers:     R00 = x                           ( Registers R01 thru R03 are to be initialized before executing "0F3" )

                             •  R01 = a
                             •  R02 = b
                             •  R03 = c
Flags: /
Subroutines: /
 
 

01  LBL "0F3"
02  STO 00    
03  CLST
04  SIGN
05  ENTER^
06  ENTER^
07  LBL 01
08  RDN
09  X<>Y
10  RCL 00    
11  *
12  RCL 01
13  R^
14  ST+ Y
15  RDN
16  /
17  RCL 02    
18  R^
19  ST+ Y
20  RDN
21  /
22  RCL 03    
23  R^
24  ST+ Y
25  ISG X
26  CLX
27  ST* Y
28  RDN
29  /
30  STO Z
31  X<>Y
32  ST+ Y
33  X#Y?
34  GTO 01    
35  END

 
( 51 bytes / SIZE 004 )
 
 

      STACK        INPUTS      OUTPUTS
           X             x    0F3(;a,b,c;x)

Example:

   2  SQRT   STO 01
   3  SQRT   STO 02
   5  SQRT   STO 03

  PI  XEQ "0F3"   >>>>   0F3(;21/2,31/2,51/2;PI) = 1.616609701         ---Execution time = 5 seconds---
 

     c)    1F2(a;b,c;x)
 

Data Registers:     R00 = x                           ( Registers R01 thru R03 are to be initialized before executing "1F2" )

                             •  R01 = a
                             •  R02 = b
                             •  R03 = c
Flags: /
Subroutines: /

-Line 16 is the unique difference between "0F3" & "1F2"
 
 

01  LBL "1F2"
02  STO 00    
03  CLST
04  SIGN
05  ENTER^
06  ENTER^
07  LBL 01
08  RDN
09  X<>Y
10  RCL 00    
11  *
12  RCL 01
13  R^
14  ST+ Y
15  RDN
16  *
17  RCL 02    
18  R^
19  ST+ Y
20  RDN
21  /
22  RCL 03    
23  R^
24  ST+ Y
25  ISG X
26  CLX
27  ST* Y
28  RDN
29  /
30  STO Z
31  X<>Y
32  ST+ Y
33  X#Y?
34  GTO 01    
35  END

 
  ( 51 bytes / SIZE 004 )
 
 

      STACK        INPUTS      OUTPUTS
           X             x    0F3(;a,b,c;x)
 
 
Example:

   2  SQRT   STO 01
   3  SQRT   STO 02
   5  SQRT   STO 03

  PI  XEQ "1F2"   >>>>   1F2(21/2;31/2,51/2;PI) = 2.767636697         ---Execution time = 8 seconds---
 

References:

[1]  Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]   http://functions.wolfram.com/