Closest harmonic quantity to an integer

I discussed within the earlier submit that the harmonic numbers H_n are by no means integers for n > 1. Within the spirit of that submit, I’d like to search out the worth of n such that H_n is closest to a given integer m.

We’ve two issues to unravel. First, how will we precisely and effectively compute harmonic numbers? For small n we will straight implement the definition. For big n, the direct strategy can be sluggish and would accumulate floating level error. However in that case we may use the asymptotic approximation

from this submit. As is commonly the case, the direct strategy will get worse as n will increase, however the asymptotic approximation will get higher as n will increase. Right here γ is the Euler-Mascheroni fixed.

The second drawback to unravel is methods to discover the worth of n in order that H_n comes closest to m with out making an attempt too many doable values of n? We are able to discard the upper order phrases above and see that n is roughly exp(m − γ).

Right 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 quantity to m
def nearest_harmonic_number(m):
    
    if m == 1:
        return 1

    guess = int(np.exp(m - gamma))    
    if H(guess) < m:
        i = guess
        whereas H(guess) < m: guess += 1 j = guess else: j = guess whereas 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 are able to use this, for instance, to search out the closest harmonic quantity to 10.

>>> nearest_harmonic_number(10)
12366
>>> H(12366)
9.99996214846655

I wrote the code with integer values of m in thoughts, however the code works effective with actual numbers. For instance, we may discover the harmonic quantity closest to √20.

>>> nearest_harmonic_number(20**0.5)
49
>>> H(49)**2
20.063280462918804

Associated posts