HP Articles Forum
[Return to the Index ]
[ Previous | Next ]
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 08 RAD 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