Computing harmonic numbers

The harmonic numbers are defined by

Harmonic numbers are sort of a discrete analog of logarithms since

As n goes to infinity, the difference between H_n and log n is Euler's constant I^3 = 0.57721" [1]

How would you compute H_n? For small n, simply use the definition. But if n is very large, there's a way to approximate H_n without having to do a large sum.

Since in the limit H_n - log n goes to I^3, a crude approximation would be

But we could do much better by adding a couple terms to the approximation above. [2] That is,

The error in the approximation above is between 0 and 1/120n⁴.

So if you used this to compute the 1000th harmonic number, the error would be less than one part in 120,000,000,000,000. Said another way, for n = 1000 the approximation differs from the exact value in the 15th significant digit, approximately the resolution of floating point numbers (i.e. IEEE 754 double precision).

And the formula is even more accurate for larger n. If we wanted to compute the millionth harmonic number, the error in our approximation would be somewhere around the 26th decimal place.

* * *

[1] See Julian Havil's excellent Gamma: Exploring Euler's Constant. It's popular-level book, but more sophisticated than most such books.

[2] There's a sequence of increasingly accurate approximations that keep adding reciprocals of even powers of n, based on truncating an asymptotic series. See Concrete Mathematics for details.