version 1.1, 2001/10/02 11:17:05 |
version 1.2, 2002/09/11 07:26:52 |
|
|
constpi(long prec) |
constpi(long prec) |
{ |
{ |
GEN p1,p2,p3,tmppi; |
GEN p1,p2,p3,tmppi; |
long l,n,n1,av1,av2; |
long l, n, n1; |
|
gpmem_t av1, av2; |
double alpha; |
double alpha; |
|
|
#define k1 545140134 |
#define k1 545140134 |
Line 56 constpi(long prec) |
|
Line 57 constpi(long prec) |
|
|
|
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 */ |
/* initialisation longueur mantisse */ |
if (prec>=4) l=4; else l=prec; |
if (prec>=4) l=4; else l=prec; |
|
|
constpi(prec); affrr(gpi,x); return x; |
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 **/ |
/** FONCTION EULER **/ |
|
|
consteuler(long prec) |
consteuler(long prec) |
{ |
{ |
GEN u,v,a,b,tmpeuler; |
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; |
if (geuler && lg(geuler) >= prec) return; |
|
|
Line 116 consteuler(long prec) |
|
Line 139 consteuler(long prec) |
|
prec++; |
prec++; |
|
|
l = prec+1; x = (long) (1 + (bit_accuracy(l) >> 2) * LOG2); |
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 */ |
n=(long)(1+3.591*x); /* z=3.591: z*[ ln(z)-1 ]=1 */ |
av2 = avma; |
av2 = avma; |
if (x < SQRTVERYBIGINT) |
if (x < SQRTVERYBIGINT) |
|
|
transc(GEN (*f)(GEN,long), GEN x, long prec) |
transc(GEN (*f)(GEN,long), GEN x, long prec) |
{ |
{ |
GEN p1,p2,y; |
GEN p1,p2,y; |
long lx,i,av,tetpil; |
long lx, i; |
|
gpmem_t av, tetpil; |
|
|
av=avma; |
av=avma; |
switch(typ(x)) |
switch(typ(x)) |
|
|
return y; |
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) |
/* INTEGER POWERING (a^|n| for integer a and integer |n| > 1) |
* |
* |
* Nonpositive powers and exponent one should be handled by gpow() and |
* Nonpositive powers and exponent one should be handled by gpow() and |
|
|
static GEN |
static GEN |
puissii(GEN a, GEN n, long s) |
puissii(GEN a, GEN n, long s) |
{ |
{ |
long av,*p,m,k,i,lim; |
gpmem_t av; |
GEN y; |
GEN y; |
|
|
if (!signe(a)) return gzero; /* a==0 */ |
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] == 1) { a = icopy(a); setsigne(a,s); return a; } |
if (n[2] == 2) return sqri(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); |
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(). |
/* return a^n as a t_REAL of precision prec. Adapted from puissii(). |
* Assume a > 0, n > 0 */ |
* Assume a > 0, n > 0 */ |
GEN |
GEN |
rpowsi(ulong a, GEN n, long prec) |
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 |
GEN |
gpowgs(GEN x, long n) |
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; |
GEN y; |
|
|
if (n == 0) return puiss0(x); |
if (n == 0) return puiss0(x); |
Line 440 gpowgs(GEN x, long n) |
|
Line 466 gpowgs(GEN x, long n) |
|
GEN |
GEN |
pow_monome(GEN x, GEN nn) |
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; |
GEN p1,y; |
if (is_bigint(nn)) err(talker,"power overflow in pow_monome"); |
if (is_bigint(nn)) err(talker,"power overflow in pow_monome"); |
n = itos(nn); m = labs(n); |
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 |
GEN |
powgi(GEN x, GEN n) |
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 (typ(n) != t_INT) err(talker,"not an integral exponent in powgi"); |
if (!sn) return puiss0(x); |
if (!sn) return puiss0(x); |
Line 511 powgi(GEN x, GEN n) |
|
Line 537 powgi(GEN x, GEN n) |
|
} |
} |
case t_PADIC: |
case t_PADIC: |
{ |
{ |
GEN p, pe; |
long e = itos(n)*valp(x), v; |
|
GEN mod, p = (GEN)x[2]; |
|
|
if (!signe(x[4])) |
if (!signe(x[4])) |
{ |
{ |
if (sn < 0) err(talker, "division by 0 p-adic in powgi"); |
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); |
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 |
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]); |
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; |
return y; |
} |
} |
case t_QFR: |
case t_QFR: |
Line 537 powgi(GEN x, GEN n) |
|
Line 564 powgi(GEN x, GEN n) |
|
case t_POL: |
case t_POL: |
if (ismonome(x)) return pow_monome(x,n); |
if (ismonome(x)) return pow_monome(x,n); |
default: |
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); |
if (sn < 0) y = ginv(y); |
return av==avma? gcopy(y): gerepileupto(av,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 |
static GEN |
ser_pui(GEN x, GEN n, long prec) |
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 |
GEN |
gpow(GEN x, GEN n, long prec) |
gpow(GEN x, GEN n, long prec) |
{ |
{ |
long av,tetpil,i,lx,tx; |
long i, lx, tx; |
|
gpmem_t av, tetpil; |
GEN y; |
GEN y; |
|
|
if (typ(n)==t_INT) return powgi(x,n); |
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); |
for (i=1; i<lx; i++) y[i]=lpui((GEN)x[i],n,prec); |
return y; |
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 (tx==t_SER) |
{ |
{ |
if (valp(x)) |
if (valp(x)) |
err(talker,"not an integer exponent for non invertible series in gpow"); |
err(talker,"not an integer exponent for non invertible series in gpow"); |
if (lg(x) == 2) return gcopy(x); /* O(1) */ |
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; |
av=avma; |
if (gcmp0(x)) |
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)); |
x = ground(gmulsg(gexpo(x),n)); |
if (is_bigint(x) || (ulong)x[2] >= (ulong)HIGHEXPOBIT) |
if (is_bigint(x) || (ulong)x[2] >= (ulong)HIGHEXPOBIT) |
err(talker,"underflow or overflow in gpow"); |
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) |
if (tx==t_INTMOD && typ(n)==t_FRAC) |
{ |
{ |
GEN p1; |
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]); |
y=cgetg(3,tx); copyifstack(x[1],y[1]); |
av=avma; |
av=avma; |
p1=mpsqrtnmod((GEN)x[2],(GEN)n[2],(GEN)x[1],NULL); |
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 |
GEN |
mpsqrt(GEN x) |
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; |
double beta; |
GEN y,p1,p2,p3; |
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++) |
for (i=1; i<=n; i++) |
{ |
{ /* execute p2 := (p2 + p1/p2)/2 */ |
l0 = l1<<1; |
l0 = l1<<1; |
if (l0<=l) |
if (l0 <= l) |
{ l2 += l1; l1=l0; } |
{ l2 += l1; l1=l0; } |
else |
else |
{ l2 += l+1-l1; l1=l+1; } |
{ 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 |
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); |
if (e & 1) err(sqrter6); |
e>>=1; setvalp(y,e); y[3] = y[2]; |
|
|
y = cgetg(5,t_PADIC); |
pp = precp(x); |
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 */ |
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 |
GEN |
gsqrt(GEN x, long prec) |
gsqrt(GEN x, long prec) |
{ |
{ |
long av,tetpil,e; |
long e; |
|
gpmem_t av, tetpil; |
GEN y,p1,p2; |
GEN y,p1,p2; |
|
|
switch(typ(x)) |
switch(typ(x)) |
Line 846 gsqrt(GEN x, long prec) |
|
Line 920 gsqrt(GEN x, long prec) |
|
return padic_sqrt(x); |
return padic_sqrt(x); |
|
|
case t_SER: |
case t_SER: |
e=valp(x); |
e = valp(x); |
if (gcmp0(x)) return zeroser(varn(x), (e+1)>>1); |
if (gcmp0(x)) return zeroser(varn(x), (e+1)>>1); |
if (e & 1) err(sqrter6); |
if (e & 1) err(sqrter6); |
av = avma; |
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 */ |
if (typ(y) == t_SER) /* generic case */ |
y[1] = evalsigne(1) | evalvalp(e>>1) | evalvarn(varn(x)); |
y[1] = evalsigne(1) | evalvalp(e>>1) | evalvarn(varn(x)); |
else /* e.g coeffs are POLMODs */ |
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); |
return transc(gsqrt,x,prec); |
} |
} |
Line 864 gsqrt(GEN x, long prec) |
|
Line 938 gsqrt(GEN x, long prec) |
|
void |
void |
gsqrtz(GEN x, GEN y) |
gsqrtz(GEN x, GEN y) |
{ |
{ |
long av=avma, prec = precision(y); |
long prec = precision(y); |
|
gpmem_t av=avma; |
|
|
if (!prec) err(infprecer,"gsqrtz"); |
if (!prec) err(infprecer,"gsqrtz"); |
gaffect(gsqrt(x,prec),y); avma=av; |
gaffect(gsqrt(x,prec),y); avma=av; |
Line 877 gsqrtz(GEN x, GEN y) |
|
Line 952 gsqrtz(GEN x, GEN y) |
|
GEN |
GEN |
rootsof1complex(GEN n, long prec) |
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 (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); |
return gerepileupto(ltop,z); |
} |
} |
/*Only the O() of y is used*/ |
/*Only the O() of y is used*/ |
GEN |
GEN |
rootsof1padic(GEN n, GEN y) |
rootsof1padic(GEN n, GEN y) |
{ |
{ |
ulong ltop=avma; |
gpmem_t ltop=avma; |
GEN z,r; |
GEN z,r; |
(void)mpsqrtnmod(gun,n,(GEN)y[2],&z); |
(void)mpsqrtnmod(gun,n,(GEN)y[2],&z); |
if (z==gzero){avma=ltop;return gzero;}/*should not happen*/ |
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); |
static GEN paexp(GEN x); |
/*compute the p^e th root of x p-adic*/ |
/*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 n,a; |
GEN p=(GEN)x[2]; |
GEN p=(GEN)x[2]; |
long v=0; |
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; |
return a; |
} |
} |
/*compute the nth root of x p-adic p prime with n*/ |
/*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 a,r; |
GEN p=(GEN)x[2]; |
GEN p=(GEN)x[2]; |
long v=0; |
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 |
else |
return gerepile(ltop,tetpil,r); |
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 p=(GEN)x[2]; |
GEN q; |
GEN q; |
long e; |
long e; |
GEN *gptr[2]; |
GEN *gptr[2]; |
if (gcmp0(x)) |
if (gcmp0(x)) |
{ |
{ |
GEN y; |
|
long m = itos(n); |
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*/ |
/*First treat the ramified part using logarithms*/ |
e=pvaluation(n,p,&q); |
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 |
GEN |
gsqrtn(GEN x, GEN n, GEN *zetan, long prec) |
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; |
long e,m; |
GEN y,z; |
GEN y,z; |
if (zetan) *zetan=gzero; |
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: |
case t_SER: |
e=valp(x);m=itos(n); |
e=valp(x);m=itos(n); |
if (gcmp0(x)) return zeroser(varn(x), (e+m-1)/m); |
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; |
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 */ |
if (typ(y) == t_SER) /* generic case */ |
y[1] = evalsigne(1) | evalvalp(e/m) | evalvarn(varn(x)); |
y[1] = evalsigne(1) | evalvalp(e/m) | evalvarn(varn(x)); |
else /* e.g coeffs are POLMODs */ |
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: |
case t_INTMOD: |
z=gzero; |
z=gzero; |
/*This is not great, but else it will generate too much trouble*/ |
/*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) |
if (zetan) |
{ |
{ |
z=cgetg(3,tx); copyifstack(x[1],z[1]); |
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) |
if (tx==t_INT && is_pm1(x) && signe(x)>0) |
y=gun; /*speed-up since there is no way to call rootsof1complex |
y=gun; /*speed-up since there is no way to call rootsof1complex |
directly from gp*/ |
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 |
else |
{ |
{ |
av=avma; |
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 |
GEN |
mpexp1(GEN x) |
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 */; |
double a,b,alpha,beta, gama = 2.0 /* optimized for SUN3 */; |
/* KB: 3.0 is better for UltraSparc */ |
/* KB: 3.0 is better for UltraSparc */ |
GEN y,p1,p2,p3,p4,unr; |
GEN y,p1,p2,p3,p4,unr; |
|
|
if (typ(x)!=t_REAL) err(typeer,"mpexp1"); |
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 */ |
l=lg(x); y=cgetr(l); av=avma; /* room for result */ |
|
|
l2 = l+1; ex = expo(x); |
l2 = l+1; ex = expo(x); |
|
|
GEN |
GEN |
mpexp(GEN x) |
mpexp(GEN x) |
{ |
{ |
long av, sx = signe(x); |
long sx = signe(x); |
|
gpmem_t av; |
GEN y; |
GEN y; |
|
|
if (!sx) return addsr(1,x); |
if (!sx) return addsr(1,x); |
if (sx<0) setsigne(x,1); |
if (sx<0) setsigne(x,1); |
av = avma; y = addsr(1, mpexp1(x)); |
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); |
return gerepileupto(av,y); |
} |
} |
|
|
static GEN |
static GEN |
paexp(GEN x) |
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; |
GEN y,r,p1; |
|
|
if (gcmp0(x)) return gaddgs(x,1); |
if (gcmp0(x)) return gaddgs(x,1); |
|
|
return gerepileupto(av, y); |
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 |
GEN |
gexp(GEN x, long prec) |
gexp(GEN x, long prec) |
{ |
{ |
long av,tetpil,ex; |
|
GEN r,y,p1,p2; |
|
|
|
switch(typ(x)) |
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); |
return transc(gexp,x,prec); |
} |
} |
Line 1298 gexp(GEN x, long prec) |
|
Line 1385 gexp(GEN x, long prec) |
|
void |
void |
gexpz(GEN x, GEN y) |
gexpz(GEN x, GEN y) |
{ |
{ |
long av=avma, prec = precision(y); |
long prec = precision(y); |
|
gpmem_t av=avma; |
|
|
if (!prec) err(infprecer,"gexpz"); |
if (!prec) err(infprecer,"gexpz"); |
gaffect(gexp(x,prec),y); avma=av; |
gaffect(gexp(x,prec),y); avma=av; |
Line 1309 gexpz(GEN x, GEN y) |
|
Line 1397 gexpz(GEN x, GEN y) |
|
/** FONCTION LOGARITHME **/ |
/** 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 |
GEN |
|
mplog2(long prec) |
|
{ |
|
GEN x = cgetr(prec); |
|
affrr(constlog2(prec), x); return x; |
|
} |
|
|
|
GEN |
mplog(GEN x) |
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 (typ(x)!=t_REAL) err(typeer,"mplog"); |
if (signe(x)<=0) err(talker,"non positive argument in 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; |
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; |
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 |
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 */ |
/* affrr needed here instead of setlg since prec may DECREASE */ |
p1 = cgetr(l2); affrr(subrr(p4,unr), p1); |
p1 = cgetr(l2); affrr(subrr(p4,unr), p1); |
|
|
p5 = addrr(p4,unr); setlg(p5,l2); |
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; |
s -= ex; |
l1 += s>>TWOPOTBITS_IN_LONG; if (l1>l2) l1=l2; |
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 |
GEN |
teich(GEN x) |
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 (typ(x)!=t_PADIC) err(talker,"not a p-adic argument in teichmuller"); |
if (!signe(x[4])) return gcopy(x); |
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 |
static GEN |
palogaux(GEN x) |
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 |
else |
{ |
{ |
long av=avma; |
|
GEN p1; |
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 |
GEN |
palog(GEN x) |
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 (!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 |
GEN |
Line 1471 log0(GEN x, long flag,long prec) |
|
Line 1614 log0(GEN x, long flag,long prec) |
|
GEN |
GEN |
glog(GEN x, long prec) |
glog(GEN x, long prec) |
{ |
{ |
long av,tetpil; |
gpmem_t av, tetpil; |
GEN y,p1,p2; |
GEN y,p1,p2; |
|
|
switch(typ(x)) |
switch(typ(x)) |
Line 1492 glog(GEN x, long prec) |
|
Line 1635 glog(GEN x, long prec) |
|
return palog(x); |
return palog(x); |
|
|
case t_SER: |
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); |
p1=gdiv(derivser(x),x); |
tetpil=avma; p1=integ(p1,varn(x)); |
tetpil=avma; p1=integ(p1,varn(x)); |
if (gcmp1((GEN)x[2])) return gerepile(av,tetpil,p1); |
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 |
void |
glogz(GEN x, GEN y) |
glogz(GEN x, GEN y) |
{ |
{ |
long av=avma, prec = precision(y); |
long prec = precision(y); |
|
gpmem_t av=avma; |
|
|
if (!prec) err(infprecer,"glogz"); |
if (!prec) err(infprecer,"glogz"); |
gaffect(glog(x,prec),y); avma=av; |
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 |
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 |
/* 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; |
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 |
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; |
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); |
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 |
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++; |
if (t) l0++; |
l1 += l0; if (l1>l2) { l0 += l2-l1; l1=l2; } |
l1 += l0; if (l1>l2) { l0 += l2-l1; l1=l2; } |
l4 += l0; |
l4 += l0; |
p3 = mulrr(p3,p2); |
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++) |
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; |
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 |
static GEN |
mpaut(GEN x) |
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)); |
return gerepileuptoleaf(av, mpsqrt(p1)); |
} |
} |
|
|
/********************************************************************/ |
/********************************************************************/ |
/** **/ |
/** COSINE **/ |
/** FONCTION COSINUS **/ |
|
/** **/ |
|
/********************************************************************/ |
/********************************************************************/ |
|
|
static GEN |
static GEN |
mpcos(GEN x) |
mpcos(GEN x) |
{ |
{ |
long mod8,av,tetpil; |
long mod8; |
|
gpmem_t av; |
GEN y,p1; |
GEN y,p1; |
|
|
if (typ(x)!=t_REAL) err(typeer,"mpcos"); |
if (typ(x)!=t_REAL) err(typeer,"mpcos"); |
if (!signe(x)) return addsr(1,x); |
if (!signe(x)) return addsr(1,x); |
|
|
av=avma; p1=mpsc1(x,&mod8); tetpil=avma; |
av = avma; p1 = mpsc1(x,&mod8); |
switch(mod8) |
switch(mod8) |
{ |
{ |
case 0: case 4: |
case 0: case 4: |
|
|
default: /* case 3: case 5: */ |
default: /* case 3: case 5: */ |
y=mpaut(p1); break; |
y=mpaut(p1); break; |
} |
} |
return gerepile(av,tetpil,y); |
return gerepileuptoleaf(av, y); |
} |
} |
|
|
GEN |
GEN |
gcos(GEN x, long prec) |
gcos(GEN x, long prec) |
{ |
{ |
long av,tetpil; |
gpmem_t av, tetpil; |
GEN r,u,v,y,p1,p2; |
GEN r,u,v,y,p1,p2; |
|
|
switch(typ(x)) |
switch(typ(x)) |
Line 1699 gcos(GEN x, long prec) |
|
Line 1836 gcos(GEN x, long prec) |
|
void |
void |
gcosz(GEN x, GEN y) |
gcosz(GEN x, GEN y) |
{ |
{ |
long av = avma, prec = precision(y); |
long prec = precision(y); |
|
gpmem_t av = avma; |
|
|
if (!prec) err(infprecer,"gcosz"); |
if (!prec) err(infprecer,"gcosz"); |
gaffect(gcos(x,prec),y); avma=av; |
gaffect(gcos(x,prec),y); avma=av; |
} |
} |
|
|
/********************************************************************/ |
/********************************************************************/ |
/** **/ |
/** SINE **/ |
/** FONCTION SINUS **/ |
|
/** **/ |
|
/********************************************************************/ |
/********************************************************************/ |
|
|
GEN |
GEN |
mpsin(GEN x) |
mpsin(GEN x) |
{ |
{ |
long mod8,av,tetpil; |
long mod8; |
|
gpmem_t av; |
GEN y,p1; |
GEN y,p1; |
|
|
if (typ(x)!=t_REAL) err(typeer,"mpsin"); |
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) |
switch(mod8) |
{ |
{ |
case 0: case 6: |
case 0: case 6: |
|
|
default: /* case 3: case 7: */ |
default: /* case 3: case 7: */ |
y=subsr(-1,p1); break; |
y=subsr(-1,p1); break; |
} |
} |
return gerepile(av,tetpil,y); |
return gerepileuptoleaf(av, y); |
} |
} |
|
|
GEN |
GEN |
gsin(GEN x, long prec) |
gsin(GEN x, long prec) |
{ |
{ |
long av,tetpil; |
gpmem_t av, tetpil; |
GEN r,u,v,y,p1,p2; |
GEN r,u,v,y,p1,p2; |
|
|
switch(typ(x)) |
switch(typ(x)) |
Line 1757 gsin(GEN x, long prec) |
|
Line 1890 gsin(GEN x, long prec) |
|
p1=gsub(p2,p1); |
p1=gsub(p2,p1); |
gsincos((GEN)x[1],&u,&v,prec); |
gsincos((GEN)x[1],&u,&v,prec); |
tetpil=avma; |
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); |
gerepilemanyvec(av,tetpil,y+1,2); |
return y; |
return y; |
|
|
Line 1777 gsin(GEN x, long prec) |
|
Line 1911 gsin(GEN x, long prec) |
|
void |
void |
gsinz(GEN x, GEN y) |
gsinz(GEN x, GEN y) |
{ |
{ |
long av=avma, prec = precision(y); |
long prec = precision(y); |
|
gpmem_t av=avma; |
|
|
if (!prec) err(infprecer,"gsinz"); |
if (!prec) err(infprecer,"gsinz"); |
gaffect(gsin(x,prec),y); avma=av; |
gaffect(gsin(x,prec),y); avma=av; |
} |
} |
|
|
/********************************************************************/ |
/********************************************************************/ |
/** **/ |
/** SINE, COSINE together **/ |
/** PROCEDURE SINUS,COSINUS **/ |
|
/** **/ |
|
/********************************************************************/ |
/********************************************************************/ |
|
|
void |
void |
mpsincos(GEN x, GEN *s, GEN *c) |
mpsincos(GEN x, GEN *s, GEN *c) |
{ |
{ |
long av,tetpil,mod8; |
long mod8; |
|
gpmem_t av, tetpil; |
GEN p1, *gptr[2]; |
GEN p1, *gptr[2]; |
|
|
if (typ(x)!=t_REAL) err(typeer,"mpsincos"); |
if (typ(x)!=t_REAL) err(typeer,"mpsincos"); |
if (!signe(x)) |
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; |
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); |
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 |
void |
gsincos(GEN x, GEN *s, GEN *c, long prec) |
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 r,u,v,u1,v1,p1,p2,p3,p4,ps,pc; |
GEN *gptr[4]; |
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); |
p3=gmul(v1,v); p4=gmul(r,u); |
gptr[0]=&p1; gptr[1]=&p2; gptr[2]=&p3; gptr[3]=&p4; |
gptr[0]=&p1; gptr[1]=&p2; gptr[2]=&p3; gptr[3]=&p4; |
gerepilemanysp(av,tetpil,gptr,4); |
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; |
return; |
|
|
case t_SER: |
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; |
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; |
pc=cgetg(ly,t_SER); *c=pc; |
ps=cgetg(lx,t_SER); *s=ps; |
ps=cgetg(lx,t_SER); *s=ps; |
pc[1] = evalsigne(1) | evalvalp(0) | evalvarn(varn(x)); |
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++) |
for (i=ex2; i<ly; i++) |
{ |
{ |
ii=i-ex; av=avma; p1=gzero; |
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)); |
p1=gadd(p1,gmulgs(gmul((GEN)x[j-ex+2],(GEN)ps[ii-j]),j)); |
tetpil=avma; |
tetpil=avma; |
pc[i]=lpile(av,tetpil,gdivgs(p1,2-i)); |
pc[i]=lpile(av,tetpil,gdivgs(p1,2-i)); |
if (ii<lx) |
if (ii<lx) |
{ |
{ |
av=avma; p1=gzero; |
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=gadd(p1,gmulgs(gmul((GEN)x[j-ex+2],(GEN)pc[i-j]),j)); |
p1=gdivgs(p1,i-2); tetpil=avma; |
p1=gdivgs(p1,i-2); tetpil=avma; |
ps[i-ex]=lpile(av,tetpil,gadd(p1,(GEN)x[i-ex])); |
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 |
static GEN |
mptan(GEN x) |
mptan(GEN x) |
{ |
{ |
long av=avma,tetpil; |
gpmem_t av=avma, tetpil; |
GEN s,c; |
GEN s,c; |
|
|
mpsincos(x,&s,&c); tetpil=avma; |
mpsincos(x,&s,&c); tetpil=avma; |
|
|
GEN |
GEN |
gtan(GEN x, long prec) |
gtan(GEN x, long prec) |
{ |
{ |
long av,tetpil; |
gpmem_t av, tetpil; |
GEN s,c; |
GEN s,c; |
|
|
switch(typ(x)) |
switch(typ(x)) |
Line 1960 gtan(GEN x, long prec) |
|
Line 2109 gtan(GEN x, long prec) |
|
void |
void |
gtanz(GEN x, GEN y) |
gtanz(GEN x, GEN y) |
{ |
{ |
long av=avma, prec = precision(y); |
long prec = precision(y); |
|
gpmem_t av=avma; |
|
|
if (!prec) err(infprecer,"gtanz"); |
if (!prec) err(infprecer,"gtanz"); |
gaffect(gtan(x,prec),y); avma=av; |
gaffect(gtan(x,prec),y); avma=av; |
Line 1969 gtanz(GEN x, GEN y) |
|
Line 2119 gtanz(GEN x, GEN y) |
|
static GEN |
static GEN |
mpcotan(GEN x) |
mpcotan(GEN x) |
{ |
{ |
long av=avma,tetpil; |
gpmem_t av=avma, tetpil; |
GEN s,c; |
GEN s,c; |
|
|
mpsincos(x,&s,&c); tetpil=avma; |
mpsincos(x,&s,&c); tetpil=avma; |
|
|
GEN |
GEN |
gcotan(GEN x, long prec) |
gcotan(GEN x, long prec) |
{ |
{ |
long av,tetpil; |
gpmem_t av, tetpil; |
GEN s,c; |
GEN s,c; |
|
|
switch(typ(x)) |
switch(typ(x)) |
Line 1988 gcotan(GEN x, long prec) |
|
Line 2138 gcotan(GEN x, long prec) |
|
return mpcotan(x); |
return mpcotan(x); |
|
|
case t_SER: |
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: |
case t_COMPLEX: |
av=avma; gsincos(x,&s,&c,prec); tetpil=avma; |
av=avma; gsincos(x,&s,&c,prec); tetpil=avma; |
return gerepile(av,tetpil,gdiv(c,s)); |
return gerepile(av,tetpil,gdiv(c,s)); |
Line 2003 gcotan(GEN x, long prec) |
|
Line 2153 gcotan(GEN x, long prec) |
|
void |
void |
gcotanz(GEN x, GEN y) |
gcotanz(GEN x, GEN y) |
{ |
{ |
long av=avma, prec = precision(y); |
long prec = precision(y); |
|
gpmem_t av=avma; |
|
|
if (!prec) err(infprecer,"gcotanz"); |
if (!prec) err(infprecer,"gcotanz"); |
gaffect(gcotan(x,prec),y); avma=av; |
gaffect(gcotan(x,prec),y); avma=av; |