OpenXM_contrib/pari-2.2/src/basemath/trans1.c - diff

Return to trans1.c CVS log

Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath

Diff for /OpenXM_contrib/pari-2.2/src/basemath/Attic/trans1.c between version 1.1 and 1.2

version 1.1, 2001/10/02 11:17:05

version 1.2, 2002/09/11 07:26:52

Line 41 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

Line 56 constpi(long prec)

Line 57 constpi(long prec)

prec++;

n=(long)(1+(prec-2)/alpha2);

n = (long)(1 + (prec-2)/alpha2);

n1=6*n-1; p1=cgetr(prec);

n1 = 6*n - 1;

p2=addsi(k2,mulss(n,k1));

p2 = addsi(k2, mulss(n,k1));

affir(p2,p1);

p1 = itor(p2, prec);

/* initialisation longueur mantisse */

if (prec>=4) l=4; else l=prec;

Line 96 mppi(long prec)

Line 97 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 **/

Line 106 void

Line 128 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;

Line 116 consteuler(long prec)

Line 139 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);

a = stor(x,l); u=mplog(a); setsigne(u,-1); affrr(u,a);

b=realun(l);

b = realun(l);

v=realun(l);

v = realun(l);

n=(long)(1+3.591*x); /* z=3.591: z*[ ln(z)-1 ]=1 */

av2 = avma;

if (x < SQRTVERYBIGINT)

Line 162 GEN

Line 185 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))

Line 254 puiss0(GEN x)

Line 278 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

Line 269 puiss0(GEN x)

Line 310 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 */

Line 283 puissii(GEN a, GEN n, long s)

Line 324 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;

av=avma; lim=stack_lim(av,1);

y = leftright_pow(a, n, NULL, &_sqri, &_muli);

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;

}

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;

gpmem_t av;

GEN y, unr = realun(prec);

GEN y;

GEN (*sq)(GEN);

sr_muldata D;

GEN (*mulsg)(long,GEN);

if (a == 1) return unr;

if (a == 1) return realun(prec);

if (a == 2) { setexpo(unr, itos(n)); return unr; }

if (a == 2) return real2n(itos(n), prec);

if (is_pm1(n)) { affsr(a, unr); return unr; }

if (is_pm1(n)) return stor(a, prec);

/* be paranoid about memory consumption */

av = avma;

av=avma; lim=stack_lim(av,1);

D.sqr = &sqri;

y = stoi(a); p = n+2; m = *p;

D.mulsg = &mulsi;

/* normalize, i.e set highest bit to 1 (we know m != 0) */

D.prec = prec;

k = 1+bfffo(m); m<<=k; k = BITS_IN_LONG-k;

D.a = a;

/* first bit is now implicit */

y = leftright_pow(stoi(a), n, (void*)&D, &_rpowsi_sqr, &_rpowsi_mul);

sq = &sqri; mulsg = &mulsi;

if (typ(y) == t_INT) y = itor(y, prec);

for (i=lgefint(n)-2;;)

return gerepileuptoleaf(av, y);

{

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);

}

GEN

gpowgs(GEN x, long n)

{

long lim,av,m,tx;

long m, tx;

static long gn[3] = {evaltyp(t_INT)|m_evallg(3), 0, 0};

gpmem_t lim, av;

static long gn[3] = {evaltyp(t_INT)|_evallg(3), 0, 0};

GEN y;

if (n == 0) return puiss0(x);

Line 440 gpowgs(GEN x, long n)

Line 466 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);

Line 462 extern GEN powrealform(GEN x, GEN n);

Line 489 extern GEN powrealform(GEN x, GEN n);

GEN

powgi(GEN x, GEN n)

{

long i,j,m,tx, sn=signe(n);

long tx, sn = signe(n);

ulong lim,av;

GEN y;

GEN y, p1;

if (typ(n) != t_INT) err(talker,"not an integral exponent in powgi");

if (!sn) return puiss0(x);

Line 511 powgi(GEN x, GEN n)

Line 537 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];

mod = (GEN)x[3]; v = ggval(n, p);

pe= (GEN)x[3]; i = ggval(n, p);

if (v == 0) mod = icopy(mod);

if (i == 0) pe = icopy(pe);

else

{

pe = mulii(pe, gpowgs(p,i));

mod = mulii(mod, gpowgs(p,v));

pe = gerepileuptoint((long)y, pe);

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[3] = (long)mod;

y[4] = (long)powmodulo((GEN)x[4], n, pe);

y[4] = (long)powmodulo((GEN)x[4], n, mod);

return y;

}

case t_QFR:

Line 537 powgi(GEN x, GEN n)

Line 564 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;

gpmem_t av = avma;

y=x; j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;

y = leftright_pow(x, n, NULL, &_sqr, &_mul);

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;

}

if (sn < 0) y = ginv(y);

return av==avma? gcopy(y): gerepileupto(av,y);

}

Line 569 powgi(GEN x, GEN n)

Line 577 powgi(GEN x, GEN n)

static GEN

ser_pui(GEN x, GEN n, long prec)

{

long av,tetpil,lx,i,j;

gpmem_t av, tetpil;

GEN p1,p2,y,alp;

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];

GEN X, Y, alp;

if (gcmp1(lead))

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;

av = avma; p1 = gzero;

y=cgetg(lx=lg(x),t_SER);

for (j=1; j<=min(i,mi); j++)

y[1] = evalsigne(1) | evalvalp(0) | evalvarn(varn(x));

y[2] = un;

for (i=3; i<lx; i++)

{

av=avma; p1=gzero;

p2 = gsubgs(gmulgs(alp,j),i);

for (j=1; j<i-1; j++)

p1 = gadd(p1, gmul(gmul(p2,(GEN)X[j]),(GEN)Y[i-j]));

{

p2 = gsubgs(gmulgs(alp,j),i-2);

p1 = gadd(p1,gmul(gmul(p2,(GEN)x[j+2]),(GEN)y[i-j]));

}

tetpil=avma; y[i]=lpile(av,tetpil,gdivgs(p1,i-2));

}

gunclone(alp); return y;

tetpil = avma; Y[i] = lpile(av,tetpil, gdivgs(p1,i));

}

av=avma; p1=gdiv(x,lead); p1[2] = un; /* in case it's inexact */

gunclone(alp); return y;

p1=gpow(p1,n,prec); p2=gpow(lead,n,prec);

tetpil=avma; return gerepile(av,tetpil,gmul(p1,p2));

}

av=avma; y=gmul(n,glog(x,prec)); tetpil=avma;

p1 = gdiv(x,lead); p1[2] = un; /* in case it's inexact */

return gerepile(av,tetpil,gexp(y,prec));

p1 = gpow(p1, n, prec);

p2 = gpow(lead,n, prec); return gmul(p1, p2);

}

GEN

gpow(GEN x, GEN n, long prec)

{

long av,tetpil,i,lx,tx;

long i, lx, tx;

gpmem_t av, tetpil;

GEN y;

if (typ(n)==t_INT) return powgi(x,n);

Line 615 gpow(GEN x, GEN n, long prec)

Line 628 gpow(GEN x, GEN n, long prec)

for (i=1; i<lx; i++) y[i]=lpui((GEN)x[i],n,prec);

return y;

}

if (tx==t_POL)

{

av=avma; y=tayl(x,gvar(x),precdl); tetpil=avma;

return gerepile(av,tetpil,gpow(y,n,prec));

}

if (tx==t_SER)

{

if (valp(x))

err(talker,"not an integer exponent for non invertible series in gpow");

if (lg(x) == 2) return gcopy(x); /* O(1) */

return ser_pui(x,n,prec);

av = avma;

return gerepileupto(av, ser_pui(x, n, prec));

}

av=avma;

if (gcmp0(x))

Line 636 gpow(GEN x, GEN n, long prec)

Line 655 gpow(GEN x, GEN n, long prec)

x = ground(gmulsg(gexpo(x),n));

if (is_bigint(x) || (ulong)x[2] >= (ulong)HIGHEXPOBIT)

err(talker,"underflow or overflow in gpow");

avma = av; y = cgetr(3);

avma = av; return realzero_bit(itos(x));

y[1] = evalexpo(itos(x));

y[2] = 0; return y;

}

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);

Line 666 gpow(GEN x, GEN n, long prec)

Line 683 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");

if (typ(x) != t_REAL) err(typeer,"mpsqrt");

s=signe(x); if (s<0) err(talker,"negative argument in mpsqrt");

s = signe(x); if (s < 0) err(talker,"negative argument in mpsqrt");

if (!s)

if (!s) return realzero_bit(expo(x) >> 1);

{

y = cgetr(3);

l = lg(x); y = cgetr(l); av0 = avma;

y[1] = evalexpo(expo(x)>>1);

y[2] = 0; return y;

}

l=lg(x); y=cgetr(l); av=avma;

p1=cgetr(l+1); affrr(x,p1);

p1 = cgetr(l+1); affrr(x,p1);

ex=expo(x); eps = ex & 1; ex >>= 1;

ex = expo(x); eps = ex & 1; ex >>= 1;

setexpo(p1,eps); setlg(p1,3);

setexpo(p1,eps); /* x = 2^(2 ex) p1 */

n=(long)(2+log((double) (l-2))/LOG2);

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));

beta=sqrt((eps+1) * (2 + p1[2]/C31));

p2 = cgetr(l+1);

p2[2]=(long)((beta-2)*C31);

p2[1] = evalexpo(0) | evalsigne(1);

if (!p2[2]) { p2[2]=HIGHBIT; setexpo(p2,1); }

p2[2] = (long)((beta-2)*C31);

for (i=3; i<=l; i++) p2[i]=0;

if (!p2[2]) { p2[2] = (long)HIGHBIT; setexpo(p2,1); }

setlg(p2,3);

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);

setlg(p1, l1+2);

divrrz(p1,p2,p3); addrrz(p2,p3,p2); setexpo(p2,expo(p2)-1);

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);

affrr(p2,y); setexpo(y, expo(y)+ex);

avma=av; return y;

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 z = mod16(x)==1? gun: stoi(3);

GEN y;

long zp;

gpmem_t av, lim;

e=valp(x); y=cgetg(5,t_PADIC); copyifstack(x[2],y[2]);

if (pp == 4) return z;

if (gcmp0(x))

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;

GEN mod;

y[3] = un;

zp = (zp<<1) - 1;

y[1] = evalvalp(e) | evalprecp(precp(x));

if (zp > pp) zp = pp;

return y;

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]);

long r = mod8(x);

if ((r!=1 && pp>=2) && (r!=5 || pp!=2)) err(sqrter5);

if (pp <= 3)

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(;;)

{

lpp=lp; lp=(lp<<1)-1;

switch(pp) {

if (lp < pp) y[3] = lshifti((GEN)y[3], lpp-1);

case 1: break;

else

case 2: if ((r&3) == 1) break;

{ lp=pp; y[3] = x[3]; }

case 3: if (r == 1) break;

setprecp(y,lp); y = gdiv(gadd(y, gdiv(x,y)), gdeux);

default: err(sqrter5);

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);

}

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]);

z = sqrt_padic(x, mod, pp, p);

if (!y[4]) err(sqrter5);

z = gerepileuptoint(av, z);

if (pp <= lp) { setprecp(y,1); return y; }

mod = icopy(mod);

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);

}

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))

Line 846 gsqrt(GEN x, long prec)

Line 920 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 = dummycopy(x); setvalp(y, 0);

y = ser_pui(x,ghalf,prec);

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)));

y = gmul(y, gpowgs(polx[varn(x)], e>>1));

return y;

return gerepileupto(av, y);

}

return transc(gsqrt,x,prec);

}

Line 864 gsqrt(GEN x, long prec)

Line 938 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;

Line 877 gsqrtz(GEN x, GEN y)

Line 952 gsqrtz(GEN x, GEN y)

GEN

rootsof1complex(GEN n, long prec)

{

ulong ltop=avma;

gpmem_t ltop = avma;

GEN z,a;

GEN z;

if (is_pm1(n)) return realun(prec);

if (lgefint(n)==3 && n[2]==2)

if (lgefint(n)==3 && n[2]==2) return stor(-1, prec);

{

a=realun(prec);

z = exp_Ir( divri(Pi2n(1, prec), n) ); /* exp(2Ipi/n) */

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);

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*/

Line 909 rootsof1padic(GEN n, GEN y)

Line 977 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;

Line 938 GEN padic_sqrtn_ram(GEN x, long e)

Line 1007 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;

Line 978 GEN padic_sqrtn_unram(GEN x, GEN n, GEN *zetan)

Line 1048 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]);

return padiczero(p, (valp(x)+m-1)/m);

y[4] = zero; e = (e+m-1)/m;

y[3] = un;

y[1] = evalvalp(e) | evalprecp(precp(x));

return y;

}

/*First treat the ramified part using logarithms*/

e=pvaluation(n,p,&q);

Line 1041 GEN padic_sqrtn(GEN x, GEN n, GEN *zetan)

Line 1108 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;

Line 1070 gsqrtn(GEN x, GEN n, GEN *zetan, long prec)

Line 1138 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 = dummycopy(x); setvalp(y, 0);

y = ser_pui(x,ginv(n),prec);

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)));

y = gmul(y, gpowgs(polx[varn(x)], e/m));

return y;

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]);

Line 1106 gsqrtn(GEN x, GEN n, GEN *zetan, long prec)

Line 1174 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;

Line 1135 gsqrtn(GEN x, GEN n, GEN *zetan, long prec)

Line 1211 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)

if (!sx) return realzero_bit(expo(x));

{

y=cgetr(3); y[1]=x[1]; y[2]=0; return y;

}

l=lg(x); y=cgetr(l); av=avma; /* room for result */

l2 = l+1; ex = expo(x);

Line 1205 mpexp1(GEN x)

Line 1279 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);

Line 1240 paexp(GEN x)

Line 1316 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<ex; i++) yd[i]=zero;

for ( ; i<ly; i++)

{

av = avma; p1 = gzero;

for (j=ex; j<=min(i,mi); j++)

p1 = gadd(p1, gmulgs(gmul((GEN)xd[j],(GEN)yd[i-j]),j));

yd[i] = lpileupto(av, gdivgs(p1,i));

}

return y;

}

av = avma; y = cgetg(lx, t_SER);

y[1] = x[1]; y[2] = zero;

for (i=3; i <lx; i++) y[i] = x[i];

p1 = gexp((GEN)x[2],prec);

y = gmul(p1, serexp(normalize(y),prec));

return gerepileupto(av, y);

}

GEN

gexp(GEN x, long prec)

{

long av,tetpil,ex;

GEN r,y,p1,p2;

switch(typ(x))

{

case t_REAL:

case t_REAL: return mpexp(x);

return mpexp(x);

case t_COMPLEX: return cxexp(x,prec);

case t_PADIC: return paexp(x);

case t_COMPLEX:

case t_SER: return serexp(x,prec);

y=cgetg(3,t_COMPLEX); av=avma;

case t_INTMOD: err(typeer,"gexp");

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;

case t_PADIC:

return paexp(x);

case t_SER:

if (gcmp0(x)) return gaddsg(1,x);

ex=valp(x); if (ex<0) err(negexper,"gexp");

if (ex)

{

long i,j,ly=lg(x)+ex;

y=cgetg(ly,t_SER);

y[1] = evalsigne(1) | evalvalp(0) | evalvarn(varn(x));

y[2] = un;

for (i=3; i<ex+2; i++) y[i]=zero;

for ( ; i<ly; i++)

{

av=avma; p1=gzero;

for (j=ex; j<i-1; j++)

p1=gadd(p1,gmulgs(gmul((GEN)x[j-ex+2],(GEN)y[i-j]),j));

tetpil=avma; y[i]=lpile(av,tetpil,gdivgs(p1,i-2));

}

return y;

}

av=avma; p1=gcopy(x); p1[2]=zero;

p2=gexp(normalize(p1),prec);

p1=gexp((GEN)x[2],prec); tetpil=avma;

return gerepile(av,tetpil,gmul(p1,p2));

case t_INTMOD:

err(typeer,"gexp");

}

return transc(gexp,x,prec);

}

Line 1298 gexp(GEN x, long prec)

Line 1385 gexp(GEN x, long prec)

void

gexpz(GEN x, GEN y)

{

long av=avma, prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"gexpz");

gaffect(gexp(x,prec),y); avma=av;

Line 1309 gexpz(GEN x, GEN y)

Line 1397 gexpz(GEN x, GEN y)

/** FONCTION LOGARITHME **/

/** **/

/********************************************************************/

/* 2 * atanh(1/3) */

GEN

constlog2(long prec)

{

static GEN glog2 = NULL;

const long _3 = 3, _9 = _3*_3;

gpmem_t av0, av;

long k, l, G;

GEN s, u, S, U, tmplog2;

if (glog2 && lg(glog2) >= 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;

gpmem_t ltop, av;

double alpha,beta,a,b;

long EX,l,l1,l2,m,n,k,ex,s;

GEN y,p1,p2,p3,p4,p5,unr;

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;

av=avma;

l = lg(x); EX = expo(x);

if (sgn > 0) p2 = divsr(1,p1); /* x<1 changer x en 1/x */

z = cgetr(l); ltop = avma;

for (m1=1; expo(p2)>0; m1++) p2 = mpsqrt(p2);

if (m1 > 1 || sgn > 0) { affrr(p2,p1); avma=av; }

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);

a = -log(alpha)/LOG2;

beta = (BITS_IN_LONG/2)*l*LOG2;

b = sqrt(BITS_IN_HALFULONG*l/3.0);

a = alpha/LOG2;

if (a <= b)

b = sqrt((BITS_IN_LONG/2)*l/3.0);

if (a<=b)

{

n=(long)(1+sqrt((BITS_IN_LONG/2)*3.0*l));

n = 1 + (long)sqrt(BITS_IN_HALFULONG*3.0*l);

m=(long)(1+b-a);

m = 1 + (long)(b-a);

l2 += m>>TWOPOTBITS_IN_LONG;

p4=cgetr(l2); affrr(p1,p4);

p4 = cgetr(l2); affrr(p1,p4);

p1 = p4; av=avma;

p1 = p4; av = avma;

for (i=1; i<=m; i++) p1 = mpsqrt(p1);

for (k=1; k<=m; k++) p1 = mpsqrt(p1);

affrr(p1,p4); avma=av;

affrr(p1,p4); avma = av;

}

else

{

n=(long)(1+beta/alpha);

n = 1 + (long)(BITS_IN_HALFULONG*l / a);

m=0; p4 = p1;

m = 0; p4 = p1;

}

unr=realun(l2);

unr = realun(l2);

p2=cgetr(l2); p3=cgetr(l2); av=avma;

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(divrr(p1,p5), y); /* = (X-1) / (X+1) ~ 0 */

affrr(mulrr(p2,p2), p3);

affrr(gsqr(y), y2); /* = y^2 */

affrr(divrs(unr,2*n+1), p4); setlg(p4,4); avma=av;

k = 2*n + 1;

affrr(divrs(unr,k), p4); setlg(p4,4); avma = av;

s=0; ex=expo(p3); l1 = 4;

s = 0; ex = expo(y2); l1 = 4;

for (i=n; i>=1; i--)

for (k -= 2; k > 0; k -= 2)

{

{ /* compute sum_i=0^n y^2i / (2i + 1), k = 2i+1 */

setlg(p3,l1); p5 = mulrr(p4,p3);

setlg(y2, l1); p5 = mulrr(p4,y2);

setlg(unr,l1); p1 = divrs(unr,2*i-1);

setlg(unr,l1); p1 = divrs(unr, k);

s -= ex;

l1 += s>>TWOPOTBITS_IN_LONG; if (l1>l2) l1=l2;

s %= BITS_IN_LONG;

s &= (BITS_IN_LONG-1);

setlg(p1,l1); setlg(p4,l1); setlg(p5,l1);

setlg(p1, l1);

p1 = addrr(p1,p5); affrr(p1,p4); avma=av;

setlg(p4, l1);

setlg(p5, l1); affrr(addrr(p1,p5), p4); avma=av;

}

setlg(p4,l2); affrr(mulrr(p2,p4), y);

setlg(p4, l2);

setexpo(y, expo(y)+m+m1);

y = mulrr(y,p4); /* = log(X)/2 */

if (sgn > 0) setsigne(y, -1);

setexpo(y, expo(y)+m+1);

avma=ltop; return y;

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;

GEN p,q,y,z,aux,p1;

long av,n,k;

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];

p = (GEN)x[2];

if (!cmpis((GEN)x[2],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);

p1 = addsi(-1, p);

if (n & 2)

z = resii(z, p);

addsiz(-1,(GEN)x[3],(GEN)y[4]);

aux = divii(addsi(-1,q),p1); n = precp(x);

else

for (k=1; k<n; k<<=1)

affsi(1,(GEN)y[4]);

z = modii(mulii(z,addsi(1,mulii(aux,addsi(-1,powmodulo(z,p1,q))))), q);

return y;

}

av = avma; p1 = addsi(-1,(GEN)x[2]);

affii(z, (GEN)y[4]); avma = av; return y;

aux = divii(addsi(-1,(GEN)x[3]),p1); n = precp(x);

for (k=1; k<n; k<<=1)

z=modii(mulii(z,addsi(1,mulii(aux,addsi(-1,powmodulo(z,p1,(GEN)x[3]))))),

(GEN)x[3]);

affii(z,(GEN)y[4]); avma=av; return y;

}

/* Let x = 1 mod p and y := (x-1)/(x+1) = 0 (p). Then

* log(x) = log(1+y) - log(1-y) = 2 \sum_{k odd} y^k / k.

* palogaux(x) returns the last sum (not multiplied by 2) */

static GEN

palogaux(GEN x)

{

long av1=avma,tetpil,k,e,pp;

long k,e,pp;

GEN y,s,y2;

GEN y,s,y2, p = (GEN)x[2];

if (egalii(gun,(GEN)x[4]))

if (egalii(gun, (GEN)x[4]))

{

y=gaddgs(x,-1);

long v = valp(x)+precp(x);

if (egalii(gdeux,(GEN)x[2]))

if (egalii(gdeux,p)) v--;

{

return padiczero(p, v);

setvalp(y,valp(y)-1);

if (!gcmp1((GEN)y[3])) y[3]=lshifti((GEN)y[3],-1);

}

return gerepilecopy(av1,y);

}

y=gdiv(gaddgs(x,-1),gaddgs(x,1)); e=valp(y); pp=e+precp(y);

y = gdiv(gaddgs(x,-1), gaddgs(x,1));

if (egalii(gdeux,(GEN)x[2])) pp--;

e = valp(y); pp = e+precp(y);

if (egalii(gdeux,p)) pp--;

else

{

long av=avma;

GEN p1;

for (p1=stoi(e); cmpsi(pp,p1)>0; pp++) p1 = mulii(p1, p);

for (p1=stoi(e); cmpsi(pp,p1)>0; pp++)

pp -= 2;

p1 = mulii(p1,(GEN)x[2]);

avma=av; pp-=2;

}

k=pp/e; if (!odd(k)) k--;

k = pp/e; if (!odd(k)) k--;

y2=gsqr(y); s=gdivgs(gun,k);

y2 = gsqr(y); s = gdivgs(gun,k);

while (k>=3)

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;

gpmem_t av = avma;

GEN p1,y;

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);

y = gsqr(x); setvalp(y,0);

affii(powmodulo((GEN)x[4],p1,(GEN)x[3]),(GEN)y[4]);

y = palogaux(y);

y=gmulgs(palogaux(y),2); tetpil=avma;

return gerepile(av,tetpil,gdiv(y,p1));

}

y=gsqr(x); setvalp(y,0); tetpil=avma;

else

return gerepile(av,tetpil,palogaux(y));

{ /* 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

Line 1471 log0(GEN x, long flag,long prec)

Line 1614 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))

Line 1492 glog(GEN x, long prec)

Line 1635 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);

Line 1509 glog(GEN x, long prec)

Line 1653 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;

Line 1517 glogz(GEN x, GEN y)

Line 1662 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

* take mmax=1518500248; on the alpha it does not seem worthwhile */

long l, l0, l1, l2, l4, i, n, m, s, t;

long l,l0,l1,l2,l4,ee,i,n,m,s,t,av;

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 (typ(x) != t_REAL) err(typeer,"mpsc1");

if (!signe(x))

l = lg(x); y = cgetr(l); av = avma;

{

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;

l++; pitemp = mppi(l+1); setexpo(pitemp,-1);

l++; pitemp = mppi(l);

p1 = addrr(x,pitemp); setexpo(pitemp,0);

setexpo(pitemp,-1);

if (expo(p1) >= bit_accuracy(min(l,lg(p1))) + 3)

p1 = addrr(x,pitemp); /* = x + Pi/4 */

err(precer,"mpsc1");

l0 = min(l, lg(p1));

p3 = divrr(p1,pitemp); p2 = mpent(p3);

if (expo(p1) >= bit_accuracy(l0) + 3) err(precer,"mpsc1");

if (signe(p2)) x = subrr(x, mulir(p2,pitemp));

p1 = cgetr(l+1); affrr(x, p1);

*ptmod8 = (signe(p1) < 0)? 4: 0;

setexpo(pitemp, 0);

if (signe(p2))

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);

long k = mod4(p1);

if (signe(p2) < 0 && mod4) mod4 = 4-mod4;

if (signe(p1) < 0 && k) k = 4-k;

*ptmod8 += mod4;

/* 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;

a = 0.5 / LOG2;

c = a+sqrt((beta+b)/LOG2);

b = 0.5 * a;

d = ((beta/c)-alpha-log(c))/LOG2;

c = a + sqrt((beta+b) / LOG2);

if (d>=0)

d = ((beta/c) - alpha - log(c)) / LOG2;

if (d >= 0)

{

m=(long)(1+d);

m = (long)(1+d);

n=(long)((1+c)/2.0);

n = (long)((1+c) / 2.0);

setexpo(p1,expo(p1)-m);

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);

p2 = realun(l2);

p4=cgetr(l2); av = avma;

x2 = cgetr(l2); av = avma;

affrr(gsqr(p1),p4); setlg(p4,4);

affrr(gsqr(x), x2);

if (n>mmax)

setlg(x2, 3);

p3 = divrs(divrs(p4,2*n+2),2*n+1);

if (n > mmax)

p3 = divrs(divrs(x2, 2*n+2), 2*n+1);

else

p3 = divrs(p4, (2*n+2)*(2*n+1));

p3 = divrs(x2, (2*n+2)*(2*n+1));

ee = -expo(p3); s=0;

s = -expo(p3);

l4 = l1 = 3 + (ee>>TWOPOTBITS_IN_LONG);

l4 = l1 = 3 + (s>>TWOPOTBITS_IN_LONG);

if (l4<=l2) { setlg(p2,l4); setlg(p4,l4); }

if (l4<=l2) { setlg(p2,l4); setlg(x2,l4); }

for (i=n; i>mmax; i--)

s = 0;

for (i=n; i>=2; i--)

{

p3 = divrs(divrs(p4,2*i),2*i-1); s -= expo(p3);

if (i > mmax)

t=s&(BITS_IN_LONG-1); l0=(s>>TWOPOTBITS_IN_LONG);

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); }

if (l4<=l2) { setlg(p2,l4); setlg(x2,l4); }

subsrz(1,p3,p2); avma=av;

subsrz(1,p3, p2); avma = av;

}

for ( ; i>=2; i--)

setlg(p2,l2); setlg(x2,l2);

{

setexpo(p2, expo(p2)-1); setsigne(p2, -signe(p2));

p3 = divrs(p4, 2*i*(2*i-1)); s -= expo(p3);

p2 = mulrr(x2,p2);

t=s&(BITS_IN_LONG-1); l0=(s>>TWOPOTBITS_IN_LONG);

/* Now p2 = sum {1<= i <=n} (-1)^i x^(2i) / (2i)! ~ cos(x) - 1 */

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);

for (i=1; i<=m; i++)

{

{ /* p2 = cos(x)-1 --> cos(2x)-1 */

p2 = mulrr(p2,addsr(2,p2)); setexpo(p2,expo(p2)+1);

p2 = mulrr(p2, addsr(2,p2)); setexpo(p2, expo(p2)+1);

}

affrr(p2,y); avma=av; return y;

}

/********************************************************************/

/* sqrt (1 - (1+x)^2) = sqrt(-x*(x+2)). Sends cos(x)-1 to |sin(x)| */

/** **/

/** FONCTION sqrt(-x*(x+2)) **/

/** **/

/********************************************************************/

static GEN

mpaut(GEN x)

{

long av = avma;

gpmem_t av = avma;

GEN p1 = mulrr(x,addsr(2,x));

GEN p1 = mulrr(x, addsr(2,x));

setsigne(p1,-signe(p1));

setsigne(p1, -signe(p1));

return gerepileuptoleaf(av, mpsqrt(p1));

}

/********************************************************************/

/** **/

/** COSINE **/

/** FONCTION COSINUS **/

/** **/

/********************************************************************/

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:

Line 1658 mpcos(GEN x)

Line 1795 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))

Line 1699 gcos(GEN x, long prec)

Line 1836 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;

}

/********************************************************************/

/** **/

/** SINE **/

/** FONCTION SINUS **/

/** **/

/********************************************************************/

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))

if (!signe(x)) return realzero_bit(expo(x));

{

y=cgetr(3); y[1]=x[1]; y[2]=0;

return y;

}

av=avma; p1=mpsc1(x,&mod8); tetpil=avma;

av = avma; p1 = mpsc1(x,&mod8);

switch(mod8)

{

case 0: case 6:

Line 1736 mpsin(GEN x)

Line 1869 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))

Line 1757 gsin(GEN x, long prec)

Line 1890 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;

Line 1777 gsin(GEN x, long prec)

Line 1911 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;

}

/********************************************************************/

/** **/

/** SINE, COSINE together **/

/** PROCEDURE SINUS,COSINUS **/

/** **/

/********************************************************************/

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];

*s = realzero_bit(expo(x));

p1[2]=0; *c=addsr(1,x); return;

*c = addsr(1,x); return;

}

av=avma; p1=mpsc1(x,&mod8); tetpil=avma;

Line 1818 mpsincos(GEN x, GEN *s, GEN *c)

Line 1952 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];

Line 1848 gsincos(GEN x, GEN *s, GEN *c, long prec)

Line 1994 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:

Line 1876 gsincos(GEN x, GEN *s, GEN *c, long prec)

Line 2023 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));

Line 1885 gsincos(GEN x, GEN *s, GEN *c, long prec)

Line 2034 gsincos(GEN x, GEN *s, GEN *c, long prec)

for (i=ex2; i<ly; i++)

{

ii=i-ex; av=avma; p1=gzero;

for (j=ex; j<ii-1; j++)

for (j=ex; j<=min(ii-2,mi); j++)

p1=gadd(p1,gmulgs(gmul((GEN)x[j-ex+2],(GEN)ps[ii-j]),j));

tetpil=avma;

pc[i]=lpile(av,tetpil,gdivgs(p1,2-i));

if (ii<lx)

{

av=avma; p1=gzero;

for (j=ex; j<=i-ex2; j++)

for (j=ex; j<=min(i-ex2,mi); j++)

p1=gadd(p1,gmulgs(gmul((GEN)x[j-ex+2],(GEN)pc[i-j]),j));

p1=gdivgs(p1,i-2); tetpil=avma;

ps[i-ex]=lpile(av,tetpil,gadd(p1,(GEN)x[i-ex]));

Line 1926 gsincos(GEN x, GEN *s, GEN *c, long prec)

Line 2075 gsincos(GEN x, GEN *s, GEN *c, long prec)

static GEN

mptan(GEN x)

{

long av=avma,tetpil;

gpmem_t av=avma, tetpil;

GEN s,c;

mpsincos(x,&s,&c); tetpil=avma;

Line 1936 mptan(GEN x)

Line 2085 mptan(GEN x)

GEN

gtan(GEN x, long prec)

{

long av,tetpil;

gpmem_t av, tetpil;

GEN s,c;

switch(typ(x))

Line 1960 gtan(GEN x, long prec)

Line 2109 gtan(GEN x, long prec)

void

gtanz(GEN x, GEN y)

{

long av=avma, prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"gtanz");

gaffect(gtan(x,prec),y); avma=av;

Line 1969 gtanz(GEN x, GEN y)

Line 2119 gtanz(GEN x, GEN y)

static GEN

mpcotan(GEN x)

{

long av=avma,tetpil;

gpmem_t av=avma, tetpil;

GEN s,c;

mpsincos(x,&s,&c); tetpil=avma;

Line 1979 mpcotan(GEN x)

Line 2129 mpcotan(GEN x)

GEN

gcotan(GEN x, long prec)

{

long av,tetpil;

gpmem_t av, tetpil;

GEN s,c;

switch(typ(x))

Line 1988 gcotan(GEN x, long prec)

Line 2138 gcotan(GEN x, long prec)

return mpcotan(x);

case t_SER:

if (gcmp0(x)) err(diver9);

if (gcmp0(x)) err(talker,"0 argument in cotan");

if (valp(x)<0) err(negexper,"gcotan"); /* fall through */

if (valp(x) < 0) err(negexper,"cotan"); /* fall through */

case t_COMPLEX:

av=avma; gsincos(x,&s,&c,prec); tetpil=avma;

return gerepile(av,tetpil,gdiv(c,s));

Line 2003 gcotan(GEN x, long prec)

Line 2153 gcotan(GEN x, long prec)

void

gcotanz(GEN x, GEN y)

{

long av=avma, prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"gcotanz");

gaffect(gcotan(x,prec),y); avma=av;

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>