Complex Hypergeometric Functions for the HP-41
Overview
1°) "The" Hypergeometric Function 2F1
2°) Generalized Hypergeometric Functions
mFp
( M-Code Routine )
3°) Application
to a few Special Functions
a) Bessel Functions:
Asymptotic Expansions
b) Fresnel Integrals:
Ascending Series & Asymptotic Series
c) Kelvin Functions:
Ascending Series & Asymptotic Series
d) Sine & Cosine
Integrals: Asymptotic Series
-The following programs deal with complex variables, but the
parameters remain real.
1°) "The" Hypergeometric Function 2F1
-This program calculates F( a,b,c; x+i.y ) for | x + i.y
| < 1
Formula: With z = x + i.y
F(a,b;c; z) = 1 + (a)1(b)1/(c)1.
z1/1! + ............. + (a)n(b)n/(c)n
. zn/n! + .......... where (a)n = a(a+1)(a+2)
...... (a+n-1) and | z | <
1
Data Registers: R00: temp ( Registers R01 R02 R03 are to be initialized before executing "2F1Z" )
• R01 = a
• R02 = b
R04 to R07: temp
• R03 = c
R08 & R09 = F(z)
Flags: /
Subroutine: "Z*Z" ( cf "Complex Functions
for the HP-41" )
| 01 LBL "2F1Z"
02 STO 04 03 X<>Y 04 STO 05 05 CLX 06 STO 00 07 STO 07 08 STO 09 09 SIGN 10 STO 06 11 STO 08 12 LBL 01 |
13 RCL 07
14 RCL 06 15 RCL 05 16 RCL 04 17 XEQ "Z*Z" 18 RCL 00 19 RCL 01 20 + 21 ST* Z 22 * 23 RCL 00 24 RCL 02 |
25 +
26 ST* Z 27 * 28 RCL 03 29 RCL 00 30 ST+ Y 31 ISG X 32 CLX 33 STO 00 34 * 35 ST/ Z 36 / |
37 STO 06
38 X<>Y 39 STO 07 40 RCL 09 41 + 42 STO 09 43 LASTX 44 - 45 ABS 46 X<>Y 47 RCL 08 48 + |
49 STO 08
50 LASTX 51 - 52 ABS 53 + 54 X#0? 55 GTO 01 56 RCL 09 57 RCL 08 58 END |
( 77 bytes / SIZE 010 )
| STACK | INPUTS | OUTPUTS |
| Y | y | B |
| X | x | A |
With F ( x+i.y) = A + i.B
Example: Compute F( 0.4 , 0.6 ; 1.7 ; 0.2 + 0.3 i )
0.4 STO 01 0.6 STO 02 1.7 STO 03
0.3 ENTER^
0.2 XEQ "2F1Z" >>>>
1.023584797
---Execution time = 54s---
X<>Y 0.049325384
-So, F( 0.4 , 0.6 ; 1.7 ; 0.2 + 0.3 i ) = 1.023584797 + 0.049325384 i
Notes:
-Execution time is reduced to 23 seconds if you replace line 17
XEQ "Z*Z" by an M-Code routine Z*Z and lines 31-32 by
X+1
-Similar programs may of course be written to calculate complex Kummer's
functions,
or generalized hypergeometric functions of complex variables.
-But an M-code routine is by far preferable:
2°) Generalized Hypergeometric Functions
mFp
( M-Code Routine )
-With z = x + i.y , HGFZ calculates:
mFp( a1,a2,....,am ; b1,b2,....,bp ; z ) = SUMk=0,1,2,..... [(a1)k(a2)k.....(am)k] / [(b1)k(b2)k.....(bp)k] . zk/k!
where (ai)k
= ai(ai+1)(ai+2) ...... (ai+k-1)
& (ai)0 = 1 ,
likewise for (bj)k
Warning: HGFZ checks that register
Rm+p exists but it does not check for alpha data !
Data Registers: R00 = unused ( Registers R01 thru Rm+p are to be initialized before executing "HGFZ" )
• R01 = a1 • R02 = a2
....................... • Rmm = am
• Rm+1 = b1 • Rm+2 = b2 ....................
• Rm+p = bp
Flags: /
Subroutines: /
09A "Z"
006 "F"
007 "G"
008 "H"
2A0 SETDEC
3C8 CLRKEY
046 C=0 S&X
C
270 RAMSLCT =
038 READDATA T
05E C=0 MS
C = | C |
070 N=C ALL
10E A=C ALL
078 READ 1(Z)
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
268 WRIT 9(Q)
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)
168 WRIT 5(M)
028 WRIT 0(T)
04E C=0 ALL
0A8 WRIT 2(Y)
1E8 WRIT 7(O)
228 WRIT 8(P)
35C PT= 12
050 LD@PT- 1
0E8 WRIT 3(X)
1A8 WRIT 6(N)
378 READ 13(c)
Loop1
03C RCR 3
070 N=C ALL
278 READ 9(Q)
Loop2
03C RCR 3
10E A=C ALL
0B0 C=N ALL
260 SETHEX
226 C=C+1 S&X
306 ?A<C S&X
0CF JC+19h=25d
070 N=C ALL
270 RAMSLCT
038 READDATA
10E A=C ALL
046 C=0 S&X
270 RAMSLCT
2A0 SETDEC
238 READ 8(P)
01D ?NCXQ
C=
060 1807
A+C
10E A=C ALL
068 WRIT 1(Z)
1B8 READ 6(N)
135 ?NCXQ
C=
060 184D
A*C
1A8 WRIT 6(N)
078 READ 1(Z)
10E A=C ALL
1F8 READ 7(O)
135 ?NCXQ
C=
060 184D
A*C
1E8 WRIT 7(O)
313 JNC -1Eh=30d
2F3 JNC -22h=34d
278 READ 9(Q)
10E A=C ALL
0B0 C=N ALL
260 SETHEX
226 C=C+1 S&X
070 N=C ALL
306 ?A<C S&X
0C7 JC+18h=24d
270 RAMSLCT
038 READDATA
10E A=C ALL
046 C=0 S&X
270 RAMSLCT
2A0 SETDEC
238 READ 8(P)
01D ?NCXQ
C=
060 1807
A+C
10E A=C ALL
068 WRIT 1(Z)
1B8 READ 6(N)
0AE A<>C ALL
261 ?NCXQ
C=
060 1898
A/C
1A8 WRIT 6(N)
1F8 READ 7(O)
10E A=C ALL
078 READ 1(Z)
261 ?NCXQ
C=
060 1898
A/C
1E8 WRIT 7(O)
313 JNC -1Eh=30d
2A0 SETDEC
00E A=0 ALL
A
35C PT=12
=
162 A=A+1@PT 1
238 READ 8(P)
01D ?NCXQ
C=
060 1807
A+C
228 WRIT 8(P)
138 READ 4(L)
10E A=C ALL
1B8 READ 6(N)
135 ?NCXQ
C=
060 184D
A*C
013 JNC +02
293 JNC -2Eh=46d
068 WRIT 1(Z)
178 READ 5(M)
10E A=C ALL
1F8 READ 7(O)
135 ?NCXQ
C=
060 184D
A*C
2BE C=-C-1 MS
10E A=C ALL
078 READ 1(Z)
01D ?NCXQ
C=
060 1807
A+C
10E A=C ALL
238 READ 8(P)
261 ?NCXQ
C=
060 1898
A/C
070 N=C ALL
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)
284 CLRF 7
36E ?A#C ALL
013 JNC+02
288 SETF 7
138 READ 4(L)
10E A=C ALL
1F8 READ 7(O)
135 ?NCXQ
C=
060 184D
A*C
068 WRIT 1(Z)
178 READ 5(M)
10E A=C ALL
1B8 READ 6(N)
135 ?NCXQ
C=
060 184D
A*C
10E A=C ALL
078 READ 1(Z)
01D ?NCXQ
C=
060 1807
A+C
10E A=C ALL
238 READ 8(P)
261 ?NCXQ
C=
060 1898
A/C
1E8 WRIT 7(O)
10E A=C ALL
0B8 READ 2(Y)
01D ?NCXQ
C=
060 1807
A+C
10E A=C ALL
0B0 C=N ALL
1A8 WRIT 6(N)
0B8 READ 2(Y)
0AE A<>C ALL
0A8 WRIT 2(Y)
3CC ?KEY
Press any key
360 ?C RTN
to stop an infinite loop
36E ?A#C ALL
217 JC -3Eh=-62d
28C ?FSET 7
207 JC -40h=-64d
345 ?NCGO
?ncgo
042 10D1
CLA
| STACK | INPUTS | OUTPUTS |
| T | m | y |
| Z | p | / |
| Y | y | B |
| X | x | A |
| L | / | x |
where mFp( x + i.y ) = A + i.B
Example: The same one: 2F1( 0.4 , 0.6 ; 1.7 ; 0.2 + 0.3 i )
0.4 STO 01 0.6 STO 02 1.7 STO 03
2 ENTER^
1 ENTER^
0.3 ENTER^
0.2 XEQ "HGFZ" >>>>
1.023584797
---Execution time = 13s---
X<>Y 0.049325384
-Thus, 2F1( 0.4 , 0.6 ; 1.7 ; 0.2 + 0.3 i ) = 1.023584797 + 0.049325384 i
Notes:
-Of course, this function may be calculated for itself.
-But several formulae may be re-written to compute kelvin functions,
asymptotic expansions...
via hypergeometric functions of complex variables.
-A few examples are given below:
3°) Application to a few Special Functions
a) Bessel Functions:
Asymptotic Expansions
Formulae:
Jn(x) = (2/(pi*x))1/2
( P.cos(x-(2n+1)pi/4) - Q.sin(x-(2n+1)pi/4) )
Yn(x) = (2/(pi*x))1/2
( P.sin(x-(2n+1)pi/4) + Q.cos(x-(2n+1)pi/4) )
With P + i Q = 2F0( (1+2n)/2
, (1-2n)/2 ;; i / (2x) )
Data Registers: R00 = x , R01-R02-R03:
temp
Flags: /
Subroutines: /
| 01 LBL "AEJY"
02 DEG 03 STO 00 04 X<>Y 05 STO 01 06 CHS 07 STO 02 08 .5 09 ST+ 01 |
10 ST+ 02
11 1/X 12 0 13 RCL 00 14 ST+ X 15 1/X 16 0 17 HGFZ 18 STO 03 |
19 CLX
20 RCL 00 21 R-D 22 RCL 01 23 90 24 * 25 - 26 STO 02 27 X<>Y |
28 P-R
29 RCL 02 30 RCL 03 31 P-R 32 ST+ T 33 RDN 34 X<>Y 35 - 36 RCL 00 |
37 PI
38 * 39 2 40 / 41 SQRT 42 ST/ Z 43 / 44 X<>Y 45 END |
( 62 bytes / SIZE 004 )
| STACK | INPUTS | OUTPUTS |
| Y | n | Yn(x) |
| X | x > 0 | Jn(x) |
Example:
PI ENTER^
11.2 XEQ "AEJY" >>>>
JPI(11.2) = 0.226432126
---Execution time = 19s---
X<>Y YPI(11.2) = -0.088603991
Notes:
-The program will not work if x is too "small" or if n is too
"large"
-With n = PI, an infinite loop occurs if x = 11.1 ( press any
key to stop the loop )
b) Fresnel Integrals:
Ascending Series & Asymptotic Series
-Fresnel integrals C(x) & S(x) may be computed by
C(x) + i S(x) = x 1F1( 1/2 ; 3/2
; i PI x2/2 )
Data Registers: R00 = x , R01-R02:
temp
Flags: /
Subroutines: /
| 01 LBL "CSX"
02 1 03 X<>Y 04 STO 00 05 X^2 06 PI 07 .5 08 STO 01 09 STO 02 10 * 11 * 12 ISG 02 13 0 14 HGFZ 15 RCL 00 16 ST* Z 17 * 18 END |
( 30 bytes / SIZE 003 )
| STACK | INPUTS | OUTPUTS |
| Y | / | S(x) |
| X | x | C(x) |
Example:
1.4 XEQ "CSX" >>>>
C(1.4) = 0.543095784
---Execution time = 10s---
X<>Y S(1.4) = 0.713525077
Notes:
-x must remain "small", say x < 2
-For x = 3 , the errors are of the order of 10 -6
and the results are meaningless with x = 4
>>> For large arguments, we can use asymptotic series which may be re-written:
C(x) + i S(x) ~ ( 1 + i ) / 2 -
( i / (x.PI) ) exp ( i PI x2 / 2 ) 2F0(
1 , 1/2 ;; -2 i / ( PI x2 ) )
Data Registers: R00 = x , R01-R02:
temp
Flags: /
Subroutines: /
| 01 LBL "AECS"
02 .5 03 STO 02 04 1/X 05 STO Z 06 X<>Y 07 STO 00 08 X^2 |
09 PI
10 * 11 / 12 CHS 13 1 14 STO 01 15 CLX 16 STO Z |
17 HGFZ
18 R-P 19 X<>Y 20 RCL 00 21 X^2 22 RCL 01 23 ASIN 24 * |
25 +
26 X<>Y 27 P-R 28 PI 29 RCL 00 30 * 31 ST/ Z 32 / |
33 CHS
34 X<>Y 35 RCL 02 36 ST+ Z 37 + 38 END |
( 53 bytes / SIZE 003 )
| STACK | INPUTS | OUTPUTS |
| Y | / | S(x) |
| X | x | C(x) |
Example:
4 XEQ "AECS" >>>> C(4)
= 0.498426033
---Execution time = 11s---
X<>Y S(4) = 0.420515754
Note:
-Divergence occurs if x = 3.9
c) Kelvin Functions:
Ascending Series & Asymptotic Series
-The Kelvin functions of the 1st kind may be expressed via Bessel functions of the 1st kind:
bern(x) + i bein(x) = Jn( x exp(3i(PI)/4) ) whence
bern(x) + i bein(x) = ( x (i-1)/sqrt(8)
)n 0F1( n+1 ; i x2/4
) / Gamma(n+1)
Data Registers: R00 = x , R01-R02-R03:
temp
Flags: /
Subroutines: "Z^X" "Z*Z" ( cf "Complex
Functions for the HP-41" ) and "1/G" ( cf "Gamma Function
for the HP-41" )
| 01 LBL "KLV1"
02 STO 00 03 X<>Y 04 STO 02 05 1 06 STO T 07 + 08 STO 01 |
09 CLX
10 2 11 / 12 X^2 13 0 14 STO T 15 HGFZ 16 STO 03 |
17 CLX
18 RCL 00 19 8 20 SQRT 21 / 22 ENTER^ 23 CHS 24 RCL 02 |
25 XEQ "Z^X"
26 R^ 27 RCL 03 28 XEQ "Z*Z" 29 STO 02 30 X<>Y 31 STO 03 32 RCL 01 |
33 XEQ "1/G"
34 ST* 02 35 ST* 03 36 RCL 03 37 RCL 02 38 END |
( 61 bytes / SIZE 004 )
| STACK | INPUTS | OUTPUTS |
| Y | n | bein(x) |
| X | x | bern(x) |
Example:
2 SQRT
PI XEQ "KLV1" >>>> -0.674095953
---Execution time = 9s---
X<>Y -1.597357212
-Thus, bersqrt(2) (PI) = -0.674095953 & beisqrt(2) (PI) = -1.597357212
Note:
-Lines 28 to 37 may be replaced by Z*Z
RCL 01 1/G ST* Z * if you
have the M-Code routines Z*Z and 1/G
>>> The asymptotic series return the 4 functions bern(x)
, bein(x) , kern(x) , kein(x). We
have:
Formulae: if x > 0
bern(x) = exp(x/sqrt(2)) (2.PI.x) -1/2
[ fn(x) cos A + gn(x) sin A ] - (1/PI) [ sin (360°.n)
kern(x) + cos (360°.n) kein(x) ]
bein(x) = exp(x/sqrt(2)) (2.PI.x) -1/2
[ fn(x) sin A - gn(x) cos A ] + (1/PI) [ cos (360°.n)
kern(x) - sin (360°.n) kein(x) ]
kern(x) = exp(-x/sqrt(2)) (PI/(2x))1/2
[ fn(-x) cos B - gn(-x) sin B ]
with A = (180°/PI) x/sqrt(2) + ( n/2 - 1/8 ) 180°
kein(x) = exp(-x/sqrt(2)) (PI/(2x))1/2
[ -fn(-x) sin B - gn(-x) cos B ]
B = A + 45°
with
fn(x) + i gn(x)
~ 2F0( (1-2n)/2 , (1+2n)/2 ;; exp(i(PI)/4)
/ ( 2x) )
fn(-x) + i gn(-x) ~
2F0(
(1-2n)/2 , (1+2n)/2 ;; - exp(i(PI)/4) / ( 2x) )
Data Registers: R00 = x , R01
to R09: temp
Flags: /
Subroutines: /
| 01 LBL "AEKV"
02 DEG 03 STO 00 04 X<>Y 05 STO 01 06 STO 03 07 CHS 08 STO 02 09 .5 10 ST+ 01 11 ST+ 02 12 1/X 13 0 14 RCL 00 15 8 16 SQRT 17 * 18 1/X 19 ENTER^ 20 HGFZ 21 STO 04 22 X<>Y 23 STO 05 |
24 2
25 0 26 LASTX 27 CHS 28 ENTER^ 29 HGFZ 30 STO 06 31 X<>Y 32 STO 07 33 RCL 00 34 2 35 SQRT 36 / 37 STO 01 38 RCL 03 39 2 40 / 41 8 42 1/X 43 + 44 PI 45 * 46 + |
47 R-D
48 STO 02 49 RCL 06 50 P-R 51 RCL 02 52 RCL 07 53 P-R 54 ST+ T 55 RDN 56 - 57 PI 58 RCL 00 59 ST+ X 60 / 61 SQRT 62 ENTER^ 63 X<> 01 64 E^X 65 ST* 01 66 / 67 ST* Z 68 * 69 STO 08 |
70 X<>Y
71 CHS 72 STO 09 73 45 74 ST- 02 75 8 76 * 77 ST* 03 78 RCL 03 79 RCL 08 80 P-R 81 RCL 03 82 RCL 09 83 P-R 84 ST+ T 85 RDN 86 - 87 STO 07 88 X<>Y 89 STO 06 90 RCL 02 91 RCL 04 92 P-R |
93 RCL 02
94 RCL 05 95 P-R 96 ST- T 97 RDN 98 + 99 RCL 01 100 ST* Z 101 * 102 RCL 06 103 - 104 X<>Y 105 RCL 07 106 + 107 PI 108 ST/ Z 109 / 110 RCL 08 111 RCL 09 112 RDN 113 X<> Z 114 END |
( 142 bytes / SIZE 010 )
| STACK | INPUTS | OUTPUTS |
| T | / | kein(x) |
| Z | / | kern(x) |
| Y | n | bein(x) |
| X | x > 0 | bern(x) |
Example:
PI ENTER^
11 XEQ "AEKV" >>>>
berPI(11) = 209.8380891
---Execution time = 19s---
RDN beiPI(11) = 17.30058118
RDN kerPI(11) = 0.0001458703398
RDN keiPI(11) = -0.0001588289186
Note:
-There is an infinite loop with n = PI & x =10.9
d) Sine & Cosine
Integrals: Asymptotic Series
Formulae:
Si(x) = pi/2 - f(x) cos x - g(x) sin x
with f(x) + i g(x) ~ (1/x) 2F0(
1 , 1 ;; i/x )
Ci(x) = f(x) sin x - g(x) cos x
Data Registers: R00 = x , R01-R02:
temp
Flags: /
Subroutines: /
| 01 LBL "AECSI"
02 DEG 03 STO 00 04 1 05 STO 01 06 STO 02 07 ST+ X |
08 0
09 RCL 00 10 1/X 11 0 12 HGFZ 13 RCL 00 14 ST/ Y |
15 ST/ Z
16 R-D 17 STO T 18 X<>Y 19 P-R 20 R^ 21 R^ |
22 P-R
23 ST- T 24 RDN 25 + 26 PI 27 2 28 / |
29 X<>Y
30 - 31 X<>Y 32 END |
( 48 bytes / SIZE 003 )
| STACK | INPUTS | OUTPUTS |
| Y | / | Si(x) |
| X | x | Ci(x) |
Example:
26.3 XEQ "AECSI" >>>>
Ci(26.3) = 0.0343063916
---Execution time = 14s---
X<>Y Si(26.3) = 1.554589810
Note:
-There is an infinite loop with x =26.2
Conclusion:
-In most cases - when it is possible - HGFZ allows to write faster &
shorter programs than using HGF+ twice.
References:
[1] Abramowitz and Stegun - "Handbook of Mathematical Functions"
- Dover Publications - ISBN 0-486-61272-4
[2] http://functions.wolfram.com/