RE: Euler-Maclaurin sumation function

[cino hilliard <hillcino368@hotmail.com>]

 
> Date: Wed, 2 Sep 2009 11:50:58 +0200
> From: Bill.Allombert@math.u-bordeaux1.fr
> To: pari-users@list.cr.yp.to
> Subject: Re: Euler-Maclaurin sumation function
>
> On Sat, Aug 01, 2009 at 11:34:01PM -0500, cino hilliard wrote:
> >
> > Hi,
> >
> > I understand the Euler-Maclaurin sumation formula is used in the Pari
> > zeta() function. I checked the code but do not quite understand it.
> >
> >
> >
> > Can someone build a Pari script that will do this. Mathematica has a
> > function Nsum that does this but I can't afford even the $295 price
> > for the home edition.
>
> ... a Pari script that will do what ?

Just as noted above an Nsum function as in Mathematica or an eulermac function
as in Maple to approximate the sums of functions. Actually, I would prefer a built-in
function eulermac(x=a,b,f(x)) which would be the  Euler-Maclaurin sumation formula 
as is also noted above to quickly approximate the sum of functions like (x-1)/log(x), x=2,n.
The sum() function in Pari is too slow for large n.
 
A poster in the primenumbers group was kind enough to write the script template
below which I modified slightly to get more output. I am able to modify the script
to do various functions. However, I was unable to generalize it into a crisp eulermac()
function as in Mathematica and Maple. 

As previously noted, a built-in, optimized function would be preferred. Yes, it
is quite a challenge but it would greatly enhance the functionality of Pari.
 
Cheers and Roebuck,
 
Cino Hilliard

 
******************************Code**************************************
\\ f(x,lx) = x/lx - 1/lx + a/lx^2 - b/lx^3; \\ modify at will
 f(x,lx) = (x-1)/lx;
 a=1; b=1; m=2; \\ chosen values
 \p60;
 terms=12;trunc=10^2;g=f(x,lx);d=[];
 bv=bernvec(terms+1);
 {
 for(n=1,2*terms-1, \\ store odd derivatives
 g=deriv(g,x)+1/x*deriv(g,lx);if(n%2,d=concat(d,g)));
 }
\\ fx(x)=f(x,log(x));
   fx(x)=f(x,log(x));
\\ dx(y)=subst(subst(d,x,y),lx,log(y));
 dx(y)=subst(subst(d,x,y),lx,log(y));
 emac(m,n) =
 {
 local(dm,dn);
 dm=dx(m);dn=dx(n);intnum(x=m,n,fx(x))+(fx(n)+fx(m))/2
 +sum(k=1,terms,bv[k+1]/((2*k)!)*(dn[k]-dm[k]));
 }
\\ s=sum(k=2,trunc-1,fx(k));
   s=sum(k=2,trunc-1,fx(k));
 sm(m,n) = \\ as requested
 {
 if(n>trunc&&trunc>m&&m>1,s-sum(k=2,m-1,fx(k))+emac(trunc,n));
 }
for(k=1,40,print([k,sm(m,sqrtint(10^k))]))
\\   for(k=1,20,print([k,sm(m, 10^k)]))