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