Inverse Harmonic and Inverse Digamma Functions (HP 49G (*) and HP 50g )

[Albert Chan]

(04-27-2016 02:53 AM)Gerson W. Barbosa Wrote:  4 P2 --> 55.0973869316 Psi --> 4.00000000039
3 P2 --> 20.5834634677 Psi --> 3.00000002131

Proceeding likewise, we get

P3:

« DUP P2 Psi OVER - - P2
»

This is good for x as low as 2, but not for x = 1:

Thanks for the link to get me here.

Regading the trouble with non-convergence, this algorithm may not be stable.

Code:

from mpmath import *
mp.pretty = True    
p1 = lambda x: exp(x) + 0.5

def p2(x):
    try: y0 = p1(x); x0 = psi(0,y0)     # y0 = inv_psi(x) guess
    except TypeError: x, x0, y0 = x     # unpack
    y0 = p1(x + (x-x0))                 # algorithm in question ...
    return x, psi(0,y0), y0             # y0 is new guess

>>> p2(1)
(1, 1.00005982215968, 3.20333496618784)
>>> p2(_) # = p3
(1, 1.00545464001878, 3.21811921983331)
>>> p2(_) # = p4
(1, 1.00011835326821, 3.20349494484173)
>>> p2(_) # = p5
(1, 1.00539674070051, 3.21796012995814)

Iterated results are oscillatory. For this case, p2 result is the best one.
Anyway, result is still bad. We might as well use Newton's method
Trivia: when x is big, psi(0,x) ≈ ln(x-0.5) → psi(1,x) ≈ 1/(x-0.5)

>>> x = p1(1)
>>> x -= (psi(0,x)-1) / psi(1,x)
>>> x -= (psi(0,x)-1) / psi(1,x)
>>> x, psi(0,x)
(3.20317146837693, 1.0)

Just for fun, I use above "perfect guess" for p2(1)

>>> p2((1, 1.0, 3.20317146837693))
(1, 1.00551381652574, 3.21828182845905)
>>> p2(_) # = p3
(1, 1.00005982215968, 3.20333496618784)
>>> p2(_) # = p4
(1, 1.00545464001878, 3.21811921983331)