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...