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
(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.5772156649015328
293596189489309783
Thanks Albert for taking some time looking at my post, your program is nice and short.
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.0946503206224769
91469435654672848+1.076674047468581
096181134287806413i
Argument with negative real part, we can use
reflection formula
\( ψ(z) = ψ(1-z) - \pi \cot(\pi\,z) \)
(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.
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.