HypergeoComplex

# 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      =
05E   C=0 MS          C = | C |
070  N=C ALL
10E  A=C ALL
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
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"
128   WRIT 4(L)
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)
03C  RCR 3
070   N=C ALL
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
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
2A0  SETDEC
01D   ?NCXQ           C=
060   1807               A+C
10E   A=C ALL
068   WRIT 1(Z)
135   ?NCXQ            C=
060   184D                A*C
1A8   WRIT 6(N)
10E   A=C ALL
135   ?NCXQ            C=
060   184D                A*C
1E8   WRIT 7(O)
313   JNC -1Eh=30d
2F3   JNC -22h=34d
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
10E   A=C ALL
046   C=0 S&X
270   RAMSLCT
2A0  SETDEC
01D   ?NCXQ           C=
060   1807               A+C
10E   A=C ALL
068   WRIT 1(Z)
0AE  A<>C ALL
261   ?NCXQ             C=
060   1898                 A/C
1A8   WRIT 6(N)
10E   A=C ALL
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
01D   ?NCXQ           C=
060   1807               A+C
228   WRIT 8(P)
10E   A=C ALL
135   ?NCXQ            C=
060   184D                A*C
013   JNC +02
293   JNC -2Eh=46d
068   WRIT 1(Z)
10E   A=C ALL
135   ?NCXQ            C=
060   184D                A*C
2BE   C=-C-1 MS
10E   A=C ALL
01D   ?NCXQ           C=
060   1807               A+C
10E   A=C ALL
261   ?NCXQ             C=
060   1898                 A/C
070   N=C ALL
10E   A=C ALL
01D   ?NCXQ           C=
060   1807               A+C
10E   A=C ALL
0AE  A<>C ALL
0E8   WRIT 3(X)
284   CLRF 7
36E   ?A#C ALL
013   JNC+02
288   SETF 7
10E   A=C ALL
135   ?NCXQ            C=
060   184D                A*C
068   WRIT 1(Z)
10E   A=C ALL
135   ?NCXQ            C=
060   184D                A*C
10E   A=C ALL
01D   ?NCXQ           C=
060   1807               A+C
10E   A=C ALL
261   ?NCXQ             C=
060   1898                 A/C
1E8   WRIT 7(O)
10E   A=C ALL
01D   ?NCXQ           C=
060   1807               A+C
10E   A=C ALL
0B0   C=N ALL
1A8   WRIT 6(N)
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/