=================================================================== RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/trans1.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/trans1.c 2001/10/02 11:17:05 1.1 +++ OpenXM_contrib/pari-2.2/src/basemath/Attic/trans1.c 2002/09/11 07:26:52 1.2 @@ -1,4 +1,4 @@ -/* $Id: trans1.c,v 1.1 2001/10/02 11:17:05 noro Exp $ +/* $Id: trans1.c,v 1.2 2002/09/11 07:26:52 noro Exp $ Copyright (C) 2000 The PARI group. @@ -41,7 +41,8 @@ void constpi(long prec) { GEN p1,p2,p3,tmppi; - long l,n,n1,av1,av2; + long l, n, n1; + gpmem_t av1, av2; double alpha; #define k1 545140134 @@ -56,10 +57,10 @@ constpi(long prec) prec++; - n=(long)(1+(prec-2)/alpha2); - n1=6*n-1; p1=cgetr(prec); - p2=addsi(k2,mulss(n,k1)); - affir(p2,p1); + n = (long)(1 + (prec-2)/alpha2); + n1 = 6*n - 1; + p2 = addsi(k2, mulss(n,k1)); + p1 = itor(p2, prec); /* initialisation longueur mantisse */ if (prec>=4) l=4; else l=prec; @@ -96,6 +97,27 @@ mppi(long prec) constpi(prec); affrr(gpi,x); return x; } +/* Pi * 2^n */ +GEN +Pi2n(long n, long prec) +{ + GEN x = mppi(prec); setexpo(x, 1+n); + return x; +} + +/* I * Pi * 2^n */ +GEN +PiI2n(long n, long prec) +{ + GEN z = cgetg(3,t_COMPLEX); + z[1] = zero; + z[2] = (long)Pi2n(n, prec); return z; +} + +/* 2I * Pi */ +GEN +PiI2(long prec) { return PiI2n(1, prec); } + /********************************************************************/ /** **/ /** FONCTION EULER **/ @@ -106,7 +128,8 @@ void consteuler(long prec) { GEN u,v,a,b,tmpeuler; - long l,n,k,x,av1,av2; + long l, n, k, x; + gpmem_t av1, av2; if (geuler && lg(geuler) >= prec) return; @@ -116,9 +139,9 @@ consteuler(long prec) prec++; l = prec+1; x = (long) (1 + (bit_accuracy(l) >> 2) * LOG2); - a=cgetr(l); affsr(x,a); u=mplog(a); setsigne(u,-1); affrr(u,a); - b=realun(l); - v=realun(l); + a = stor(x,l); u=mplog(a); setsigne(u,-1); affrr(u,a); + b = realun(l); + v = realun(l); n=(long)(1+3.591*x); /* z=3.591: z*[ ln(z)-1 ]=1 */ av2 = avma; if (x < SQRTVERYBIGINT) @@ -162,7 +185,8 @@ GEN transc(GEN (*f)(GEN,long), GEN x, long prec) { GEN p1,p2,y; - long lx,i,av,tetpil; + long lx, i; + gpmem_t av, tetpil; av=avma; switch(typ(x)) @@ -254,6 +278,23 @@ puiss0(GEN x) return y; } +static GEN +_sqr(void *data /* ignored */, GEN x) { + (void)data; return gsqr(x); +} +static GEN +_mul(void *data /* ignored */, GEN x, GEN y) { + (void)data; return gmul(x,y); +} +static GEN +_sqri(void *data /* ignored */, GEN x) { + (void)data; return sqri(x); +} +static GEN +_muli(void *data /* ignored */, GEN x, GEN y) { + (void)data; return mulii(x,y); +} + /* INTEGER POWERING (a^|n| for integer a and integer |n| > 1) * * Nonpositive powers and exponent one should be handled by gpow() and @@ -269,7 +310,7 @@ puiss0(GEN x) static GEN puissii(GEN a, GEN n, long s) { - long av,*p,m,k,i,lim; + gpmem_t av; GEN y; if (!signe(a)) return gzero; /* a==0 */ @@ -283,80 +324,65 @@ puissii(GEN a, GEN n, long s) if (n[2] == 1) { a = icopy(a); setsigne(a,s); return a; } if (n[2] == 2) return sqri(a); } - /* be paranoid about memory consumption */ - av=avma; lim=stack_lim(av,1); - y = a; p = n+2; m = *p; - /* normalize, i.e set highest bit to 1 (we know m != 0) */ - k = 1+bfffo(m); m<<=k; k = BITS_IN_LONG-k; - /* first bit is now implicit */ - for (i=lgefint(n)-2;;) - { - for (; k; m<<=1,k--) - { - y = sqri(y); - if (m < 0) y = mulii(y,a); /* first bit is set: multiply by base */ - if (low_stack(lim, stack_lim(av,1))) - { - if (DEBUGMEM>1) err(warnmem,"puissii"); - y = gerepileuptoint(av,y); - } - } - if (--i == 0) break; - m = *++p; k = BITS_IN_LONG; - } + av = avma; + y = leftright_pow(a, n, NULL, &_sqri, &_muli); setsigne(y,s); return gerepileuptoint(av,y); } +typedef struct { + long prec, a; + GEN (*sqr)(GEN); + GEN (*mulsg)(long,GEN); +} sr_muldata; + +static GEN +_rpowsi_mul(void *data, GEN x, GEN y/*unused*/) +{ + sr_muldata *D = (sr_muldata *)data; + return D->mulsg(D->a, x); +} + +static GEN +_rpowsi_sqr(void *data, GEN x) +{ + sr_muldata *D = (sr_muldata *)data; + if (typ(x) == t_INT && lgefint(x) >= D->prec) + { /* switch to t_REAL */ + D->sqr = &gsqr; + D->mulsg = &mulsr; + x = itor(x, D->prec); + } + return D->sqr(x); +} + /* return a^n as a t_REAL of precision prec. Adapted from puissii(). * Assume a > 0, n > 0 */ GEN rpowsi(ulong a, GEN n, long prec) { - long av,*p,m,k,i,lim; - GEN y, unr = realun(prec); - GEN (*sq)(GEN); - GEN (*mulsg)(long,GEN); + gpmem_t av; + GEN y; + sr_muldata D; - if (a == 1) return unr; - if (a == 2) { setexpo(unr, itos(n)); return unr; } - if (is_pm1(n)) { affsr(a, unr); return unr; } - /* be paranoid about memory consumption */ - av=avma; lim=stack_lim(av,1); - y = stoi(a); p = n+2; m = *p; - /* normalize, i.e set highest bit to 1 (we know m != 0) */ - k = 1+bfffo(m); m<<=k; k = BITS_IN_LONG-k; - /* first bit is now implicit */ - sq = &sqri; mulsg = &mulsi; - for (i=lgefint(n)-2;;) - { - for (; k; m<<=1,k--) - { - y = sq(y); - if (m < 0) y = mulsg(a,y); - if (lgefint(y) >= prec && typ(y) == t_INT) /* switch to t_REAL */ - { - sq = &gsqr; mulsg = &mulsr; - affir(y, unr); y = unr; - } - - if (low_stack(lim, stack_lim(av,1))) - { - if (DEBUGMEM>1) err(warnmem,"rpuisssi"); - y = gerepileuptoleaf(av,y); - } - } - if (--i == 0) break; - m = *++p, k = BITS_IN_LONG; - } - if (typ(y) == t_INT) { affir(y, unr); y = unr ; } - return gerepileuptoleaf(av,y); + if (a == 1) return realun(prec); + if (a == 2) return real2n(itos(n), prec); + if (is_pm1(n)) return stor(a, prec); + av = avma; + D.sqr = &sqri; + D.mulsg = &mulsi; + D.prec = prec; + D.a = a; + y = leftright_pow(stoi(a), n, (void*)&D, &_rpowsi_sqr, &_rpowsi_mul); + if (typ(y) == t_INT) y = itor(y, prec); + return gerepileuptoleaf(av, y); } GEN gpowgs(GEN x, long n) { - long lim,av,m,tx; - static long gn[3] = {evaltyp(t_INT)|m_evallg(3), 0, 0}; + long m, tx; + gpmem_t lim, av; + static long gn[3] = {evaltyp(t_INT)|_evallg(3), 0, 0}; GEN y; if (n == 0) return puiss0(x); @@ -440,7 +466,8 @@ gpowgs(GEN x, long n) GEN pow_monome(GEN x, GEN nn) { - long n,m,i,j,dd,tetpil, av = avma; + long n, m, i, j, dd; + gpmem_t tetpil, av = avma; GEN p1,y; if (is_bigint(nn)) err(talker,"power overflow in pow_monome"); n = itos(nn); m = labs(n); @@ -462,9 +489,8 @@ extern GEN powrealform(GEN x, GEN n); GEN powgi(GEN x, GEN n) { - long i,j,m,tx, sn=signe(n); - ulong lim,av; - GEN y, p1; + long tx, sn = signe(n); + GEN y; if (typ(n) != t_INT) err(talker,"not an integral exponent in powgi"); if (!sn) return puiss0(x); @@ -511,25 +537,26 @@ powgi(GEN x, GEN n) } case t_PADIC: { - GEN p, pe; + long e = itos(n)*valp(x), v; + GEN mod, p = (GEN)x[2]; + if (!signe(x[4])) { if (sn < 0) err(talker, "division by 0 p-adic in powgi"); - y=gcopy(x); setvalp(y, itos(n)*valp(x)); return y; + return padiczero(p, e); } y = cgetg(5,t_PADIC); - p = (GEN)x[2]; - pe= (GEN)x[3]; i = ggval(n, p); - if (i == 0) pe = icopy(pe); + mod = (GEN)x[3]; v = ggval(n, p); + if (v == 0) mod = icopy(mod); else { - pe = mulii(pe, gpowgs(p,i)); - pe = gerepileuptoint((long)y, pe); + mod = mulii(mod, gpowgs(p,v)); + mod = gerepileuptoint((gpmem_t)y, mod); } - y[1] = evalprecp(precp(x)+i) | evalvalp(itos(n) * valp(x)); + y[1] = evalprecp(precp(x)+v) | evalvalp(e); icopyifstack(p, y[2]); - y[3] = (long)pe; - y[4] = (long)powmodulo((GEN)x[4], n, pe); + y[3] = (long)mod; + y[4] = (long)powmodulo((GEN)x[4], n, mod); return y; } case t_QFR: @@ -537,31 +564,12 @@ powgi(GEN x, GEN n) case t_POL: if (ismonome(x)) return pow_monome(x,n); default: - av=avma; lim=stack_lim(av,1); - p1 = n+2; m = *p1; - y=x; j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j; - for (i=lgefint(n)-2;;) - { - for (; j; m<<=1,j--) - { - y=gsqr(y); - if (low_stack(lim, stack_lim(av,1))) - { - if(DEBUGMEM>1) err(warnmem,"[1]: powgi"); - y = gerepileupto(av, y); - } - if (m<0) y=gmul(y,x); - if (low_stack(lim, stack_lim(av,1))) - { - if(DEBUGMEM>1) err(warnmem,"[2]: powgi"); - y = gerepileupto(av, y); - } - } - if (--i == 0) break; - m = *++p1, j = BITS_IN_LONG; - } + { + gpmem_t av = avma; + y = leftright_pow(x, n, NULL, &_sqr, &_mul); if (sn < 0) y = ginv(y); return av==avma? gcopy(y): gerepileupto(av,y); + } } } @@ -569,42 +577,47 @@ powgi(GEN x, GEN n) static GEN ser_pui(GEN x, GEN n, long prec) { - long av,tetpil,lx,i,j; - GEN p1,p2,y,alp; + gpmem_t av, tetpil; + GEN y, p1, p2, lead; - if (gvar9(n) > varn(x)) + if (gvar9(n) <= varn(x)) return gexp(gmul(n, glog(x,prec)), prec); + lead = (GEN)x[2]; + if (gcmp1(lead)) { - GEN lead = (GEN)x[2]; - if (gcmp1(lead)) + GEN X, Y, alp; + long lx, mi, i, j, d; + + av = avma; alp = gclone(gadd(n,gun)); avma = av; + lx = lg(x); + y = cgetg(lx,t_SER); + y[1] = evalsigne(1) | evalvalp(0) | evalvarn(varn(x)); + X = x+2; + Y = y+2; + + d = mi = lx-3; while (mi>=1 && gcmp0((GEN)X[mi])) mi--; + Y[0] = un; + for (i=1; i<=d; i++) { - av=avma; alp = gclone(gadd(n,gun)); avma=av; - y=cgetg(lx=lg(x),t_SER); - y[1] = evalsigne(1) | evalvalp(0) | evalvarn(varn(x)); - y[2] = un; - for (i=3; i= (ulong)HIGHEXPOBIT) err(talker,"underflow or overflow in gpow"); - avma = av; y = cgetr(3); - y[1] = evalexpo(itos(x)); - y[2] = 0; return y; + avma = av; return realzero_bit(itos(x)); } if (tx==t_INTMOD && typ(n)==t_FRAC) { GEN p1; - if (!isprime((GEN)x[1])) err(talker,"modulus must be prime in gpow"); + if (!BSW_psp((GEN)x[1])) err(talker,"modulus must be prime in gpow"); y=cgetg(3,tx); copyifstack(x[1],y[1]); av=avma; p1=mpsqrtnmod((GEN)x[2],(GEN)n[2],(GEN)x[1],NULL); @@ -666,124 +683,181 @@ gpow(GEN x, GEN n, long prec) GEN mpsqrt(GEN x) { - long l,l0,l1,l2,s,eps,n,i,ex,av; + gpmem_t av, av0; + long l,l0,l1,l2,s,eps,n,i,ex; double beta; GEN y,p1,p2,p3; - if (typ(x)!=t_REAL) err(typeer,"mpsqrt"); - s=signe(x); if (s<0) err(talker,"negative argument in mpsqrt"); - if (!s) - { - y = cgetr(3); - y[1] = evalexpo(expo(x)>>1); - y[2] = 0; return y; - } - l=lg(x); y=cgetr(l); av=avma; + if (typ(x) != t_REAL) err(typeer,"mpsqrt"); + s = signe(x); if (s < 0) err(talker,"negative argument in mpsqrt"); + if (!s) return realzero_bit(expo(x) >> 1); + + l = lg(x); y = cgetr(l); av0 = avma; - p1=cgetr(l+1); affrr(x,p1); - ex=expo(x); eps = ex & 1; ex >>= 1; - setexpo(p1,eps); setlg(p1,3); + p1 = cgetr(l+1); affrr(x,p1); + ex = expo(x); eps = ex & 1; ex >>= 1; + setexpo(p1,eps); /* x = 2^(2 ex) p1 */ - n=(long)(2+log((double) (l-2))/LOG2); - p2=cgetr(l+1); p2[1]=evalexpo(0) | evalsigne(1); - beta=sqrt((eps+1) * (2 + p1[2]/C31)); - p2[2]=(long)((beta-2)*C31); - if (!p2[2]) { p2[2]=HIGHBIT; setexpo(p2,1); } - for (i=3; i<=l; i++) p2[i]=0; - setlg(p2,3); + n = (long)(2 + log((double) (l-2))/LOG2); + beta = sqrt((eps+1) * (2 + p1[2]/C31)); + p2 = cgetr(l+1); + p2[1] = evalexpo(0) | evalsigne(1); + p2[2] = (long)((beta-2)*C31); + if (!p2[2]) { p2[2] = (long)HIGHBIT; setexpo(p2,1); } + for (i=3; i<=l; i++) p2[i] = 0; - p3=cgetr(l+1); l-=2; l1=1; l2=3; + p3 = cgetr(l+1); l -= 2; l1 = 1; l2 = 3; + av = avma; for (i=1; i<=n; i++) - { + { /* execute p2 := (p2 + p1/p2)/2 */ l0 = l1<<1; - if (l0<=l) + if (l0 <= l) { l2 += l1; l1=l0; } else { l2 += l+1-l1; l1=l+1; } - setlg(p3,l1+2); setlg(p1,l1+2); setlg(p2,l2); - divrrz(p1,p2,p3); addrrz(p2,p3,p2); setexpo(p2,expo(p2)-1); + setlg(p1, l1+2); + setlg(p3, l1+2); + setlg(p2, l2); + affrr(divrr(p1,p2), p3); + affrr(addrr(p2,p3), p2); setexpo(p2, expo(p2)-1); + avma = av; } - affrr(p2,y); setexpo(y,expo(y)+ex); - avma=av; return y; + affrr(p2,y); setexpo(y, expo(y)+ex); + avma = av0; return y; } +/* O(p^e) */ +GEN +padiczero(GEN p, long e) +{ + GEN y = cgetg(5,t_PADIC); + y[4] = zero; + y[3] = un; + copyifstack(p,y[2]); + y[1] = evalvalp(e) | evalprecp(0); + return y; +} + +/* assume x unit, precp(x) = pp > 3 */ static GEN -padic_sqrt(GEN x) +sqrt_2adic(GEN x, long pp) { - long av = avma, av2,lim,e,r,lpp,lp,pp; - GEN y; + GEN z = mod16(x)==1? gun: stoi(3); + long zp; + gpmem_t av, lim; - e=valp(x); y=cgetg(5,t_PADIC); copyifstack(x[2],y[2]); - if (gcmp0(x)) + if (pp == 4) return z; + zp = 3; /* number of correct bits in z (compared to sqrt(x)) */ + + av = avma; lim = stack_lim(av,2); + for(;;) { - y[4] = zero; e = (e+1)>>1; - y[3] = un; - y[1] = evalvalp(e) | evalprecp(precp(x)); - return y; + GEN mod; + zp = (zp<<1) - 1; + if (zp > pp) zp = pp; + mod = shifti(gun, zp); + z = addii(z, resmod2n(mulii(x, mpinvmod(z,mod)), zp)); + z = shifti(z, -1); /* (z + x/z) / 2 */ + if (pp == zp) return z; + + if (zp < pp) zp--; + + if (low_stack(lim,stack_lim(av,2))) + { + if (DEBUGMEM > 1) err(warnmem,"padic_sqrt"); + z = gerepileuptoint(av,z); + } } +} + +/* x unit defined modulo modx = p^pp, p != 2, pp > 0 */ +static GEN +sqrt_padic(GEN x, GEN modx, long pp, GEN p) +{ + GEN mod, z = mpsqrtmod(x, p); + long zp = 1; + gpmem_t av, lim; + + if (!z) err(sqrter5); + if (pp <= zp) return z; + + av = avma; lim = stack_lim(av,2); + mod = p; + for(;;) + { /* could use the hensel_lift_accel idea. Not really worth it */ + GEN inv2; + zp <<= 1; + if (zp < pp) mod = sqri(mod); else { zp = pp; mod = modx; } + inv2 = shifti(mod, -1); /* = (mod + 1)/2 = 1/2 */ + z = addii(z, resii(mulii(x, mpinvmod(z,mod)), mod)); + z = mulii(z, inv2); + z = modii(z, mod); /* (z + x/z) / 2 */ + if (pp <= zp) return z; + + if (low_stack(lim,stack_lim(av,2))) + { + GEN *gptr[2]; gptr[0]=&z; gptr[1]=&mod; + if (DEBUGMEM>1) err(warnmem,"padic_sqrt"); + gerepilemany(av,gptr,2); + } + } +} + +static GEN +padic_sqrt(GEN x) +{ + long pp, e = valp(x); + GEN z,y,mod, p = (GEN)x[2]; + gpmem_t av; + + if (gcmp0(x)) return padiczero(p, (e+1) >> 1); if (e & 1) err(sqrter6); - e>>=1; setvalp(y,e); y[3] = y[2]; + + y = cgetg(5,t_PADIC); pp = precp(x); - if (egalii(gdeux, (GEN)x[2])) + mod = (GEN)x[3]; + x = (GEN)x[4]; /* lift to int */ + e >>= 1; av = avma; + if (egalii(gdeux, p)) { - lp=3; y[4] = un; r = mod8((GEN)x[4]); - if ((r!=1 && pp>=2) && (r!=5 || pp!=2)) err(sqrter5); - if (pp <= lp) { setprecp(y,1); return y; } - - x = dummycopy(x); setvalp(x,0); setvalp(y,0); - av2=avma; lim = stack_lim(av2,2); - y[3] = lstoi(8); - for(;;) + long r = mod8(x); + if (pp <= 3) { - lpp=lp; lp=(lp<<1)-1; - if (lp < pp) y[3] = lshifti((GEN)y[3], lpp-1); - else - { lp=pp; y[3] = x[3]; } - setprecp(y,lp); y = gdiv(gadd(y, gdiv(x,y)), gdeux); - if (lp < pp) lp--; - if (cmpii((GEN)y[4], (GEN)y[3]) >= 0) - y[4] = lsubii((GEN)y[4], (GEN)y[3]); - - if (pp <= lp) break; - if (low_stack(lim,stack_lim(av2,2))) - { - if (DEBUGMEM > 1) err(warnmem,"padic_sqrt"); - y = gerepileupto(av2,y); + switch(pp) { + case 1: break; + case 2: if ((r&3) == 1) break; + case 3: if (r == 1) break; + default: err(sqrter5); } + z = gun; + pp = 1; } - y = gcopy(y); + else + { + if (r != 1) err(sqrter5); + z = sqrt_2adic(x, pp); + z = gerepileuptoint(av, z); + pp--; + } + mod = shifti(gun, pp); } else /* p != 2 */ { - lp=1; y[4] = (long)mpsqrtmod((GEN)x[4],(GEN)x[2]); - if (!y[4]) err(sqrter5); - if (pp <= lp) { setprecp(y,1); return y; } - - x = dummycopy(x); setvalp(x,0); setvalp(y,0); - av2=avma; lim = stack_lim(av2,2); - for(;;) - { - lp<<=1; - if (lp < pp) y[3] = lsqri((GEN)y[3]); - else - { lp=pp; y[3] = x[3]; } - setprecp(y,lp); y = gdiv(gadd(y, gdiv(x,y)), gdeux); - - if (pp <= lp) break; - if (low_stack(lim,stack_lim(av2,2))) - { - if (DEBUGMEM > 1) err(warnmem,"padic_sqrt"); - y = gerepileupto(av2,y); - } - } + z = sqrt_padic(x, mod, pp, p); + z = gerepileuptoint(av, z); + mod = icopy(mod); } - setvalp(y,e); return gerepileupto(av,y); + y[1] = evalprecp(pp) | evalvalp(e); + copyifstack(p,y[2]); + y[3] = (long)mod; + y[4] = (long)z; return y; } GEN gsqrt(GEN x, long prec) { - long av,tetpil,e; + long e; + gpmem_t av, tetpil; GEN y,p1,p2; switch(typ(x)) @@ -846,17 +920,17 @@ gsqrt(GEN x, long prec) return padic_sqrt(x); case t_SER: - e=valp(x); + e = valp(x); if (gcmp0(x)) return zeroser(varn(x), (e+1)>>1); if (e & 1) err(sqrter6); av = avma; - /* trick: ser_pui assumes valp(x) = 0 */ - y = ser_pui(x,ghalf,prec); + y = dummycopy(x); setvalp(y, 0); + y = ser_pui(y, ghalf, prec); if (typ(y) == t_SER) /* generic case */ y[1] = evalsigne(1) | evalvalp(e>>1) | evalvarn(varn(x)); else /* e.g coeffs are POLMODs */ - y = gerepileupto(av, gmul(y, gpowgs(polx[varn(x)], e>>1))); - return y; + y = gmul(y, gpowgs(polx[varn(x)], e>>1)); + return gerepileupto(av, y); } return transc(gsqrt,x,prec); } @@ -864,7 +938,8 @@ gsqrt(GEN x, long prec) void gsqrtz(GEN x, GEN y) { - long av=avma, prec = precision(y); + long prec = precision(y); + gpmem_t av=avma; if (!prec) err(infprecer,"gsqrtz"); gaffect(gsqrt(x,prec),y); avma=av; @@ -877,26 +952,19 @@ gsqrtz(GEN x, GEN y) GEN rootsof1complex(GEN n, long prec) { - ulong ltop=avma; - GEN z,a; + gpmem_t ltop = avma; + GEN z; if (is_pm1(n)) return realun(prec); - if (lgefint(n)==3 && n[2]==2) - { - a=realun(prec); - setsigne(a,-1); - return a; - } - a=mppi(prec); setexpo(a,2); /* a=2*pi */ - a=divri(a,n); - z = cgetg(3,t_COMPLEX); - gsincos(a,(GEN*)(z+2),(GEN*)(z+1),prec); + if (lgefint(n)==3 && n[2]==2) return stor(-1, prec); + + z = exp_Ir( divri(Pi2n(1, prec), n) ); /* exp(2Ipi/n) */ return gerepileupto(ltop,z); } /*Only the O() of y is used*/ GEN rootsof1padic(GEN n, GEN y) { - ulong ltop=avma; + gpmem_t ltop=avma; GEN z,r; (void)mpsqrtnmod(gun,n,(GEN)y[2],&z); if (z==gzero){avma=ltop;return gzero;}/*should not happen*/ @@ -909,9 +977,10 @@ rootsof1padic(GEN n, GEN y) } static GEN paexp(GEN x); /*compute the p^e th root of x p-adic*/ -GEN padic_sqrtn_ram(GEN x, long e) +GEN +padic_sqrtn_ram(GEN x, long e) { - ulong ltop=avma; + gpmem_t ltop=avma; GEN n,a; GEN p=(GEN)x[2]; long v=0; @@ -938,9 +1007,10 @@ GEN padic_sqrtn_ram(GEN x, long e) return a; } /*compute the nth root of x p-adic p prime with n*/ -GEN padic_sqrtn_unram(GEN x, GEN n, GEN *zetan) +GEN +padic_sqrtn_unram(GEN x, GEN n, GEN *zetan) { - ulong ltop=avma,tetpil; + gpmem_t ltop=avma, tetpil; GEN a,r; GEN p=(GEN)x[2]; long v=0; @@ -978,22 +1048,19 @@ GEN padic_sqrtn_unram(GEN x, GEN n, GEN *zetan) else return gerepile(ltop,tetpil,r); } -GEN padic_sqrtn(GEN x, GEN n, GEN *zetan) + +GEN +padic_sqrtn(GEN x, GEN n, GEN *zetan) { - ulong ltop=avma,tetpil; + gpmem_t ltop=avma, tetpil; GEN p=(GEN)x[2]; GEN q; long e; GEN *gptr[2]; if (gcmp0(x)) { - GEN y; long m = itos(n); - e=valp(x); y=cgetg(5,t_PADIC); copyifstack(x[2],y[2]); - y[4] = zero; e = (e+m-1)/m; - y[3] = un; - y[1] = evalvalp(e) | evalprecp(precp(x)); - return y; + return padiczero(p, (valp(x)+m-1)/m); } /*First treat the ramified part using logarithms*/ e=pvaluation(n,p,&q); @@ -1041,7 +1108,8 @@ GEN padic_sqrtn(GEN x, GEN n, GEN *zetan) GEN gsqrtn(GEN x, GEN n, GEN *zetan, long prec) { - long av,tetpil,i,lx,tx; + long i, lx, tx; + gpmem_t av, tetpil; long e,m; GEN y,z; if (zetan) *zetan=gzero; @@ -1070,19 +1138,19 @@ gsqrtn(GEN x, GEN n, GEN *zetan, long prec) case t_SER: e=valp(x);m=itos(n); if (gcmp0(x)) return zeroser(varn(x), (e+m-1)/m); - if (e % m) err(talker,"incorrect valuation in gsqrt"); + if (e % m) err(talker,"incorrect valuation in gsqrtn"); av = avma; - /* trick: ser_pui assumes valp(x) = 0 */ - y = ser_pui(x,ginv(n),prec); + y = dummycopy(x); setvalp(y, 0); + y = ser_pui(y, ginv(n), prec); if (typ(y) == t_SER) /* generic case */ y[1] = evalsigne(1) | evalvalp(e/m) | evalvarn(varn(x)); else /* e.g coeffs are POLMODs */ - y = gerepileupto(av, gmul(y, gpowgs(polx[varn(x)], e/m))); - return y; + y = gmul(y, gpowgs(polx[varn(x)], e/m)); + return gerepileupto(av, y); case t_INTMOD: z=gzero; /*This is not great, but else it will generate too much trouble*/ - if (!isprime((GEN)x[1])) err(talker,"modulus must be prime in gsqrtn"); + if (!BSW_psp((GEN)x[1])) err(talker,"modulus must be prime in gsqrtn"); if (zetan) { z=cgetg(3,tx); copyifstack(x[1],z[1]); @@ -1106,6 +1174,14 @@ gsqrtn(GEN x, GEN n, GEN *zetan, long prec) if (tx==t_INT && is_pm1(x) && signe(x)>0) y=gun; /*speed-up since there is no way to call rootsof1complex directly from gp*/ + else if (gcmp0(x)) + { + if (gsigne(n) < 0) err(gdiver2); + if (isinexactreal(x)) + y = realzero_bit( itos( gfloor(gdivsg(gexpo(x), n)) ) ); + else + y = realzero(prec); + } else { av=avma; @@ -1135,16 +1211,14 @@ gsqrtn(GEN x, GEN n, GEN *zetan, long prec) GEN mpexp1(GEN x) { - long l,l1,l2,i,n,m,ex,s,av,av2, sx = signe(x); + long l, l1, l2, i, n, m, ex, s, sx = signe(x); + gpmem_t av, av2; double a,b,alpha,beta, gama = 2.0 /* optimized for SUN3 */; /* KB: 3.0 is better for UltraSparc */ GEN y,p1,p2,p3,p4,unr; if (typ(x)!=t_REAL) err(typeer,"mpexp1"); - if (!sx) - { - y=cgetr(3); y[1]=x[1]; y[2]=0; return y; - } + if (!sx) return realzero_bit(expo(x)); l=lg(x); y=cgetr(l); av=avma; /* room for result */ l2 = l+1; ex = expo(x); @@ -1205,20 +1279,22 @@ mpexp1(GEN x) GEN mpexp(GEN x) { - long av, sx = signe(x); + long sx = signe(x); + gpmem_t av; GEN y; if (!sx) return addsr(1,x); if (sx<0) setsigne(x,1); av = avma; y = addsr(1, mpexp1(x)); - if (sx<0) { y = divsr(1,y); setsigne(x,-1); } + if (sx<0) { y = ginv(y); setsigne(x,-1); } return gerepileupto(av,y); } static GEN paexp(GEN x) { - long k,av, e = valp(x), pp = precp(x), n = e + pp; + long k, e = valp(x), pp = precp(x), n = e + pp; + gpmem_t av; GEN y,r,p1; if (gcmp0(x)) return gaddgs(x,1); @@ -1240,57 +1316,68 @@ paexp(GEN x) return gerepileupto(av, y); } +static GEN +cxexp(GEN x, long prec) +{ + GEN r,p1,p2, y = cgetg(3,t_COMPLEX); + gpmem_t av = avma, tetpil; + r=gexp((GEN)x[1],prec); + gsincos((GEN)x[2],&p2,&p1,prec); + tetpil = avma; + y[1] = lmul(r,p1); + y[2] = lmul(r,p2); + gerepilemanyvec(av,tetpil,y+1,2); + return y; +} + +static GEN +serexp(GEN x, long prec) +{ + gpmem_t av; + long i,j,lx,ly,ex,mi; + GEN p1,y,xd,yd; + + ex = valp(x); + if (ex < 0) err(negexper,"gexp"); + if (gcmp0(x)) return gaddsg(1,x); + lx = lg(x); + if (ex) + { + ly = lx+ex; y = cgetg(ly,t_SER); + mi = lx-1; while (mi>=3 && gcmp0((GEN)x[mi])) mi--; + mi += ex-2; + y[1] = evalsigne(1) | evalvalp(0) | evalvarn(varn(x)); + /* zd[i] = coeff of X^i in z */ + xd = x+2-ex; yd = y+2; ly -= 2; + yd[0] = un; + for (i=1; i= prec) return glog2; + + tmplog2 = newbloc(prec); + *tmplog2 = evaltyp(t_REAL) | evallg(prec); + av0 = avma; + l = prec+1; G = bit_accuracy(l+1); + + s = S = divrs(realun(l), _3); + u = U = mpcopy(s); + av = avma; + for (k = 3; ; k += 2) + { + u = divrs(u, _9); + if (bit_accuracy(l) - expo(u) > G) { + l--; if (l <= 2) break; + setlg(U,l); + affrr(s,S); s = S; + affrr(u,U); u = U; avma = av; + } + s = addrr(s, divrs(u,k)); + } + setexpo(s, -1); affrr(s, tmplog2); + gunclone(glog2); glog2 = tmplog2; + avma = av0; return glog2; +} + GEN +mplog2(long prec) +{ + GEN x = cgetr(prec); + affrr(constlog2(prec), x); return x; +} + +GEN mplog(GEN x) { - long l,l1,l2,m,m1,n,i,ex,s,sgn,ltop,av; - double alpha,beta,a,b; - GEN y,p1,p2,p3,p4,p5,unr; + gpmem_t ltop, av; + long EX,l,l1,l2,m,n,k,ex,s; + double alpha,a,b; + GEN z,p1,y,y2,p4,p5,unr; if (typ(x)!=t_REAL) err(typeer,"mplog"); if (signe(x)<=0) err(talker,"non positive argument in mplog"); - l=lg(x); sgn = cmpsr(1,x); if (!sgn) return realzero(l); - y=cgetr(l); ltop=avma; - l2 = l+1; p2=p1=cgetr(l2); affrr(x,p1); - av=avma; - if (sgn > 0) p2 = divsr(1,p1); /* x<1 changer x en 1/x */ - for (m1=1; expo(p2)>0; m1++) p2 = mpsqrt(p2); - if (m1 > 1 || sgn > 0) { affrr(p2,p1); avma=av; } + av = avma; + l = lg(x); EX = expo(x); + z = cgetr(l); ltop = avma; + l2 = l+1; y=p1=cgetr(l2); affrr(x,p1); + setexpo(p1, 0); + if (gcmp1(p1)) { + if (!EX) { avma = av; return realzero(l); } + affrr(mulsr(EX, mplog2(l)), z); + avma = ltop; return z; + } + /* 1 < p1 < 2 */ + av = avma; l -= 2; alpha = 1+p1[2]/C31; if (!alpha) alpha = 0.00000001; - l -= 2; alpha = -log(alpha); - beta = (BITS_IN_LONG/2)*l*LOG2; - a = alpha/LOG2; - b = sqrt((BITS_IN_LONG/2)*l/3.0); - if (a<=b) + a = -log(alpha)/LOG2; + b = sqrt(BITS_IN_HALFULONG*l/3.0); + if (a <= b) { - n=(long)(1+sqrt((BITS_IN_LONG/2)*3.0*l)); - m=(long)(1+b-a); + n = 1 + (long)sqrt(BITS_IN_HALFULONG*3.0*l); + m = 1 + (long)(b-a); l2 += m>>TWOPOTBITS_IN_LONG; - p4=cgetr(l2); affrr(p1,p4); - p1 = p4; av=avma; - for (i=1; i<=m; i++) p1 = mpsqrt(p1); - affrr(p1,p4); avma=av; + p4 = cgetr(l2); affrr(p1,p4); + p1 = p4; av = avma; + for (k=1; k<=m; k++) p1 = mpsqrt(p1); + affrr(p1,p4); avma = av; } else { - n=(long)(1+beta/alpha); - m=0; p4 = p1; + n = 1 + (long)(BITS_IN_HALFULONG*l / a); + m = 0; p4 = p1; } - unr=realun(l2); - p2=cgetr(l2); p3=cgetr(l2); av=avma; + unr = realun(l2); + y = cgetr(l2); + y2 = cgetr(l2); av = avma; + /* want to compute log(X), X ~ 1 (X = p4) */ + /* y = (X-1)/(X+1). log(X) = log(1+y) - log(1-y) = 2 \sum_{k odd} y^k / * k */ + /* affrr needed here instead of setlg since prec may DECREASE */ p1 = cgetr(l2); affrr(subrr(p4,unr), p1); p5 = addrr(p4,unr); setlg(p5,l2); - affrr(divrr(p1,p5), p2); - affrr(mulrr(p2,p2), p3); - affrr(divrs(unr,2*n+1), p4); setlg(p4,4); avma=av; + affrr(divrr(p1,p5), y); /* = (X-1) / (X+1) ~ 0 */ + affrr(gsqr(y), y2); /* = y^2 */ + k = 2*n + 1; + affrr(divrs(unr,k), p4); setlg(p4,4); avma = av; - s=0; ex=expo(p3); l1 = 4; - for (i=n; i>=1; i--) - { - setlg(p3,l1); p5 = mulrr(p4,p3); - setlg(unr,l1); p1 = divrs(unr,2*i-1); + s = 0; ex = expo(y2); l1 = 4; + for (k -= 2; k > 0; k -= 2) + { /* compute sum_i=0^n y^2i / (2i + 1), k = 2i+1 */ + setlg(y2, l1); p5 = mulrr(p4,y2); + setlg(unr,l1); p1 = divrs(unr, k); s -= ex; l1 += s>>TWOPOTBITS_IN_LONG; if (l1>l2) l1=l2; - s %= BITS_IN_LONG; - setlg(p1,l1); setlg(p4,l1); setlg(p5,l1); - p1 = addrr(p1,p5); affrr(p1,p4); avma=av; + s &= (BITS_IN_LONG-1); + setlg(p1, l1); + setlg(p4, l1); + setlg(p5, l1); affrr(addrr(p1,p5), p4); avma=av; } - setlg(p4,l2); affrr(mulrr(p2,p4), y); - setexpo(y, expo(y)+m+m1); - if (sgn > 0) setsigne(y, -1); - avma=ltop; return y; + setlg(p4, l2); + y = mulrr(y,p4); /* = log(X)/2 */ + setexpo(y, expo(y)+m+1); + if (EX) y = addrr(y, mulsr(EX, mplog2(l2))); + affrr(y, z); avma = ltop; return z; } GEN teich(GEN x) { - GEN aux,y,z,p1; - long av,n,k; + GEN p,q,y,z,aux,p1; + long n, k; + gpmem_t av; if (typ(x)!=t_PADIC) err(talker,"not a p-adic argument in teichmuller"); if (!signe(x[4])) return gcopy(x); - y=cgetp(x); z=(GEN)x[4]; - if (!cmpis((GEN)x[2],2)) + p = (GEN)x[2]; + q = (GEN)x[3]; + z = (GEN)x[4]; + y = cgetp(x); av = avma; + if (egalii(p, gdeux)) + z = (mod4(z) & 2)? addsi(-1,q): gun; + else { - n = mod4(z); - if (n & 2) - addsiz(-1,(GEN)x[3],(GEN)y[4]); - else - affsi(1,(GEN)y[4]); - return y; + p1 = addsi(-1, p); + z = resii(z, p); + aux = divii(addsi(-1,q),p1); n = precp(x); + for (k=1; k0; pp++) - p1 = mulii(p1,(GEN)x[2]); - avma=av; pp-=2; + for (p1=stoi(e); cmpsi(pp,p1)>0; pp++) p1 = mulii(p1, p); + pp -= 2; } - k=pp/e; if (!odd(k)) k--; - y2=gsqr(y); s=gdivgs(gun,k); - while (k>=3) + k = pp/e; if (!odd(k)) k--; + y2 = gsqr(y); s = gdivgs(gun,k); + while (k > 2) { - k-=2; s=gadd(gmul(y2,s),gdivgs(gun,k)); + k -= 2; s = gadd(gmul(y2,s), gdivgs(gun,k)); } - tetpil=avma; return gerepile(av1,tetpil,gmul(s,y)); + return gmul(s,y); } GEN palog(GEN x) { - long av=avma,tetpil; - GEN p1,y; + gpmem_t av = avma; + GEN y, p = (GEN)x[2]; if (!signe(x[4])) err(talker,"zero argument in palog"); - if (cmpis((GEN)x[2],2)) + if (egalii(p, gdeux)) { - y=cgetp(x); p1=gsubgs((GEN)x[2],1); - affii(powmodulo((GEN)x[4],p1,(GEN)x[3]),(GEN)y[4]); - y=gmulgs(palogaux(y),2); tetpil=avma; - return gerepile(av,tetpil,gdiv(y,p1)); + y = gsqr(x); setvalp(y,0); + y = palogaux(y); } - y=gsqr(x); setvalp(y,0); tetpil=avma; - return gerepile(av,tetpil,palogaux(y)); + else + { /* compute log(x^(p-1)) / (p-1) */ + GEN mod = (GEN)x[3], p1 = subis(p,1); + y = cgetp(x); + y[4] = (long)powmodulo((GEN)x[4], p1, mod); + p1 = divii(subis(mod,1), p1); /* 1/(1-p) */ + y = gmul(palogaux(y), mulis(p1, -2)); + } + return gerepileupto(av,y); } GEN @@ -1471,7 +1614,7 @@ log0(GEN x, long flag,long prec) GEN glog(GEN x, long prec) { - long av,tetpil; + gpmem_t av, tetpil; GEN y,p1,p2; switch(typ(x)) @@ -1492,7 +1635,8 @@ glog(GEN x, long prec) return palog(x); case t_SER: - av=avma; if (valp(x)) err(negexper,"glog"); + av=avma; + if (valp(x) || gcmp0(x)) err(talker,"log is not analytic at 0"); p1=gdiv(derivser(x),x); tetpil=avma; p1=integ(p1,varn(x)); if (gcmp1((GEN)x[2])) return gerepile(av,tetpil,p1); @@ -1509,7 +1653,8 @@ glog(GEN x, long prec) void glogz(GEN x, GEN y) { - long av=avma, prec = precision(y); + long prec = precision(y); + gpmem_t av=avma; if (!prec) err(infprecer,"glogz"); gaffect(glog(x,prec),y); avma=av; @@ -1517,136 +1662,128 @@ glogz(GEN x, GEN y) /********************************************************************/ /** **/ -/** FONCTION SINCOS-1 **/ +/** SINE, COSINE **/ /** **/ /********************************************************************/ +/* Reduce x0 mod Pi/2 to x in [-Pi/4, Pi/4]. Return cos(x)-1 */ static GEN -mpsc1(GEN x, long *ptmod8) +mpsc1(GEN x0, long *ptmod8) { - const long mmax = 23169; + const long mmax = 23169; /* largest m such that (2m+2)(2m+1) < 2^31 */ /* on a 64-bit machine with true 128 bit/64 bit division, one could - * take mmax=1518500248; on the alpha it does not seem worthwhile - */ - long l,l0,l1,l2,l4,ee,i,n,m,s,t,av; + * take mmax=1518500248; on the alpha it does not seem worthwhile */ + long l, l0, l1, l2, l4, i, n, m, s, t; + gpmem_t av; double alpha,beta,a,b,c,d; - GEN y,p1,p2,p3,p4,pitemp; + GEN y,p1,p2,p3,x2,pitemp, x = x0; - if (typ(x)!=t_REAL) err(typeer,"mpsc1"); - if (!signe(x)) - { - y=cgetr(3); - y[1] = evalexpo((expo(x)<<1) - 1); - y[2] = 0; *ptmod8=0; return y; - } - l=lg(x); y=cgetr(l); av=avma; + if (typ(x) != t_REAL) err(typeer,"mpsc1"); + l = lg(x); y = cgetr(l); av = avma; - l++; pitemp = mppi(l+1); setexpo(pitemp,-1); - p1 = addrr(x,pitemp); setexpo(pitemp,0); - if (expo(p1) >= bit_accuracy(min(l,lg(p1))) + 3) - err(precer,"mpsc1"); - p3 = divrr(p1,pitemp); p2 = mpent(p3); - if (signe(p2)) x = subrr(x, mulir(p2,pitemp)); - p1 = cgetr(l+1); affrr(x, p1); + l++; pitemp = mppi(l); + setexpo(pitemp,-1); + p1 = addrr(x,pitemp); /* = x + Pi/4 */ + l0 = min(l, lg(p1)); + if (expo(p1) >= bit_accuracy(l0) + 3) err(precer,"mpsc1"); - *ptmod8 = (signe(p1) < 0)? 4: 0; - if (signe(p2)) + setexpo(pitemp, 0); + p1 = mpent( divrr(p1,pitemp) ); /* round ( x / (Pi/2) ) */ + if (signe(p1)) x = subrr(x, mulir(p1,pitemp)); /* x mod Pi/2 */ + p2 = cgetr(l); affrr(x, p2); x = p2; + + *ptmod8 = (signe(x) < 0)? 4: 0; + if (signe(p1)) { - long mod4 = mod4(p2); - if (signe(p2) < 0 && mod4) mod4 = 4-mod4; - *ptmod8 += mod4; + long k = mod4(p1); + if (signe(p1) < 0 && k) k = 4-k; + /* x = x0 - k*Pi/2 (mod 2Pi) */ + *ptmod8 += k; } - if (gcmp0(p1)) alpha=1000000.0; + + if (gcmp0(x)) alpha = 1000000.0; else { - m=expo(p1); alpha=(m<-1022)? -1-m*LOG2: -1-log(fabs(rtodbl(p1))); + long e = expo(x); + alpha = -1 - ((e<-1022)? e*LOG2: log(fabs(rtodbl(x)))); } beta = 5 + bit_accuracy(l)*LOG2; - a=0.5/LOG2; b=0.5*a; - c = a+sqrt((beta+b)/LOG2); - d = ((beta/c)-alpha-log(c))/LOG2; - if (d>=0) + a = 0.5 / LOG2; + b = 0.5 * a; + c = a + sqrt((beta+b) / LOG2); + d = ((beta/c) - alpha - log(c)) / LOG2; + if (d >= 0) { - m=(long)(1+d); - n=(long)((1+c)/2.0); - setexpo(p1,expo(p1)-m); + m = (long)(1+d); + n = (long)((1+c) / 2.0); + setexpo(x, expo(x)-m); } - else { m=0; n=(long)((1+beta/alpha)/2.0); } + else { m = 0; n = (long)((2+beta/alpha) / 2.0); } l2=l+1+(m>>TWOPOTBITS_IN_LONG); - p2=realun(l2); setlg(p2,4); - p4=cgetr(l2); av = avma; - affrr(gsqr(p1),p4); setlg(p4,4); - - if (n>mmax) - p3 = divrs(divrs(p4,2*n+2),2*n+1); + p2 = realun(l2); + x2 = cgetr(l2); av = avma; + affrr(gsqr(x), x2); + + setlg(x2, 3); + if (n > mmax) + p3 = divrs(divrs(x2, 2*n+2), 2*n+1); else - p3 = divrs(p4, (2*n+2)*(2*n+1)); - ee = -expo(p3); s=0; - l4 = l1 = 3 + (ee>>TWOPOTBITS_IN_LONG); - if (l4<=l2) { setlg(p2,l4); setlg(p4,l4); } - for (i=n; i>mmax; i--) + p3 = divrs(x2, (2*n+2)*(2*n+1)); + s = -expo(p3); + l4 = l1 = 3 + (s>>TWOPOTBITS_IN_LONG); + if (l4<=l2) { setlg(p2,l4); setlg(x2,l4); } + s = 0; + for (i=n; i>=2; i--) { - p3 = divrs(divrs(p4,2*i),2*i-1); s -= expo(p3); - t=s&(BITS_IN_LONG-1); l0=(s>>TWOPOTBITS_IN_LONG); + if (i > mmax) + p3 = divrs(divrs(x2, 2*i), 2*i-1); + else + p3 = divrs(x2, 2*i*(2*i-1)); + s -= expo(p3); + t = s & (BITS_IN_LONG-1); l0 = (s>>TWOPOTBITS_IN_LONG); if (t) l0++; l1 += l0; if (l1>l2) { l0 += l2-l1; l1=l2; } l4 += l0; p3 = mulrr(p3,p2); - if (l4<=l2) { setlg(p2,l4); setlg(p4,l4); } - subsrz(1,p3,p2); avma=av; + if (l4<=l2) { setlg(p2,l4); setlg(x2,l4); } + subsrz(1,p3, p2); avma = av; } - for ( ; i>=2; i--) - { - p3 = divrs(p4, 2*i*(2*i-1)); s -= expo(p3); - t=s&(BITS_IN_LONG-1); l0=(s>>TWOPOTBITS_IN_LONG); - if (t) l0++; - l1 += l0; if (l1>l2) { l0 += l2-l1; l1=l2; } - l4 += l0; - p3 = mulrr(p3,p2); - if (l4<=l2) { setlg(p2,l4); setlg(p4,l4); } - subsrz(1,p3,p2); avma=av; - } - if (l4<=l2) { setlg(p2,l4); setlg(p4,l4); } - setexpo(p4,expo(p4)-1); setsigne(p4, -signe(p4)); - p2 = mulrr(p4,p2); + setlg(p2,l2); setlg(x2,l2); + setexpo(p2, expo(p2)-1); setsigne(p2, -signe(p2)); + p2 = mulrr(x2,p2); + /* Now p2 = sum {1<= i <=n} (-1)^i x^(2i) / (2i)! ~ cos(x) - 1 */ for (i=1; i<=m; i++) - { - p2 = mulrr(p2,addsr(2,p2)); setexpo(p2,expo(p2)+1); + { /* p2 = cos(x)-1 --> cos(2x)-1 */ + p2 = mulrr(p2, addsr(2,p2)); setexpo(p2, expo(p2)+1); } affrr(p2,y); avma=av; return y; } -/********************************************************************/ -/** **/ -/** FONCTION sqrt(-x*(x+2)) **/ -/** **/ -/********************************************************************/ - +/* sqrt (1 - (1+x)^2) = sqrt(-x*(x+2)). Sends cos(x)-1 to |sin(x)| */ static GEN mpaut(GEN x) { - long av = avma; - GEN p1 = mulrr(x,addsr(2,x)); - setsigne(p1,-signe(p1)); + gpmem_t av = avma; + GEN p1 = mulrr(x, addsr(2,x)); + setsigne(p1, -signe(p1)); return gerepileuptoleaf(av, mpsqrt(p1)); } /********************************************************************/ -/** **/ -/** FONCTION COSINUS **/ -/** **/ +/** COSINE **/ /********************************************************************/ static GEN mpcos(GEN x) { - long mod8,av,tetpil; + long mod8; + gpmem_t av; GEN y,p1; if (typ(x)!=t_REAL) err(typeer,"mpcos"); if (!signe(x)) return addsr(1,x); - av=avma; p1=mpsc1(x,&mod8); tetpil=avma; + av = avma; p1 = mpsc1(x,&mod8); switch(mod8) { case 0: case 4: @@ -1658,13 +1795,13 @@ mpcos(GEN x) default: /* case 3: case 5: */ y=mpaut(p1); break; } - return gerepile(av,tetpil,y); + return gerepileuptoleaf(av, y); } GEN gcos(GEN x, long prec) { - long av,tetpil; + gpmem_t av, tetpil; GEN r,u,v,y,p1,p2; switch(typ(x)) @@ -1699,32 +1836,28 @@ gcos(GEN x, long prec) void gcosz(GEN x, GEN y) { - long av = avma, prec = precision(y); + long prec = precision(y); + gpmem_t av = avma; if (!prec) err(infprecer,"gcosz"); gaffect(gcos(x,prec),y); avma=av; } /********************************************************************/ -/** **/ -/** FONCTION SINUS **/ -/** **/ +/** SINE **/ /********************************************************************/ GEN mpsin(GEN x) { - long mod8,av,tetpil; + long mod8; + gpmem_t av; GEN y,p1; if (typ(x)!=t_REAL) err(typeer,"mpsin"); - if (!signe(x)) - { - y=cgetr(3); y[1]=x[1]; y[2]=0; - return y; - } + if (!signe(x)) return realzero_bit(expo(x)); - av=avma; p1=mpsc1(x,&mod8); tetpil=avma; + av = avma; p1 = mpsc1(x,&mod8); switch(mod8) { case 0: case 6: @@ -1736,13 +1869,13 @@ mpsin(GEN x) default: /* case 3: case 7: */ y=subsr(-1,p1); break; } - return gerepile(av,tetpil,y); + return gerepileuptoleaf(av, y); } GEN gsin(GEN x, long prec) { - long av,tetpil; + gpmem_t av, tetpil; GEN r,u,v,y,p1,p2; switch(typ(x)) @@ -1757,7 +1890,8 @@ gsin(GEN x, long prec) p1=gsub(p2,p1); gsincos((GEN)x[1],&u,&v,prec); tetpil=avma; - y[1]=lmul(p2,u); y[2]=lmul(p1,v); + y[1]=lmul(p2,u); + y[2]=lmul(p1,v); gerepilemanyvec(av,tetpil,y+1,2); return y; @@ -1777,29 +1911,29 @@ gsin(GEN x, long prec) void gsinz(GEN x, GEN y) { - long av=avma, prec = precision(y); + long prec = precision(y); + gpmem_t av=avma; if (!prec) err(infprecer,"gsinz"); gaffect(gsin(x,prec),y); avma=av; } /********************************************************************/ -/** **/ -/** PROCEDURE SINUS,COSINUS **/ -/** **/ +/** SINE, COSINE together **/ /********************************************************************/ void mpsincos(GEN x, GEN *s, GEN *c) { - long av,tetpil,mod8; + long mod8; + gpmem_t av, tetpil; GEN p1, *gptr[2]; if (typ(x)!=t_REAL) err(typeer,"mpsincos"); if (!signe(x)) { - p1=cgetr(3); *s=p1; p1[1]=x[1]; - p1[2]=0; *c=addsr(1,x); return; + *s = realzero_bit(expo(x)); + *c = addsr(1,x); return; } av=avma; p1=mpsc1(x,&mod8); tetpil=avma; @@ -1818,10 +1952,22 @@ mpsincos(GEN x, GEN *s, GEN *c) gerepilemanysp(av,tetpil,gptr,2); } +/* return exp(ix), x a t_REAL */ +GEN +exp_Ir(GEN x) +{ + gpmem_t av = avma; + GEN v = cgetg(3,t_COMPLEX); + mpsincos(x, (GEN*)(v+2), (GEN*)(v+1)); + if (!signe(x)) return gerepilecopy(av, (GEN)v[1]); + return v; +} + void gsincos(GEN x, GEN *s, GEN *c, long prec) { - long av,tetpil,ii,i,j,ex,ex2,lx,ly; + long ii, i, j, ex, ex2, lx, ly, mi; + gpmem_t av, tetpil; GEN r,u,v,u1,v1,p1,p2,p3,p4,ps,pc; GEN *gptr[4]; @@ -1848,7 +1994,8 @@ gsincos(GEN x, GEN *s, GEN *c, long prec) p3=gmul(v1,v); p4=gmul(r,u); gptr[0]=&p1; gptr[1]=&p2; gptr[2]=&p3; gptr[3]=&p4; gerepilemanysp(av,tetpil,gptr,4); - ps[1]=(long)p1; ps[2]=(long)p2; pc[1]=(long)p3; pc[2]=(long)p4; + ps[1]=(long)p1; pc[1]=(long)p3; + ps[2]=(long)p2; pc[2]=(long)p4; return; case t_SER: @@ -1876,6 +2023,8 @@ gsincos(GEN x, GEN *s, GEN *c, long prec) } ly=lx+2*ex; + mi = lx-1; while (mi>=3 && gcmp0((GEN)x[mi])) mi--; + mi += ex-2; pc=cgetg(ly,t_SER); *c=pc; ps=cgetg(lx,t_SER); *s=ps; pc[1] = evalsigne(1) | evalvalp(0) | evalvarn(varn(x)); @@ -1885,14 +2034,14 @@ gsincos(GEN x, GEN *s, GEN *c, long prec) for (i=ex2; i