The previous post includes code for solving the equation
Hn = m
i.e. finding the value of n for which the nth harmonic number is the closest to m. It works well for small values of m. It works for large m in the sense that the solution is very close to m, but it’s not necessarily the best solution.
For example, set m = 100. The code returns
n = 15092688622113830917200248731913020965388288
and indeed for that value of n,
Hn − 100 ≈ 3 × 10−15
and that’s as much as we could hope for with IEEE 754 floats.
The approximation
n = exp(m −γ)
is very good for large values of m. Using Mathematica we can find the exact value of n.
f[n_] := Log[n] + EulerGamma + 1/(2 n) - 1/(12 n^2) n = Floor[Exp[100 - EulerGamma]]; N[f[n], 50] 100.00000000000000000000000000000000000000000000900 N[f[n - 1], 50] 99.999999999999999999999999999999999999999999942747
So
n = 15092688622113788323693563264538101449859497
A similar process can find the solution to
Hn = 1000
is
n = 110611511026604935641074705584421138393028001852577373936470952377218354575172401275457597579044729873152469512963401398362087144972181770571895264066114088968182356842977823764462179821981744448731785408629116321919957856034605877855212667092287520105386027668843119590555646814038787297694678647529533718769401069269427475868793531944696435696745559289326610132208504257721469829210704462876574915362273129090049477919400226313586033
In case you’re wondering whether my function for computing harmonic numbers is accurate enough, it’s actually overkill, with error O(1/120n4).