numerical summation of series, rational case.

/* sumninf(a(n),p)  computes an approximate value of 
sum(a(n),n,1,inf) with a(n) a rational function.

It accelerate the convergence rate of a(n): 
there is b(n) such that lim a(n)/b(n) = 1 then

a(n) - b(n) converges to zero more quickly than a(n) and 
sum(b(n),n,1,inf)  known.

Repeating this procedure we speed up the covergence rate,
 p is the number of times we speed up a(n).

Finally sum(a(n),1,inf) = sum(a(n),1,20) + sum(a(n),20,inf)
 aprox
sum(a(n),n,1,20)+ sum(b(n),n,20,inf)
because r(n) =  a(n) - b(n) converges very quickly to zero
 this procedure  gives a good approximation.  */

--source code -------------------------------------------

:lisp (defun principal(ex) 
(let* ((num (cadr ex)) (den (cddr ex)) (c1) (c2 (third den))
 (deg) (g2 (second den)))
(if (atom  num)  (setq c1  num g1 0) 
(setq c1  (third num) g1 (second num)))
(setq deg (- g2 g1))  (list '(MLIST SIMP) c1 c2 deg))) ;

g1(ex):= block([l],l:?principal(rat(ex)), ex -
l[1]/l[2]*product(1/(n+i),i,0,l[3]-1));

matchdeclare(x_,integerp);

defrule(uno,1/(n+x_),sum(-1/k,k,1,x_+20));


sumninf(ex0,p):=block([ex],ex:rat(ex0),
	for i:1 thru p do ex:g1(ex),
	sum(ex0,n,1,20)+
            apply1(partfrac(ex0 - ex,n),uno));
---------------------------------------------------

/* examples 


(%i87) sum(1/n^2,n,1,inf),simpsum,numer;
(%o87) 1.644934066848226

(%i73) sumninf(1/n^2,14),float;
(%o73) 1.644934066846179


(%i86) sum(1/n^4,n,1,inf),simpsum,numer;
(%o86) 1.082323233711138

(%i83) sumninf(1/n^4,15),float;
(%o83) 1.082323233717282

(%i94) y:rat(1/n^2);
(%o94) 1/n^2
(%i96) g1(%);
(%o96) 2/(n^4+3*n^3+2*n^2)
(%i97) g1(%);
(%o97) 6/(n^5+6*n^4+11*n^3+6*n^2)
(%i98) g1(%);
(%o98) 24/(n^6+10*n^5+35*n^4+50*n^3+24*n^2)
(%i99) partfrac(%-y,n);
(%o99) -1/(4*(n+4))+4/(3*(n+3))-
3/(n+2)+4/(n+1)-25/(12*n)
*/