I mentioned in the previous post that the harmonic numbers Hn are never integers for n > 1. In the spirit of that post, I’d like to find the value of n such that Hn is closest to a given integer m.
We have two problems to solve. First, how do we accurately and efficiently compute harmonic numbers? For small n we can directly implement the definition. For large n, the direct approach would be slow and would accumulate floating point error. But in that case we could use the asymptotic approximation
from this post. As is often the case, the direct approach gets worse as n increases, but the asymptotic approximation gets better as n increases. Here γ is the Euler-Mascheroni constant.
The second problem to solve is how to find the value of n so that Hn comes closest to m without trying too many possible values of n? We can discard the higher order terms above and see that n is roughly exp(m − γ).
Here’s the code.
import numpy as np
gamma = 0.57721566490153286
def H(n):
if n < 1000:
return sum([1/k for k in range(1, n+1)])
else:
n = float(n)
return np.log(n) + gamma + 1/(2*n) - 1/(12*n**3)
# return n such that H_n is closest harmonic number to m
def nearest_harmonic_number(m):
if m == 1:
return 1
guess = int(np.exp(m - gamma))
if H(guess) < m:
i = guess
while H(guess) < m: guess += 1 j = guess else: j = guess while H(guess) > m:
guess -= 1
i = guess
x = np.array([abs(H(k) - m) for k in range(i, j+1)])
return i + np.argmin(x)
We can use this, for example, to find the closest harmonic number to 10.
>>> nearest_harmonic_number(10) 12366 >>> H(12366) 9.99996214846655
I wrote the code with integer values of m in mind, but the code works fine with real numbers. For example, we could find the harmonic number closest to √20.
>>> nearest_harmonic_number(20**0.5) 49 >>> H(49)**2 20.063280462918804