The Museum of HP Calculators

# Hypergeometric Functions for the HP-41

This program is Copyright © 2006 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

1°) Hypergeometric Functions  ( slightly improved )
a) Program#1
b) Program#2
2°) Generalized Hypergeometric Functions   ( new )

1°) Hypergeometric Funtions

-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)!  F(a-c+1;b-c+1;2-c;x)

-This may be used in some of the 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 Functions for the HP-41"

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               ( lines 100 to 131 may be replaced by  RCL 00   XEQ "HGF" )
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)

where   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  .....

2°) Generalized Hypergeometric Funtions

-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] . xn/n!

where (ai)k = ai(ai+1)(ai+2) ...... (ai+k-1)   &  (ai)0 = 1    ,    likewise for  (bi)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 )

Example:    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  =  Hypergeometric Function
1F1  =  Kummer's Function

Reference:        Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4