Post Reply 
Digamma Function for the free42
04-22-2022, 02:20 PM
Post: #1
Digamma Function for the free42
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 
Find all posts by this user
Quote this message in a reply
04-23-2022, 02:21 AM
Post: #2
RE: Digamma Function for the free42
(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
Find all posts by this user
Quote this message in a reply
04-23-2022, 05:06 PM
Post: #3
RE: Digamma Function for the free42
Thanks Albert for taking some time looking at my post, your program is nice and short.
Find all posts by this user
Quote this message in a reply
04-23-2022, 05:56 PM
Post: #4
RE: Digamma Function for the free42
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) \)
Find all posts by this user
Quote this message in a reply
04-23-2022, 06:07 PM
Post: #5
RE: Digamma Function for the free42
(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.
Visit this user's website Find all posts by this user
Quote this message in a reply
04-24-2022, 04:46 PM
Post: #6
RE: Digamma Function for the free42
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.
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)