The Museum of HP Calculators

# Gamma , Digamma , Polygamma & Beta Functions for the HP-41

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