=================================================================== RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/trans3.c,v retrieving revision 1.1 retrieving revision 1.2 diff -u -p -r1.1 -r1.2 --- OpenXM_contrib/pari-2.2/src/basemath/Attic/trans3.c 2001/10/02 11:17:06 1.1 +++ OpenXM_contrib/pari-2.2/src/basemath/Attic/trans3.c 2002/09/11 07:26:53 1.2 @@ -1,4 +1,4 @@ -/* $Id: trans3.c,v 1.1 2001/10/02 11:17:06 noro Exp $ +/* $Id: trans3.c,v 1.2 2002/09/11 07:26:53 noro Exp $ Copyright (C) 2000 The PARI group. @@ -20,6 +20,11 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston, /** **/ /********************************************************************/ #include "pari.h" +extern GEN rpowsi(ulong a, GEN n, long prec); +extern GEN divrs2_safe(GEN x, long i); +extern void dcxlog(double s, double t, double *a, double *b); +extern double dnorm(double s, double t); +extern GEN trans_fix_arg(long *prec, GEN *s0, GEN *sig, gpmem_t *av, GEN *res); GEN cgetc(long l) @@ -36,34 +41,351 @@ fix(GEN x, long l) y = cgetr(l); gaffect(x, y); return y; } +long +lgcx(GEN z) +{ + long tz = typ(z); + + switch(tz) + { + case t_INT: case t_FRAC: case t_RFRAC: case t_QUAD: return BIGINT; + case t_REAL: return lg(z); + case t_COMPLEX: return min(lgcx((GEN)z[1]),lgcx((GEN)z[2])); + default: err(typeer,"lgcx"); + } + return 0; +} + +GEN +setlgcx(GEN z, long l) +{ + long tz = typ(z); + GEN p1; + + switch(tz) + { + case t_INT: case t_FRAC: case t_RFRAC: case t_QUAD: return z; + case t_REAL: p1 = cgetr(l); affrr(z,p1); return p1; + case t_COMPLEX: p1 = cgetc(l); gaffect(z,p1); return p1; + default: err(typeer,"setlgcx"); return gzero; /* not reached */ + } +} + +/* force z to be of type real/complex */ + +GEN +setlgcx2(GEN z, long l) +{ + long tz = typ(z); + GEN p1; + + switch(tz) + { + case t_INT: case t_FRAC: case t_RFRAC: case t_QUAD: case t_REAL: + p1 = cgetr(l); gaffect(z,p1); return p1; + case t_COMPLEX: p1 = cgetc(l); gaffect(z,p1); return p1; + default: err(typeer,"setlgcx"); return gzero; /* not reached */ + } +} + +/* a exporter ou ca existe deja ? */ + +long +isint(GEN n, long *ptk) +{ + long tn=typ(n); + GEN p1,p2; + + switch(tn) + { + case t_INT: + *ptk = itos(n); return 1; + case t_REAL: + p1 = gfloor(n); + if (gcmp(p1,n)==0) {*ptk = itos(p1); return 1;} + else return 0; + case t_FRAC: case t_FRACN: + p1 = dvmdii((GEN)n[1],(GEN)n[2],&p2); + if (!signe(p2)) {*ptk = itos(p1); return 1;} + else return 0; + case t_COMPLEX: + if (gcmp0((GEN)n[2])) return isint((GEN)n[1],ptk); + else return 0; + case t_QUAD: + if (gcmp0((GEN)n[3])) return isint((GEN)n[2],ptk); + else return 0; + default: err(typeer,"isint"); return 0; /* not reached */ + } +} + +double +norml1(GEN n, long prec) +{ + long tn=typ(n); + gpmem_t av; + double res; + + switch(tn) + { + case t_INT: case t_REAL: case t_FRAC: case t_FRACN: + return gtodouble(gabs(n,prec)); + case t_COMPLEX: + return norml1((GEN)n[1],prec)+norml1((GEN)n[2],prec); + case t_QUAD: + av = avma; res = norml1(gmul(n,realun(prec)),prec); avma = av; + return res; + default: err(typeer,"norml1"); return 0.; /* not reached */ + } +} + /***********************************************************************/ /** **/ -/** K BESSEL FUNCTION **/ +/** BESSEL FUNCTIONS **/ /** **/ /***********************************************************************/ +/* computes sum_{k=0}^m n!*((-1)^flag*z^2/4)^k/(k!*(k+n)!) */ + +static GEN +_jbessel(GEN n, GEN z, long flag, long m) +{ + long k, limit; + gpmem_t av; + GEN p1,s; + + p1 = gmul2n(gsqr(z),-2); if (flag & 1) p1 = gneg(p1); + if (typ(z) == t_SER) + { + long v = valp(z); + k = lg(p1)-2 - v; + if (v < 0) err(negexper,"jbessel"); + if (k <= 0) return gadd(gun, zeroser(varn(z), 2*v)); + p1 = gprec(p1, k); + } + s = gun; + av = avma; limit = stack_lim(av,1); + for (k=m; k>=1; k--) + { + s = gadd(gun,gdiv(gdivgs(gmul(p1,s),k),gaddsg(k,n))); + if (low_stack(limit,stack_lim(av,1))) + { + if (DEBUGMEM>1) err(warnmem,"jbessel"); + s = gerepilecopy(av, s); + } + } + return s; +} + +static GEN +jbesselintern(GEN n, GEN z, long flag, long prec) +{ + long tz=typ(z), i, lz, lim, k=-1, ki, precnew; + gpmem_t av=avma, tetpil; + double B,N,L,x; + GEN p1,p2,y,znew,nnew; + + switch(tz) + { + case t_REAL: case t_COMPLEX: + i = precision(z); if (i) prec = i; + if (isint(setlgcx2(n,prec),&ki)) + { + k = labs(ki); + p2 = gdiv(gpowgs(gmul2n(z,-1),k),mpfact(k)); + if ((flag&1) && (ki<0) && (k&1)) p2 = gneg(p2); + } + else p2 = gdiv(gpow(gmul2n(z,-1),n,prec),ggamma(gaddgs(n,1),prec)); + if (gcmp0(z)) return gerepilecopy(av, p2); + else + { + x = gtodouble(gabs(z,prec)); + L = x*1.3591409; + B = bit_accuracy(prec)*LOG2/(2*L); + N = 1 + B; +/* 3 Newton iterations are enough except in pathological cases */ + N = (N + B)/(log(N)+1); N = (N + B)/(log(N)+1); N = (N + B)/(log(N)+1); + lim = max((long)(L*N),2); + precnew = prec; + if (x >= 1.0) precnew += 1 + (long)(x/(LOG2*BITS_IN_LONG)); + znew = setlgcx(z,precnew); + if (k >= 0) p1 = setlgcx(_jbessel(stoi(k),znew,flag,lim),prec); + else + { + i = precision(n); + nnew = (i && (i < precnew)) ? setlgcx(n,precnew) : n; + p1 = setlgcx(_jbessel(nnew,znew,flag,lim),prec); + } + tetpil = avma; return gerepile(av,tetpil,gmul(p2,p1)); + } + + case t_SER: + if (isint(setlgcx2(n,prec),&ki)) + { + k = labs(ki); + p1 = _jbessel(stoi(k),z,flag,lg(z)-2); + } + else p1 = _jbessel(n,z,flag,lg(z)-2); + return gerepilecopy(av,p1); + + case t_VEC: case t_COL: case t_MAT: + lz=lg(z); y=cgetg(lz,typ(z)); + for (i=1; i=0) l = linit; + else l = lgcx(z) - 1 + ((-2*k*gz)>>TWOPOTBITS_IN_LONG); + if (l>prec) prec = l; + prec += (-gz)>>TWOPOTBITS_IN_LONG; + z = setlgcx(z,prec); + p1 = _jbesselh(k,z,prec); + p1 = gmul(gsqrt(gdiv(gmul2n(z,1),mppi(prec)),prec),p1); + tetpil = avma; return gerepile(av,tetpil,setlgcx(p1,linit)); + + case t_SER: + if (gcmp0(z)) return gpowgs(z,k); + av = avma; + if (valp(z) < 0) err(negexper,"jbesselh"); + z = gprec(z, lg(z)-2 + (2*k+1)*valp(z)); + p1 = gdiv(_jbesselh(k,z,prec),gpowgs(z,k)); + for (i=1; i<=k; i++) p1 = gmulgs(p1,2*i+1); + return gerepilecopy(av,p1); + + case t_VEC: case t_COL: case t_MAT: + lz=lg(z); y=cgetg(lz,typ(z)); + for (i=1; i>TWOPOTBITS_IN_LONG); + lnew = l + rab; prec += rab; + } + else lnew = l; + yfin=cgetr(l); l1=lnew+1; + av=avma; x = fix(gx, lnew); y=cgetr(lnew); u=cgetr(l1); v=cgetr(l1); c=cgetr(l1); d=cgetr(l1); e=cgetr(l1); f=cgetr(l1); nu2=gmulgs(gsqr(nu),-4); @@ -88,7 +410,7 @@ kbessel(GEN nu, GEN gx, long prec) } gmulz(s,zf,u); t=gmul2n(t,-1); gdivgsz(gadd(gmul(t,zf),gmul(u,nu)),-n2,v); avma=av1; - affsr(n2,q=cgetr(l1)); + q = stor(n2, l1); r=gmul2n(x,1); av1=avma; for(;;) { @@ -131,9 +453,257 @@ kbessel(GEN nu, GEN gx, long prec) } gmulz(s,zf,y); } - gdivz(y,mpexp(x),y); avma=av; return y; + gdivz(y,mpexp(x),yfin); + avma=av; return yfin; } +/* computes sum_{k=0}^m ((-1)^flag*z^2/4)^k/(k!*(k+n)!)*(H(k)+H(k+n)) + +sum_{k=0}^{n-1}((-1)^(flag+1)*z^2/4)^(k-n)*(n-k-1)!/k!. Ici n + doit etre un long. Attention, contrairement a _jbessel, pas de n! devant. + Quand flag > 1, calculer exactement les H(k) et factorielles */ + +static GEN +_kbessel(long n, GEN z, long flag, long m, long prec) +{ + long k, limit; + gpmem_t av; + GEN p1,p2,p3,s,*tabh; + + p1 = gmul2n(gsqr(z),-2); if (flag & 1) p1 = gneg(p1); + if (typ(z) == t_SER) + { + long v = valp(z); + k = lg(p1)-2 - v; + if (v < 0) err(negexper,"kbessel"); + if (k <= 0) return gadd(gun, zeroser(varn(z), 2*v)); + p1 = gprec(p1, k); + } + tabh = (GEN*)cgetg(m+n+2,t_VEC); tabh[1] = gzero; + if (flag <= 1) + { + s = realun(prec); tabh[2] = s; + for (k=2; k<=m+n; k++) + { + s = divrs(addsr(1,mulsr(k,s)),k); tabh[k+1] = s; + } + } + else + { + s = gun; tabh[2] = s; + for (k=2; k<=m+n; k++) + { + s = gdivgs(gaddsg(1,gmulsg(k,s)),k); tabh[k+1] = s; + } + } + s = gadd(tabh[m+1],tabh[m+n+1]); + av = avma; limit = stack_lim(av,1); + for (k=m; k>=1; k--) + { + s = gadd(gadd(tabh[k],tabh[k+n]),gdiv(gmul(p1,s),mulss(k,k+n))); + if (low_stack(limit,stack_lim(av,1))) + { + if (DEBUGMEM>1) err(warnmem,"kbessel"); + s = gerepilecopy(av, s); + } + } + p3 = (flag <= 1) ? mpfactr(n,prec) : mpfact(n); + s = gdiv(s,p3); + if (n) + { + p1 = gneg(ginv(p1)); + p2 = gmulsg(n,gdiv(p1,p3)); + for (k=n-1; k>=0; k--) + { + s = gadd(s,p2); + p2 = gmul(p2,gmul(mulss(k,n-k),p1)); + } + } + return s; +} + +static GEN +kbesselintern(GEN n, GEN z, long flag, long prec) +{ + long tz=typ(z), i, k, ki, lz, lim, precnew, fl, fl2, ex; + gpmem_t av=avma, tetpil; + double B,N,L,x,rab; + GEN p1,p2,y,p3,znew,nnew,pplus,pmoins,s,c; + + fl = (flag & 1) == 0; + switch(tz) + { + case t_REAL: case t_COMPLEX: + if (gcmp0(z)) err(talker,"zero argument in a k/n bessel function"); + i = precision(z); if (i) prec = i; + x = gtodouble(gabs(z,prec)); +/* Experimentally optimal on a PIII 500 Mhz. Your optimum may vary. */ + if ((x > (8*(prec-2)+norml1(n,prec)-3)) && !flag) return kbessel(n,z,prec); + precnew = prec; + if (x >= 1.0) + { + rab = x/(LOG2*BITS_IN_LONG); if (fl) rab *= 2; + precnew += 1 + (long)rab; + } + znew = setlgcx(z,precnew); + if (isint(setlgcx2(n,precnew),&ki)) + { + k = labs(ki); + L = x*1.3591409; + B = (bit_accuracy(prec)*LOG2)/(2*L); + if (fl) B += 0.367879; + N = 1 + B; +/* 3 Newton iterations are enough except in pathological cases */ + N = (N + B)/(log(N)+1); N = (N + B)/(log(N)+1); N = (N + B)/(log(N)+1); + lim = max((long)(L*N),2); + p1 = _kbessel(k,znew,flag,lim,precnew); + p1 = gmul(gpowgs(gmul2n(znew,-1),k),p1); + p2 = gadd(mpeuler(precnew),glog(gmul2n(znew,-1),precnew)); + p3 = jbesselintern(stoi(k),znew,flag,precnew); + p2 = gsub(gmul2n(p1,-1),gmul(p2,p3)); + p2 = setlgcx(p2,prec); + if (fl == 0) {p1 = mppi(prec); p2 = gmul2n(gdiv(p2,p1),1);} + fl = (fl && (k&1)) || (!fl && (ki>=0 || (k&1)==0)); + tetpil = avma; return gerepile(av,tetpil,fl ? gneg(p2) : gcopy(p2)); + } + else + { + i = precision(n); + nnew = (i && (i < precnew)) ? setlgcx(n,precnew) : n; + p2 = mppi(precnew); gsincos(gmul(nnew,p2),&s,&c,precnew); + ex = gexpo(s); + if (ex < 0) + { + rab = (-ex)/(LOG2*BITS_IN_LONG); if (fl) rab *= 2; + precnew += 1 + (long)rab; + } + nnew = (i && (i < precnew)) ? setlgcx(n,precnew) : n; + znew = setlgcx(znew,precnew); + p2 = mppi(precnew); gsincos(gmul(nnew,p2),&s,&c,precnew); + pplus = jbesselintern(nnew,znew,flag,precnew); + pmoins = jbesselintern(gneg(nnew),znew,flag,precnew); + if (fl) p1 = gmul(gsub(pmoins,pplus),gdiv(p2,gmul2n(s,1))); + else p1 = gdiv(gsub(gmul(c,pplus),pmoins),s); + tetpil = avma; return gerepile(av,tetpil,setlgcx(p1,prec)); + } + + case t_SER: + if (isint(setlgcx2(n,prec),&ki)) + { + k = labs(ki); + p1 = _kbessel(k,z,flag+2,lg(z)-2,prec); + return gerepilecopy(av,p1); + } + else + { + fl2 = isint(setlgcx2(gmul2n(n,1),prec),&ki); + if (!fl2) + err(talker,"cannot give a power series result in k/n bessel function"); + else + { + k = labs(ki); n = gmul2n(stoi(k),-1); + fl2 = (k&3)==1; + pmoins = jbesselintern(gneg(n),z,flag,prec); + if (fl) + { + pplus = jbesselintern(n,z,flag,prec); + p2 = gpowgs(z,-k); if (fl2 == 0) p2 = gneg(p2); + p3 = gmul2n(divii(mpfact(k + 1),mpfact((k + 1) >> 1)),-(k + 1)); + p3 = gdivgs(gmul2n(gsqr(p3),1),k); + p2 = gmul(p2,p3); + p1 = gsub(pplus,gmul(p2,pmoins)); + } + else p1 = pmoins; + tetpil = avma; + return gerepile(av,tetpil,fl2 ? gneg(p1) : gcopy(p1)); + } + } + + case t_VEC: case t_COL: case t_MAT: + lz=lg(z); y=cgetg(lz,typ(z)); + for (i=1; i= -1) p1=ghalf; @@ -225,7 +796,7 @@ hyperu(GEN a, GEN b, GEN gx, long prec) GEN kbessel2(GEN nu, GEN x, long prec) { - long av = avma,tetpil; + gpmem_t av = avma, tetpil; GEN p1,p2,x2,a,pitemp; if (typ(x)==t_REAL) prec = lg(x); @@ -242,7 +813,7 @@ GEN incgam(GEN s, GEN x, long prec) { GEN p1,z = cgetr(prec); - long av = avma; + gpmem_t av = avma; if (gcmp0(x)) return ggamma(s,prec); if (typ(x)!=t_REAL) { gaffect(x,z); x=z; } @@ -273,7 +844,8 @@ GEN incgam1(GEN a, GEN x, long prec) { GEN p2,p3,y, z = cgetr(prec); - long av=avma, av1, l,n,i; + long l, n, i; + gpmem_t av=avma, av1; double m,mx; if (typ(x) != t_REAL) { gaffect(x,z); x=z; } @@ -296,7 +868,8 @@ GEN incgam2(GEN a, GEN x, long prec) { GEN b,p1,p2,p3,y, z = cgetr(prec); - long av = avma,av1,l,n,i; + long l, n, i; + gpmem_t av = avma, av1; double m,mx; if (typ(x) != t_REAL) { gaffect(x,z); x=z; } @@ -328,7 +901,8 @@ GEN incgam3(GEN s, GEN x, long prec) { GEN b,p1,p2,p3,y, z = cgetr(prec); - long av = avma, av1,l,n,i; + long l, n, i; + gpmem_t av = avma, av1; if (typ(x) != t_REAL) { gaffect(x,z); x=z; } l=lg(x); n = -bit_accuracy(l)-1; @@ -364,7 +938,7 @@ GEN incgam4(GEN a, GEN x, GEN g, long prec) { GEN p1, z = cgetr(prec); - long av = avma; + gpmem_t av = avma; if (typ(x) != t_REAL) { gaffect(x,z); x=z; } if (gcmp(subrs(x,1),a) > 0 || gsigne(greal(a)) <= 0) @@ -384,7 +958,8 @@ incgam0(GEN a, GEN x, GEN z,long prec) GEN eint1(GEN x, long prec) { - long av = avma,tetpil,l,n,i; + long l, n, i; + gpmem_t av = avma, tetpil; GEN p1,p2,p3,p4,run,y; if (typ(x) != t_REAL) { gaffect(x,p1=cgetr(prec)); x = p1;} @@ -454,7 +1029,8 @@ eint1(GEN x, long prec) GEN veceint1(GEN C, GEN nmax, long prec) { - long av,av1,k,n,nstop,i,cd,nmin,G,a, chkpoint; + long k, n, nstop, i, cd, nmin, G, a, chkpoint; + gpmem_t av, av1; GEN y,e1,e2,F0,F,M2,f,den,minvn,mcn,p1,vdiff,ap,unr,zeror,deninit; if (!nmax) return eint1(C,prec); @@ -527,7 +1103,7 @@ veceint1(GEN C, GEN nmax, long prec) GEN gerfc(GEN x, long prec) { - long av=avma; + gpmem_t av=avma; GEN p1,p2; if (typ(x)!=t_REAL) { p1=cgetr(prec); gaffect(x,p1); x=p1; } @@ -548,7 +1124,8 @@ gerfc(GEN x, long prec) static GEN czeta(GEN s, long prec) { - long av,n,p,n1,flag1,flag2,i,i2; + long n, p, n1, flag1, flag2, i, i2; + gpmem_t av; double st,sp,sn,ssig,ns,alpha,beta,maxbeta,xinf; GEN y,z,res,sig,ms,p1,p2,p3,p31,pitemp; @@ -677,7 +1254,8 @@ next_bin(GEN y, long n, long k) static GEN izeta(long k, long prec) { - long av=avma,av2,kk,n,li, limit; + long kk, n, li; + gpmem_t av=avma, av2, limit; GEN y,p1,pi2,qn,z,q,binom; /* treat trivial cases */ @@ -689,7 +1267,7 @@ izeta(long k, long prec) return gerepileuptoleaf(av, divrs(y,k-1)); } if (k > bit_accuracy(prec)+1) return realun(prec); - pi2 = mppi(prec); setexpo(pi2,2); /* 2Pi */ + pi2 = Pi2n(1, prec); if ((k&1) == 0) { p1 = mulrr(gpuigs(pi2,k),absr(bernreal(k,prec))); @@ -715,18 +1293,17 @@ izeta(long k, long prec) y = mulrr(divrr(gpuigs(pi2,k),mpfactr(kk,prec)),y); av2 = avma; limit = stack_lim(av2,1); - qn = gsqr(q); z = divsr(1,addrs(q,-1)); + qn = gsqr(q); z = ginv( addrs(q,-1) ); for (n=2; ; n++) { - p1 = divsr(1,mulir(gpuigs(stoi(n),k),addrs(qn,-1))); + p1 = ginv( mulir(gpuigs(stoi(n),k),addrs(qn,-1)) ); z = addrr(z,p1); if (expo(p1)< li) break; qn = mulrr(qn,q); if (low_stack(limit,stack_lim(av2,1))) { - GEN *gptr[2]; gptr[0]=&z; gptr[1]=&qn; if (DEBUGMEM>1) err(warnmem,"izeta"); - gerepilemany(av2,gptr,2); + gerepileall(av2,2, &z, &qn); } } setexpo(z,expo(z)+1); @@ -757,9 +1334,8 @@ izeta(long k, long prec) qn=mulrr(qn,q); if (low_stack(limit,stack_lim(av2,1))) { - GEN *gptr[2]; gptr[0]=&z; gptr[1]=&qn; if (DEBUGMEM>1) err(warnmem,"izeta"); - gerepilemany(av2,gptr,2); + gerepileall(av2,2, &z, &qn); } } setexpo(z,expo(z)+1); @@ -780,13 +1356,13 @@ pows(long x, long n) /* return n^-s, n > 1 odd. tab[q] := q^-s, q prime power */ static GEN -n_s(long n, GEN *tab) +n_s(ulong n, GEN *tab) { byteptr d = diffptr + 2; GEN x = NULL; long p, e; - for (p = 3; n > 1; p += *d++) + for (p = 3; n > 1; ) { e = svaluation(n, p, &n); if (e) @@ -794,6 +1370,7 @@ n_s(long n, GEN *tab) GEN y = tab[pows(p,e)]; if (!x) x = y; else x = gmul(x,y); } + NEXT_PRIME_VIADIFF_CHECK(p,d); } return x; } @@ -804,7 +1381,7 @@ GEN czeta(GEN s0, long prec); static GEN izeta(long k, long prec) { - long av = avma; + gpmem_t av = avma; GEN y,p1,pi2; /* treat trivial cases */ @@ -827,12 +1404,6 @@ izeta(long k, long prec) return czeta(stoi(k), prec); } -extern GEN rpowsi(ulong a, GEN n, long prec); -extern GEN divrs2_safe(GEN x, long i); -extern void dcxlog(double s, double t, double *a, double *b); -extern double dnorm(double s, double t); -extern GEN trans_fix_arg(long *prec, GEN *s0, GEN *sig, long *av, GEN *res); - /* s0 a t_INT, t_REAL or t_COMPLEX. * If a t_INT, assume it's not a trivial case (i.e we have s0 > 1, odd) * */ @@ -841,11 +1412,12 @@ czeta(GEN s0, long prec) { GEN s, u, a, y, res, tes, sig, invn2, p1, unr; GEN sim, ms, s1, s2, s3, s4, s5, *tab, tabn; - long p, i, sqn, nn, lim, lim2, ct, av, av2 = avma, avlim; + long p, i, sqn, nn, lim, lim2, ct; + gpmem_t av, av2 = avma, avlim; int funeq = 0; byteptr d; - if (DEBUGLEVEL) timer2(); + if (DEBUGLEVEL>2) (void)timer2(); s = trans_fix_arg(&prec,&s0,&sig,&av,&res); if (gcmp0(s)) { y = gneg(ghalf); goto END; } if (signe(sig) <= 0 || expo(sig) < -1) @@ -888,7 +1460,7 @@ czeta(GEN s0, long prec) lim = (long) ceil(l); if (lim < 2) lim = 2; l2 = (lim+ssig/2.-.25); nn = 1 + (long)ceil( sqrt(l2*l2 + st*st/4) * la / PI ); - if (DEBUGLEVEL) fprintferr("lim, nn: [%ld, %ld]\n",lim,nn); + if (DEBUGLEVEL>2) fprintferr("lim, nn: [%ld, %ld]\n",lim,nn); if ((ulong)nn >= maxprime()) err(primer1); } prec++; unr = realun(prec); /* one extra word of precision */ @@ -897,43 +1469,48 @@ czeta(GEN s0, long prec) d = diffptr + 1; if (typ(s0) == t_INT) { /* no explog for 1/p^s */ - for (p=2; p < nn; p += *d++) + for (p=2; p < nn;) { tab[p] = divrr(unr, rpowsi(p, s0, prec)); + NEXT_PRIME_VIADIFF_CHECK(p,d); + } a = divrr(unr, rpowsi(nn, s0, prec)); } else { /* general case */ ms = gneg(s); p1 = cgetr(prec); - for (p=2; p < nn; p += *d++) + for (p=2; p < nn;) { affsr(p, p1); tab[p] = gexp(gmul(ms, mplog(p1)), prec); + NEXT_PRIME_VIADIFF_CHECK(p,d); } affsr(nn,p1); a = gexp(gmul(ms, mplog(p1)), prec); } sqn = (long)sqrt(nn-1.); d = diffptr + 2; /* fill in odd prime powers */ - for (p=3; p <= sqn; p += *d++) + for (p=3; p <= sqn; ) { ulong oldq = p, q = p*p; while (q<(ulong)nn) { tab[q] = gmul(tab[p], tab[oldq]); oldq = q; q *= p; } + NEXT_PRIME_VIADIFF_CHECK(p,d); } - if (DEBUGLEVEL) msgtimer("tab[q^-s] from 1 to N-1"); + if (DEBUGLEVEL>2) msgtimer("tab[q^-s] from 1 to N-1"); tabn = cgetg(nn, t_VECSMALL); ct = 0; for (i = nn-1; i; i>>=1) tabn[++ct] = (i-1)>>1; sim = y = unr; for (i=ct; i > 1; i--) { - long j, av2 = avma; + long j; + gpmem_t av2 = avma; for (j=tabn[i]+1; j<=tabn[i-1]; j++) sim = gadd(sim, n_s(2*j+1, tab)); sim = gerepileupto(av2, sim); y = gadd(sim, gmul(tab[2],y)); } y = gadd(y, gmul2n(a,-1)); - if (DEBUGLEVEL) msgtimer("sum from 1 to N-1"); + if (DEBUGLEVEL>2) msgtimer("sum from 1 to N-1"); invn2 = divrs(unr, nn*nn); lim2 = lim<<1; tes = bernreal(lim2, prec); @@ -976,7 +1553,7 @@ czeta(GEN s0, long prec) u = gmul(gmul(tes,invn2), gmul2n(s2, -1)); tes = gmulsg(nn, gaddsg(1, u)); } - if (DEBUGLEVEL) msgtimer("Bernoulli sum"); + if (DEBUGLEVEL>2) msgtimer("Bernoulli sum"); /* y += tes a / (s-1) */ y = gadd(y, gmul(tes, gdiv(a, gsub(s, unr)))); @@ -984,7 +1561,8 @@ END: if (funeq) { GEN pitemp = mppi(prec); setexpo(pitemp,2); /* 2Pi */ - y = gmul(gmul(y, ggamma(s,prec)), gpui(pitemp,gneg(s),prec)); + y = gmul(gmul(y, ggamma(gprec_w(s,prec-1),prec)), + gpui(pitemp,gneg(s),prec)); setexpo(pitemp,0); /* Pi/2 */ y = gmul2n(gmul(y, gcos(gmul(pitemp,s),prec)), 1); } @@ -1020,7 +1598,8 @@ gzeta(GEN x, long prec) void gzetaz(GEN x, GEN y) { - long av=avma, prec = precision(y); + long prec = precision(y); + gpmem_t av=avma; if (!prec) err(infprecer,"gzetaz"); gaffect(gzeta(x,prec),y); avma=av; @@ -1038,7 +1617,8 @@ gzetaz(GEN x, GEN y) static GEN cxpolylog(long m, GEN x, long prec) { - long av=avma,li,i,n,bern_upto; + long li, i, n, bern_upto; + gpmem_t av=avma; GEN p1,z,h,q,s; if (gcmp1(x)) return izeta(m,prec); @@ -1070,7 +1650,8 @@ cxpolylog(long m, GEN x, long prec) GEN polylog(long m, GEN x, long prec) { - long av,av1,limpile,l,e,i,G,sx; + long l, e, i, G, sx; + gpmem_t av, av1, limpile; GEN X,z,p1,p2,n,y,logx; if (m<0) err(talker,"negative index in polylog"); @@ -1142,7 +1723,8 @@ dilog(GEN x, long prec) GEN polylogd0(long m, GEN x, long flag, long prec) { - long k,l,av,fl,m2; + long k, l, fl, m2; + gpmem_t av; GEN p1,p2,p3,y; m2=m&1; av=avma; @@ -1188,7 +1770,8 @@ polylogdold(long m, GEN x, long prec) GEN polylogp(long m, GEN x, long prec) { - long k,l,av,fl,m2; + long k, l, fl, m2; + gpmem_t av; GEN p1,y; m2=m&1; av=avma; @@ -1234,7 +1817,8 @@ polylogp(long m, GEN x, long prec) GEN gpolylog(long m, GEN x, long prec) { - long i,lx,av=avma,v,n; + long i, lx, v, n; + gpmem_t av=avma; GEN y,p1,p2; if (m<=0) @@ -1272,6 +1856,7 @@ gpolylog(long m, GEN x, long prec) p1=glog(gsub(gun,x),prec); return gerepileupto(av, gneg(p1)); } + if (gcmp0(x)) return gcopy(x); if (valp(x)<=0) err(impl,"polylog around a!=0"); v=varn(x); n=(lg(x)-2)/valp(x); y=ggrando(polx[v],lg(x)-2); for (i=n; i>=1; i--) @@ -1294,7 +1879,8 @@ gpolylog(long m, GEN x, long prec) void gpolylogz(long m, GEN x, GEN y) { - long av=avma, prec = precision(y); + long prec = precision(y); + gpmem_t av=avma; if (!prec) err(infprecer,"gpolylogz"); gaffect(gpolylog(m,x,prec),y); avma=av; @@ -1322,18 +1908,15 @@ qq(GEN x, long prec) if (tx==t_PADIC) return x; if (is_scalar_t(tx)) { - long l=precision(x); - GEN p1,p2,q; - - if (tx != t_COMPLEX || gcmp((GEN)x[2],gzero)<=0) + long l = precision(x); + if (!l) l = prec; + if (tx != t_COMPLEX || gsigne((GEN)x[2]) <= 0) err(talker,"argument must belong to upper half-plane"); - if (!l) l=prec; p1=mppi(l); setexpo(p1,2); - p2=cgetg(3,t_COMPLEX); p2[1]=zero; p2[2]=(long)p1; /* 2*I*Pi */ - q=gmul(x,p2); return gexp(q,prec); + return gexp(gmul(x, PiI2(l)), l); /* e(x) */ } - if (tx != t_POL && tx != t_SER) - err(talker,"bad argument for modular function"); + if (tx == t_SER) return x; + if (tx != t_POL) err(talker,"bad argument for modular function"); return tayl(x,gvar(x),precdl); } @@ -1347,8 +1930,7 @@ inteta(GEN q) y=gun; qn=gun; ps=gun; if (tx==t_PADIC) { - if (valp(q) <= 0) - err(talker,"non-positive valuation in inteta"); + if (valp(q) <= 0) err(talker,"non-positive valuation in eta"); for(;;) { p1=gneg_i(gmul(ps,gmul(q,gsqr(qn)))); @@ -1358,31 +1940,31 @@ inteta(GEN q) } else { - long l = 0, v = 0; /* gcc -Wall */ - long av = avma, lim = stack_lim(av,3); + long l, v = 0; /* gcc -Wall */ + gpmem_t av = avma, lim = stack_lim(av, 3); if (is_scalar_t(tx)) l = -bit_accuracy(precision(q)); else { - v=gvar(q); tx = 0; - if (valp(q) <= 0) - err(talker,"non-positive valuation in inteta"); + v = gvar(q); l = lg(q)-2; tx = 0; + if (valp(q) <= 0) err(talker,"non-positive valuation in eta"); } for(;;) { - p1=gneg_i(gmul(ps,gmul(q,gsqr(qn)))); - y=gadd(y,p1); qn=gmul(qn,q); ps=gmul(p1,qn); - y=gadd(y,ps); + p1 = gneg_i(gmul(ps,gmul(q,gsqr(qn)))); + /* qn = q^n + * ps = (-1)^n q^(n(3n+1)/2) + * p1 = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */ + y = gadd(y,p1); qn = gmul(qn,q); ps = gmul(p1,qn); + y = gadd(y,ps); if (tx) { if (gexpo(ps)-gexpo(y) < l) return y; } else - { if (gval(ps,v)>=precdl) return y; } + { if (gval(ps,v) >= l) return y; } if (low_stack(lim, stack_lim(av,3))) { - GEN *gptr[3]; - if(DEBUGMEM>1) err(warnmem,"inteta"); - gptr[0]=&y; gptr[1]=&qn; gptr[2]=&ps; - gerepilemany(av,gptr,3); + if(DEBUGMEM>1) err(warnmem,"eta"); + gerepileall(av, 3, &y, &qn, &ps); } } } @@ -1391,7 +1973,7 @@ inteta(GEN q) GEN eta(GEN x, long prec) { - long av = avma; + gpmem_t av = avma; GEN q = qq(x,prec); return gerepileupto(av,inteta(q)); } @@ -1400,15 +1982,15 @@ eta(GEN x, long prec) GEN trueeta(GEN x, long prec) { - long tx=typ(x), av=avma, tetpil,l; + long tx=typ(x), l; + gpmem_t av=avma, tetpil; GEN p1,p2,q,q24,n,z,m,unapprox; if (!is_scalar_t(tx)) err(typeer,"trueeta"); if (tx != t_COMPLEX || gsigne((GEN)x[2])<=0) err(talker,"argument must belong to upper half-plane"); l=precision(x); if (l) prec=l; - p1=mppi(prec); setexpo(p1,2); - p2=cgetg(3,t_COMPLEX); p2[1]=zero; p2[2]=(long)p1; /* 2*I*Pi */ + p2 = PiI2(prec); z=gexp(gdivgs(p2,24),prec); /* exp(2*I*Pi/24) */ unapprox=gsub(gun,gpuigs(stoi(10),-8)); m=gun; @@ -1434,7 +2016,8 @@ eta0(GEN x, long flag,long prec) GEN jell(GEN x, long prec) { - long av=avma,tetpil,tx = typ(x); + long tx = typ(x); + gpmem_t av=avma, tetpil; GEN p1; if (!is_scalar_t(tx) || tx == t_PADIC) @@ -1456,7 +2039,7 @@ jell(GEN x, long prec) GEN wf2(GEN x, long prec) { - long av=avma,tetpil; + gpmem_t av=avma, tetpil; GEN p1,p2; p1=gsqrt(gdeux,prec); @@ -1468,7 +2051,7 @@ wf2(GEN x, long prec) GEN wf1(GEN x, long prec) { - long av=avma,tetpil; + gpmem_t av=avma, tetpil; GEN p1,p2; p1=trueeta(gmul2n(x,-1),prec); p2=trueeta(x,prec); @@ -1479,7 +2062,7 @@ wf1(GEN x, long prec) GEN wf(GEN x, long prec) { - long av=avma,tetpil; + gpmem_t av=avma, tetpil; GEN p1,p2; p1=gdiv(trueeta(gmul2n(gaddgs(x,1),-1),prec),trueeta(x,prec)); @@ -1506,7 +2089,7 @@ sagm(GEN x, long prec) { GEN p1,a,b,a1,b1,y; long l,ep; - ulong av; + gpmem_t av; if (gcmp0(x)) return gcopy(x); switch(typ(x)) @@ -1570,7 +2153,8 @@ sagm(GEN x, long prec) GEN agm(GEN x, GEN y, long prec) { - long av,tetpil, ty=typ(y); + long ty=typ(y); + gpmem_t av, tetpil; GEN z; if (is_matvec_t(ty)) @@ -1587,7 +2171,8 @@ agm(GEN x, GEN y, long prec) GEN logagm(GEN q) { - long av=avma,prec=lg(q),tetpil,s,n,lim; + long prec=lg(q), s, n, lim; + gpmem_t av=avma, tetpil; GEN y,q4,q1; if (typ(q)!=t_REAL) err(typeer,"logagm"); @@ -1606,7 +2191,7 @@ logagm(GEN q) GEN glogagm(GEN x, long prec) { - long av,tetpil; + gpmem_t av, tetpil; GEN y,p1,p2; switch(typ(x)) @@ -1644,9 +2229,13 @@ glogagm(GEN x, long prec) GEN theta(GEN q, GEN z, long prec) { - long av=avma,tetpil,l,n; + long l, n; + gpmem_t av=avma, tetpil; GEN ps,qn,qnold,y,zy,lq,ps2,p1,k,zold; + if (!is_scalar_t(typ(q)) || !is_scalar_t(typ(z))) + err(impl,"theta for non-scalar types"); + l=precision(q); if (l) prec=l; p1=realun(prec); z=gmul(p1,z); if (!l) q=gmul(p1,q); @@ -1681,7 +2270,8 @@ theta(GEN q, GEN z, long prec) GEN thetanullk(GEN q, long k, long prec) { - long av=avma,tetpil,l,n; + long l, n; + gpmem_t av=avma, tetpil; GEN p1,ps,qn,y,ps2; l=precision(q); @@ -1703,65 +2293,3 @@ thetanullk(GEN q, long k, long prec) return gerepile(av,tetpil,gmul(p1,y)); } -GEN -jbesselh(GEN n, GEN z, long prec) -{ - long av,tetpil,k,l,i,lz; - GEN y,p1,p2,zinv,p0,s,c; - - if (typ(n)!=t_INT) err(talker,"not an integer index in jbesselh"); - k=itos(n); if (k<0) err(impl,"ybesselh"); - - switch(typ(z)) - { - case t_REAL: case t_COMPLEX: - if (gcmp0(z)) return gzero; - av=avma; zinv=ginv(z); - l=precision(z); if (l>prec) prec=l; - gsincos(z,&s,&c,prec); - p1=gmul(zinv,s); - if (k) - { - p0=p1; p1=gmul(zinv,gsub(p0,c)); - for (i=2; i<=k; i++) - { - p2=gsub(gmul(gmulsg(2*i-1,zinv),p1),p0); - p0=p1; p1=p2; - } - } - p2=gsqrt(gdiv(gmul2n(z,1),mppi(prec)),prec); - tetpil=avma; return gerepile(av,tetpil,gmul(p2,p1)); - - case t_VEC: case t_COL: case t_MAT: - lz=lg(z); y=cgetg(lz,typ(z)); - for (i=1; i