hp41programs

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      =
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/