## 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...