+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: Not HP Calculators (/forum-7.html)
+--- Forum: Not quite HP Calculators - but related (/forum-8.html)
+--- Thread: Digamma Function for the free42 (/thread-18293.html)
Digamma Function for the free42 - Juan14 - 04-22-2022 02:20 PM
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
RE: Digamma Function for the free42 - Albert Chan - 04-23-2022 02:21 AM
(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 ENDFor euler_gamma constant, it matched 16 digits
1 XEQ "Psi" +/- → 0.5772156649015328293596189489309783
RE: Digamma Function for the free42 - Juan14 - 04-23-2022 05:06 PM
Thanks Albert for taking some time looking at my post, your program is nice and short.
RE: Digamma Function for the free42 - Albert Chan - 04-23-2022 05:56 PM
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 ENDExample, Ψ(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) \)
RE: Digamma Function for the free42 - Thomas Okken - 04-23-2022 06:07 PM
(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.
RE: Digamma Function for the free42 - Juan14 - 04-24-2022 04:46 PM
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.