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

04272016, 02:53 AM
(This post was last modified: 05082016 09:04 PM by Gerson W. Barbosa.)
Post: #1




Inverse Harmonic and Inverse Digamma Functions (HP 49G (*) and HP 50g )
.
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(EulerMascheroni 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 EulerMascheroni 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. 

04272016, 04:57 PM
Post: #2




RE: Inverse Harmonic and Inverse Digamma Functions (HP 49G (*) and HP 50g )
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 

04272016, 06:37 PM
(This post was last modified: 04272016 06:41 PM by Gerson W. Barbosa.)
Post: #3




RE: Inverse Harmonic and Inverse Digamma Functions (HP 49G (*) and HP 50g )
(04272016 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. 

08282020, 12:25 AM
Post: #4




RE: Inverse Harmonic and Inverse Digamma Functions (HP 49G (*) and HP 50g )
(04272016 02:53 AM)Gerson W. Barbosa Wrote: 4 P2 > 55.0973869316 Psi > 4.00000000039 Thanks for the link to get me here. Regading the trouble with nonconvergence, this algorithm may not be stable. Code: from mpmath import * >>> 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(x0.5) → psi(1,x) ≈ 1/(x0.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) 

« Next Oldest  Next Newest »

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