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
would be natural to jump to the conclusion of equation (1)...until 
realizing that the summation diverges.