Thanks Egan and Walter for your interest!
By analyzing the difference between the actual harmonic number and the well known asymptotic approximation, Hn = ln(n) + EulerGamma, I obtained the following approximation:
1 1 1 1
SUM(k=1,n,1/k) ~ ln(n) + EulerGamma + ----- - -------- + --------- - --------
2 4 6
2 n 12 n 120 n 252 n
Nothing new here though, that's just equation (13) in this mathworld article:
http://mathworld.wolfram.com/HarmonicNumber.html
In order to obtain the inverse function one could think of just isolating n in the expression above. The problem is this cannot be done, at least with usual algebra. The best I was able to do was isolating n in h ~ EulerGamma + ln(n) + 1/(2n) and getting
1
n ~ - --------------------------------
/ 1 \
2 W | - ------------------- |
| h - EulerGamma |
\ 2 e /
where W() is Lambert's W function and EulerGamma is the Euler-Mascheroni constant, 0.577215664901532860606512090082...
Then I added x as a correction term to the exponent of e, isolated it, which was not hard to do, and built up a table of the differences in relation to the actual values and tried some curve fitting. But I did not obtain anything worthwile.
Later I found a better approximation formula in http://groups.google.com/group/sci.math/browse_frm/thread/f2a9c879857e2ac9:
n ~ 1/2*(1/Sqrt(-3*W(-1/(12*e^(2*(h - EulerGamma))))) - 1)
By using the same technique I had tried earlier, I eventually obtained
n ~ 1/2*(1/Sqrt(-3*W(-1/(12*e^(2*(h - EulerGamma + 1/(150/11*e^(4*h))))))) - 1)
The actual constants I obtained were 3.99989985 and 13.63604051 which are very close to 4 and 150/11. This makes me believe the exact constants might have a true mathematical meaning, but I may be wrong.
The approximation works nicely as it can be seen in the table below:
h n (actual n)
------------------------------------------------
1.00 1.00022409778 (1)
1.50 2.00002207370 (2)
2.45 6.00000020420 (6)
2.93 10.0108471094 (10.0108470927)
3.50 18.0907442948 (18.0907442943)
4.00 30.1532900556 (30.1532900558)
For h greater than 11, the asymptotic approximation n ~ e^(h - EulerGamma) - 1/2 is used in the program.
The correction term I added will prevent the formula to be evaluated for h greater than 289 or so on the hp 33s due to overflow, but that is not necessary as we have seen. If such were the case, the following interesting property I have observed could be used:
n(h+k)
lim -------- = e^k
h -> inf n(h)
Regards,
Gerson.
------------------------------------------------------
Harmonic Numbers Program (hp 33s)
A0001 LBL A A0016 252
A0002 LN A0017 /
A0003 LASTx A0018 120
A0004 2 A0019 1/x
A0005 * A0020 -
A0006 1/x A0021 *
A0007 + A0022 12
A0008 STO A A0023 1/x
A0009 LASTx A0024 +
A0010 x^2 A0025 *
A0011 4 A0026 STO- A
A0012 * A0027 RCL A
A0013 ENTER A0028 5.77215664902E-1
A0014 ENTER A0029 +
A0015 ENTER A0030 RTN
LENGTH: 162
CHKSUM: 9920
Edited: 25 July 2008, 12:30 p.m.