Digamma Function for the free42

[Albert Chan]

(08-28-2020 09:26 PM)Albert Chan Wrote:  \(\qquad\qquad\exp( \psi(x+1/2)) = x
+\frac{1}{24 \cdot x
+\Large\frac{1}{\frac{10}{37} \cdot x
+\frac{1}{\frac{689976}{74381} \cdot x \;
+\; \cdots}}} \)

"Last" constant may be adjusted to fit better. 74381/689976 * 3.7 ≈ .3989

>>> from mpmath import *
>>> mp.dps = 34
>>> def psi0(x): x-=mpf(.5); return log(x + 1/(24*x + 37/(10*x+3986/(1000*x))))
...
>>> for x in range(31, 40): print '%d\t%+g' % (x, psi(0,x) - psi0(x))
...
31      -5.06859e-16
32      -2.81682e-16
33      -1.30994e-16
34      -3.12469e-17
35      +3.35827e-17
36      +7.44676e-17
37      +9.89633e-17
38      +1.12294e-16
39      +1.18083e-16

Sweet spot is to get ψ(x+33), then adjust back down, ψ(x) = ψ(x+1) - 1/x

Code:

00 { 58-Byte Prgm }
01▸LBL "Psi"
02 1
03 -
04 STO 11
05 33.5
06 +
07 ENTER
08 ENTER
09 1/X
10 .3986
11 ×
12 +
13 3.7
14 X<>Y
15 ÷
16 X<>Y
17 24
18 ×
19 +
20 1/X
21 +
22 LN
23 33
24 X<>Y
25▸LBL 00
26 RCL 11
27 RCL+ ST Z
28 1/X
29 -
30 DSE ST Y
31 GTO 00
32 END

For euler_gamma constant, it matched 16 digits

1 XEQ "Psi" +/-      → 0.5772156649015328293596189489309783