The Museum of HP Calculators
The Museum of HP Calculators
This program is Copyright © 2006-2007 by Jean-Marc Baillard and is used here by permission.
This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.
Overview
1°) Gamma Function
a) Program#1
b) Program#2 ( more accurate results
)
c) Reciprocal of the Gamma Function
d) Reciprocal of the Gamma Function (
program#2 )
e) Logarithm of the Gamma Function
f) Gamma Function in the complex plane
2°) Digamma Function
a) Real Arguments
b) Complex Arguments
3°) Polygamma Functions
4°) Incomplete Gamma Function
a) Via Kummer's Function, program#1
b) Via Kummer's Function, program#2
(
new )
c) A Continued Fraction
5°) Beta Function
a) Program#1
b) Positive Integer Arguments
6°) Incomplete Beta Function
a) Via Hypergeometric Functions
b) A Continued Fraction
1°) Gamma Function
a) Program#1
-The asymptotic formula e-x xx-1/2 (2pi)1/2
exp ( 1 + 1/(12x) -1/(360x3) ) is used if x > 16
-Otherwise, an integer n is added to x such that x+n >
16 and this formula is then applied to x+n
with the relation: Gam(x+n) = (x+n-1) (x+n-2) ...... (x+1)
x Gam(x)
-Gam(x) is defined if x # 0 ; -1 ; -2 ; .......
-Synthetic register M may be replaced by any unused register.
Data Registers: /
Flags: /
Subroutines: /
01 LBL "GAM"
02 FRC
03 X#0?
04 GTO 00
05 LASTX
06 1
07 -
08 FACT
09 RTN
10 LBL 00
11 16
12 LASTX
13 STO M
14 GTO 02
15 LBL 01
16 1
17 +
18 ST* M
19 LBL 02
20 X<Y?
21 GTO 01
22 ENTER^
23 X^2
24 SIGN
25 LASTX
26 30
27 *
28 1/X
29 -
30 X<>Y
31 12
32 *
33 /
34 E^X
35 0
36 X<> M
37 /
38 X<>Y
39 1
40 E^X
41 /
42 R^
43 Y^X
44 *
45 X<>Y
46 PI
47 *
48 ST+ X
49 SQRT
50 *
51 END
( 69 bytes / SIZE 000 )
| STACK | INPUT | OUTPUT |
| X | x | Gam(x) |
Example: Compute Gam(PI) and Gam (-6.14)
PI XEQ "GAM" yields
2.288037783 ( in 5 seconds )
-6.14 R/S
------ -0.007872567 ( in 7 seconds
)
b) Program#2
-In order to avoid a loss of significant digits if x < 10 ,
the following routine uses more terms of the asymptotic series.
Data Registers: /
Flags: /
Subroutines: /
01 LBL "GAM+"
02 FRC
lines 02 to 10 take advantage of the built in FACT function
03 X#0?
but they can be deleted after replacing line 14 by X<>Y
04 GTO 00
05 LASTX
06 1
07 -
08 FACT
09 RTN
10 LBL 00
11 5
12 STO M
13 ST/ M
14 LASTX
15 LBL 01
16 X>Y?
17 GTO 02
18 ST* M
19 1
20 +
21 GTO 01
22 LBL 02
23 ENTER^
24 ENTER^
25 X^2
26 99
27 *
28 1/X
29 140
30 1/X
31 -
32 X<>Y
33 X^2
34 /
35 105
36 1/X
37 +
38 X<>Y
39 X^2
40 /
41 30
42 1/X
43 -
44 X<>Y
45 X^2
46 /
47 1
48 +
49 12
50 /
51 X<>Y
52 ST/ Y
53 -
54 E^X
55 X<>Y
56 R^
57 .5
58 ST- Y
these lines avoid an OUT OF RANGE for x > 57
59 *
60 Y^X
61 ST* Y
62 *
63 PI
64 ST+ X
65 SQRT
66 *
67 0
68 X<> M
69 /
70 END
( 98 bytes / SIZE 000 )
| STACK | INPUT | OUTPUT |
| X | x | Gam(x) |
Example: Compute Gam(PI) and Gam (-6.14)
PI XEQ "GAM+" yields
2.288037795
( in 4 seconds )
-6.14 R/S
------ -0.007872567231 ( in
6 seconds ) ( the last 2 decimals should be 20
)
c) Reciprocal
of the Gamma Function
-It's sometimes worthwhile to use a program which computes 1/Gam(x)
because this function is defined for any x-value.
-Thus, it may avoid tests when x = 0 ; -1 ; -2 ; -3 ; ......
Data Registers: /
Flags: /
Subroutines: /
01 LBL "1/G"
02 STO M
03 16
04 X<>Y
05 GTO 02
06 LBL 01
07 1
08 +
09 ST* M
10 LBL 02
11 X<Y?
12 GTO 01
13 ENTER^
14 ENTER^
15 X^2
16 30
17 *
18 1/X
19 1
20 -
21 X<>Y
22 12
23 *
24 /
25 E^X
26 0
27 X<> M
28 *
29 X<>Y
30 1
31 E^X
32 /
33 R^
34 Y^X
35 /
36 X<>Y
37 PI
38 *
39 ST+ X
40 SQRT
41 /
42 END
( 59 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| X | x | 1 / gam(x) |
Example: PI XEQ "1/G"
>>>> 1/Gam(pi) = 0.437055720 ( in 5 seconds
)
-3 R/S
>>>> 1/Gam(-3) = 0
d) Reciprocal
of the Gamma Function ( Program#2 )
- "1/G+" uses the same method as "GAM+" for more accurate results
Data Registers: /
Flags: /
Subroutines: /
01 LBL "1/G+"
02 5
03 STO M
04 ST/ M
05 X<>Y
06 LBL 01
07 X>Y?
08 GTO 02
09 ST* M
10 1
11 +
12 GTO 01
13 LBL 02
14 ENTER^
15 ENTER^
16 X^2
17 99
18 *
19 1/X
20 140
21 1/X
22 -
23 X<>Y
24 X^2
25 /
26 105
27 1/X
28 +
29 X<>Y
30 X^2
31 /
32 30
33 1/X
34 -
35 X<>Y
36 X^2
37 /
38 1
39 +
40 12
41 CHS
42 /
43 X<>Y
44 ST/ Y
45 +
46 E^X
47 0
48 X<> M
49 *
50 X<>Y
51 R^
52 CHS
53 .5
54 ST+ Y
55 *
56 Y^X
57 ST* Y
58 *
59 PI
60 ST+ X
61 SQRT
62 /
63 END
( 89 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| X | x | 1 / gam(x) |
Example: PI XEQ "1/G+"
>>>> 1/Gam(pi) = 0.4370557174 ( in 4 seconds
)
-3 R/S
>>>> 1/Gam(-3) = 0
( in 5 seconds )
e) Logarithm of
the Gamma Function ( large arguments )
-This program evaluates log | Gam(x) | and it's especially useful if n > 69
Data Registers: /
Flags: /
Subroutines: /
01 LBL "LOGAM"
02 ENTER^
03 ABS
04 LOG
05 STO M
06 16
07 RCL Z
08 LBL 01
09 1
10 +
11 STO Z
12 ABS
13 LOG
14 ST+ M
15 X<> Z
16 X<Y?
17 GTO 01
18 ENTER^
19 X^2
20 SIGN
21 LASTX
22 30
23 *
24 1/X
25 -
26 X<>Y
27 12
28 *
29 /
30 X<>Y
31 -
32 10
33 LN
34 /
35 X<>Y
36 ENTER^
37 LOG
38 *
39 +
40 X<>Y
41 PI
42 *
43 ST+ X
44 SQRT
45 LOG
46 +
47 0
48 X<> M
49 -
50 END
( 72 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| X | x | log | gam(x) | |
Example: 1000 XEQ "LOGAM"
yields log | gam(1000) | = 2564.604644
whence gam(1000) = 4.02387 102564
-If you prefer ln | gam(x) | instead of log
| gam(x) | , delete lines 32-33-34 and replace all the LOG
with LN
f) Gamma Function in
the complex plane
-This program employs the same asymptotic formula as "GAM"
Data Registers: R00 to R05: temp
Flags: /
Subroutine: /
01 LBL "GAMZ"
02 RAD
03 STO 01
04 X<>Y
05 STO 02
06 1
07 STO 03
08 CLX
09 STO 04
10 LBL 01
11 16
12 RCL 01
13 X>Y?
14 GTO 02
15 RCL 02
16 X<>Y
17 R-P
18 ST* 03
19 X<>Y
20 ST+ 04
21 1
22 ST+ 01
23 GTO 01
24 LBL 02
25 RCL 02
26 RCL 01
27 R-P
28 2
29 CHS
30 ST* Z
31 Y^X
32 360
33 CHS
34 /
35 P-R
36 12
37 1/X
38 +
39 R-P
40 RCL 02
41 RCL 01
42 R-P
43 RDN
44 ST- Z
45 X<> T
46 /
47 P-R
48 X<>Y
49 RCL 02
50 -
51 STO 00
52 X<>Y
53 RCL 01
54 -
55 STO 05
56 RCL 02
57 RCL 01
58 .5
59 -
60 R-P
61 RCL 02
62 RCL 01
63 R-P
64 LN
65 R-P
66 RDN
67 ST+ Z
68 X<> T
69 *
70 P-R
71 X<>Y
72 RCL 00
73 +
74 RCL 04
75 -
76 X<>Y
77 RCL 05
78 +
79 E^X
80 RCL 03
81 /
82 PI
83 ST+ X
84 SQRT
85 *
86 P-R
87 DEG
88 END
( 113 bytes / SIZE 006 )
| STACK | INPUTS | OUTPUTS |
| Y | y | B |
| X | x | A |
where Gam(x+iy) = A + iB
Example: Compute Gam(1+2i)
2 ENTER^
1 XEQ "GAMZ" yields
0.151904003
X<>Y
0.019804881 thus,
Gam(1+2i) = 0.151904003 + 0.019804881. i
( in 21 seconds )
Note: -If you prefer to obtain Ln(Gam(x+iy))
, simply replace line 86 by LN
-For instance,
Ln(Gam(2+3i)) = -1.876078781 + 0.129646321 i
2°) Digamma Function
a) Real
Arguments
-The Digamma ( or Psi ) function is defined by Psi(x) =
Gam'(x)/Gam(x) where Gam'(x) is the first derivative
of Gam(x)
-The asymptotic expansion Psi(x) = ln x - 1/(2x) -1/(12x2)
+ 1/(120x4) - 1/(252x6) + 1/(240x8)
is used for x > 8
together with the property: Psi(x+1) = Psi(x) + 1/x
Data Registers: /
Flags: /
Subroutines: /
-Synthetic registers M & N may be replaced by R00 & R01
-In this case, delete line 41 and replace line 02 with
0 STO 00 RDN
01 LBL "PSI"
02 CLA
03 8
04 X<>Y
05 LBL 01
06 1/X
07 ST+ M
08 X<> L
09 1
10 +
11 X<Y?
12 GTO 01
13 1/X
14 STO N
15 X^2
16 ENTER^
17 STO Z
18 20
19 /
20 21
21 1/X
22 -
23 *
24 .1
25 +
26 *
27 1
28 -
29 *
30 12
31 /
32 RCL N
33 LN
34 LASTX
35 2
36 /
37 +
38 -
39 RCL M
40 -
41 CLA
42 END
( 61 bytes / SIZE 000 )
| STACK | INPUT | OUTPUT |
| X | x | Psi(x) |
Example: Compute Digamma(-1.6) & Digamma(1)
-1.6 XEQ "PSI" >>> Digam(-1.6)
= Psi(-1.6) = -0.269717876 ( in 4 seconds )
1
R/S >>> Psi(1) =
-0.577215665 = the opposite of the Euler's constant ( 0.5772156649
)
b) Complex
Arguments
Data Registers: R00 thru R05: temp
Flags: /
Subroutines: /
01 LBL "PSIZ"
02 STO 01
03 X<>Y
04 STO 02
05 CLX
06 STO 03
07 STO 04
08 LBL 01
09 10
10 RCL 01
11 X>Y?
12 GTO 02
13 RCL 02
14 X<>Y
15 R-P
16 1/X
17 X<>Y
18 CHS
19 X<>Y
20 P-R
21 ST+ 03
22 X<>Y
23 ST+ 04
24 1
25 ST+ 01
26 GTO 01
27 LBL 02
28 RCL 02
29 RCL 01
30 R-P
31 X^2
32 STO 00
33 1/X
34 X<>Y
35 ST+ X
36 STO 05
37 CHS
38 X<>Y
39 21
40 /
41 P-R
42 .1
43 -
44 R-P
45 RCL 00
46 /
47 X<>Y
48 RCL 05
49 -
50 X<>Y
51 P-R
52 1
53 +
54 R-P
55 12
56 /
57 RCL 02
58 RCL 01
59 R-P
60 RDN
61 ST- Z
62 X<> T
63 /
64 P-R
65 .5
66 +
67 R-P
68 RCL 02
69 RCL 01
70 R-P
71 RDN
72 ST- Z
73 X<> T
74 /
75 P-R
76 RCL 02
77 RCL 01
78 RAD
79 R-P
80 LN
81 DEG
82 R^
83 ST- Z
84 X<> T
85 -
86 RCL 04
87 ST- Z
88 X<> 03
89 -
90 END
( 118 bytes / SIZE 006 )
| STACK | INPUTS | OUTPUTS |
| Y | y | b |
| X | x | a |
with Psi(x+i.y) = a + i.b
Example: 0.7 ENTER^
1.6 XEQ "PSIZ" >>>> 0.276737830
X<>Y 0.546421305 ( in 23
seconds )
Whence Psi(1.6+0.7.i) = 0.276737830
+ i. 0.546421305
3°) Polygamma Functions
-The following program calculates the nth derivative of the Psi function
( n = 0 ; 1 ; 2 ; ..... )
- Psi'(x) = the Trigamma function , Psi''(x) = the Tetragamma
function ... etc ...
-The asymptotic expansion of the Psi-function is derived n times and
the recurrence relation Psi(n)(x+1) = Psi(n)(x)
+ (-1)nn! x-n-1 is used.
Data Registers: /
Flags: /
Subroutines: /
-Synthetic registers M N O may be replaced by R00 R01 R02
-In this case, delete line 95 and replace line 02 with
0 STO 00 RDN
01 LBL "PSIN"
02 CLA
03 X<>Y
04 STO O
05 8
06 +
07 X<>Y
08 STO N
09 LBL 01
10 RCL O
11 1
12 ST+ N
13 +
14 CHS
15 Y^X
16 ST+ M
17 CLX
18 RCL N
19 X<Y?
20 GTO 01
21 1/X
22 X^2
23 STO Y
24 RCL O
25 7
26 +
27 FACT
28 7
29 FACT
30 /
31 *
32 20
33 /
34 RCL O
35 5
36 +
37 FACT
38 2520
39 /
40 -
41 *
42 RCL O
43 3
44 +
45 FACT
46 60
47 /
48 +
49 *
50 RCL O
51 1
52 +
53 FACT
54 -
55 *
56 12
57 /
58 RCL N
59 RCL O
60 Y^X
61 ST/ Y
62 RCL N
63 *
64 RCL O
65 FACT
66 ST* M
67 X<>Y
68 ST+ X
69 /
70 -
71 RCL O
72 X=0?
73 GTO 02
74 1
75 -
76 FACT
77 RCL N
78 RCL O
79 Y^X
80 /
81 -
82 GTO 03
83 LBL 02
84 X<> N
85 LN
86 +
87 LBL 03
88 RCL M
89 -
90 1
91 CHS
92 RCL O
93 Y^X
94 *
95 CLA
96 END
( 138 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| Y | n | / |
| X | x | Psi(n)(x) |
( if n = 0 we have the Digamma ( or Psi- ) function
if n = 1 ------------ Trigamma function
if n = 2 ------------ Tetragamma function and so on ... )
Examples: Calculate Digamma(-1.6) Trigamma(-1.6) Tetragamma(-1.6) Pentagamma(-1.6)
0 ENTER^
-1.6 XEQ "PSIN" >>> Digam(-1.6) = Psi(-1.6)
= -0.269717877
1 ENTER^
-1.6 R/S
>>> Trigam(-1.6) = 10.44375936
2 ENTER^
-1.6 R/S
>>> Tetragam(-1.6) = -22.49158811
3 ENTER^
-1.6 R/S
>>> Pentagam(-1.6) = 283.4070827
... etc ... ( execution time is about 13 seconds for these
examples )
4°) Incomplete Gamma Function
a) Via Kummer's
Function, program#1
-The incomplete gamma function is defind by igam(a,x)
= §0x e -t ta-1
dt ( § denotes integrals )
but this program uses the following relation:
Formula: igam(a,x) =
(xa/a).exp(-x).M(1,a+1,x)
( M = Kummer's Function )
Data Registers:
R00 = x , R01 = 1 , R02 = 1+a
, R03 = a
Flags: /
Subroutine: "KUM" ( cf "Kummer's Function
for the HP-41" )
01 LBL "IGAM"
02 X<>Y
03 STO 03
04 1
05 STO 01
06 +
07 STO 02
08 X<>Y
09 XEQ "KUM"
10 RCL 00
11 E^X
12 /
13 RCL 00
14 RCL 03
15 ST/ Z
16 Y^X
17 *
18 END
( 32 bytes / SIZE 004 )
| STACK | INPUTS | OUTPUTS |
| Y | a | / |
| X | x | igam(a,x) |
Examples:
3 ENTER^
4 XEQ "IGAM" >>>>
igam(3,4) = 1.523793389 ( in 13 seconds )
1.2 ENTER^
1.7 R/S
>>>> igam( 1.2 , 1.7 ) = 0.697290898
( in 10s )
-Another "incomplete gamma function" is also defined by P(a,x)
= igam(a,x)/GAMMA(a)
-When x tends to infiniy, igam(a,x) tends to GAMMA(a)
b) Via Kummer's
Function, program#2
-"IGAM" may produce inaccurate results if x is a large negative number.
-The following routine is preferable in this case.
Formula: igam(a,x) =
(xa/a).M(a,a+1,-x)
( M = Kummer's Function )
Data Registers:
R00 = -x , R01 = a , R02 = 1+a
Flags: /
Subroutine: "KUM" ( cf "Kummer's Function
for the HP-41" )
01 LBL "IGAM2"
02 X<>Y
03 STO 01
04 1
05 +
06 STO 02
07 X<>Y
08 CHS
09 XEQ "KUM"
10 RCL 00
11 CHS
12 RCL 01
13 ST/ Z
14 Y^X
15 *
16 END
( 31 bytes / SIZE 003 )
| STACK | INPUTS | OUTPUTS |
| Y | a | / |
| X | x | igam(a,x) |
Example:
3
ENTER^
-20 XEQ "IGAM2" >>>>
igam(3,-20) = -1.756298007 10 -11 ( in 31
seconds )
-All the digits are correct.
-With the same parameters, "IGAM" gives -1.756454406 10
-11
c) A continued
Fraction
-We may also use the continued fraction hereunder to compute GAM(a,x) = GAMMA(a) - igam(a,x)
1 1.(a-1) 2.(a-2)
GAM(a,x) = e -x xa [ --------
-------- --------- .................... ]
x+1-a+ x+3-a+ x+5-a+
Data Registers:
R00 thru R08 are used by "CFR" R01 = x , R09 = a
Flags: /
Subroutine: "CFR" ( cf "Continued
Fractions for the HP-41" )
01 LBL "IGF"
02 STO 01
03 X<>Y
04 STO 09
05 CLX
06 X<>Y
07 "T"
( the same character as line 19 )
08 ASTO 00
09 CLA
10 XEQ "CFR"
11 RCL 01
12 E^X
13 /
14 RCL 01
15 RCL 09
16 Y^X
17 *
18 RTN
19 LBL "T"
20 RCL 01
21 RCL 02
22 ST+ X
23 DSE X
24 +
25 RCL 09
26 ST- Y
27 RCL 02
28 1
29 -
30 -
31 X=0?
32 RTN
33 LASTX
34 *
35 X=0?
36 SIGN
37 END
( 58 bytes / SIZE 010 )
| STACK | INPUTS | OUTPUTS |
| Y | a | / |
| X | x | GAM(a,x) |
Example:
PI 7 XEQ "IGF" >>>> GAM(PI,7) = 0.07985329094 ( in 15 seconds )
-This method is especially useful for "large" x-values.
5°) Beta Function
a) Program#1
-The Beta function is defined by the integral: B(x,y) =
§01 tx-1 (1-t)y-1
dt
-But we have also B(x,y) = B(y,x) = GAM(x).GAM(y)/GAM(x+y)
Data Registers: /
Flags: /
Subroutine: "GAM+" or "GAM"
( Synthetic registers N , O may be replaced by any standard registers )
01 LBL "BETA"
02 STO N
03 X<>Y
04 STO O
05 XEQ "GAM+"
06 X<> N
07 ST+ O
08 XEQ "GAM+"
09 ST* N
10 RCL O
11 XEQ "GAM+"
12 ST/ N
13 X<> N
14 CLA
15 END
( 47 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| Y | y | / |
| X | x | B(x,y) |
Example: 1 E^X PI
XEQ"BETA" >>>> B(e,Pi) = 0.03789029877 ( in 11
seconds )
-An OUT OF RANGE will occur if x+y is greater than about
70.9
-So it's often better to use "LOGAM" provided all the gammas are positive.
b) Positive
Integer Arguments
-If the arguments are both positive integers m , n , we can avoid any
subroutine
because in this case, B(m,n) = (m-1)!(n-1)!/(m+n-1)!
01 LBL "BMN"
02 X>Y?
lines 02-03 is useful to reduce execution time
03 X<>Y
04 ST+ Y
05 1
06 GTO 02
07 LBL 01
08 RCL Y
09 *
10 R^
11 /
12 LBL 02
13 DSE Z
14 DSE Y
15 GTO 01
16 RCL Z
17 /
18 END
( 33 bytes / SIZE 000 )
| STACK | INPUTS | OUTPUTS |
| Y | m | / |
| X | n | B(m,n) |
Example: 100 ENTER^
200 XEQ "BMN" >>>> B(100,200) = 3.607285453
10 -84 ( in 32 seconds )
-This elementary method also avoids many overflows
-This program does not check the arguments, and negative m , n
may produce an infinite loop: add TEXT0 or
STO X or another NOP after line 13
6°) Incomplete Beta Function
a) Via Hypergeometric
Functions
-The incomplete Beta function is defined by Bx(a,b)
= §0x ta-1 (1-t)b-1
dt ( a , b > 0 )
-We have also Bx(a,b) = (xa/a) (1-x)b
F(1,a+b,a+1,x) where F is the hypergeometric function
and this formula is used by the program hereafter.
Data Registers: R00 = x , R01 = 1 ,
R02 = a+b , R03 = a+1 , R04 = a
Flags: /
Subroutine: "HGF" ( cf "Hypergeometric
Functions for the HP-41" )
01 LBL "IBETA"
02 STO 00
03 RDN
04 STO 02
05 X<>Y
06 ST+ 02
07 STO 03
08 R^
09 X<>Y
10 Y^X
11 LASTX
12 /
13 1
14 STO 01
15 ST+ 03
16 RCL 00
17 -
18 R^
19 Y^X
20 *
21 STO 04
22 RCL 00
23 XEQ "HGF"
24 RCL 04
25 *
26 END
( 42 bytes / SIZE 005 )
| STACK | INPUTS | OUTPUTS |
| Z | a | / |
| Y | b | / |
| X | x | Bx(a,b) |
Examples:
PI
1 E^X
0.7 XEQ "IBETA" >>>> B0.7(Pi,e)
= 0.02962304602 ( in 57 seconds )
21 ENTER^
40 ENTER^
0.4 R/S
>>>> B0.4(21,40) = 4.898975630
10 -18 ( in 46s )
-There is a simpler formula to compute Bx(a,b) = (xa/a)
F(a,1-b,a+1,x) but it would yield meaningless results if b
is large ( example2 above )
-Otherwise, the convergence is often faster and the program is shorter.
-If x is very close to 1 , you can also use the relation
Bx(a,b) = B(a,b) - B1-x(b,a)
-Another "Incomplete Beta Function" is defined by Ix(a,b)
= Bx(a,b) / B(a,b)
( use "BETA" to compute Ix(a,b) )
b) A continued Fraction
1 d1 d2
-We have another relation to calculate the incomplete beta function:
Bx(a,b) = (xa/a) (1-x)b [
---- ---- ---- ........ ]
1+ 1+ 1+
where d2m = m(b-m).x/[(a+2m)(a+2m-1)]
and d2m+1 = -(a+m)(a+b+m).x/[(a+2m)(a+2m+1)]
Data Registers:
R00 thru R08 are used by "CFR" R01 = x , R09 = a , R10 = b
Flags: /
Subroutine: "CFR" ( cf "Continued
Fractions for the HP-41" )
01 LBL "IBETA2"
02 RDN
03 STO 10
04 RDN
05 STO 09
06 1
07 R^
08 "T"
the same character as line 25
09 ASTO 00
10 CLA
11 XEQ "CFR"
12 1/X
13 RCL 01
14 RCL 09
15 ST/ Z
16 Y^X
17 *
18 1
19 RCL 01
20 -
21 RCL 10
22 Y^X
23 *
24 RTN
25 LBL "T"
26 RCL 02
27 2
28 /
29 ENTER^
30 INT
31 X=Y?
32 GTO 01
33 RCL 09
34 +
35 STO Y
36 RCL 10
37 +
38 *
39 CHS
40 GTO 02
41 LBL 01
42 RCL 10
43 RCL Y
44 -
45 *
46 LBL 02
47 RCL 01
48 *
49 RCL 02
50 RCL 09
51 +
52 RCL X
53 1
54 -
55 *
56 /
57 1
58 X<>Y
59 END
( 86 bytes / SIZE 011 )
| STACK | INPUTS | OUTPUTS |
| Z | a | / |
| Y | b | / |
| X | x | Bx(a,b) |
Example:
1 E^X PI
0.4 XEQ "IBETA2" >>>> B0.4(e,Pi)
= 0.01476755416 ( in 20 seconds )
-Like "IBETA" , if x is close to 1 , you can use the relation
Bx(a,b) = B(a,b) - B1-x(b,a)
References:
Abramowitz and Stegun , "Handbook of Mathematical Functions" - Dover
Publications - ISBN 0-486-61272-4
Press , Vetterling , Teukolsky , Flannery , "Numerical Recipes in C++"
- Cambridge University Press - ISBN 0-521-75033-4
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall