=================================================================== RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/trans2.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/trans2.c 2001/10/02 11:17:06 1.1 +++ OpenXM_contrib/pari-2.2/src/basemath/Attic/trans2.c 2002/09/11 07:26:53 1.2 @@ -1,4 +1,4 @@ -/* $Id: trans2.c,v 1.1 2001/10/02 11:17:06 noro Exp $ +/* $Id: trans2.c,v 1.2 2002/09/11 07:26:53 noro Exp $ Copyright (C) 2000 The PARI group. @@ -15,7 +15,7 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston, /********************************************************************/ /** **/ -/** TRANSCENDENTAL FONCTIONS **/ +/** TRANSCENDENTAL FUNCTIONS **/ /** (part 2) **/ /** **/ /********************************************************************/ @@ -33,15 +33,13 @@ static GEN mpach(GEN x); static GEN mpatan(GEN x) { - long l,l1,l2,n,m,u,i,av0,av,lp,e,sx,s; + long l, l1, l2, n, m, u, i, lp, e, sx, s; + gpmem_t av0, av; double alpha,beta,gama=1.0,delta,fi; GEN y,p1,p2,p3,p4,p5,unr; sx=signe(x); - if (!sx) - { - y=cgetr(3); y[1]=x[1]; y[2]=0; return y; - } + if (!sx) return realzero_bit(expo(x)); l = lp = lg(x); if (sx<0) setsigne(x,1); u = cmprs(x,1); @@ -126,7 +124,7 @@ mpatan(GEN x) GEN gatan(GEN x, long prec) { - long av,tetpil; + gpmem_t av, tetpil; GEN y,p1; switch(typ(x)) @@ -161,7 +159,8 @@ gatan(GEN x, long prec) void gatanz(GEN x, GEN y) { - long av=avma,prec = precision(y); + long prec = precision(y); + gpmem_t av=avma; if (!prec) err(infprecer,"gatanz"); gaffect(gatan(x,prec),y); avma=av; @@ -177,7 +176,8 @@ gatanz(GEN x, GEN y) static GEN mpasin(GEN x) { - long l,u,v,av; + long l, u, v; + gpmem_t av; GEN y,p1; u=cmprs(x,1); v=cmpsr(-1,x); @@ -199,13 +199,14 @@ mpasin(GEN x) GEN gasin(GEN x, long prec) { - long av,tetpil,l,sx; + long l, sx; + gpmem_t av, tetpil; GEN y,p1; switch(typ(x)) { case t_REAL: sx=signe(x); - if (!sx) { y=cgetr(3); y[1]=x[1]; y[2]=0; return y; } + if (!sx) return realzero_bit(expo(x)); if (sx<0) setsigne(x,1); if (cmpsr(1,x)>=0) { setsigne(x,sx); return mpasin(x); } @@ -248,7 +249,8 @@ gasin(GEN x, long prec) void gasinz(GEN x, GEN y) { - long av=avma, prec = precision(y); + long prec = precision(y); + gpmem_t av=avma; if (!prec) err(infprecer,"gasinz"); gaffect(gasin(x,prec),y); avma=av; @@ -264,7 +266,8 @@ gasinz(GEN x, GEN y) static GEN mpacos(GEN x) { - long l,u,v,sx,av; + long l, u, v, sx; + gpmem_t av; GEN y,p1,p2; u=cmprs(x,1); v=cmpsr(-1,x); sx = signe(x); @@ -274,12 +277,7 @@ mpacos(GEN x) y=mppi(2-l); setexpo(y,0); return y; } l=lg(x); - if (!u) - { - y = cgetr(3); - y[1] = evalexpo(-(bit_accuracy(l)>>1)); - y[2] = 0; return y; - } + if (!u) return realzero_bit( -(bit_accuracy(l)>>1) ); if (!v) return mppi(l); y=cgetr(l); av=avma; @@ -309,7 +307,8 @@ mpacos(GEN x) GEN gacos(GEN x, long prec) { - long av,tetpil,l,sx; + long l, sx; + gpmem_t av, tetpil; GEN y,p1; switch(typ(x)) @@ -355,7 +354,8 @@ gacos(GEN x, long prec) void gacosz(GEN x, GEN y) { - long av=avma, prec = precision(y); + long prec = precision(y); + gpmem_t av=avma; if (!prec) err(infprecer,"gacosz"); gaffect(gacos(x,prec),y); avma=av; @@ -377,11 +377,7 @@ mparg(GEN x, GEN y) sx=signe(x); sy=signe(y); if (!sy) { - if (sx>0) - { - theta=cgetr(3); theta[1]=y[1]-expo(x); - theta[2]=0; return theta; - } + if (sx>0) return realzero_bit(expo(y) - expo(x)); return mppi(lg(x)); } prec = lg(y); if (prec=0) return mpach(x); + if (gcmpgs(x,1) >= 0) return mpach(x); - y=cgetg(3,t_COMPLEX); - if (gcmpgs(x,-1)>=0) + y = cgetg(3,t_COMPLEX); + if (gcmpgs(x,-1) >= 0) { - y[2]=lmpacos(x); y[1]=zero; - return y; + y[1] = zero; + y[2] = lmpacos(x); return y; } - av=avma; p1=mpach(gneg_i(x)); tetpil=avma; - y[1]=lpile(av,tetpil,gneg(p1)); - y[2]=lmppi(lg(x)); - return y; + av = avma; + y[1] = lpileupto(av, gneg(mpach(gneg_i(x)))); + y[2] = lmppi(lg(x)); return y; case t_COMPLEX: - av=avma; p1=gaddsg(-1,gsqr(x)); - p1=gadd(x,gsqrt(p1,prec)); tetpil=avma; - y=glog(p1,prec); - if (signe(y[2])<0) { tetpil=avma; y=gneg(y); } - return gerepile(av,tetpil,y); + av = avma; p1 = gaddsg(-1,gsqr(x)); + p1 = gadd(x, gsqrt(p1,prec)); /* x + sqrt(x^2-1) */ + y = glog(p1,prec); + if (signe(y[2]) < 0) y = gneg(y); + return gerepileupto(av, y); case t_SER: - av=avma; if (valp(x)<0) err(negexper,"gach"); - p1=gdiv(derivser(x),gsqrt(gsubgs(gsqr(x),1),prec)); - y=integ(p1,varn(x)); - if (!valp(x) && gcmp1((GEN)x[2])) return gerepileupto(av,y); - if (valp(x)) + av = avma; v = valp(x); + if (v < 0) err(negexper,"gach"); + if (gcmp0(x)) { - p1=cgetg(3,t_COMPLEX); p1[1]=zero; p1[2]=lmppi(prec); - setexpo(p1[2],0); + if (!v) return gcopy(x); + return gerepileupto(av, gadd(x, PiI2n(prec,-1))); } - else p1=gach((GEN)x[2],prec); - tetpil=avma; return gerepile(av,tetpil,gadd(p1,y)); + p1 = gdiv(derivser(x), gsqrt(gsubgs(gsqr(x),1),prec)); + y = integ(p1, varn(x)); + if (v) + p1 = PiI2n(prec, -1); /* I Pi/2 */ + else + { + p1 = (GEN)x[2]; + if (gcmp1(p1)) return gerepileupto(av,y); + p1 = gach(p1, prec); + } + return gerepileupto(av, gadd(p1,y)); case t_INTMOD: case t_PADIC: err(typeer,"gach"); @@ -778,7 +775,8 @@ gach(GEN x, long prec) void gachz(GEN x, GEN y) { - long av=avma,prec = precision(y); + long prec = precision(y); + gpmem_t av=avma; if (!prec) err(infprecer,"gachz"); gaffect(gach(x,prec),y); avma=av; @@ -790,17 +788,14 @@ gachz(GEN x, GEN y) /** **/ /********************************************************************/ +/* |x| < 1 */ static GEN mpath(GEN x) { - long av; + gpmem_t av; GEN y,p1; - if (!signe(x)) - { - y=cgetr(3); y[1]=x[1]; y[2]=0; - return y; - } + if (!signe(x)) return realzero_bit(expo(x)); y=cgetr(lg(x)); av=avma; p1 = addrs(divsr(2,subsr(1,x)),-1); affrr(mplog(p1), y); avma=av; @@ -810,7 +805,7 @@ mpath(GEN x) GEN gath(GEN x, long prec) { - long av,tetpil; + gpmem_t av, tetpil; GEN y,p1; switch(typ(x)) @@ -847,7 +842,8 @@ gath(GEN x, long prec) void gathz(GEN x, GEN y) { - long av=avma, prec = precision(y); + long prec = precision(y); + gpmem_t av=avma; if (!prec) err(infprecer,"gathz"); gaffect(gath(x,prec),y); avma=av; @@ -863,7 +859,8 @@ gathz(GEN x, GEN y) void mpbern(long nb, long prec) { - long n,m,i,j,d,d1,d2,l,av,av2,code0; + long n, m, i, j, d, d1, d2, l, code0; + gpmem_t av, av2; GEN p1,p2, B; if (nb < 0) nb = 0; @@ -888,7 +885,7 @@ mpbern(long nb, long prec) { fprintferr("caching Bernoulli numbers 2*%ld to 2*%ld, prec = %ld\n", i,nb,prec); - timer2(); + (void)timer2(); } p2 = p1; av2=avma; @@ -926,7 +923,7 @@ bernreal(long n, long prec) { GEN B; - if (n==1) { B=cgetr(prec); affsr(-1,B); setexpo(B,-1); return B; } + if (n==1) { B = stor(-1, prec); setexpo(B,-1); return B; } if (n<0 || n&1) return gzero; n >>= 1; mpbern(n+1,prec); B=cgetr(prec); affrr(bern(n),B); return B; @@ -936,22 +933,26 @@ bernreal(long n, long prec) static GEN bernfracspec(long k) { - long n,av,lim, K = k+1; - GEN s,c,N,h; + ulong n, K = k+1; + gpmem_t av, lim; + GEN s, c, N, b; - c = N = stoi(K); s = gun; h = gzero; + c = N = utoi(K); s = gun; b = gzero; av = avma; lim = stack_lim(av,2); - for (n=2; ; n++) + for (n=2; ; n++) /* n <= k+1 */ { - c = gdivgs(gmulgs(c,n-k-2),n); - h = gadd(h, gdivgs(gmul(c,s), n)); - if (n==K) return gerepileupto(av,h); - N[2] = n; s = addii(s, gpuigs(N,k)); + c = diviiexact(muliu(c,k+2-n), utoi(n)); + if (n & 1) setsigne(c, 1); else setsigne(c, -1); + /* c = (-1)^(n-1) binomial(k+1, n), s = 1^k + ... + (n-1)^k */ + + b = gadd(b, gdivgs(mulii(c,s), n)); + if (n == K) return gerepileupto(av, b); + + N[2] = n; s = addii(s, gpowgs(N,k)); if (low_stack(lim, stack_lim(av,2))) { - GEN *gptr[3]; gptr[0]=&c; gptr[1]=&h; gptr[2]=&s; if (DEBUGMEM>1) err(warnmem,"bernfrac"); - gerepilemany(av,gptr,3); + gerepileall(av,3, &c,&b,&s); } } } @@ -960,45 +961,35 @@ GEN bernfrac(long k) { if (!k) return gun; - if (k==1) return gneg(ghalf); - if (k<0 || k&1) return gzero; + if (k == 1) return gneg(ghalf); + if (k < 0 || k & 1) return gzero; return bernfracspec(k); } -static GEN -bernvec2(long k) -{ - GEN B = cgetg(k+2,t_VEC); - long i; - - for (i=1; i<=k; i++) - B[i+1]=(long)bernfracspec(i<<1); - B[1]=un; return B; -} - /* mpbern as exact fractions */ GEN bernvec(long nb) { - long n,m,i,j,d1,d2,av,tetpil; - GEN p1,y; + GEN y = cgetg(nb+2,t_VEC); + long n, i; - if (nb < 300) return bernvec2(nb); + if (nb > 46340 && BITS_IN_LONG == 32) err(impl, "bernvec for n > 46340"); - y=cgetg(nb+2,t_VEC); y[1]=un; - for (i=1; i<=nb; i++) - { - av=avma; p1=gzero; - n=8; m=5; d1=i-1; d2=2*i-3; - for (j=d1; j>0; j--) - { - p1 = gmulsg(n*m, gadd(p1,(GEN)y[j+1])); - p1 = gdivgs(p1, d1*d2); - n+=4; m+=2; d1--; d2-=2; + y[1] = un; + for (n = 1; n <= nb; n++) + { /* compute y[n+1] = B_{2n} */ + gpmem_t av = avma; + GEN b = gmul2n(stoi(1-2*n), -1); /* 1 + (2n+1)B_1 = (1-2n) /2 */ + GEN c = gun; + ulong u1 = 2*n + 1, u2 = n, d1 = 1, d2 = 1; + + for (i = 1; i < n; i++) + { + c = diviiexact(muliu(c, u1*u2), utoi(d1*d2)); /* = binomial(2n+1, 2*i) */ + b = gadd(b, gmul(c, (GEN)y[i+1])); + u1 -= 2; u2--; d1++; d2 += 2; } - p1 = gsubsg(1,gdivgs(gaddsg(1,p1),2*i+1)); - tetpil=avma; p1 = gmul2n(p1,-2*i); - y[i+1] = lpile(av,tetpil,p1); + y[n+1] = lpileupto(av, gdivgs(b, -(1+2*n))); } return y; } @@ -1013,7 +1004,8 @@ static GEN mpgamma(GEN x) { GEN y,p1,p2,p3,p4,p5,p6,p7,p71,pitemp; - long l,l1,l2,u,i,k,e,s,sx,n,p,av,av1; + long l, l1, l2, u, i, k, e, s, sx, n, p; + gpmem_t av, av1; double alpha,beta,dk; sx=signe(x); if (!sx) err(gamer2); @@ -1087,7 +1079,8 @@ static GEN cxgamma(GEN x, long prec) { GEN y,p1,p2,p3,p4,p5,p6,p7,p71,pitemp; - long l,l1,l2,u,i,k,e,s,n,p,av,av1; + long l, l1, l2, u, i, k, e, s, n, p; + gpmem_t av, av1; double alpha,beta,dk; if (gcmp0((GEN)x[2])) return ggamma((GEN)x[1],prec); @@ -1200,7 +1193,7 @@ double dnorm(double s, double t) { return s*s + t*t; } GEN -trans_fix_arg(long *prec, GEN *s0, GEN *sig, long *av, GEN *res) +trans_fix_arg(long *prec, GEN *s0, GEN *sig, gpmem_t *av, GEN *res) { GEN s, p1; long l; @@ -1232,10 +1225,11 @@ GEN gammanew(GEN s0, long la, long prec) { GEN s, u, a, y, res, tes, sig, invn2, p1, unr, nnx, pitemp; - long i, nn, lim, av, av2, avlim; + long i, lim, nn; + gpmem_t av, av2, avlim; int funeq = 0; - if (DEBUGLEVEL) timer2(); + if (DEBUGLEVEL) (void)timer2(); s = trans_fix_arg(&prec,&s0,&sig,&av,&res); if (signe(sig) <= 0 || expo(sig) < -1) @@ -1278,7 +1272,6 @@ gammanew(GEN s0, long la, long prec) nn = 1; if (DEBUGLEVEL) fprintferr("lim, nn: [%ld, %ld]\n",lim,nn); - if ((ulong)nn >= maxprime()) err(primer1); } prec++; unr = realun(prec); @@ -1338,15 +1331,53 @@ gammanew(GEN s0, long la, long prec) GEN ggamma(GEN x, long prec) { + gpmem_t av=avma, tetpil; + long m, ma; + GEN p1,p2,p3; + switch(typ(x)) { case t_INT: if (signe(x)<=0) err(gamer2); - return transc(ggamma,x,prec); + if (cmpis(x,481177) > 0) err(talker,"argument too large in ggamma"); +/* heuristic */ + if (cmpis(x,350 + 70*(prec-2)) > 0) return transc(ggamma,x,prec); + p2 = cgetr(prec); av = avma; + p1 = mpfact(itos(x) - 1); + affir(p1,p2); avma = av; + return p2; case t_REAL: return mpgamma(x); + case t_FRAC: + if (cmpis((GEN)x[2],2) == 0) + { + p2 = cgetr(prec); av = avma; + if (cmpis(mpabs((GEN)x[1]),962354) > 0) + err(talker,"argument too large in ggamma"); +/* heuristic */ + if (cmpis((GEN)x[1],200 + 50*(prec-2)) > 0) + return transc(ggamma,x,prec); + p1 = gsqrt(mppi(prec),prec); + m = itos((GEN)x[1]) - 1; ma = labs(m); + p3 = gmul2n(divii(mpfact(ma),mpfact(ma/2)),-ma); + if (m >= 0) affrr(gmul(p3,p1),p2); + else + { + affrr(gdiv(p1,p3),p2); + if ((m&3) == 2) setsigne(p2,-1); + } + avma = av; + return p2; + } + else return transc(ggamma,x,prec); + + case t_FRACN: + p1 = gred(x); + tetpil = avma; + return gerepile(av,tetpil,ggamma(p1,prec)); + case t_COMPLEX: return cxgamma(x,prec); @@ -1364,7 +1395,8 @@ ggamma(GEN x, long prec) void ggammaz(GEN x, GEN y) { - long av=avma, prec = precision(y); + long prec = precision(y); + gpmem_t av=avma; if (!prec) err(infprecer,"ggammaz"); gaffect(ggamma(x,prec),y); avma=av; @@ -1374,7 +1406,8 @@ static GEN mplngamma(GEN x) { GEN z,y,p1,p2,p3,p4,p5,p6,p7,p71,pitemp; - long l,l1,l2,u,i,k,e,f,s,sx,n,p,av,av1; + long l, l1, l2, u, i, k, e, f, s, sx, n, p; + gpmem_t av, av1; double alpha,beta,dk; sx=signe(x); if (!sx) err(talker,"zero argument in mplngamma"); @@ -1466,14 +1499,15 @@ mplngamma(GEN x) } /* t_REAL result */ y[3] = y[0]; y += 3; - affrr(p4,y); avma = (long)y; return y; + affrr(p4,y); avma = (gpmem_t)y; return y; } static GEN cxlngamma(GEN x, long prec) { GEN y,p1,p2,p3,p4,p5,p6,p7,p71,pitemp; - long l,l1,l2,flag,i,k,e,s,n,p,av,av1; + long l, l1, l2, flag, i, k, e, s, n, p; + gpmem_t av, av1; double alpha,beta,dk; if (gcmp0((GEN)x[2])) return glngamma((GEN)x[1],prec); @@ -1581,14 +1615,21 @@ cxlngamma(GEN x, long prec) GEN glngamma(GEN x, long prec) { - long i,av,n; - GEN y,p1; + long i, n; + gpmem_t av; + GEN y,p1,p2; switch(typ(x)) { case t_INT: if (signe(x)<=0) err(gamer2); - return transc(glngamma,x,prec); + p2 = cgetr(prec); av = avma; +/* heuristic */ + if (cmpis(x,200 + 50*(prec-2)) > 0) + return transc(glngamma,x,prec); + p1 = glog(mpfact(itos(x) - 1),prec); + affrr(p1,p2); avma = av; + return p2; case t_REAL: return mplngamma(x); @@ -1618,7 +1659,8 @@ glngamma(GEN x, long prec) void glngammaz(GEN x, GEN y) { - long av=avma, prec = precision(y); + long prec = precision(y); + gpmem_t av=avma; if (!prec) err(infprecer,"glngammaz"); gaffect(glngamma(x,prec),y); avma=av; @@ -1633,7 +1675,8 @@ glngammaz(GEN x, GEN y) static GEN mpgamd(long x, long prec) { - long i,j,a,l,av; + long i, j, a, l; + gpmem_t av; GEN y,p1,p2; a = labs(x); @@ -1660,7 +1703,7 @@ mpgamd(long x, long prec) void mpgamdz(long s, GEN y) { - long av=avma; + gpmem_t av=avma; affrr(mpgamd(s,lg(y)),y); avma=av; } @@ -1668,7 +1711,7 @@ mpgamdz(long s, GEN y) GEN ggamd(GEN x, long prec) { - long av,tetpil; + gpmem_t av, tetpil; switch(typ(x)) { @@ -1690,7 +1733,8 @@ ggamd(GEN x, long prec) void ggamdz(GEN x, GEN y) { - long av=avma, prec = precision(y); + long prec = precision(y); + gpmem_t av=avma; if (!prec) err(infprecer,"ggamdz"); gaffect(ggamd(x,prec),y); avma=av; @@ -1702,18 +1746,78 @@ ggamdz(GEN x, GEN y) /** **/ /********************************************************************/ -#if 1 +GEN +psinew(GEN s0, long prec) +{ + gpmem_t av; + GEN sum,z,a,res,tes,in2,sig,s,unr; + long lim,nn,k; + const long la = 3; + + + if (DEBUGLEVEL>2) (void)timer2(); + s = trans_fix_arg(&prec,&s0,&sig,&av,&res); + if (signe(sig) <= 0) + { + GEN pi = mppi(prec); + z = gsub(psinew(gsubsg(1,s), prec), gmul(pi, gcotan(gmul(pi,s), prec))); + gaffect(z, res); avma = av; return res; + } + + { + double ssig = rtodbl(sig); + double st = rtodbl(gimag(s)); + double l; + { + double rlog, ilog; /* log (s-gamma) */ + dcxlog(ssig - rtodbl(mpeuler(3)), st, &rlog,&ilog); + l = dnorm(rlog,ilog); + } + if (l < 0.000001) l = 0.000001; + l = log(l) / 2.; + lim = 2 + (long)ceil((pariC2*(prec-2) - l) / (2*(1+log((double)la)))); + if (lim < 2) lim = 2; + + l = (2*lim-1)*la / (2.*PI); + l = l*l - st*st; + if (l < 0.) l = 0.; + nn = (long)ceil( sqrt(l) - ssig ); + if (nn < 1) nn = 1; + if (DEBUGLEVEL>2) fprintferr("lim, nn: [%ld, %ld]\n",lim,nn); + } + prec++; unr = realun(prec); /* one extra word of precision */ + + a = gdiv(unr, gaddgs(s, nn)); /* 1 / (s+n) */ + sum = gmul2n(a,-1); + for (k = 0; k < nn; k++) + sum = gadd(sum, gdiv(unr, gaddgs(s, k))); + z = gsub(glog(gaddgs(s, nn), prec), sum); + if (DEBUGLEVEL>2) msgtimer("sum from 0 to N-1"); + + tes = divrs(bernreal(2*lim, prec), 2*lim); + in2 = gsqr(a); + for (k=2*lim-2; k>=2; k-=2) + { + tes = gadd(gmul(in2,tes), divrs(bernreal(k, prec), k)); + } + if (DEBUGLEVEL>2) msgtimer("Bernoulli sum"); + z = gsub(z, gmul(in2,tes)); + gaffect(z, res); avma = av; return res; +} + +#if 0 static GEN mppsi(GEN z) /* version p=2 */ { - long l,n,k,x,xx,av,av1,tetpil; + long l, n, k, x, xx; + gpmem_t av, av1, tetpil; GEN zk,u,v,a,b; av=avma; l=lg(z); x=(long)(1 + (bit_accuracy(l)>>1)*LOG2 + 1.58*rtodbl(absr(z))); if (expo(z)>=15 || x>46340) err(impl,"psi(x) for x>=29000"); xx=x*x; n=(long)(1+3.591*x); - affsr(x,a=cgetr(l)); + a = stor(x, l); a = mplog(a); gaffect(a,u=cgetr(l)); gaffsg(1,b=cgetr(l)); @@ -1727,7 +1831,6 @@ mppsi(GEN z) /* version p=2 */ } tetpil=avma; return gerepile(av,tetpil,divrr(u,v)); } -#else static GEN /* by Manfred Radimersky */ mppsi(GEN z) { @@ -1750,13 +1853,12 @@ mppsi(GEN z) gerepile(head, tail, subrr(s, a)); } - a = cgetr(len); - affsr(0, a); + a = stor(0, len); x = cgetr(len); affrr(z, x); tail = avma; while(cmprs(x, num >> 2) < 0) { - gaddz(a, divsr(1, x), a); + gaddz(a, ginv(x), a); gaddgsz(x, 1, x); avma = tail; } @@ -1764,14 +1866,14 @@ mppsi(GEN z) s = mplog(x); gsubz(s, a, s); b = gmul2n(x, 1); - gsubz(s, divsr(1, b), s); + gsubz(s, ginv(b), s); mpbern(num, len); affsr(-1, a); gdivgsz(a, 2, a); f = mulrr(x, x); - f = divsr(1, f); + f = ginv(f); k = 1; do { gmulz(a, f, a); @@ -1783,13 +1885,11 @@ mppsi(GEN z) return gerepilecopy(head, s); } -#endif - -#if 1 static GEN cxpsi(GEN z, long prec) /* version p=2 */ { - long l,n,k,x,xx,av,av1,tetpil; + long l, n, k, x, xx; + gpmem_t av, av1, tetpil; GEN zk,u,v,a,b,p1; if (gcmp0((GEN)z[2])) return gpsi((GEN)z[1],prec); @@ -1810,7 +1910,6 @@ cxpsi(GEN z, long prec) /* version p=2 */ } tetpil=avma; return gerepile(av,tetpil,gdiv(u,v)); } -#else GEN cxpsi(GEN z, long prec) /* by Manfred Radimersky */ { @@ -1877,12 +1976,9 @@ gpsi(GEN x, long prec) { switch(typ(x)) { - case t_REAL: - return mppsi(x); + case t_REAL: case t_COMPLEX: + return psinew(x,prec); - case t_COMPLEX: - return cxpsi(x,prec); - case t_INTMOD: case t_PADIC: err(typeer,"gpsi"); case t_SER: @@ -1894,7 +1990,8 @@ gpsi(GEN x, long prec) void gpsiz(GEN x, GEN y) { - long av=avma, prec = precision(y); + long prec = precision(y); + gpmem_t av=avma; if (!prec) err(infprecer,"gpsiz"); gaffect(gpsi(x,prec),y); avma=av;