Digamma Function for the free42

[Juan14]

The next program uses the formula provided by Albert Chan here. The program shifts the argument value 35 units to improve precision and then by the recurrence formula ψ(x+1)=ψ(x)+1/x calculates the desired value. The values I obtained seem to be accurate in the first 15 digits. It would be interesting if some one can try it in the original hp 42s.

PHP Code:

00 { 73-Byte Prgm }
01 LBL “Psi”
02 STO 11
03 34.5
04 +
05 74381
06 689976
07 RCL× ST Z
08 ÷
09 10
10 37
11 ÷
12 RCL× ST Z
13 +
14 1/X
15 24
16 RCL× ST Z
17 +
18 1/X
19 +
20 LN
21 35
22 X<>Y
23 LBL 00
24 -1
25 RCL+ ST Z
26 RCL+ 11
27 1/X
28 -
29 DSE ST Y
30 GTO 00
32 END 


[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

[Juan14]

Thanks Albert for taking some time looking at my post, your program is nice and short.

[Albert Chan]

Previous Psi code seems unable to cope with complex argument.
Storing complex numbers in memory 11 cause a "Invalid Type" error. (why ?)

This version do it all in the stack.

Code:

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

Example, Ψ(1+i)

1+i XEQ "Psi"      → 0.094650320622476991469435654672848+1.076674047468581096181134287806413i

Argument with negative real part, we can use reflection formula

\( ψ(z) = ψ(1-z) - \pi \cot(\pi\,z) \)

[Thomas Okken]

(04-23-2022 05:56 PM)Albert Chan Wrote:  Storing complex numbers in memory 11 cause a "Invalid Type" error. (why ?)

The numbered registers are cells in a matrix named REGS. If REGS is a real matrix, the numbered registers can hold real numbers or strings, and if REGS is a complex matrix, the numbered registers can contain complex numbers (and only complex numbers!).

This is all HP-42S behavior. I have an option on my to-do list to allow using a list for REGS, which would mean you could store values of any type. No ETA on that feature yet, though.

[Juan14]

Thanks Albert, the new improvement loos neat, I didn’t even consider complex numbers, I used the register 11 just in case someone wanted to try it in the hp 42s, I was using a local variable (LSTO “x”) in my free42 and also FUNC 11 at the start of the program so the program behaves like a build in function..
I always thought the hp 50g should have a list that acts like the registers in the RPN calculators with the arithmetic functions, that implementation in the free42 would be awesome.