The Museum of HP Calculators

HP Articles Forum

Fast DiGamma Function for Positive and Negative Real Arguments

Posted by Steven Thomas Smith on 26 Sept 2009, 6:06 a.m.

This code computes the logarithmic derivative of the Gamma function, digamma(x) = Gamma'(x)/Gamma(x). It works for all real arguments, positive and negative, and is robust for very small arguments. It simply uses the reflection formula

digamma(x+1) = digamma(x) + 1/x

to move to a region where the asymptotic expansion for digamma yields full precision. When x is negative, it uses the analytic continuation term. It is based on Jean-Marc Baillard's PSI, which uses the asymptotic expansion for positive arguments.

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