Re: Inverse Of Harmonic Numbers
David W. Cantrell
This (presumably new) result provides a simple and accurate way to
approximate that inverse. I also correct an old error of mine.
--------------------------------------------------------
For background, readers may wish to look at the previous part of this
thread "Inverse Of Harmonic Numbers", started by Leroy Quet on 2001 Nov. 11,
<http://groups.google.com/group/sci.math/browse_frm/thread/f2a9c879857e2ac9>.
As before, let InvH(x) denote the desired inverse and g Euler's gamma
constant. I had noted that InvH(x) can be approximated well by
e^(x - g) - 1/2 when x >> 0. At that time, however, I had not known that
that approximation had already appeared in the literature; see, for example,
R. P. Boas, Jr. and J. W. Wrench, Jr. "Partial sums of the harmonic series"
_Amer. Math. Monthly_ 78 (1971) 864-870.
--------------------------------------------------------
The (presumably new) result:
With u = e^(x - g), an asymptotic series for InvH(x) has the form
-1/2 + Sum(j=0, oo, c_j u^(1-2j)). The list of coefficients c_j is
{1, -1/24, 3/640, -1525/580608, 615881/199065600, -3058641/504627200,
38800188510523/2191186722816000, -3213747182969063/44497945755648000,
100462329712125/255806104666112,
-43865443313064357090353257/15953645581139831685120000,
4543042335221166932765440567147/188420950968830433165312000000,
-103986681387361620043171941/401521614736326656000000, ...},
or approximately
{1, -0.0416667, 0.0046875, -0.00262656, 0.00309386, -0.00606119,
0.0177074, -0.0722224, 0.392728, -2.74956, 24.1111, -258.982, ...},
which can be obtained by series reversion. Mathematica code for the above
list of the first twelve coefficients is
n = 12; coeffs = InverseSeries[
Exp[Series[HarmonicNumber[x - 1/2], {x, Infinity, 2n - 1}] - EulerGamma],
u][[3]]; Table[coeffs[[2i - 1]], {i, 1, n}]
I suspect that this asymptotic series for InvH(x) is new. (But if anyone
has encountered it before, a reference would be much appreciated!) One
reason I suspect that it's new is that the coefficients c_j are not
presently found in OEIS. (I plan to submit them to OEIS later.) BTW,
speaking of OEIS, readers might be interested in
<http://www.research.att.com/~njas/sequences/A002387>, in particular, Dean
Hickerson's simple heuristic argument that, with probability > 97%, the
least integer k such that H(k) > n, a positive integer, is always given by
k = floor(e^(n - g) + 1/2). Of course, one can now use more terms of the
asymptotic series to get an expression for k which has an even higher
probability of always being correct. For example, using the same heuristic
argument and employing two more terms of the series, we can get
k = floor(u + 1/2 - 1/(24u) + 3/(640u^3)) where u = e^(n - g), with
probability > 99.995%.
Comments, references, etc. are welcomed.
-----------------------------------------------
I need to take this opportunity to correct an error I made in the old part
of the thread. On 2001 Nov. 12, I had said, in "responding" to myself
> Perhaps it's worth noting that, if one is interested in inverting just
> the _discrete_ H(n), then we can get a precise result by merely rounding
> to the nearest integer. In other words, for all positive integers n,
>
> Round( e^(H(n)-g) ) = n
>
>> But as x increases without bound, the error (not just the
>> relative error) approaches 0. As examples,
>>
>> e^(H(25)-g) - 1/2 = 25.0016 and
>>
>> e^(H(1000)-g) - 1/2 = 1000.000042
It is true that "inverting just the _discrete_ H(n), then we can get a
precise result by merely rounding to the nearest integer." But, as the
quoted numerical examples of mine make obvious, I _should_ then have said
| In other words, for all positive integers n,
|
| Round( e^(H(n)-g) - 1/2 ) = n
I must suppose that's what I had intended to say, although a correct
alternative would have been to have said that
Floor( e^(H(n)-g) ) = n for all positive integers n.
-----------------------------------------------
David W. Cantrell