Hello, ratsum.mac computes an approximation of sum(a(n),n,n0,inf) with a(n)
a rational function. This is the code.
/* ratsum.mac v1.0 */
/* ratsum computes an approximation of sum(an,n,n0,inf) with a_n a rational
function.
The algorithm uses two facts:
(1) the series b_k = 1/( n(n+1)(n+2)...(n+k)) is telescopic and its summ is
1/k *k!.
(2) for any rational a_n we take b_k such that lim a(n)/b_k = c #0, c
constant. Then
a(n) - c *b(n) converges to zero more quickly than a(n),
write r_1(n) = a(n) - c*b(n),
iterating this process p times we obtain
a_n = c_1*b_1 + c_2*b_2 + ... c_s*b_s + r_p(n)
for n big, the term r_p decreaces very quickly (the difference between the
degree of numerator
and degree of denominator increase in each step) hence
sum(a_n) is approx a_1 + ... + a_terms + sum c_i b_i,n,terms+1,inf)
parameters:
p:17 is the number of times the acceleration process is done.
terms:30 is the partial summation of a_n used by the algorithm,
both parameters can be increased without problems to increase the
accuracy of the algorithm.
*/
principal(ex,n):= block([num,den,h1,h2,c1,c2],num:ratnum(ex),
den:ratdenom(ex),
h1:hipow(num,n),h2:hipow(den,n),
[ratcoef(num,n,h1)/ratcoef(den,n,h2),h2 - h1]);
ratsum(ratio,n,start):=block([ex,l,k,sb,terms,p],p:17,terms:30,
ex:rat(ratio),sb:0,
for i:1 thru p do
(l:principal(rat(ex),n),
k:l[2]-1,
ex:ex - l[1]*product(1/(n+i),i,0,k),
sb:sb + l[1]/k *product(1/i,i,terms+start,terms+start+k-1)),
sum(ratio,n,start,terms+start-1) + sb);
/* examples
(%i45) ratsum(1/(n*(n+1)),n,1);
(%o45) 1
(%i46) ratsum(1/(n*(n+1)*(n+2)),n,5);
(%o46) 1/60
(%i51) ratsum(1/n^2,n,1);
(%o51) 1696150783182270964333069992550871/
1031136029927435887508487994560000
(%i52) float(%);
(%o52) 1.644934066848226
(%i53) float(%pi^2/6);
(%o53) 1.644934066848226
*/