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) Gamma & Digamma Functions
     a-3) 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
 

Last update:   "GPS" in paragraph 2°)-a-2)
 
 

1°) Gamma Function
 

        a) Program#1
 

-The asymptotic formula  e^-x x^x-1/2 (2pi)^1/2 exp ( 1 + 1/(12x) -1/(360x³) ) 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.x² ) + 1/( 105.x⁴ ) - 1/( 140.x⁶ ) + 1/( 99.x⁸ ) - .... ] 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 ²⁰⁰  ( 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  10²⁵⁶⁴
 

-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/(12x²) + 1/(120x⁴) - 1/(252x⁶) + 1/(240x⁸)  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)  Gamma & Digamma Functions
 

-"GPS" computes gamma & digamma functions ( almost ) simultaneously.
-Asymptotic expansions are used in both cases.
 

Data Registers:    R00 & R01: temp
Flags:  /
Subroutines: /
   
 

01 LBL "GPS"  02 SIGN  03 1  04 STO 00            05 CLST  06 8  07 LASTX  08 LBL 01  09 ST* 00  10 1/X  11 ST- Z  12 X<> L  13 1  14 +  15 X<Y?  16 GTO 01  17 STO 01

18 1/X  19 LASTX            20 LN  21 R^  22 +  23 RCL Y  24 2  25 /  26 -  27 X<>Y  28 X^2  29 5  30 %  31 21  32 1/X  33 -  34 RCL Y

35 *  36 .1  37 +  38 RCL Y  39 *  40 1  41 -  42 *  43 12  44 /  45 +  46 X<> 00  47 RCL 01            48 RCL 01  49 210  50 *  51 53

52 X<>Y  53 /  54 RCL Y  55 +  56 30  57 *  58 1/X  59 RCL Y            60 +  61 12  62 *  63 1/X  64 RCL Y  65 -  66 E^X  67 PI

68 ST+ X  69 SQRT  70 R^  71 /  72 *  73 X<>Y  74 R^  75 .5  76 ST- Y  77 *  78 Y^X  79 ST* Y  80 *  81 RCL 00            82 X<>Y  83 END

 
       ( 114 bytes / SIZE 002 )
 
 

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

Examples:

  PI  XEQ "GPS"   >>>>   Gam(Pi) = 2.288037795                                 ---Execution time = 6s---
                             X<>Y    Psi(Pi)  = 0.977213308

  7.28  CHS   R/S   >>>>   Gam(-7.28) = 0.0004577130753                  ---Execution time = 9s---
                              X<>Y    Psi(-7.28)  = 4.651194213
 

            a-3)  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⁽ⁿ⁾(x+1) = Psi⁽ⁿ⁾(x) + (-1)ⁿn! 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) =  §₀^x  e ^-t t^a-1 dt    ( § denotes integrals )
  but this program uses the following relation:

Formula:      igam(a,x) =  (x^a/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) =  (x^a/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 x^a
                          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) = §₀¹  t^x-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   B_x(a,b) = §₀^x  t^a-1 (1-t)^b-1 dt       ( a , b > 0 )
-We have also  B_x(a,b) = (x^a/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"  >>>>  B_0.7(Pi,e) =  0.02962304602   ( in 57 seconds )

   21   ENTER^
   40   ENTER^
  0.4      R/S        >>>>   B_0.4(21,40) =  4.898975630 10 ^-18   ( in 46s )

-There is a simpler formula to compute B_x(a,b) = (x^a/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  B_x(a,b) = B(a,b) -  B_1-x(b,a)
-Another "Incomplete Beta Function" is defined by  I_x(a,b) = B_x(a,b) / B(a,b)             ( use "BETA" to compute I_x(a,b)  )
 

        b) A continued Fraction
 

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

    B_x(a,b) =  (x^a/a) (1-x)^b  [ 1 / ( 1 + d₁ / ( 1 + d₂ / ( 1 + ... ))) ]

     where  d_2m = m(b-m).x/[(a+2m)(a+2m-1)]    and   d_2m+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"  >>>>  B_0.4(e,Pi) =  0.01476755416   ( in 20 seconds )

-Like "IBETA" , if  x is close to 1 , you can use the relation  B_x(a,b) = B(a,b) -  B_1-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