# Finding how many terms of the harmonic series must be summed to exceed x?

The harmonic series is the sum

1 + 1/2 + 1/3 + 1/4 + 1/5 + 1/6 + ... + 1/n + ...

It is known that this sum diverges, meaning (informally) that the sum is infinite and (more formally) that for any real number x, there there is some number n such t that the sum of the first n terms of the harmonic series is greater than x. For example, given x = 3, we have that

1 + 1/2 + 1/3 + ... + 1/11 = 83711/27720 ≈ 3.02

So eleven terms must be summed together to exceed 3.

Consider the following question

Given an integer x, find the smallest value of n such that the sum of the first n terms of the harmonic series exceeds x.

Clearly we can compute this by just adding in more and more terms of the harmonic series, but this seems like it could be painfully slow. The best bound I'm aware of on the number of terms necessary is 2O(n), which uses the fact that

1 + 1/2 + 1/3 + 1/4 + 1/5 + 1/6 + 1/7 + 1/8 + ...

is greater than

1 + (1/2) + (1/4 + 1/4) + (1/8 + 1/8 + 1/8 + 1/8) + ...

which is in turn

1 + 1/2 + 1/2 + 1/2 + ...

where each new 1/2 takes twice as many terms as the previous to accrue. This means that the brute-force solution is likely to be completely infeasible for any reasonably large choice of x.

Is there way to calculate the harmonic series which requires less operations than the brute-force solution?

-

## migrated from stackoverflow.comJan 22 '12 at 18:06

This question came from our site for professional and enthusiast programmers.

How exact does the answer need to be? Is it acceptable to limit x to an integer? What the maximum x you need to be able to handle? – David Schwartz Jan 21 '12 at 22:06
@DavidSchwartz- I'm mostly interested in the case where x is an integer, and there is no upper bound on x - let's assume that we're using BigIntegers on a machine with unbounded memory. Also, I would like an absolute exact answer if at all possible. – templatetypedef Jan 21 '12 at 22:10
See oeis.org/A002387 for some information. – DSM Jan 21 '12 at 22:17
If an approximate answer will do e^(n - .57721 - 1/(2n)) is pretty close. – David Schwartz Jan 21 '12 at 22:19

The DiGamma function (the derivative of the logarithm of the Gamma function) is directly related to the harmonic numbers: ψ(n) = Hn-1 - γ, where γ is Euler's constant (0.577...).

You can use one of the approximations in the Wikipedia article to compute an approximate value for Hn, and then use that in a standard root finding algorithm like bisection to locate the solution.

-

If you replace the sum by an integral you get the logarithm function so that $\ln(n)$ is a first approximation.

In fact the Euler $\gamma$ constant ($.577215664901532860606512090$) may be defined by the following formula : $\displaystyle \gamma=\lim_{n \to \infty} \left(H_n-\ln(n+1/2)\right)$

From this you may deduce the equivalence as $n \to \infty$ : $$H_n \thicksim \gamma + \ln(n+1/2)$$
(for $n=10^6$ we get about 14 digits of precision)

And revert this (David Schwartz proposed a similar idea) to get the value $n$ required to get a sum $s$ : $$n(s) \approx e^{s-\gamma} -\frac12$$

The first integer to cross the $s$ should be given by $\lfloor e^{s-\gamma}+\frac12\rfloor\;$ ('should be' because of the little error made on $H_n$ compensated by the low probability of people testing values much higher than 20 :-)).

Example : the sum will cross the value $20$ for $n$ evaluated at $\rm floor(\rm exp(20-gamma)+0.5)= \rm round(\rm exp(20-gamma))= 272400600$ and indeed (this is not a proof!) :
$H_{272400599}=19.9999999977123$
$H_{272400600}=20.0000000013833$

-

-
Interesting thanks! It seems that Benoit Cloitre proposed in 2002 a more precise conjecture "for $n\gt 1$, $a(n) = \lfloor e^{n-\gamma}+\frac12\rfloor$" (or a(n)= round(exp(n-gamma)) ). – Raymond Manzoni Jan 22 '12 at 21:00

I think you meant O(2^(2n)) = O(4^n) instead of O(2^n). The "right" bound is O(exp(n)), which holds because H(m) >= ln (m + 1).

    i+1
1    /
- >= | (1/x) dx = ln (i+1) - ln i
i    /
i

m        m
---      ---
\   1    \
>  - >=  >  (ln (i+1) - ln i) = ln (m+1) - ln 1 = ln (m+1).
/   i    /
---      ---
i=1      i=1


The main obstacle to a provably fast exact algorithm is the usual table maker's dilemma: we don't have a good handle on exactly how close H(m) can be to an integer, so it's not clear a priori how much precision is needed for the approximate methods. There's a similar issue with many optimization problems involving Euclidean distances not being in NP: we don't know how to test the sign of a sum of square roots efficiently.

-

The $n^{th}$ harmonic number, $H_n$ has an asymptotic expansion of the form:

$\hspace{2cm} \displaystyle H_n \sim \ln{n}+\gamma+\frac{1}{2n}-\sum_{k=1}^\infty \frac{B_{2k}}{2k n^{2k}}=\ln{n}+\gamma+\frac{1}{2n}-\frac{1}{12n^2}+\frac{1}{120n^4}-...$

-