cino hilliard on Wed, 02 Sep 2009 16:14:24 +0200

[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

 RE: Euler-Maclaurin sumation function

 > 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 scriptto 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, itis 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)]))