Recurrences For Harmonic Sums

Let h(n) denote the nth partial sum of the harmonic series, and define
H[k] = h(2^k).  It's not hard to show that H[n] approaches a linear
recurring sequence satisfying

                H[k] - 2 H[k-1] + H[k-2]  =  0

as k increases.  The error of this recurrence decreases by a factor of 
two at each step.  Correcting for this error, we arrive at the third-
order recurrence

            2 H[k] - 5 H[k-1] + 4 H[k-2] - H[k-3]  =  0

Eliminating the main error term from this relation leads to a fourth-
order recurrence, and so on.  The coefficients of these recurrences 
up through order six are tabulated below (without signs):

(2)                1      2      1
(4)                2      5      4      1
(16)               8     22     21      8      1
(64)             128    360    358    149     24     1
(256)           8192  23168  23272   9894   1685    88    1

The numbers in parentheses on the left indicate the factor by which 
the error of the respective recurrence decreases at each step.  If 
c[i,j] denotes the coefficient in the ith row and the jth column, 
and r_i denotes the "error factor" of the ith row, then the rule 
of formation for these coefficients can be expressed as

        c[i,j]  =  (r_(i-1)) c[i-1,j] + c[i-1,j-1]              

Note that the sum of the elements of each row, with alternating signs, 
is zero.  Also, the "error factor" increases by a multiple of 4 at 
each level except the first, so we have r_1 = 2 and r_i = 4^(i-1) 
for i>1.

This same approach works for other geometric arguments, e.g., if we 
define G[k] = h(3^k), we can derive a sequence of recurrences with the
coefficients listed below (without signs):

(3)              1      2      1
(9)              3      7      5      1
(81)            27     66     52     14     1

The rule of formation is the same, except that in this case we have 
r_1 = 3 and r_i = 9^(i-1) for i>1.

To generalize this a bit, for any fixed integer m>1 define 

                        H[k] = h(m^k)

Since h(n) increases logrithmically, it's natural to expect that the 
above recurrence is a good starting point for any m, because

    ln(m^k)  +  A ln(m^(k-1))  +  B ln(m^(k-2)) 

             =  k ln(m)  +  A (k-1) ln(m)  +  B (k-2) ln(m)

Setting this equal to zero and dividing by ln(m) gives

        k + A (k-1) + B (k-2)   =   k(1+A+B) - (A+2B)  =  0
 
which implies the two conditions

             1 + A + B = 0            A + 2B = 0

so we must have A=-2 and B=1.  However, the third order recurrence 
for h(m^k) gives only two constraints on three coefficients, i.e.,

           1 + A + B + C = 0           A + 2B + 3C  =  0

and similarly for higher order recurrences.  Thus, all recurrences for
ln(m^k) of order greater than 2 are underspecified.  Anyway, taking the
order 2 recurrence we see that the characteristic polynomial splits 
into linear factors

                   x^2 - 2x + 1  =  (x-1)^2

It turns out that the characteristic polynomials of all the asymptotic
recurrences for H[k] described previously also split into linear factors.
Specifically, for m=2 we have

                       2x^3 - 5x^2 + 4x - 1 = (2x-1)(x-1)^2
  
              8x^4 - 22x^3 + 21x^2 - 8x + 1 = (4x-1)(2x-1)(x-1)^2

128x^5 - 360x^4 + 358x^3 - 149x^2 + 24x - 1 = (16x-1)(4x-1)(2x-1)(x-1)^2

and so on.  Similarly for m=3 we have

                      3x^3 - 7x^2 + 5x - 1  =  (3x-1)(x-1)^2

            27x^4 - 66x^3 + 52x^2 - 14x + 1 = (9x-1)(3x-1)(x-1)^2

and so on.  In general, the characteristic polynomial of the asymptotic
Jth-order recurrence (J>3) for H[k] with any fixed m is

                                 J-3
      P(x)  =  (mx-1) (x-1)^2  PROD  (m^(2t) x - 1)
                                 t=1

so the characteristic roots (eigen values) are

          1, 1, 1/m, 1/m^2, 1/m^4, 1/m^6,..., 1/m^(2(J-3))

This allows us to give an explicit formula for the terms of the 
Jth-order linear recurring sequence of values that approximate H[k].  
Let Q be a JxJ matrix with the components

                Q[i,1] = 1              for i=1 to J
                Q[i,2] = i              for i=1 to J
                Q[i,3] = 1/m^i          for i=1 to J
                Q[i,j] = 1/m^2i(j-3)    for i=1 to J and j=4 to J

and let R be a column vector with the components

                       R[i] = H[i-1]     for i=1 to J

Then the asymptotic recurrence for H[k] has the following closed-form
solution using the exact values H[0] through H[J-1] as initial values:

   H[k-1]  =  [ 1  k  1/m^k  1/m^2k  1/m^4k ... 1/m^(2k(J-3)) ] Q^-1 R

Clearly, as k increases, all but the first two terms of this expression
go to zero, so for large values of k we have approximately

                    H[k-1] = C[1] + C[2]*k

where C[1] and C[2] are the first two components of the column vector 
C = Q^-1 R.  The difference between H[k-1] and h(2^(k-1)) is therefore
approximated by

         C[1] + C[2] k  - (k-1) ln(2)

                   =   C[1] + ln(2)   +   k ( C[2] - ln(2) )

For high order recurrence this value becomes virtually independent of 
k, which suggests that C[2] approaches ln(2).  Making this substitution, 
we would expect to find that

                   gamma ~=  C[1] + C[2]

To illustrate, with m=2 the closed-form solution for the 6th-order 
recurrence gives 
                     C[1] + C[2]  =  0.5772175....

which is approximately equal to gamma = 0.5772156...