[seqfan] Attempt at resummation of almost-geometric series A372908 A372909

[Richard J. Mathar]

There are numeric constants defined in the style of almost-geometric series,
Sum_{k>=0} (a^k+1)/(b^k+1) for integer pairs a and b, 0<a<b.
If the "+1" in the numerator and denominator were absent, this
are geometric series 
  Sum_{k>=0) (a/b)^k = 1/(1-a/b).
This is related with constants like A372907, A372908, A372909 (apart from
signs).
One idea of series acceleration is to expand the denominator
in the generalized series in a geometric series, to swap the order of summations, 
to pack two adjacent terms of the new alternating series into a single one...:

sum_{k>=0} (a^k+1)/(b^k+1)
= sum_{k>=0} (a^k+1)/b^k * 1/(1+1/b^k) 
= sum_{k>=0} (a^k+1)/b^k sum_{j>=0} (-1)^j b^(-j*k)
= sum_{j>=0} (-1)^j \sum_{k >=0} (a^k+1)/b^k b^(-j*k)
= sum_{j>=0} (-1)^j \sum_{k >=0} a^k/b^k b^(-j*k)
  + sum_{j>=0} (-1)^j \sum_{k >=0} 1/b^k b^(-j*k)
= sum_{j>=0} (-1)^j \sum_{k >=0} a^k/b^k b^(-j*k)
  + sum_{j>=0} (-1)^j \sum_{k >=0} 1/b^k b^(-j*k)
= sum_{j>=0} (-1)^j b^(1+j)/[b^(1+j)-a]
  + sum_{j>=0} (-1)^j b^(1+j)/[b^(1+j)-1]
= sum_{j>=0} (-1)^j b^(1+j)*[2*b^(1+j)-1-a] / [ {b^(1+j)-a}*{b^(1+j)-1} ]
= sum_{j=0,2,4,..} b^(1+j)*[2*b^(1+j)-1-a] / [ {b^(1+j)-a}*{b^(1+j)-1} ]
  - sum_{j=1,3,5,..} b^(1+j)*[2*b^(1+j)-1-a] / [ {b^(1+j)-a}*{b^(1+j)-1} ]
= sum_{j>=0} b^(1+2j)*[2*b^(1+2j)-1-a] / [ {b^(1+2j)-a}*{b^(1+2j)-1} ]
  - sum_{j>=0} b^(2+2j)*[2*b^(2+2j)-1-a] / [ {b^(2+2j)-a}*{b^(2+2j)-1} ]
= sum_{j>=0} 
  b^(2j+1)*(b-1)*[ (1+a)*b^(4j+3) -2a*b^(2j+2)-2*a*b^(2j+1)+a(a+1)] 
    / [ {b^(2+2j)-a}*{b^(2+2j)-1}*{b^(1+2j)-a}*{b^(1+2j)-1} ]
The acceleration effect is that the leading powers in the
numerators of that j-sum are b^(6j+5) whereas they are b^(8j+6) in
the denominators, so convergence is ~ b^(-2j-1).

My question is: if I implement this in Maple, the result is
one less than the naive direct summation of the k-sum. So where in this
calculation did a constant term "1" vanish?

# compute sum_{k>=0) (a^k+1)/(b^k+1)
# by resumption, taking Nterms in the new sum
A := proc(a,b,Nterms)
    local val,j ;
    val :=0.0 ;
    for j from 0 to Nterms-1 do
        val := val+ b^(2*j+1)*(b-1)*( 
            (1+a)*b^(4*j+3) 
            -2*a*b^(2*j+2)
            -2*a*b^(2*j+1)
            +a*(a+1))
            / (b^(2+2*j)-a)
        /(b^(2+2*j)-1)
        /(b^(1+2*j)-a)
        /(b^(1+2*j)-1) ;
        val := evalf(val) ;
    end do:
    val ;
end proc:

# Evaluate convergence as with A372909, a=2, b=4:
# Should be 2.16249.. but is 1.16249... with the program:
Digits := 20;
A(2,4,10) ;
A(2,4,20) ;
A(2,4,30) ;

--
Richard