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) ₀F₁(;a;x)
    b) ₀F₃(;a,b,c;x)
    c) ₁F₂(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 Lim_{t 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)₁(b)₁/(c)₁. x¹/1! + ............. + (a)_n(b)_n/(c)_n . xⁿ/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:

Lim_{t tends to c} ( F(a,b,t,x) )/Gam(t) = [ (a)_1-c(b)_1-c/(1-c)! ] x^1-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) = Lim_{t -> 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 = Lim_t->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) = Lim_{t -> 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 = Lim_t->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) = Lim_{t -> 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 )

_mF_p( a₁,a₂,....,a_m; b₁,b₂,....,b_p; x ) = SUM_{k=0,1,2,.....} [(a₁)_k(a₂)_k.....(a_m)_k] / [(b₁)_k(b₂)_k.....(b_p)_k] . x^k/k!

where (a_i)_k = a_i(a_i+1)(a_i+2) ...... (a_i+k-1) & (a_i)₀ = 1 , likewise for (b_j)_k

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

• R01 = a₁ • R02 = a₂ ....................... • Rmm = a_m
• Rm+1 = b₁ • Rm+2 = b₂ .................... • Rm+p = b_p
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	_mF_p( a₁,a₂,....,a_m; b₁,b₂,....,b_p; x )
X	x	_mF_p( a₁,a₂,....,a_m; b₁,b₂,....,b_p; x )

Examples: Calculate ₃F₄( 1 , 4 , 7; 2 , 3 , 6 , 5; Pi ) and ₄F₃( 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" >>>>   ₃F₄( 1 , 4 , 7; 2 , 3 , 6 , 5; Pi ) = 1.631019643 ( in 27 seconds )

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

Remarks:

₂F₁ = "The" hypergeometric Function
₁F₁ = Kummer's Function

b) M-Code Routine

-This M-code routine calculates the sum SUM_{k=k0,k0+1,k0+2,.....} [(a₁)_k(a₂)_k.....(a_m)_k] / [(b₁)_k(b₂)_k.....(b_p)_k] . x^k/k!
-Here, k₀ is not necessarily 0
-@HGF checks that register Rm+p does exist ( NONEXISTENT is displayed otherwise )

-It takes   k₀   in synthetic register O       and     u_k0 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
08F JC +17d=11h
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
0B8   READ 2(Y)
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
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
138   READ 4(L)
13D C=
060 AB*C
0A8 WRIT 2(Y)
0F8   READ 3(X)
025   C=
060 AB+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
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 u_k
X	x	_mF_p( a₁,a₂,....,a_m; b₁,b₂,....,b_p; 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" >>>> ₃F₄( 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    >>>> ₄F₃( 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 b_j 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( a₁,a₂,....,a_m; b₁,b₂,....,b_p; x ) = SUM_{k=0,1,2,.....} [ (a₁)_k(a₂)_k.....(a_m)_k] / [Gam(k+b₁) Gam(k+b₂).....Gam(k+b_p)] . x^k/k!

-If no b_j is a negative integer, we simply divide F by the product Gam(b₁) Gam(b₂).....Gam(b_p)
-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 = a₁ • R02 = a₂ ....................... • Rmm = a_m
• Rm+1 = b₁ • Rm+2 = b₂ .................... • Rm+p = b_p
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) = _mF_p( a₁,a₂,....,a_m; b₁,b₂,....,b_p; x ) if Y = +p hypergeometric function

and f(x) = _mF^~_p( a₁,a₂,....,a_m; b₁,b₂,....,b_p; 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+" >>>>   ₃F₄( 1 , 4 , 7; 2 , 3 , 6 , 5; Pi ) = 1.631019643 ( in 6 seconds )

    3   ENTER^
-4   ENTER^               ~
   PI     R/S      >>>>   ₃F₄( 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     >>>> ₂F₂( 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:

_mF_p( a₁,a₂,....,a_m; b₁,b₂,....,b_p; x ) if Y-input > 0
_mF^~_p( a₁,a₂,....,a_m; b₁,b₂,....,b_p; 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 = a₁    • R02 = a₂   .......................   • Rmm = a_m
                                    • Rm+1 = b₁ • Rm+2 = b₂ ....................   • Rm+p = b_p
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 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) = _mF_p( a₁,a₂,....,a_m; b₁,b₂,....,b_p; x ) if Y = + p hypergeometric function

and f(x) = _mF^~_p( a₁,a₂,....,a_m; b₁,b₂,....,b_p; 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+" >>>>   ₃F₄( 1 , 4 , 7; 2 , 3 , 6 , 5; Pi ) = 1.631019643 ( in 5 seconds )

    3   ENTER^
-4   ENTER^                         ~
   PI XEQ "HGF+"   >>>>   ₃F₄( 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+"   >>>>   ₂F₂( 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

   ₀F₁ is useful to compute Bessel and Airy functions
   ₀F₃ facilitates the calculation of Kelvin functions.
   ₁F₂ 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) ₀F₁(;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	₀F₁(;a;x)

Example:

2 SQRT STO 01

PI XEQ "0F1" >>>> ₀F₁(;2^1/2;PI) = 5.198957195

b) ₀F₃(;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	₀F₃(;a,b,c;x)

Example:

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

PI XEQ "0F3" >>>> ₀F₃(;2^1/2,3^1/2,5^1/2;PI) = 1.616609701 ---Execution time = 5 seconds---

c) ₁F₂(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	₀F₃(;a,b,c;x)

Example:

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

PI XEQ "1F2" >>>> ₁F₂(2^1/2;3^1/2,5^1/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/

hp41programs

Hypergeometric Functions for the HP-41