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

[Gerson W. Barbosa]

.
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:



InvPsi will accept arguments greater or equal -1. About 0.05 s, 0.10 sor 0.50 s, depending on the arguments.