+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: Inverse Harmonic and Inverse Digamma Functions (HP 49G (*) and HP 50g ) (/thread-6157.html)
Inverse Harmonic and Inverse Digamma Functions (HP 49G (*) and HP 50g ) - Gerson W. Barbosa - 04-27-2016 02:53 AM
.
InvH:
« .577215664902 - InvPsi 1. -
»
InvPsi:
« DUPDUP DUP InvP InvP DUP InvP InvP DUPDUP InvP InvP DUP InvP DUP EXP .5 + Psi OVER - - EXP .5 +
»
InvP:
« DUP EXP .5 + Psi OVER - - EXP .5 + Psi OVER - -
»
(*) On the HP 49G, replace Psi with 0 PSI.
Examples:
4.012007 InvH --> 30.523595226 (in less than half a second); this is a very belated solution to one of the famous Valentin's challenges (#3 here).
.577215664902 InvH --> 0.46163214497 ; InvH(Euler-Mascheroni constant) = xmin (local minimum of the continuous factorial function).
4 pi 2 / - 3 ENTER 2 LN * - InvH --> 0.24999999999 ; one of the special values for fractional arguments examples here.
1.5 InvH --> 2. ; H2
1 Psi --> -0.577215664902 InvPsi --> 1.00000000009
2 Psi --> 0.422784335098 InvPsi --> 2.
Background:
Let H(x) be the continuous function associated with Harmonic Numbers. Then,
\[H(x)=\gamma +\psi \left ( x+1 \right )\]
\[\psi \left ( x+1 \right )=H(x)-\gamma\]
\[x+1 =\psi^{-1} \left ( H(x)-\gamma \right )\]
\[x+1 =\psi^{-1} \left ( H(x)-\gamma \right )\]
\[x =\psi^{-1} \left ( H(x)-\gamma \right )-1\]
or
\[H^{-1}(x) =\psi^{-1} \left ( x-\gamma \right )-1\]
That is, in order to obtain the inverse of H(x) we only need the Inverse Digamma Function and the Euler-Mascheroni constant. No problem with the constant, but the Inverse Digamma might be a problem since Digamma is not easily invertible.
A rough approximation is
\[\psi^{-1} \left ( x\right )=e^{x}+\frac{1}{2}\]
The equivalent HP 50g program is
P1:
« EXP .5 +
»
However, this is good only for x >= 10, not good enough for our purposes:
10 P1 --> 22026.9657948 Psi --> 10.0000000001.
But,
9 P1 --> 8103.58392758 Psi --> 9.00000000064
8 P1 --> 2981.45798704 Psi --> 8.00000000469
7 P1 --> 1097.13315843 Psi --> 7.00000003465
So, let's try to improve the accuracy a bit:
P2:
« DUP P1 Psi OVER - - P1
»
10 P2 --> 22026.9657926 Psi --> 9.99999999999
9 P2 --> 8103.58392239 Psi --> 8.99999999999
8 P2 --> 2981.45797306 Psi --> 8.
7 P2 --> 1097.13312043 Psi --> 7.
6 P2 --> 403.928690211 Psi --> 6.
5 P2 --> 148.912878357 Psi --> 5.00000000001
But,
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:
2 P3 --> 7.88342863120 Psi --> 2.
1 P3 --> 3.20317150637 Psi --> 1.0000000139
A couple more steps suffice for x around -0.6 and greater, which is good enough for our purposes. That's what the InvPsi program above does, albeit in an inelegant way. Also, this is just an intuitive and somewhat inneficient approach. Better methods suggestions are welcome.
Edited to fix a couple of typos.
Edited again to include a printout of my current directory:
![[Image: 26830083161_611c747a6c_b.jpg]](https://c2.staticflickr.com/8/7670/26830083161_611c747a6c_b.jpg)
InvPsi will accept arguments greater or equal -1. About 0.05 s, 0.10 sor 0.50 s, depending on the arguments.
RE: Inverse Harmonic and Inverse Digamma Functions (HP 49G (*) and HP 50g ) - John Keith - 04-27-2016 04:57 PM
Cool, thanks! I was just trying to figure out an inverse harmonic number function the other day and getting nowhere. Neither Mathworld nor Wikipedia seemed to have much to say on the subject.
John
RE: Inverse Harmonic and Inverse Digamma Functions (HP 49G (*) and HP 50g ) - Gerson W. Barbosa - 04-27-2016 06:37 PM
(04-27-2016 04:57 PM)John Keith Wrote: Cool, thanks! I was just trying to figure out an inverse harmonic number function the other day and getting nowhere. Neither Mathworld nor Wikipedia seemed to have much to say on the subject.
Thanks for your interest!
If the arguments are always plain harmonic numbers, that is, the ones obtained from the discrete function, then the inverse function can be implemented simply as
InvHn:
« .577215664902 - EXP IP »
Examples:
137 ENTER 60 / InvHn --> 5. ; 137/60 = H(5)
10 Hx --> 2.92896825397 ; H(10)
InvHn --> 10.
2E10 --> 24.2962137754 ; H(2x10^10)
InvHn --> 20000000000.
where Hx is the program for the continuous function:
Hx:
« 1. + Psi .577215664902 + »
Regards,
Gerson.
RE: Inverse Harmonic and Inverse Digamma Functions (HP 49G (*) and HP 50g ) - Albert Chan - 08-28-2020 12:25 AM
(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)