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  X0?   62  GTO 03   63  RCL 03   64  XY   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
38D  ?NCXQ
008   02E3
106   A=C S&X
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"
128  WRIT 4(L)
0E8  WRIT 3(X)
2A0  SETDEC
03C  RCR 3
268   WRIT 9(Q)
0B0  C=N ALL
03C  RCR 3
10E  A=C ALL
260   SETHEX
226   C=C+1 S&X
306   ?A<C S&X
08F   JC +17d=11h
268   WRIT 9(Q)
270   RAMSLCT
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
2A0  SETDEC
01D  ?NCXQ             C=
060   1807                A+C
13D  C=
060   AB*C
0A8   WRIT 2(Y)
353   JNC -22d=-16h
333   JNC -26d=-1Ah
0B0   C=N ALL
10E   A=C ALL
260   SETHEX
226   C=C+1 S&X
268   WRIT 9(Q)
306   ?A<C S&X
08F   JC +17d=11h
270   RAMSLCT
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
2A0  SETDEC
01D  ?NCXQ             C=
060   1807                 A+C
10E   A=C ALL
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
01D  ?NCXQ             C=
060   1807                 A+C
1E8   WRIT 7(O)
10E   A=C ALL
0AE  A<>C ALL
261   ?NCXQ             C=
060   1898                 A/C
13D  C=
060   AB*C
0A8  WRIT 2(Y)
025   C=
060   AB+C
10E   A=C ALL
0AE  A<>C ALL
0E8   WRIT 3(X)
3CC  ?KEY
360   ?C RTN
36E   ?A#C ALL
267   JC-52d=-34h                          replace this word by     277   JC -50d=-32h    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  XY 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  XY   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
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
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
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"
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
03C   RCR 3
1E8   WRIT 7(O)
260    SETHEX
10E   A=C ALL
226   C=C+1 S&X
1E8   WRIT 7(O)
306   ?A<C S&X
147   JC+28h=+40d
270   RAMSLCT
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
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.
10E   A=C ALL
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
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
0AE  A<>C ALL
0E8   WRIT 3(X)
3C4  ST=0                 C
045   ?NCXQ            =
06C  1B11               A^C
10E   A=C ALL
135   ?NCXQ           C=
060   184D              A*C
1A8  WRIT 6(N)
10E   A=C ALL
261   ?NCXQ                 C=
060   1898                     A/C
1A8  WRIT 6(N)
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)
03C  RCR 3
268   WRIT 9(Q)
03C   RCR 3
10E   A=C ALL
260   SETHEX
226   C=C+1 S&X
306   ?A<C S&X
097   JC+12h=18d
268   WRIT 9(Q)
270   RAMSLCT
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
2A0  SETDEC
01D   ?NCXQ             C=
060   1807                  A+C
10E   A=C ALL
135   ?NCXQ           C=
060   184D              A*C
1A8  WRIT 6(N)
34B   JNC-17h=-23d
2BB  JNC-29h=-41d
10E   A=C ALL
260   SETHEX
226   C=C+1 S&X
268   WRIT 9(Q)
306   ?A<C S&X
0DF  JC+1Bh=+27d
270   RAMSLCT
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
2A0  SETDEC
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
10E   A=C ALL
261   ?NCXQ                 C=
060   1898                     A/C
1A8  WRIT 6(N)
2FB  JNC-21h=-33d
2A0  SETDEC
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)
0A8  WRIT 2(Y)
0E8  WRIT 3(X)                     The rest of the routine is similar to @HGF
03C  RCR 3
268   WRIT 9(Q)
03C  RCR 3
10E  A=C ALL
260   SETHEX
226   C=C+1 S&X
306   ?A<C S&X
097   JC +18d=12h
268   WRIT 9(Q)
270   RAMSLCT
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
2A0  SETDEC
01D  ?NCXQ             C=
060   1807                 A+C
10E   A=C ALL
135   ?NCXQ            C=
060   184D                A*C
0A8   WRIT 2(Y)
34B   JNC -17h=-23d
32B   JNC -1Bh=-27d
10E   A=C ALL
260   SETHEX
226   C=C+1 S&X
268   WRIT 9(Q)
306   ?A<C S&X
08F   JC +17d=11h
270   RAMSLCT
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
2A0  SETDEC
01D  ?NCXQ             C=
060   1807                 A+C
10E   A=C ALL
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
01D  ?NCXQ             C=
060   1807                 A+C
068   WRIT 1(Z)
10E   A=C ALL
0AE  A<>C ALL
261   ?NCXQ             C=
060   1898                 A/C
10E   A=C ALL
135   ?NCXQ            C=
060   184D                A*C
0A8  WRIT 2(Y)
10E   A=C ALL
01D  ?NCXQ             C=
060   1807                 A+C
10E   A=C ALL
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/