## Bernoulli Numbers and Harmonic Series

```Let H(n) denote the sum of the first n terms of the harmonic series,
i.e.,
1     1     1           1
H(n) = 1 + --- + --- + --- + ... + ---
2     3     4           n

It's well known that

H(n) -> log(n) + gamma       as         n -> inf

where gamma = 0.57721... is Euler's constant and log is the natural
logarithm.  Interestingly, it appears that the discrepancy between
H(n) and log(n)+gamma can be computed very directly from the Bernoulli
numbers.  This relation not only allows us to determine H(n) very
precisely, it can also be used to compute the value of gamma from a
known value of H(n).

For any integer n>1 it appears that

inf   B_k
H(n) = log(n) + gamma  -  SUM   -----               (1)
k=1  k n^k

where B_k denotes the kth Bernoulli number, i.e.,

B_0 = 1     B_1 = -1/2     B_2 = 1/6     B_3 = 0    B_4 = -1/30

and so on.  The summation _appears_ to rapidly converge, so that just
ten non-zero terms are needed to give

H(256) = 6.12434496281728...

Conversely, given (for example) the fact that H(8) = 2283/840, we can
easily compute

2283                    inf   B_k
gamma  =   ----  -  3 log(2)  +   SUM   -----             (2)
840                    k=1  k n^k

Again taking just the first ten non-zero terms of the summation, this
formula gives
gamma = 0.577215664901533....

which is correct to 15 decimal places.  Of course, the value of log(2)
is easily computed using the formula

/ 1 + 1/3 \         /  1      1       1         \
log(2) = log( --------- )  =  2 (  --- + ----- + ----- + ...  )
\ 1 - 1/3 /         \ 1*3   3*3^3   5*3^5       /

Anyway, more digits of gamma can be determined simply by using a larger
value of n (e.g., n=16) with the same number of terms in the Bernoulli
summation.  However...

When I said the summation _appears_ to converge very rapidly, I was
referring to the fact that the terms initially become extremely small
as k increases.  This is because the denominator k n^k increases rapidly
whereas the numerator B_k actually gets smaller...for a while.  But
in the long run the values of the Bernoulli numbers increase at a
super-exponential rate, as can be seen from the inequality

2(2k)!
| B_(2k) |   >   -----------                   (3)
(2 PI)^(2k)

Therefore, in the long run the numerators of the terms in the summation
rise faster than the denominators, so the summations in equations (1)
and (2) actually diverge.  This strikes me as a very interesting example
of asymptotic divergent series.

The reader may have noticed that (1) is given directly by the Euler-
Maclaurin summation formula.  I didn't initially recognize this as
an EM summation because of the round-about path by which I arrived at
it.  For any positive integers b and N you can construct a Nth order
linear recurring sequence approximating the values of H(b^k),
k=0,1,2,...  For example, setting  h(k) = H(2^k) we have the
(approximate) 3rd order recurrence

2h(k) - 5h(k-1) +  4h(k-2) -  h(k-3)   =  0

Taking the exact initial values h(0)=1, h(1)=3/2, and h(2)=25/12, the
closed-form solution of this recurrence gives the (approximate) values

2        2        1
h(k)  =  --- k +  ---  +  -----
3        3      3*2^k

Of course, this is exact only for the three initial values (k=0,1,2),
so let's call this function h3(k).  Similarly the fourth order
"harmonic recurrence" gives the function

218       110       17           64
h4(k)  =   --- k  +  ---  +  ------  -  ----------
315       189     35*2^k     945*2^(2k)

and so on.  Obviously the first j values of hj(k) are the exact values
h(k) of the harmonic series.  Going to higher and higher order formulas
it quickly becomes apparent that the coefficient of k is approaching
ln(2), and the constant coefficient is approaching gamma.  For example,
the constant coefficient of h7(k) is

330229884852951965355239887485346
---------------------------------  =  0.5772156...
572108310987828197293898319807675

Also, the coefficients of the subsequent terms approach B_n/n, so it