hp41programs

Gamma2

Gamma Functions(2) for the HP-41


Overview
 

 1°)  Multivariate Gamma Function
 2°)  Multivariate Digamma Function
 3°)  Q-Gamma Function
 4°)  Elliptic Gamma Function ( focal program & M-code routine )
 5°)  Hyperbolic Gamma Function
 

-A few programs that generalize the usual Gamma function.
 

1°)  Multivariate Gamma Function
 

-This function may be computed by

   Gm (a) = PI m(m-1)/4   Prodj=1,2,...,m  Gam ( a - ( m-1)/2 )    where m is a positive integer and  Gam is the usual Gamma function.
 

Data Registers:   R00 = m
                              R01 = a   R02 = Gm(a)   R03: temp
Flags: /
Subroutine:   GAM ( M-Code Routine )  It may be replace by a focal program "GAM" ( cf "Gamma Function for the HP41" )
                                                                   In this case, shift R00 to R03 in the listing below to avoid registers coflicts...
 
 

 01 LBL "MVG"
 02 STO 01
 03 X<>Y
 04 STO 00
 05 STO 03
 06 1
 07 STO 02
 08 LBL 01
 09 RCL 01
 10 RCL 03
 11 1
 12 -
 13 2
 14 /
 15 -
 16 GAM
 17 ST* 02
 18 DSE 03
 19 GTO 01
 20 PI
 21 RCL 00
 22 RCL 00
 23 1
 24 -
 25 *
 26 4
 27 /
 28 Y^X
 29 RCL 02
 30 *
 31 STO 02
 32 END

 
( 44 bytes / SIZE 004 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             m            /
           X             a         Gm(a)

 
Example:    Evaluate   G4 ( 1.23 )

     4     ENTER^
   1.23  XEQ "MVG"   >>>>   G4 ( 1.23 ) = - 650.1892181                   ---Execution time = 9s---

Note:

"MVG" will run slower if you use a focal "GAM" subroutine.
 

2°)  Multivariate Digamma Function
 

-Similarly, the multivariate digamma function is defined by:

    Psim (a) = SUMj=1,2,...,m  Psi ( a + (1-m)/2 )
 

Data Registers:   R00 = m
                              R01 = a   R02 = Psim(a)
Flags: /
Subroutine:   PSI ( M-Code Routine )     It may be replace by a focal program "PSI" ( cf "Gamma Function for the HP41" )
                                                                   In this case, shift R00 to R02 in the listing below to avoid registers coflicts...
 
 

 01 LBL "MVP"
 02 STO 01
 03 X<>Y
 04 STO 00
 05 CLX
 06 STO 02
 07 LBL 01
 08 RCL 01
 09 1
 10 RCL 00
 11 -
 12 2
 13 /
 14 +
 15 PSI
 16 ST+ 02
 17 DSE 00
 18 GTO 01
 19 RCL 02
 20 END

 
( 32 bytes / SIZE 003 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             m            /
           X             a         Gm(a)

 
Example:    Evaluate   G4 ( PI )

   4    ENTER^
  PI   XEQ "MVP"   >>>>   G4 ( PI ) = 2.418961631                    ---Execution time = 7s---

Note:

-Of course, "MVP" will run slower if you use a focal "PSI" subroutine.
 

3°)  Q-Gamma Function
 

-This function is defined by an infinite product:
 

  Gq ( x )  =  ( 1 - q )1-x  ¶n=0,1,2,.... ( 1 - qn+1 ) / ( 1 - qn+x )
 

Data Registers:   R00 = q
                              R01 = x   R02 = qx   R03 = Gq ( x )        R04: temp
Flags: /
Subroutines: /
 
 

 01 LBL "QGAM"
 02 STO 01
 03 X<>Y
 04 STO 00
 05 X<>Y
 06 Y^X
 07 STO 02
 08 1
 09 STO 04
 10 RCL 00         
 11 -
 12 1
 13 RCL 01
 14 -
 15 Y^X
 16 STO 03 
 17 LBL 01
 18 RCL 02         
 19 RCL 00
 20 RCL 04
 21 ST* Z
 22 *
 23 STO 04         
 24 X<>Y
 25 1
 26 ST- Z
 27 -
 28 /
 29 RCL 03
 30 *
 31 STO 03         
 32 LASTX
 33 X#Y?
 34 GTO 01
 35 END

 
     ( 47 bytes / SIZE 005 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             q        Gq ( x )
           X             x        Gq ( x )

 
Example:    Evaluate   G0.3 ( 3.14 )

     0.3   ENTER^
    3.14  XEQ "QGAM"  >>>>   G0.3 ( 3.14 ) = 1.358251666                  ---Execution time = 12s---

Note:

-We assume 0 < q < 1
-As q tends to 1 , Gq ( x ) tends to the usual Gamma(x)
 

4°)  Elliptic Gamma Function
 

  Gell ( x , y , z )  is defined by the double infinite product    ¶m=0,1,2,....  ¶n=0,1,2,....   ( 1 - ym+1 zn+1 / x ) / ( 1 - ym zn  x )
 
 

Data Registers:   R00 = x
                              R01 = a   R02 = b   R03 = Gell ( x , y , z )        R04-R05: temp
Flags: /
Subroutines: /
 
 

 01 LBL "ELG"
 02 STO 00    
 03 RDN
 04 STO 01
 05 X<>Y
 06 STO 02
 07 1
 08 STO 03
 09 STO 04
 10 STO 05
 11 LBL 01
 12 RCL 01
 13 RCL 02
 14 *
 15 RCL 04    
 16 RCL 05
 17 *
 18 ST* Y
 19 RCL 00 
 20 X^2
 21 *
 22 RCL 00
 23 ST- Z
 24 -
 25 /
 26 RCL 03    
 27 *
 28 STO 03
 29 RCL 01 
 30 ST* 04
 31 X<> L
 32 X#Y?
 33 GTO 01
 34 RCL 01
 35 STO 04
 36 RCL 05    
 37 RCL 02
 38 ST* Z
 39 *
 40 STO 05 
 41 ST* Y
 42 RCL 00
 43 X^2
 44 *
 45 RCL 00
 46 ST- Z
 47 -
 48 /
 49 RCL 03    
 50 *
 51 STO 03 
 52 LASTX
 53 X#Y?
 54 GTO 01
 55 END

 
      ( 72 bytes / SIZE 006 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             z             /
           Y             y   Gell ( x , y , z )
           X             x   Gell ( x , y , z )

 
Example:    Evaluate   Gell ( 12 , 0.1 , 0.2 )

    0.2  ENTER^
    0.1  ENTER^
     12  XEQ "ELG"  >>>>   -1.176280598                             ---Execution time = 56s---
 

Note:

-The product will diverge for many arguments.
-Even if it converges, execution time may be prohibitive without a good emulator...

-With V41 you'll get:    Gell ( 0.8 , 0.6 , 0.7 ) = 31.41284778         and      Gell ( 1.28 , 0.6 , 0.7 ) = -74963.68445
-With free42, it yields   Gell ( 0.8 , 0.6 , 0.7 ) = 31.4128471558     and      Gell ( 1.28 , 0.6 , 0.7 ) = -74963.6802183

-You can also add  VIEW X  after line 31 or 56  to follow the convergence.

-Here is an M-Code version of the same program:
-There is no check for alpha data and the alpha register is cleared
 

087  "G"
00C  "L"
005  "E"
2A0  SETDEC
0F8  C=X
128   L=C
04E  C
35C  =
050   1
168  M=C
1A8  N=C
0E8  X=C
1B8  C=N    Loop
10E  A=C ALL
178  C=M
135  C=
060  A*C
070  N=C ALL    ( cpu register N )
138  C=L
13D  C=
060  AB*C
009  C=
060  AB-1
0F0  C<>N ALL
10E  A=C ALL
078  C=Z
135  C=
060  A*C
0B8  C=Y
13D  C=
060  AB*C
138  C=L
269  C=
060  AB/C
009  C=
060  AB-1
0B0  C=N ALL
269  C=
060  AB/C
0F8  C=X
13D  C=
060  AB*C
10E  A=C ALL
0F8  C=X
0AE  A<>C ALL
0E8  X=C
36E  ?A#C ALL
043  JNC+08
178  C=M
10E  A=C ALL
0B8  C=Y
135  C=
060  A*C
168  M=C ALL
2B3  JNC-42d             goto loop
0B8  C=Y
168  M=C ALL
078  C=Z
10E  A=C ALL
1B8  C=N
135  C=
060  A*C
1A8  N=C
070  N=C ALL
138  C=L
13D  C=
060  AB*C
009  C=
060  AB-1
0F0  C<>N ALL
10E  A=C ALL
078  C=Z
135  C=
060  A*C
0B8  C=Y
13D  C=
060  AB*C
138  C=L
269  C=
060  AB/C
009  C=
060  AB-1
0B0  C=N ALL
269  C=
060  AB/C
0F8  C=X
13D  C=
060  AB*C
10E  A=C ALL
0F8  C=X
0AE  A<>C ALL
0E8  X=C
36E  ?A#C ALL
2CF  JC-39d
345  ?NCGO
042  CLA

( 96 words )
 
 

      STACK        INPUTS      OUTPUTS
           T             T             T
           Z             z             z
           Y             y             y
           X             x   Gell ( x , y , z )
           L             L             x

 
Example:    Evaluate   Gell ( 0.1 , 0.2 , 0.3 )

    0.3  ENTER^
    0.2  ENTER^
    0.1  XEQ "ELG"   >>>>     Gell ( 0.1 , 0.2 , 0.3 ) = 0.291691838               ---Execution time = 48s---

Note:

-Press ENTER^  ON   to stop a possible infinite loop or if execution time is too long.
 

5°)  Hyperbolic Gamma Function
 

   Ghyp ( a , b , x )  =  Exp ( i.J )   where   J = §0infinity  [ Sin ( 2x.y ) / ( 2 Sinh ( a.y ) ) / Sinh ( b.y )  - x / ( a.b.y ) ]  dy / y

-So, Ghyp ( a , b , x )  is a complex number.

 "HYG" evaluates the integral by Gauss-Legendre 3 point formula.
 
 

Data Registers:   R00 = function name

                              R01 = 0       R03 = nbre of subintervals    R05 to R08 are used by "GL3"                 R10 = x   R11 = a   R12 = b
                              R02 = Pi/2   R04 = Integral                      R09 is unused
Flag:  F22
Subroutines:  "GL3"  ( cf "Numerical Integration for the HP41" )   or another integrator...
 
 
 

 01 LBL "HYG"
 02 "A^B^X"
 03 PROMPT
 04 STO 10
 05 RDN
 06 STO 12
 07 X<>Y
 08 STO 11
 09 RAD
 10 CLX
 11 STO 01
 12 PI
 13 2
 14 /
 15 STO 02
 16 CF 22
 17 "N=?"
 18 PROMPT
 19 FC?C 22
 20 1
 21 STO 03
 22 "T"
 23 ASTO 00
 24 LBL 10
 25 XEQ "GL3"
 26 RTN
 27 "N=?"
 28 PROMPT
 29 FC?C 22
 30 GTO 00
 31 STO 03
 32 GTO 10
 33 LBL "T"
 34 TAN
 35 ENTER
 36 ENTER
 37 ENTER
 38 ST+ X
 39 RCL 10      
 40 *
 41 SIN
 42 ST+ X
 43 X<>Y
 44 RCL 11
 45 *
 46 E^X-1
 47 LASTX
 48 CHS
 49 E^X-1
 50 -
 51 /
 52 X<>Y
 53 RCL 12      
 54 *
 55 E^X-1
 56 LASTX
 57 CHS
 58 E^X-1
 59 -
 60 /
 61 RCL 10
 62 RCL 11
 63 /
 64 RCL 12      
 65 /
 66 R^
 67 /
 68 -
 69 X<>Y
 70 ST/ Y
 71 X^2
 72 ENTER
 73 SIGN
 74 +
 75 *
 76 RTN
 77 LBL 00
 78 RCL 04      
 79 1
 80 P-R
 81 DEG
 82 END

 
    ( 119 bytes / SIZE 013 )
 
 

      STACK        INPUTS      OUTPUT'        ............    OUTPUT''''     FINALLY
           Z             a             /             /            /            /
           Y             b             /             /            /  Im Ghyp( a,b,x )
           X             x             J             J            J  Re Ghyp( a,b,x )

 
Example:    Evaluate   Ghyp ( 1 , 2 , 3 )    in other words,  a = 1   b = 2   x = 3

    XEQ "HYG"   >>>>    "A^B^X"

    1  ENTER^
    2  ENTER^
    3    R/S           >>>>    "N=?"                    enter an N-value ( for instance 2 ) or the HP41 will take N = 1  if you press R/S without any digit entry.

    2    R/S           >>>>    -7.383610
          R/S           >>>>    "N=?"                     if you choose N = 8

    8    R/S           >>>>    -7.395757
          R/S           >>>>    "N=?"

   16   R/S           >>>>    -7.395793
          R/S           >>>>    "N=?"                    simply press R/S without any digit entry if you think the precision is acceptable

          R/S           >>>>     0.442325
                          X<>Y    -0.896855

-So   Ghyp ( 1 , 2 , 3 )  =  0.442325 - 0.896855 i

Note:

-Since the HP41 cannot deal with number > E100 , you can get an OUTOFRANGE if n is too large.
-In the example above, this happens with n = 21

-So the precision may be limited...