hp41programs

Gamma

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


Overview
 

 1°) Gamma Function

   a) Program#1
   b) Program#2 ( more accurate results )
   c) A Continued Fraction

     c-1) Focal Program
     c-2) M-Code Routine

   d) Lanczos Formula

     d-1) Focal Program
     d-2) M-Code Routine

   e)  Reciprocal of the Gamma Function
   f)  Reciprocal of the Gamma Function - more accurate results + an M-Code Routine
  g)  Logarithm of the Gamma Function
  h)  Gamma Function in the complex plane

     h-1)  Program#1
     h-2)  Program#2
     h-3)  Program#3
     h-4)  Lngamma Function

 2°) Digamma Function

   a) Real Arguments

     a-1) Focal Programs
     a-2) M-Code Routine

   b) Complex Arguments

     b-1) Program#1
     b-2) Program#2

 3°) Polygamma Functions

 4°) Incomplete Gamma Function

   a) Via Kummer's Function, program#1
   b) Via Kummer's Function, program#2
   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.
 

    Gam(x) ~  [ x^(x-1/2) ] exp [ -x + ( 1/12.x ) ( 1 - ( 1/30.x2 ) + 1/( 105.x4 ) - 1/( 140.x6 ) + 1/( 99.x8 ) - .... ] SQRT(2.Pi)
 
 

Data Registers:  /
Flags: /
Subroutines: /

-Lines 02 to 10 take advantage of the built in FACT function
 but they can be deleted after replacing line 14 by  X<>Y
 
 

01  LBL "GAM+"
02  FRC 
03  X#0? 
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
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)  A Continued Fraction
 

           c-1)  Focal Program
 

-Alternatively, we can also use the formula:

      Gam(x) = [ x^(x-1/2) ] sqrt(2.Pi) exp [ -x + ( 1/12 )/( x + ( 1/30 )/( x + ( 53/210)/( x + (195/371)/( x + ...  )))) ]

-The accuracy is the same as "GAM+"
 

Data Registers: /
Flags: /
Subroutines: /

-The last lines produce  Gam(n) = 9.999999999 E99  if  n is a non-positive integer
-Lines 53 to 57 may be replaced by   X=0?   MAXR    if you have written an M-code routine that places  9.999999999 E99  in register X
 
 

01  LBL "GAM3"
02  1
03  X<>Y
04  5
05  X<>Y
06  LBL 01 
07  ST* Z
08  1
09  +
10  X<Y?
11  GTO 01 
12  STO Y
13  371
14  *
15  195
16  X<>Y
17  /
18  RCL Y          
19  +
20  210
21  *
22  53
23  X<>Y 
24  /
25  RCL Y
26  +
27  30
28  *
29  1/X
30  RCL Y          
31  +
32  12
33  *
34  1/X
35  RCL Y 
36  -
37  E^X
38  PI
39  ST+ X 
40  SQRT
41  R^
42  X#0?
43  /
44  *
45  X<>Y           
46  R^
47  .5
48  ST- Y
49  *
50  Y^X
51  ST* Y           
52  *
53  X#0?
54  RTN
55  DEG 
56  90
57  TAN 
58  END 

 
   ( 88 bytes / SIZE 000 )
 
 

      STACK         INPUT      OUTPUT
           X              x        Gam(x)

 
Examples:

      PI   XEQ "GAM3"   >>>>      Gam(PI)    =  2.288037794                          ( in 3.5 seconds )
   -7.28  RTN    R/S       >>>>    Gam(-7.28) =  0.0004577130755                  ( in 5.6 seconds )
     -4     RTN    R/S       >>>>      Gam(-4)    =   infinity = 9.999999999 E99
 

-"GAM3" uses neither data registers nor synthetic registers. "GAM+" may be re-written in the same way!
-Note that Gam(n) is computed without any roundoff error for n = 1 , 2 , 3 , 4 , 5 , 6 , 8 , 11
-"GAM+" & "GAM3"  may be improved by using the reflection formula to reduce execution time for large negative arguments ( see below ).
 

           c-2)  M-Code Routine
 

-This routine uses the same continued fraction when x > 0
-For negative arguments, the reflection formula:  Gam(x) = -PI/[ x Gam(-x) sin ( 180° x ) ]   is also employed.
 

08D  "M"                                 The name of this routine is "GAM"
001   "A"
007   "G"
2A0   SETDEC                        if you want to check for alpha data, replace this line by the 2 words  361  ?NCXQ
0F8   READ3(X)                                                                                                                             050  14D8
128   WRIT4(L)                      x is saved in L-register ( like a built-in function )
2EE   ?C#0  ALL
09B   JNC+19d
084   CLRF5              C
0ED  ?NCXQ             =
064   193B               frc(C)
2EE  ?C#0  ALL
09F   JC+ 19d
0F8   READ3(X)
2FE   ?C#0 MS
057   JC+ 10d
00E   A=0  ALL       A
1BE  A=A-1  MS
35C  PT=12             =
162   A=A+1 @PT  -1
01D  ?NCXQ         C=
060   1807              A+C
3B1   ?NCXQ         C=
060   18EC            fact(C)        if x is a positive integer,  Gam(x) = (x-1)!  is calculated by the built-in function FACT
02B   JNC+05
04E   C=0 ALL        C
27A   C=C-1 M       =
130    LDI S&X
099   099              9.999999999 E99              if x = 0 or is a negative integer, this routine returns  9.999999999 E99
0E8   WRIT3(X)
3E0   RTN
05E   C=0 MS       C= |C|
268   WRIT9(Q)
0F8   READ3(X)
05E   C=0 MS       C= |C|
070   N=C  ALL
10E   A=C ALL
04E   C=0 ALL        C
2BE  C=-C-1 MS
35C  PT=12             =
150   LD@PT- 5     -5
01D  ?NCXQ          C=
060   1807              A+C
104   CLRF8
2FE   ?C#0 MS
053   JNC+10d
108   SETF8
04E   C=0 ALL        C
35C  PT=12             =
150   LD@PT- 5      5
10E   A=C ALL
278   READ9(Q)
01D  ?NCXQ          C=
060   1807              A+C
070   N=C  ALL
0B8  READ2(Y)
0F0  C<>N ALL
10E   A=C ALL
0A8  WRIT2(Y)
04E   C=0 ALL        C
2BE  C=-C-1 MS
266  C=C-1 S&X    =
35C  PT=12
150   LD@PT- 5   -0.5
01D  ?NCXQ          C=
060   1807              A+C
0E8   WRIT3(X)
10E   A=C ALL
0B8   READ2(Y)
0AE  A<>C  ALL
3C4   ST=0            C
045   ?NCXQ         =
06C  1B11            A^C
0F0  C<>N ALL
10E   A=C ALL
0B8  READ2(Y)
268   WRIT9(Q)
0AE  A<>C  ALL
0A8  WRIT2(Y)
04E   C=0 ALL       C
35C  PT=12
0D0  LD@PT- 3
1D0  LD@PT- 7     =
050   LD@PT- 1
130   LDI S&X
002   002               371
135   ?NCXQ       C=
060   184D          A*C
10E   A=C ALL
04E   C=0 ALL       C
35C  PT=12
050  LD@PT- 1
250  LD@PT- 9     =
150  LD@PT- 5
130   LDI S&X
002   002               195
0AE  A<>C  ALL
261   ?NCXQ        C=
060   1898            A/C
10E   A=C ALL
278   READ9(Q)
01D  ?NCXQ          C=
060   1807              A+C
10E   A=C ALL
04E   C=0 ALL       C
35C  PT=12
090  LD@PT- 2
050  LD@PT- 1     =
130   LDI S&X
002   002               210
135   ?NCXQ       C=
060   184D          A*C
10E   A=C ALL
04E   C=0 ALL       C
35C  PT=12
150  LD@PT- 5      =
0D0  LD@PT- 3
226  C=C+1 S&X  53
0AE  A<>C  ALL
261   ?NCXQ        C=
060   1898            A/C
10E   A=C ALL
278   READ9(Q)
01D  ?NCXQ          C=
060   1807              A+C
10E   A=C ALL
04E   C=0 ALL       C
35C  PT=12            =
0D0  LD@PT- 3
226  C=C+1 S&X  30
135   ?NCXQ       C=
060   184D          A*C
22D  ?NCXQ       C=
060   188D           1/C
10E   A=C ALL
278   READ9(Q)
01D  ?NCXQ          C=
060   1807              A+C
10E   A=C ALL
04E   C=0 ALL       C
35C  PT=12
050  LD@PT- 1      =
090  LD@PT- 2
226  C=C+1 S&X  12
135   ?NCXQ       C=
060   184D          A*C
22D  ?NCXQ       C=
060   188D           1/C
10E   A=C ALL
278   READ9(Q)
2BE  C=-C-1 MS   C=-C
01D  ?CXQ             C=
061   1807              A+C
044   CLRF4           C
029   ?NCXQ          =
068   1A0A          exp(C)
10E   A=C ALL
0B0  C=N ALL
135   ?NCXQ         C=
060   184D            A*C
070   N=C ALL
138   READ4(L)
05E   C=0 MS       C= |C|
0E8   WRIT3(X)
10C   ?FSET8
0A3   JNC+20d
0B0   C=N ALL
10E   A=C ALL
0F8   READ3(X)
261  ?NCXQ          C=
060   1898              A/C
070   N=C ALL
0F8   READ3(X)
00E   A=0 ALL        A
35C  PT=12             =
162   A=A+1@PT   1
01D  ?NCXQ          C=
060   1807              A+C
0E8   WRIT3(X)
10E   A=C ALL
278   READ9(Q)
37A  ?A#C M
01B  JNC +03
31A  ?A<C M
377   JC -18d
0B0  C=N ALL
10E   A=C ALL
04E   C=0 ALL        C
35C  PT=12
090   LD@PT- 2
150   LD@PT- 5
010   LD@PT- 0
190   LD@PT- 6
190   LD@PT- 6     =
090   LD@PT- 2
210   LD@PT- 8
090   LD@PT- 2
1D0  LD@PT- 7
150   LD@PT- 5  sqrt(2xPI)
135   ?NCXQ         C=
060   184D            A*C
0E8   WRIT3(X)
10E   A=C ALL
138   READ4(L)
2FE  ?C#0 MS
3A0  ?NC RTN                                      the routine stops here if x > 0 is not an integer
135   ?NCXQ         C=
060   184D            A*C
268   WRIT9(Q)
138   READ4(L)
05E   C=0 MS       C= |C|
0E8   WRIT3(X)
130   LDI S&X       C
080   080                =
0C5  ?NCXQ
028   0A31          rnd0(X)                    C is rounded to the nearest integer without changing the display settings
2A0  SETDEC
0E8   WRIT3(X)
10E   A=C ALL
04E   C=0 ALL      C
35C  PT=12           =
090   LD@PT- 2   2
044   CLRF4         C
070   N=C ALL    =
171   ?NCXQ
064   195C         AmodC
2EE   ?C#0 ALL
027   JC +04
278   READ9(Q)
2BE  C=-C-1 MS   C=-C
268   WRIT9(Q)
0F8   READ3(X)
10E   A=C ALL
278   READ9(Q)
0E8   WRIT3(X)
138   READ4(L)
01D  ?NCXQ          C=
060   1807              A+C
10E   A=C ALL
04E   C=0 ALL         C
35C   PT=12
050   LD@PT- 1       =
210   LD@PT- 8
130   LDI S&X
002   002                 180
135   ?NCXQ         C=
060   184D            A*C
070   N=C ALL
3C4  ST=0               C
144   CLRF6
284   CLRF7           =
229   ?NCXQ
048   128A           sin(C°)                    without changing the trig mode
10E   A=C ALL
0F8   READ3(X)
135   ?NCXQ         C=
060   184D            A*C
10E   A=C ALL
285   ?NCXQ            C
064   19A1
1EE   C=C+C ALL
1EE   C=C+C ALL    =
3CE  RSHFC ALL
23A   C=C+1 M
046   C=0 S&X         PI
0AE  A<>C ALL
261  ?NCXQ          C=
060   1898              A/C
0E8   WRIT3(X)
3E0   RTN
 
 

      STACK        INPUTS      OUTPUTS
           T             T            T
           Z             Z            Z
           Y             Y            Y
           X             x        Gam(x)
           L             L            x

 
Examples:

      PI    XEQ "GAM"  returns    2.288037797
    -0.2  XEQ "GAM"   ------    -5.821148569
     41   XEQ "GAM"    ------    8.159152832 E47
     -4    XEQ "GAM"   ------     9.999999999 E99
 

-Execution time varies between 0.05 s and 2.2 seconds.
-The longest execution time occurs when  -1 < x < 0
-Synthetic register Q is altered.
-See also the note after the second M-Code routine ( which is slightly longer ).
 

       d)  Lanczos Formula
 

           d-1)  Focal Program
 

-If x > 0   Gam(x) =  ( x + 4.5 ) x-1/2 exp(-(x+4.5))  sqrt(2.PI) [ 1 + a/x + b/(x+1) + c/(x+2) + d/(x+3) + e/(x+4) + f/(x+5) + ....  ]

        with             a =  76.1800917295               c =  24.0140982408              e =  0.00120865097387
                           b = -86.5053203294              d = -1.23173957245              f = -0.000005395239385

-In fact, the "1" before  a/x  is  1.00000000019...

-If x < 0   the reflection formula   Gam(x) = -PI/[ x Gam(-x) sin ( 180° x ) ]   is used
 

Data Registers: /
Flag:  F09
Subroutines: /
 
 

  01  LBL "GAM++"
  02  DEG
  03  FRC
  04  X#0?
  05  GTO 02 
  06  LASTX
  07  X>0?
  08  GTO 01
  09  90
  10  TAN
  11  RTN
  12  LBL 01 
  13  1
  14  -
  15  FACT
  16  RTN
  17  LBL 02
  18  LASTX
  19  CF 09
  20  X<0?
  21  SF 09
  22  ABS
  23  185348
  24  RCL Y 
  25  5
  26  +
  27  *
  28  1/X
  29  CHS
  30  X<>Y 
  31  4
  32  +
  33  827.36871 
  34  *
  35  1/X
  36  +
  37  X<>Y
  38  3
  39  +
  40  1.231739572  
  41  X<>Y
  42  /
  43  -
  44  X<>Y
  45  2
  46  +
  47  24.01409824
  48  X<>Y
  49  /
  50  +
  51  X<>Y 
  52  1
  53  +
  54  86.50532033 
  55  X<>Y
  56  /
  57  -
  58  76.18009173  
  59  R^
  60  /
  61  +
  62  1
  63  +
  64  PI
  65  ST+ X
  66  SQRT
  67  *
  68  X<>Y
  69  4.5
  70  +
  71  E^X
  72  ST/ Y
  73  CLX
  74  .5
  75  ST- T
  76  ST* T 
  77  X<> L 
  78  R^ 
  79  Y^X
  80  ST* Y
  81  *
  82  FC?C 09        
  83  RTN
  84  *
  85  X<>Y
  86  FIX 0
  87  RND 
  88  2
  89  MOD 
  90  1 
  91  CHS 
  92  X<>Y 
  93  Y^X
  94  * 
  95  X<>Y 
  96  RND
  97  R^
  98  - 
  99  PI 
100  R-D
101  *
102  SIN
103  *
104  PI
105  X<>Y             
106  /
107  FIX 9
108  END

 
   ( 189 bytes / SIZE 000 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             /           | x |
           X             x        Gam(x)

  | x |  is saved in  Y- Z- T-registers if x is not an integer

Examples:

      PI    XEQ "GAM++"  yields     2.288037795               ( in 5 seconds )
    70.9  XEQ "GAM++"  ------     7.827382822 E99        ---------------
  -12.9  XEQ "GAM++"  ------    -2.117236215 E-09      ( in 6 seconds )
     -4    XEQ "GAM++"  ------     9.999999999 E99        likewise for all the non-positive integers

-This routine is much faster than "GAM+" for large negative arguments.
 

           d-2)  M-Code Routine
 

-The Lanczos approximation is used in this M-code routine.
-I've used simplified - and somewhat misleading - mnemonics!
-For instance,   0F8  C=X   should be
                        0F8  READ 3(X)
 
 

08D  "M"                                 The name of this routine is "GAM"
001   "A"
007   "G"
0F8   C=X                               first executable word
361   chk alpha data
050   & xeq SETDEC
128   L=C                               x is saved in L-register ( like a built-in function )
2EE   ?C#0
09B   JNC+19d
084   C
0ED  =
064   frc(C)
2EE  ?C#0
09F   JC+ 19d
0F8   C=X
2FE   ?C#0 MS
057   JC+ 10d
00E   A
1BE   =
35C
162   -1
01D  C=
060   A+C
3B1   C=
060   fact(C)                                if x is a positive integer,  Gam(x) = (x-1)!  is calculated by the built-in function FACT
02B   JNC+05
04E   C
27A   =
130
099   9.999999999 E99             if x = 0 or is a negative integer, this routine returns  9.999999999 E99
0E8   X=C
3E0   RTN
0F8   C=X
104   cf 8
2FE  ?C#0 MS
013   JNC+02
108   sf 8
05E   C= | C |
268   Q=C
10E   A=C
0B8   C=Y
070   N=C                                  the CPU-register N is used ( but the "synthetic" register N is not )
04E   C
35C   =
110
150   4.5
01D   C=
060    A+C
0A8   Y=C
278    C=Q
10E    A=C
04E    C
2BE    =
266
35C
150    -0.5
01D    C=
060    A+C
0E8   X=C
10E   A=C
0B8   C=Y
0AE   A<>C
3C4   C
045    =
06C   A^C
0F0   C<>N
10E   A=C
0B8   C=Y
0AE   A<>C
0A8   Y=C
0AE   A<>C
029   C=
068   exp(C)
10E   A=C
0B0   C=N
0AE  A<>C
261   C=
060   A/C
268   Q=C
138   C=L
05E   C=| C |
0E8   X=C
10E   A=C
04E   C
35C   =
150   5
01D   C=
060   A+C
10E   A=C
04E   C= coefficient f
2BE
35C
150
0D0
250
150
090
21C
250
250
110
0AE  A<>C
261   C=
060   A/C
070   N=C
0F8   C=X
10E   A=C
04E   C
35C   =
110   4
01D   C=
060   A+C
10E   A=C
04E   C= coefficient e
266
266
266
35C
050
090
010
210
190
150
050
0AE   A<>C
261   C=
060   A/C
10E   A=C
0B0   C=N
01D   C=
060   A+C
070   N=C
0F8   C=X
10E   A=C
04E   C
35C   =
0D0   3
01D   C=
060   A+C
10E   A=C
04E   C= coefficient d
2BE
35C
050
090
0D0
050
1D0
0D0
250
150
1D0
090
0AE   A<>C
261   C=
060   A/C
10E   A=C
0B0   C=N
01D   C=
060   A+C
070   N=C
0F8   C=X
10E   A=C
04E   C
35C   =
090    2
01D   C=
060   A+C
10E   A=C
04E   C= coefficient c
35C
090
110
010
050
110
010
250
210
090
110
226
0AE   A<>C
261   C=
060   A/C
10E   A=C
0B0   C=N
01D   C=
060   A+C
070   N=C
0F8   C=X
00E   A
35C   =
162    1
01D   C=
060   A+C
10E   A=C
04E   C= coefficient b
2BE
35C
210
190
150
010
150
0D0
090
010
0D0
0D0
226
0AE   A<>C
261   C=
060   A/C
10E   A=C
0B0   C=N
01D   C=
060   A+C
070   N=C
04E   C= coefficient a
35C
1D0
190
050
210
15C
250
050
1D0
0D0
226
10E   A=C
0F8   C=X
261   C=
060   A/C
10E   A=C
0B0   C=N
01D  C=
060   A+C
00E   A
35C   =
162    1
01D   C=
060   A+C
10E   A=C
04E   C=sqrt(2xPI)
35C
090
150
010
190
190
090
210
090
1D0
150
135   C=
060   A*C
10E   A=C
278   C=Q
135   C=
060   A*C
0E8   X=C
10C   fs? 8
3A0   ?NC RTN              the program stops here if x > 0 is not an integer
10E   A=C
138   C=L
135   C=
060   A*C
268   Q=C
138   C=L
05E   C= | C |
0E8   X=C
130   C
080   =
0C5
028    rnd0(X)                  C is rounded to the nearest integer without changing the display settings
2A0   SETDEC
0E8   X=C
10E   A=C
04E   C
35C   =
090    2
044   C
070   =
171
064   A mod C
2EE  ?C#0?
027   JC+04
278   C=Q
2BE   C= -C
268   Q=C
0F8   C=X
10E   A=C
278   C=Q
0E8   X=C
138   C=L
01D   C=
060   A+C
10E   A=C
04E   C=180
35C
050
210
130
002
135   C=
060   A*C
070   C=sin(C)  in degrees without changing the trig mode
3C4
144
284
229
048
10E   A=C
0F8   C=X
135   C=
060   A*C
10E   A=C
285   C=PI
064
1EE
1EE
3CE
23A
046
0AE   A<>C
261   C=
060   A/C
0E8   X=C
3E0   RTN
 
 

      STACK        INPUTS      OUTPUTS
           T             T            T
           Z             Z            Z
           Y             Y            Y
           X             x        Gam(x)
           L             L            x

 
Examples:

      PI    XEQ "GAM"  gives      2.288037795               ( in 1.4 second )
    70.9  XEQ "GAM"  ------     7.827382826 E99        ---------------
  -12.9  XEQ "GAM"  ------    -2.117236215 E-09      ( in 2 second )          execution time is longer for negative arguments: we have to use the reflection formula!
     41   XEQ "GAM"  ------     8.159152832 E47        ( in 0.2 second )
     -4    XEQ "GAM"  ------     9.999999999 E99        ( in 0.05 second )     likewise for all the non-positive integers

Note:

-This routine neither checks for underflow nor for overflow, for example:

    84.7  XEQ "GAM"  produces  8.761082139  E-75   apparently meaningless!

 but if you key in  E90  /   it yields   8.761082139  E35   so   Gam(84.7) = 8.761082139  E125        ( the "9" should be a "4" )

-If the argument is a positive integer, the built-in FACT is used and it works in the same way,  for instance:

    84  XEQ "GAM"  gives  3.945523970  E-76  which is in fact  3.945523970  E124

-Of course, the display is ambiguous and I let you judge whether it is worthwhile or not
-The maximum integer is 100
-The maximum fractional argument is of the order of  99.6 ,  beyond these values, the result is really meaningless!!
-In fact, ( x + 4.5 ) x-1/2 must remain smaller than  10 200  ( for the fractional arguments )

-Otherwise, this function works like a built-in function: it saves x in L-register and registers Y Z T are unchanged
  the "alpha registers"  M N O P are unchanged too. However, the "scratch" register Q is altered ( like do the built-in  Y^X  SIN  ... etc ... )
 

       e) 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)

 
Examples:

   PI  XEQ "1/G"  >>>>  1/Gam(pi)  =  0.437055720  ( in 5 seconds )
  -3       R/S         >>>>   1/Gam(-3) = 0
 

       f) Reciprocal of the Gamma Function - more accurate results + an M-Code routine
 

- "1/G+" uses the same continued fraction as "GAM3" for more accurate results
 

Data Registers: /
Flags: /
Subroutines: /
 
 

 01  LBL "1/G+"
 02  1
 03  X<>Y
 04  5
 05  X<>Y
 06  LBL 01 
 07  ST* Z
 08  1
 09  +
 10  X<Y?
 11  GTO 01
 12  STO Y
 13  371
 14  *
 15  195
 16  X<>Y
 17  /
 18  RCL Y      
 19  +
 20  210
 21  *
 22  53
 23  X<>Y
 24  /
 25  RCL Y 
 26  +
 27  30
 28  *
 29  1/X
 30  RCL Y       
 31  +
 32  12
 33  *
 34  1/X
 35  RCL Y 
 36  X<>Y
 37  -
 38  E^X
 39  PI
 40  ST+ X        
 41  SQRT
 42  R^
 43  X<>Y
 44  /
 45  *
 46  X<>Y
 47  ENTER^    
 48  CHS
 49  .5
 50  ST+ Y 
 51  *
 52  Y^X
 53  ST* Y 
 54  *
 55 END

 
   ( 84 bytes / SIZE 000 )
 
 

      STACK        INPUTS      OUTPUTS
           X              x       1 / gam(x)

 
Examples:

    PI  XEQ "1/G+"  >>>>     1/Gam(pi)    =  0.4370557174          ( in 4 seconds )
   -3        R/S          >>>>      1/Gam(-3)   =  0                               ( in 5 seconds )
 -41.7     R/S          >>>>   1/Gam(-41.7) =  1.176053851 E50   ( in 12.5 seconds )
 

-The following M-Code routine employs the same continued fraction.
-In this listing, I've used simplified mnemonics but I hope it is legible however.
 

087  "G"
02F  "/"
031  "1"
0F8  C=X
361   chk alpha data
050   & xeq SETDEC
128  L=C
2EE  ?C#0
3A0  ?NC RTN
3C8  CLRKEY
3C4  ST=0
04E  C=0
0E8  X=C
138  C=L
084  C
0ED  =
064  frc(C)
2EE  ?C#0
02F  JC+05
138  C=L
2FE  ?C<0
360  ?C RTN
04E  C=0
268  Q=C
138  C=L
070  N=C
10E  A=C
04E  C
35C  =
190  6
0E8  X=C
2BE  C=-C
01D  ?CXQ
061  C=A+C
2FE  ?C<0
05B  JNC+11d
0F8  C=X
10E  A=C
278  C=Q
35C  PT=12
2FE  ?C<0
013  JNC+02
162  A=A+1 @PT
01D  C=
060  A+C
070  N=C
0B0  C=N
10E  A=C
04E  C
35C
0D0  =
1D0
050  3.71
135  C=
060  A*C
239  C=
060  1/AB
04E  C
35C
050  =
250
150  1.95
13D  C=
060  AB*C
0B0  C=N
025  C=
060  AB+C
04E  C
35C
090  =
050
226  21
13D  C=
060  AB*C
239  C=
060  1/AB
04E  C
35C  =
150
0D0  5.3
13D  C=
060  AB*C
0B0  C=N
025  C=
060  AB+C
04E  C
35C  =
0D0
226  30
13D  C=
060  AB*C
239  C=
060  1/AB
0B0  C=N
025  C=
060  AB+C
04E  C
35C
050  =
090
226  12
13D  C=
060  AB*C
239  C=
060  1/AB
1BE  A=A-1 MS
0B0  C=N
025  C=
060  AB+C
089  AB
064  STO Q+
04E  C
35C
150  =
266  .5
10E  A=C
0B0  C=N
2BE  C=-C
01D  C=
061  A+C
0E8  X=C
0B0  C=N
084  C
115  =
06C  LN(C)
0F8  C=X
13D  C=
060  AB*C
0D1  RCL
064  Q+
031  AB+
060  CM
044  C
035  =
068  e^AB
089  AB
064  STO Q+
138  C=L
0E8  X=C
10E  A=C
0B0  C=N
36E  ?A#C
08B  JNC+17d
3CC  ?KEY
360  ?C RTN
0A9  AB
064  <>Q+
0F8  C=X
13D  C=
060  AB*C
089  AB
064  STO Q+
0F8  C=X
00E  A
35C  =
162  1
01D  C=
060  A+C
363  JNC-20d
271  C
064  =
1EE
1EE  2 x PI ( 13 digits )
0EE  B<>C
00E  A=0
305  C=
060  sqrt(AB)
0A9  AB
064  <>Q+
0D1  RCL
064  Q+
275  C=
060  AB/CM
0E8  X=C
3E0  RTN

( 175 words )
 
 

      STACK        INPUTS      OUTPUTS
           X              x       1 / gam(x)

 
Examples:

With  x =  PI      it yields  0.4370557172
With  x = -41.7  it yields  1.176053852 E50  in 3.4 seconds
With  x = -3       it yields  0

-No reflection formula is used, so execution time may be prohibitive if x is a large negative fractional number
-Press any key to stop the loop in this case.
-This routine neither checks for underflow nor for overflow
 

       g) 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
 

      h) Gamma Function in the complex plane
 

            h-1)  Program#1
 

-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
 

            h-2)  Program#2
 

 "GAMZ+" employs the same continued fraction as "GAM3"
 

Data Registers:  R00 thru R03: temp
Flags: /
Subroutines:  "E^Z"  &  "Z^Z"  ( cf "Complex Functions for the HP-41" )

-The M-Code routines  1/Z , Z*Z , Z/Z                                   ( cf "A few M-Code Routines for the HP-41" )
  may be replaced by XEQ "1/Z" ,  XEQ "Z*Z" , XEQ "Z/Z"      ( cf "Complex Functions for the HP-41" )
 
 

01  LBL "GAMZ+"
02  STO 00 
03  X<>Y 
04  STO 03
05  CLX 
06  STO 02 
07  SIGN 
08  STO 01 
09  LBL 01
10  RCL 03
11  RCL 00
12  RCL 02
13  RCL 01
14  Z*Z
15  STO 01
16  X<>Y
17  STO 02
18  CLX
19  SIGN
20  ST+ 00
21  RCL 00 
22  5
23  X>Y?
24  GTO 01 
25  RCL 03 
26  RCL 00           
27  371
28  ST* Z
29  *
30  1/Z
31  195
32  ST* Z
33  *
34  X<>Y
35  RCL 03
36  +
37  X<>Y
38  RCL 00
39  +
40  210
41  ST* Z
42  *
43  1/Z
44  53
45  ST* Z
46  *
47  X<>Y
48  RCL 03          
49  +
50  X<>Y
51  RCL 00 
52  +
53  30
54  ST* Z
55  *
56  1/Z
57  X<>Y
58  RCL 03
59  +
60  X<>Y
61  RCL 00
62  +
63  12
64  ST* Z
65  *
66  1/Z
67  X<>Y
68  RCL 03          
69  -
70  X<>Y
71  RCL 00 
72  -
73  XEQ "E^Z"
74  PI
75  ST+ X
76  SQRT
77  ST* Z
78  *
79  X<> 00
80  X<>Y
81  X<> 03 
82  RCL Y
83  .5
84  -
85  X<>Y
86  STO T
87  X<>Y
88  XEQ "Z^Z"     
89  RCL 03 
90  RCL 00 
91  Z*Z
92  RCL 02
93  RCL 01
94  Z/Z
95  END

 
    ( 144 bytes / SIZE 004 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             y            B
           X             x             A

      where  Gam(x+iy) = A + iB

Example:

    -1.6   ENTER^
     2.4   XEQ "GAMZ+"   >>>>     0.248941260                         ---Execution time = 8s---
                                        X<>Y    -0.633220122

Whence,    Gam ( 2.4 - 1.6 i ) = 0.248941260 - 0.633220122 i

Note:

-This program may be modified to become a SIZE 000 program, which is useful for "bionic" functions:
 

            h-3)  Program#3
 

-Instead of data registers R00-R01-R02-R03, this variant uses synthetic registers  M  N  O  P
-The M-code routines  Z*Z  1/Z  Z/Z  are also employed
 

Data Registers: /
Flags: /
Subroutines:    1/Z , Z*Z , Z/Z                         ( cf "A few M-Code Routines for the HP-41" )
 
 

  01  LBL "GAMZ"
  02  RAD
  03  X<>Y
  04  STO O
  05  X<>Y
  06  0
  07  STO N
  08  SIGN
  09  STO M
  10  LBL 01
  11  CLX
  12  RCL N
  13  RCL M
  14  Z*Z
  15  STO M
  16  RDN
  17  STO N
  18  CLX
  19  SIGN
  20  +
  21  5
  22  X>Y?
  23  GTO 01 
  24  RDN
  25  STO Z
  26  371
  27  ST* Z
  28  *
  29  1/Z
  30  195
  31  ST* Z
  32  *
  33  X<>Y
  34  RCL O         
  35  +
  36  X<>Y
  37  R^
  38  +
  39  210
  40  ST* Z
  41  *
  42  1/Z
  43  53
  44  ST* Z
  45  *
  46  X<>Y
  47  RCL O 
  48  +
  49  X<>Y
  50  R^
  51  +
  52  30
  53  ST* Z
  54  *
  55  1/Z
  56  X<>Y
  57  RCL O         
  58  +
  59  X<>Y
  60  R^
  61  +
  62  12
  63  ST* Z
  64  *
  65  1/Z
  66  X<>Y
  67  RCL O 
  68  -
  69  X<>Y
  70  R^
  71  -
  72  E^X
  73  P-R
  74  PI
  75  ST+ X
  76  SQRT
  77  ST* Z
  78  *
  79  STO P         
  80  X<> T
  81  X<>Y
  82  X<> O
  83  RCL Y
  84  PI
  85  SIGN
  86  ST+ X
  87  1/X
  88  -
  89  RCL Y
  90  R^
  91  R-P
  92  LN
  93  Z*Z
  94  E^X
  95  P-R
  96  RCL O         
  97  RCL P
  98  Z*Z
  99  RCL N
100  RCL M
101  Z/Z
102  DEG
103  CLA
104  END

 
     ( 161 bytes / SIZE 000 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             y             B
           X             x             A

      where  Gam(x+iy) = A + iB

Example:

    -1.6   ENTER^
     2.4   XEQ "GAMZ"   >>>>     0.248941260                         ---Execution time = 8s---
                                      X<>Y    -0.633220122

Whence,    Gam ( 2.4 - 1.6 i ) = 0.248941260 - 0.633220122 i
 

            h-4)  LnGamma Function
 

"LNGZ" computes  Lngamma(z) = (z-1/2) Ln(z) - z + (1/2) Ln (2.PI)  + ( 1/12 )/( z + ( 1/30 )/( z + ( 53/210)/( z + (195/371)/( z + ...  ))))
 

Data Registers: /
Flags: /
Subroutines:    1/Z   Z*Z   Z+Z   Z-Z   X+1                   ( cf "A few M-Code Routines for the HP-41" )
 
 

 01  LBL "LNGZ"
 02  CLA
 03  RAD
 04  X<>Y
 05  STO O
 06  X<>Y
 07  LBL 01
 08  STO Z
 09  5
 10  X<Y?
 11  GTO 02
 12  RDN
 13  R-P
 14  LN
 15  ST+ M
 16  RDN
 17  ST+ N
 18  CLX
 19  RCL O
 20  X<>Y
 21  X+1
 22  GTO 01       
 23  LBL 02
 24  CLX
 25  371
 26  ST* Z
 27  *
 28  1/Z
 29  195
 30  ST* Z
 31  *
 32  X<>Y
 33  RCL O 
 34  +
 35  X<>Y
 36  R^
 37  +
 38  210
 39  ST* Z
 40  *
 41  1/Z
 42  53
 43  ST* Z
 44  *
 45  X<>Y
 46  RCL O        
 47  +
 48  X<>Y
 49  R^
 50  +
 51  30
 52  ST* Z
 53  *
 54  1/Z
 55  X<>Y
 56  RCL O
 57  +
 58  X<>Y
 59  R^ 
 60  +
 61  12
 62  ST* Z
 63  *
 64  1/Z
 65  X<>Y
 66  RCL O        
 67  -
 68  X<>Y
 69  R^
 70  -
 71  PI
 72  ST+ X
 73  SQRT
 74  LN
 75  +
 76  STO P
 77  X<> T
 78  X<>Y
 79  X<> O
 80  RCL Y
 81  PI
 82  SIGN
 83  ST+ X
 84  1/X
 85  -
 86  RCL Y        
 87  R^
 88  R-P
 89  LN
 90  Z*Z
 91  RCL O 
 92  RCL P
 93  Z+Z
 94  RCL N
 95  RCL M
 96  Z-Z
 97  DEG
 98  CLA
 99  END

 
          ( 153 bytes / SIZE 000 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             y             B
           X             x             A

      where  Lngamma (x+i.y) = A + i.B

Example:

        4  ENTER^
        3  XEQ "LNGZ"  >>>>  -1.756626790                          ---Execution time = 7s---
                                   X<>Y   4.742664442

Whence  Lngamma (3+4.i) =  -1.756626790 + 4.742664442 i

Notes:

-Remark that  Ln ( Gamma ( 3+4.i ) )  = -1.7566... - 1.5405... i  is different.
-So, LnGamma is not always the same as Ln ( Gamma ), though the real parts are always equal.
-Unlike Ln ( Gamma ) , LnGamma has a single branch cut: the negative real semi-axis.
 

2°) Digamma Function
 

        a) Real Arguments
 

             a-1)  Focal Programs
 

-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: /
 
 

 01  LBL "PSI"    
 02  SIGN
 03  CLST
 04  8
 05  LASTX
 06  LBL 01 
 07  1/X
 08  ST- Z
 09  X<> L
 10  1
 11  +
 12  X<Y?
 13  GTO 01      
 14  1/X
 15  LAST
 16  LN
 17  R^
 18  +
 19  RCL Y        
 20  2
 21  /
 22  -
 23  X<>Y 
 24  X^2
 25  5
 26  %
 27  21
 28  1/X
 29  -
 30  RCL Y        
 31  *
 32  .1
 33  +
 34  RCL Y 
 35  *
 36  1
 37  -
 38  *
 39  12
 40  /
 41  +
 42  END           

 
      ( 59 bytes / SIZE 000 )
 
 

      STACK         INPUT      OUTPUT
           X              x        Psi(x)

 
Examples:

     PI  XEQ "PSI"   >>>>   Psi(Pi) = 0.977213308       ( in 3 seconds )
      1         R/S        >>>>   Psi(1) = -0.577215665 = the opposite of the Euler's constant  ( 0.5772156649 )
   -7.28     R/S        >>>>   Psi(-7.28) = 4.651194213    ( in 6 seconds )
 

Note:

-Long execution time is to be expected for large negative numbers.
-The version hereafter employs the reflection formula  Psi(-x) = 1/x + Psi(x) + Pi/Tan ( 180° x )  for negative arguments.
 

Data Registers:  /
Flag:   F09
Subroutines: /

-Synthetic registers M , N , O may be replaced by R00 , R01 , R02
-In this case, delete line 65 and replace line 02 with    0   STO 00   RDN
 
 

01  LBL "PSI"   
02  CLA
03  CF 09
04  X<0?
05  SF 09
06  ABS
07  STO O 
08  8
09  X<>Y
10  LBL 01 
11  1/X
12  ST+ M
13  X<> L
14  1
15  +
16  X<Y?
17  GTO 01      
18  1/X
19  STO N 
20  X^2
21  STO Y 
22  5
23  %
24  21
25  1/X
26  -
27  *
28  .1
29  +
30  *
31  1
32  -
33  *
34  12
35  /
36  RCL N 
37  LN
38  LASTX      
39  2
40  /
41  +
42  -
43  RCL M
44  -
45  FC?C 09    
46  GTO 02
47  RCL O
48  1/X
49  +
50  PI
51  RCL O 
52  ENTER^
53  FIX 0
54  RND 
55  -
56  PI
57  R-D
58  *
59  DEG
60  TAN
61  /
62  +
63  FIX 9
64  LBL 02      
65  CLA
66  END

 
      ( 93 bytes / SIZE 000 )
 
 

      STACK         INPUT      OUTPUT
           X              x        Psi(x)

 
Examples:

     PI  XEQ "PSI"   >>>>   Psi(Pi) = 0.977213308       ( in 3 seconds )
      1         R/S        >>>>   Psi(1) = -0.577215665 = the opposite of the Euler's constant  ( 0.5772156649 )
   -7.28     R/S        >>>>   Psi(-7.28) = 4.651194216    ( in 3 seconds )
 -1234.5   R/S        >>>>   Psi(-1234.5) = 7.118826276    ( in 3 seconds )

-The reflection formula is quite interesting to reduce execution time for large negative arguments.
 

             a-2)  M-Code Routine
 

-This M-Code routine employs the same method as the 1st version of the focal program "PSI"
-There is no check for alpha data or for overflow.
 

089  "I"
013  "S"
010  "P"
2A0  SETDEC
3C8  CLRKEY
0F8  C=X
128   L=C
070   N=C ALL
084   C
0ED  =
064  frc(C)
2FE  ?C<0
01B  JNC+03
001  C=
060  AB+1
04E  C
35C  =
210   8
025  C=
060  AB+C
0F0  C<>N ALL
2BE  C=-C
025   C=
061  AB+C
2FE  ?C<0
01B  JNC+03
138  C=L
070  N=C ALL
0B0  C=N ALL
10E  A=C ALL
135  C=
060  A*C
239  C=
060  1/AB
268  Q=C
04E  C
35C
090  =
226  20
269  C=
060  AB/C
2BE  C=-C
0E8  X=C
04E  C
35C
090  =
050
226  21
22D  C=
060  1/C
0F8  C=X
025  C=
060  AB+C
278  C=Q
13D  C=
060  AB*C
04E  C
2BE
266  =
35C
050  -0.1
025  C=
060  AB+C
278  C=Q
13D  C=
060  AB*C
001  C=
060  AB+1
0B0  C=N ALL
269  C=
060  AB/C
04E  C
35C  =
190  6
269  C=
060  AB/C
001  C=
060  AB+1
0B0  C=N ALL
269  C=
060  AB/C
04E  C
35C
2BE  =
090  -2
269  C=
060  AB/C
089  AB
064  STO Q+
138  C=L
0E8  X=C                   LOOP
3CC  ?KEY
360  ?C RTN
10E  A=C ALL
0B0  C=N ALL
36E  ?A#C ALL
093  JNC+18d
0AE A<>C ALL
2BE  C=-C
22D  C=
061  1/C
0D1  RCL
064   Q+
031  C=
060  AB+CM
089  AB
064  STO Q+
0F8  C=X
00E  A
35C  =
162  1
01D  C=
060  A+C
34B  JNC-23d       GOTO LOOP
3C4  C
115  =
06C  Ln(C)
0D1  RCL
064  Q+
031  C=
060  AB+CM
0E8  X=C
3E0  RTN

 ( 123 words )
 
 

      STACK         INPUT      OUTPUT
           X              x         Psi(x)

 
Examples:

     PI   XEQ "PSI"   >>>>   Psi(Pi) = 0.9772133081       ---Execution time = 1.4s---
      1    XEQ "PSI"   >>>>   Psi(1) = -0.5772156649      ---Execution time = 1.2s---
  -7.28  XEQ "PSI"  >>>>   Psi(-7.28) = 4.651194215    ---Execution time = 2.4s---
  -41.7  XEQ "PSI"  >>>>   Psi(-41.7) = 1.459942947    ---Execution time = 5.2s---

Notes:

-The results are more accurate because we have employed 13-digit routines.
-However, the reflection formula is not used.
-So, the execution time may be too long for large negative arguments:
-Press any key to stop the loop in that case.
 

        b) Complex Arguments
 

            b-1)  Program#1
 

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 +  0.546421305  i
 

            b-2)  Program#2
 

-This variant does not use any data register.
-However, it employs several M-Code routines
 

Data Registers:  /
Flags: /
Subroutines:  Z^2  1/Z  Z/Z  Z+Z  Z-Z  X+1
 
 

 01  LBL "PSIZ"
 02  STO M
 03  X<>Y
 04  STO N
 05  X<>Y
 06  STO O
 07  10
 08  X<=Y?
 09  GTO 00
 10  X<>Y
 11  FRC
 12  X<0?
 13  X+1 
 14  +
 15  STO O     
 16  LBL 00
 17  RCL N
 18  RCL O
 19  Z^2
 20  21
 21  ST* Z
 22  *
 23  1/Z
 24  .1
 25  -
 26  RCL N
 27  RCL O
 28  Z^2
 29  Z/Z
 30  X+1
 31  RCL N
 32  RCL O     
 33  Z/Z
 34  6
 35  ST/ Z
 36  /
 37  X+1
 38  RCL N
 39  RCL O
 40  Z/Z
 41  2
 42  CHS
 43  ST/ Z
 44  /
 45  RCL N
 46  RCL O     
 47  RAD
 48  R-P
 49  LN
 50  Z+Z
 51  LBL 01
 52  RCL M
 53  RCL O
 54  X=Y?
 55  GTO 00
 56  CLX
 57  RCL N
 58  X<>Y
 59  1/Z
 60  Z-Z
 61  1
 62  ST+ M
 63  RDN
 64  GTO 01     
 65  LBL 00
 66  R^
 67  R^
 68  DEG
 69  CLA
 70  END

 
    ( 118 bytes / SIZE 000 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             y             b
           X             x             a

     with  Psi(x+i.y) = a + i.b

Examples:

   0.7  ENTER^
   1.6  XEQ "PSIZ"  >>>>  0.276737830    X<>Y    0.546421305                       ---Execution time = 9s---

  Whence   Psi( 1.6+0.7 i ) = 0.276737830 +  0.546421305  i

-Likewise, we get:      Psi( -12.4 + 1.3 i ) = 2.563568864 + 3.039764738 i                       ---Execution time = 16s---
 

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  we have the Trigamma function
      if n = 2  we have the 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" 
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)!

-Lines 02-03 are useful to reduce execution time
 
 

 01  LBL "BMN"
 02  X>Y? 
 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
 

-We have another relation to calculate the incomplete beta function:

    Bx(a,b) =  (xa/a) (1-x)b  [ 1 / ( 1 + d1 / ( 1 + d2 / ( 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"
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:

[1]  Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]  Press , Vetterling , Teukolsky , Flannery , "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4