/********************************************************************/
/** **/
/** TRANSCENDENTAL FONCTIONS **/
/** **/
/********************************************************************/
/* $Id: trans1.c,v 1.2 1999/09/23 17:50:57 karim Exp $ */
#include "pari.h"
#ifdef LONG_IS_64BIT
# define C31 (9223372036854775808.0) /* VERYBIGINT * 1.0 */
# define SQRTVERYBIGINT 3037000500 /* ceil(sqrt(VERYBIGINT)) */
# define CBRTVERYBIGINT 2097152 /* ceil(cbrt(VERYBIGINT)) */
#else
# define C31 (2147483648.0)
# define SQRTVERYBIGINT 46341
# define CBRTVERYBIGINT 1291
#endif
/********************************************************************/
/** **/
/** FONCTION PI **/
/** **/
/********************************************************************/
void
constpi(long prec)
{
GEN p1,p2,p3,tmppi;
long l,n,n1,av1,av2;
double alpha;
#define k1 545140134
#define k2 13591409
#define k3 640320
#define alpha2 (1.4722004/(BYTES_IN_LONG/4)) /* 3*log(k3/12)/(32*log(2)) */
if (gpi && lg(gpi) >= prec) return;
av1 = avma; tmppi = newbloc(prec);
*tmppi = evaltyp(t_REAL) | evallg(prec);
n=(long)(1+(prec-2)/alpha2);
n1=6*n-1; p1=cgetr(prec);
p2=addsi(k2,mulss(n,k1));
affir(p2,p1);
/* initialisation longueur mantisse */
if (prec>=4) l=4; else l=prec;
setlg(p1,l); alpha=l;
av2 = avma;
while (n)
{
if (n < CBRTVERYBIGINT)
p3 = divrs(mulsr(n1-4,mulsr(n1*(n1-2),p1)),n*n*n);
else
{
if (n1 < SQRTVERYBIGINT)
p3 = divrs(divrs(mulsr(n1-4,mulsr(n1*(n1-2),p1)),n*n),n);
else
p3 = divrs(divrs(mulsr(n1-4,mulsr(n1,mulsr(n1-2,p1))),n*n),n);
}
p3 = divrs(divrs(p3,100100025), 327843840);
subisz(p2,k1,p2); subirz(p2,p3,p1);
alpha += alpha2; l = (long)(1+alpha); /* nouvelle longueur mantisse */
if (l>prec) l=prec;
setlg(p1,l); avma = av2;
n--; n1-=6;
}
p1 = divsr(53360,p1);
mulrrz(p1,gsqrt(stoi(k3),prec), tmppi);
gunclone(gpi); avma = av1; gpi = tmppi;
}
GEN
mppi(long prec)
{
GEN x = cgetr(prec);
constpi(prec+1); affrr(gpi,x); return x;
}
/********************************************************************/
/** **/
/** FONCTION EULER **/
/** **/
/********************************************************************/
void
consteuler(long prec)
{
GEN u,v,a,b,tmpeuler;
long l,n,k,x,av1,av2;
if (geuler && lg(geuler) >= prec) return;
av1 = avma; tmpeuler = newbloc(prec);
*tmpeuler = evaltyp(t_REAL) | evallg(prec);
l = prec+1; x = (long) (1 + (bit_accuracy(l) >> 2) * LOG2);
a=cgetr(l); affsr(x,a); u=mplog(a); setsigne(u,-1); affrr(u,a);
b=realun(l);
v=realun(l);
n=(long)(1+3.591*x); /* z=3.591: z*[ ln(z)-1 ]=1 */
av2 = avma;
if (x < SQRTVERYBIGINT)
{
long xx = x*x;
for (k=1; k<=n; k++)
{
divrsz(mulsr(xx,b),k*k, b);
divrsz(addrr(divrs(mulsr(xx,a),k),b),k, a);
addrrz(u,a,u); addrrz(v,b,v); avma = av2;
}
}
else
{
GEN xx = mulss(x,x);
for (k=1; k<=n; k++)
{
divrsz(mulir(xx,b),k*k, b);
divrsz(addrr(divrs(mulir(xx,a),k),b),k, a);
addrrz(u,a,u); addrrz(v,b,v); avma = av2;
}
}
divrrz(u,v,tmpeuler);
gunclone(geuler); avma = av1; geuler = tmpeuler;
}
GEN
mpeuler(long prec)
{
GEN x=cgetr(prec);
consteuler(prec+1); affrr(geuler,x); return x;
}
/********************************************************************/
/** **/
/** CONVERSION DE TYPES POUR FONCTIONS TRANSCENDANTES **/
/** **/
/********************************************************************/
GEN
transc(GEN (*f)(GEN,long), GEN x, long prec)
{
GEN p1,p2,y;
long lx,i,av,tetpil;
av=avma;
switch(typ(x))
{
case t_INT: case t_FRAC: case t_FRACN:
p1=cgetr(prec); gaffect(x,p1); tetpil=avma;
return gerepile(av,tetpil,f(p1,prec));
case t_COMPLEX: case t_QUAD:
av=avma; p1=gmul(x,realun(prec)); tetpil=avma;
return gerepile(av,tetpil,f(p1,prec));
case t_POL: case t_RFRAC: case t_RFRACN:
p1=tayl(x,gvar(x),precdl); tetpil=avma;
return gerepile(av,tetpil,f(p1,prec));
case t_VEC: case t_COL: case t_MAT:
lx=lg(x); y=cgetg(lx,typ(x));
for (i=1; i<lx; i++)
y[i] = (long) f((GEN)x[i],prec);
return y;
case t_POLMOD:
av=avma; p1=roots((GEN)x[1],prec); lx=lg(p1); p2=cgetg(lx,t_COL);
for (i=1; i<lx; i++) p2[i]=lpoleval((GEN)x[2],(GEN)p1[i]);
tetpil=avma; y=cgetg(lx,t_COL);
for (i=1; i<lx; i++)
y[i] = (long)f((GEN)p2[i],prec);
return gerepile(av,tetpil,y);
case t_QFR: case t_QFI:
err(talker,"quadratic forms cannot be used in transcendental functions");
}
return f(x,prec);
}
/*******************************************************************/
/* */
/* POWERING */
/* */
/*******************************************************************/
GEN real_unit_form(GEN x);
GEN imag_unit_form(GEN x);
static GEN
puiss0(GEN x)
{
long lx,i;
GEN y;
switch(typ(x))
{
case t_INT: case t_REAL: case t_FRAC: case t_FRACN: case t_PADIC:
return gun;
case t_INTMOD:
y=cgetg(3,t_INTMOD); copyifstack(x[1],y[1]); y[2]=un; break;
case t_COMPLEX:
y=cgetg(3,t_COMPLEX); y[1]=un; y[2]=zero; break;
case t_QUAD:
y=cgetg(4,t_QUAD); copyifstack(x[1],y[1]);
y[2]=un; y[3]=zero; break;
case t_POLMOD:
y=cgetg(3,t_POLMOD); copyifstack(x[1],y[1]);
y[2]=lpolun[varn(x[1])]; break;
case t_POL: case t_SER: case t_RFRAC: case t_RFRACN:
return polun[gvar(x)];
case t_MAT:
lx=lg(x); if (lx==1) return cgetg(1,t_MAT);
if (lx != (lg(x[1]))) err(mattype1,"gpowgs");
y = idmat(lx-1);
for (i=1; i<lx; i++) coeff(y,i,i) = lpuigs(gcoeff(x,i,i),0);
break;
case t_QFR: return real_unit_form(x);
case t_QFI: return imag_unit_form(x);
case t_VEC: case t_COL:
err(typeer,"gpowgs");
}
return y;
}
/* INTEGER POWERING (a^n for integer a and integer n > 1)
*
* Nonpositive powers and exponent one should be handled by gpow() and
* gpowgs() in trans1.c, which at the time of this writing is the only place
* where the following is [slated to be] used.
*
* Use left shift binary algorithm (RS is wasteful: multiplies big numbers,
* with LS one of them is the base, hence small). If result is nonzero, its
* sign is set to s (= +-1) regardless of what it would have been. This makes
* life easier for gpow()/gpowgs(), which otherwise might do a
* setsigne(gun,-1)... -- GN1998May04
*/
static GEN
puissii(GEN a, GEN n, long s)
{
long av,*p,man,k,nb,lim;
GEN y;
if (!signe(a)) return gzero; /* a==0 */
if (lgefint(a)==3)
{ /* easy if |a| < 3 */
if (a[2] == 1) return (s>0)? gun: negi(gun);
if (a[2] == 2) { a = shifti(gun, itos(n)); setsigne(a,s); return a; }
}
if (lgefint(n)==3)
{ /* or if |n| < 3 */
if (n[2] == 1) { a = icopy(a); setsigne(a,s); return a; }
if (n[2] == 2) return sqri(a);
}
/* be paranoid about memory consumption */
av=avma; lim=stack_lim(av,1);
y = a; p = n+2; man = *p;
/* normalize, i.e set highest bit to 1 (we know man != 0) */
k = 1+bfffo(man); man<<=k; k = BITS_IN_LONG-k;
/* first bit is now implicit */
for (nb=lgefint(n)-2;;)
{
for (; k; man<<=1,k--)
{
y = sqri(y);
if (low_stack(lim, stack_lim(av,1)))
{
if (DEBUGMEM>1) err(warnmem,"[1]: puissii");
y = gerepileupto(av,y);
}
if (man < 0)
{ /* first bit is set: multiply by the base */
y = mulii(y,a);
if (low_stack(lim, stack_lim(av,1)))
{
if (DEBUGMEM>1) err(warnmem,"[2]: puissii");
y = gerepileupto(av,y);
}
}
}
nb--; if (nb == 0) break;
man = *++p, k = BITS_IN_LONG;
}
setsigne(y,s); return gerepileupto(av,y);
}
GEN
gpowgs(GEN x, long n)
{
long lim,av,m,tx;
static long gn[3] = {evaltyp(t_INT)|m_evallg(3), 0, 0};
GEN y;
if (n == 0) return puiss0(x);
if (n == 1) return gcopy(x);
if (n ==-1) return ginv(x);
if (n>0)
{ gn[1] = evalsigne( 1) | evallgefint(3); gn[2]= n; }
else
{ gn[1] = evalsigne(-1) | evallgefint(3); gn[2]=-n; }
/* If gpowgs were only ever called from gpow, the switch wouldn't be needed.
* In fact, under current word and bit field sizes, an integer power with
* multiword exponent will always overflow. But it seems easier to call
* puissii|powmodulo() with a mock-up GEN as 2nd argument than to write
* separate versions for a long exponent. Note that n = HIGHBIT is an
* invalid argument. --GN
*/
switch((tx=typ(x)))
{
case t_INT:
{
long sx=signe(x), sr = (sx<0 && (n&1))? -1: 1;
if (n>0) return puissii(x,(GEN)gn,sr);
if (!sx) err(talker, "division by zero in gpowgs");
if (is_pm1(x)) return (sr < 0)? icopy(x): gun;
/* n<0, |x|>1 */
y=cgetg(3,t_FRAC); setsigne(gn,1);
y[1]=(sr>0)? un: lnegi(gun);
y[2]=(long)puissii(x,(GEN)gn,1); /* force denominator > 0 */
return y;
}
case t_INTMOD:
y=cgetg(3,tx); copyifstack(x[1],y[1]);
y[2]=(long)powmodulo((GEN)(x[2]),(GEN)gn,(GEN)(x[1]));
return y;
case t_FRAC: case t_FRACN:
{
long sr = ((n&1) && (signe(x[1])!=signe(x[2])))? -1: 1;
if (n<0)
{
if (!signe((GEN)x[1]))
err(talker, "division by zero fraction in gpowgs");
setsigne((GEN)gn,1);
if (is_pm1((GEN)x[1])) /* +-1/x[2] inverts to an integer */
return puissii((GEN)x[2],(GEN)gn,sr);
y=cgetg(3,tx);
y[1]=(long)puissii((GEN)x[2],(GEN)gn,sr);
y[2]=(long)puissii((GEN)x[1],(GEN)gn,1);
return y;
}
else /* n > 0 */
{
y=cgetg(3,tx);
if (!signe((GEN)x[1])) { y[1]=zero; y[2]=un; return y; }
y[1]=(long)puissii((GEN)x[1],(GEN)gn,sr);
y[2]=(long)puissii((GEN)x[2],(GEN)gn,1);
return y;
}
}
case t_POL: case t_POLMOD:
return powgi(x,gn);
case t_RFRAC: case t_RFRACN:
{
av=avma; y=cgetg(3,tx); m = labs(n);
y[1]=lpuigs((GEN)x[1],m);
y[2]=lpuigs((GEN)x[2],m);
if (n<0) y=ginv(y); /* let ginv worry about normalizations */
return gerepileupto(av,y);
}
default:
m = labs(n);
av=avma; y=NULL; lim=stack_lim(av,1);
for (; m>1; m>>=1)
{
if (m&1) y = y? gmul(y,x): x;
x=gsqr(x);
if (low_stack(lim, stack_lim(av,1)))
{
GEN *gptr[2]; gptr[0]=&x; gptr[1]=&y;
if(DEBUGMEM>1) err(warnmem,"[3]: gpowgs");
gerepilemany(av,gptr,y? 2: 1);
}
}
y = y? gmul(y,x): x;
if (n<=0) y=ginv(y);
return gerepileupto(av,y);
}
}
GEN
pow_monome(GEN x, GEN nn)
{
long n,m,i,j,dd,tetpil, av = avma;
GEN p1,y;
if (is_bigint(nn)) err(talker,"power overflow in pow_monome");
n = itos(nn); m = labs(n);
j=lgef(x); dd=(j-3)*m+3;
p1=cgetg(dd,t_POL); m = labs(n);
p1[1] = evalsigne(1) | evallgef(dd) | evalvarn(varn(x));
for (i=2; i<dd-1; i++) p1[i]=zero;
p1[i]=lpuigs((GEN)x[j-1],m);
if (n>0) return p1;
tetpil=avma; y=cgetg(3,t_RFRAC);
y[1]=(long)denom((GEN)p1[i]);
y[2]=lmul(p1,(GEN)y[1]); return gerepile(av,tetpil,y);
}
/* n is assumed to be an integer */
GEN
powgi(GEN x, GEN n)
{
long lim,av,i,j,m,tx, sn=signe(n);
GEN y, p1;
if (typ(n) != t_INT) err(talker,"not an integral exponent in powgi");
if (!sn) return puiss0(x);
switch(tx=typ(x))
{/* catch some types here, instead of trying gpowgs() first, because of
* the simpler call interface to puissii() and powmodulo() -- even though
* for integer/rational bases other than 0,+-1 and non-wordsized
* exponents the result will be certain to overflow. --GN
*/
case t_INT:
{
long sx=signe(x), sr = (sx<0 && mod2(n))? -1: 1;
if (sn>0) return puissii(x,n,sr);
if (!sx) err(talker, "division by zero in powgi");
if (is_pm1(x)) return (sr < 0)? icopy(x): gun;
/* n<0, |x|>1 */
y=cgetg(3,t_FRAC); setsigne(n,1); /* temporarily replace n by abs(n) */
y[1]=(sr>0)? un: lnegi(gun);
y[2]=(long)puissii(x,n,1);
setsigne(n,-1); return y;
}
case t_INTMOD:
y=cgetg(3,tx); copyifstack(x[1],y[1]);
y[2]=(long)powmodulo((GEN)x[2],n,(GEN)x[1]);
return y;
case t_FRAC: case t_FRACN:
{
long sr = (mod2(n) && (signe(x[1])!=signe(x[2]))) ? -1 : 1;
if (sn<0)
{
if (!signe((GEN)x[1]))
err(talker, "division by zero fraction in powgi");
setsigne(n,1);
if (is_pm1((GEN)x[1])) /* +-1/x[2] inverts to an integer */
y=puissii((GEN)x[2],n,sr);
else
{
y=cgetg(3,tx);
y[1]=(long)puissii((GEN)x[2],n,sr);
y[2]=(long)puissii((GEN)x[1],n,1);
}
setsigne(n,-1); return y;
}
else /* n > 0 */
{
y=cgetg(3,tx);
if (!signe((GEN)x[1])) { y[1]=zero; y[2]=un; return y; }
y[1]=(long)puissii((GEN)x[1],n,sr);
y[2]=(long)puissii((GEN)x[2],n,1);
return y;
}
}
case t_POL:
if (ismonome(x)) return pow_monome(x,n);
default:
av=avma; lim=stack_lim(av,1);
p1 = n+2; m = *p1;
y=x; j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
for (i=lgefint(n)-2;;)
{
for (; j; m<<=1,j--)
{
y=gsqr(y);
if (low_stack(lim, stack_lim(av,1)))
{
if(DEBUGMEM>1) err(warnmem,"[1]: powgi");
y = gerepileupto(av, y);
}
if (m<0) y=gmul(y,x);
if (low_stack(lim, stack_lim(av,1)))
{
if(DEBUGMEM>1) err(warnmem,"[2]: powgi");
y = gerepileupto(av, y);
}
}
if (--i == 0) break;
m = *++p1, j = BITS_IN_LONG;
}
if (sn<0) y = ginv(y);
return av==avma? gcopy(y): gerepileupto(av,y);
}
}
/* we suppose n != 0, and valp(x) = 0 */
static GEN
ser_pui(GEN x, GEN n, long prec)
{
long av,tetpil,lx,i,j;
GEN p1,p2,y,alp;
if (gvar9(n) > varn(x))
{
GEN lead = (GEN)x[2];
if (gcmp1(lead))
{
av=avma; alp = gclone(gadd(n,gun)); avma=av;
y=cgetg(lx=lg(x),t_SER);
y[1] = evalsigne(1) | evalvalp(0) | evalvarn(varn(x));
y[2] = un;
for (i=3; i<lx; i++)
{
av=avma; p1=gzero;
for (j=1; j<i-1; 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;
}
av=avma; p1=gdiv(x,lead); p1[2] = un; /* in case it's inexact */
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;
return gerepile(av,tetpil,gexp(y,prec));
}
GEN
gpow(GEN x, GEN n, long prec)
{
long av,tetpil,i,lx,tx;
GEN y;
if (typ(n)==t_INT) return powgi(x,n);
tx = typ(x);
if (is_matvec_t(tx))
{
lx=lg(x); y=cgetg(lx,tx);
for (i=1; i<lx; i++) y[i]=lpui((GEN)x[i],n,prec);
return y;
}
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;
if (gcmp0(x))
{
n = greal(n);
if (gsigne(n)<=0) err(talker,"zero to a non positive exponent in gpow");
i = precision(x); if (!i) return gcopy(x);
x=ground(gmulsg(gexpo(x),n));
if (lgefint(x)<=3)
{
i = itos(x);
if (labs(i) < (long)HIGHEXPOBIT)
{
avma=av; y=cgetr(3); y[1]=evalexpo(i); y[2]=0;
return y;
}
}
err(talker,"underflow or overflow in gpow");
}
i = (long) precision(n); if (i) prec=i;
y=gmul(n,glog(x,prec)); tetpil=avma;
return gerepile(av,tetpil,gexp(y,prec));
}
/********************************************************************/
/** **/
/** RACINE CARREE **/
/** **/
/********************************************************************/
GEN
mpsqrt(GEN x)
{
long l,l0,l1,l2,s,eps,n,i,ex,av;
double beta;
GEN y,p1,p2,p3;
if (typ(x)!=t_REAL) err(typeer,"mpsqrt");
s=signe(x); if (s<0) err(talker,"negative argument in mpsqrt");
if (!s)
{
y = cgetr(3);
y[1] = evalexpo(expo(x)>>1);
y[2] = 0; return y;
}
l=lg(x); y=cgetr(l); av=avma;
p1=cgetr(l+1); affrr(x,p1);
ex=expo(x); eps = ex & 1; ex >>= 1;
setexpo(p1,eps); setlg(p1,3);
n=(long)(2+log((double) (l-2))/LOG2);
p2=cgetr(l+1); p2[1]=evalexpo(0) | evalsigne(1);
beta=sqrt((eps+1) * (2 + p1[2]/C31));
p2[2]=(long)((beta-2)*C31);
if (!p2[2]) { p2[2]=HIGHBIT; setexpo(p2,1); }
for (i=3; i<=l; i++) p2[i]=0;
setlg(p2,3);
p3=cgetr(l+1); l-=2; l1=1; l2=3;
for (i=1; i<=n; i++)
{
l0 = l1<<1;
if (l0<=l)
{ l2 += l1; l1=l0; }
else
{ l2 += l+1-l1; l1=l+1; }
setlg(p3,l1+2); setlg(p1,l1+2); setlg(p2,l2);
divrrz(p1,p2,p3); addrrz(p2,p3,p2); setexpo(p2,expo(p2)-1);
}
affrr(p2,y); setexpo(y,expo(y)+ex);
avma=av; return y;
}
static GEN
padic_sqrt(GEN x)
{
long av = avma, av2,lim,e,r,lpp,lp,pp;
GEN y;
e=valp(x); y=cgetg(5,t_PADIC); copyifstack(x[2],y[2]);
if (gcmp0(x))
{
y[4] = zero; e = (e+1)>>1;
y[3] = lpuigs((GEN)x[2],e);
y[1] = evalvalp(e) | evalprecp(precp(x));
return y;
}
if (e & 1) err(sqrter6);
e>>=1; setvalp(y,e); y[3] = y[2];
pp = precp(x);
if (egalii(gdeux, (GEN)x[2]))
{
lp=3; y[4] = un; r = mod8((GEN)x[4]);
if ((r!=1 && pp>=2) && (r!=5 || pp!=2)) err(sqrter5);
if (pp <= lp) { setprecp(y,1); return y; }
x = dummycopy(x); setvalp(x,0); setvalp(y,0);
av2=avma; lim = stack_lim(av2,2);
y[3] = lstoi(8);
for(;;)
{
lpp=lp; lp=(lp<<1)-1;
if (lp < pp) y[3] = lshifti((GEN)y[3], lpp-1);
else
{ lp=pp; y[3] = x[3]; }
setprecp(y,lp); y = gdiv(gadd(y, gdiv(x,y)), gdeux);
if (lp < pp) lp--;
if (cmpii((GEN)y[4], (GEN)y[3]) >= 0)
y[4] = lsubii((GEN)y[4], (GEN)y[3]);
if (pp <= lp) break;
if (low_stack(lim,stack_lim(av2,2)))
{
if (DEBUGMEM > 1) err(warnmem,"padic_sqrt");
y = gerepileupto(av2,y);
}
}
y = gcopy(y);
}
else /* p != 2 */
{
lp=1; y[4] = (long)mpsqrtmod((GEN)x[4],(GEN)x[2]);
if (!y[4]) err(sqrter5);
if (pp <= lp) { setprecp(y,1); return y; }
x = dummycopy(x); setvalp(x,0); setvalp(y,0);
av2=avma; lim = stack_lim(av2,2);
for(;;)
{
lp<<=1;
if (lp < pp) y[3] = lsqri((GEN)y[3]);
else
{ lp=pp; y[3] = x[3]; }
setprecp(y,lp); y = gdiv(gadd(y, gdiv(x,y)), gdeux);
if (pp <= lp) break;
if (low_stack(lim,stack_lim(av2,2)))
{
if (DEBUGMEM > 1) err(warnmem,"padic_sqrt");
y = gerepileupto(av2,y);
}
}
}
setvalp(y,e); return gerepileupto(av,y);
}
GEN
gsqrt(GEN x, long prec)
{
long av,tetpil,e;
GEN y,p1,p2;
switch(typ(x))
{
case t_REAL:
if (signe(x)>=0) return mpsqrt(x);
y=cgetg(3,t_COMPLEX); y[1]=zero;
setsigne(x,1); y[2]=lmpsqrt(x); setsigne(x,-1);
return y;
case t_INTMOD:
y=cgetg(3,t_INTMOD); copyifstack(x[1],y[1]);
p1 = mpsqrtmod((GEN)x[2],(GEN)y[1]);
if (!p1) err(sqrter5);
y[2] = (long)p1; return y;
case t_COMPLEX:
y=cgetg(3,t_COMPLEX); av=avma;
if (gcmp0((GEN)x[2]))
{
long tx=typ(x[1]);
if ((is_intreal_t(tx) || is_frac_t(tx)) && gsigne((GEN)x[1]) < 0)
{
y[1]=zero; p1=gneg_i((GEN)x[1]); tetpil=avma;
y[2]=lpile(av,tetpil,gsqrt(p1,prec));
return y;
}
y[1]=lsqrt((GEN)x[1],prec);
y[2]=zero; return y;
}
p1=gsqr((GEN)x[1]); p2=gsqr((GEN)x[2]);
p1=gsqrt(gadd(p1,p2),prec);
if (gcmp0(p1))
{
y[1]=lsqrt(p1,prec);
y[2]=lcopy((GEN)y[1]);
return y;
}
if (gsigne((GEN)x[1])<0)
{
p2=gsub(p1,(GEN)x[1]); p1=gmul2n(p2,-1);
y[2]=lsqrt(p1,prec); y[1]=ldiv((GEN)x[2],gmul2n((GEN)y[2],1));
tetpil=avma;
y = (gsigne((GEN)x[2]) > 0)? gcopy(y): gneg(y);
return gerepile(av,tetpil,y);
}
p1=gmul2n(gadd(p1,(GEN)x[1]),-1);
tetpil=avma; y[1]=lpile(av,tetpil,gsqrt(p1,prec));
av=avma; p1=gmul2n((GEN)y[1],1); tetpil=avma;
y[2]=lpile(av,tetpil,gdiv((GEN)x[2],p1));
return y;
case t_PADIC:
return padic_sqrt(x);
case t_SER:
e=valp(x);
if (gcmp0(x)) return zeroser(varn(x), (e+1)>>1);
if (e & 1) err(sqrter6);
/* trick: ser_pui assumes valp(x) = 0 */
y = ser_pui(x,ghalf,prec);
y[1] = evalsigne(1) | evalvalp(e>>1) | evalvarn(varn(x));
return y;
}
return transc(gsqrt,x,prec);
}
void
gsqrtz(GEN x, GEN y)
{
long av=avma, prec = precision(y);
if (!prec) err(infprecer,"gsqrtz");
gaffect(gsqrt(x,prec),y); avma=av;
}
/********************************************************************/
/** **/
/** FONCTION EXPONENTIELLE-1 **/
/** **/
/********************************************************************/
#ifdef LONG_IS_64BIT
# define EXMAX 46
#else
# define EXMAX 22
#endif
GEN
mpexp1(GEN x)
{
long l,l1,l2,i,n,m,ex,s,av,av2, sx = signe(x);
double a,b,alpha,beta, gama = 2.0 /* optimized for SUN3 */;
/* KB: 3.0 is better for UltraSparc */
GEN y,p1,p2,p3,p4,unr;
if (typ(x)!=t_REAL) err(typeer,"mpexp1");
if (!sx)
{
y=cgetr(3); y[1]=x[1]; y[2]=0; return y;
}
l=lg(x); y=cgetr(l); av=avma; /* room for result */
l2 = l+1; ex = expo(x);
if (ex > EXMAX) err(talker,"exponent too large in exp");
alpha = -1-log(2+x[2]/C31)-ex*LOG2;
beta = 5 + bit_accuracy(l)*LOG2;
a = sqrt(beta/(gama*LOG2));
b = (alpha + 0.5*log(beta*gama/LOG2))/LOG2;
if (a>=b)
{
n=(long)(1+sqrt(beta*gama/LOG2));
m=(long)(1+a-b);
l2 += m>>TWOPOTBITS_IN_LONG;
}
else
{
n=(long)(1+beta/alpha);
m=0;
}
unr=realun(l2);
p2=rcopy(unr); setlg(p2,4);
p4=cgetr(l2); affrr(x,p4); setsigne(p4,1);
if (m) setexpo(p4,ex-m);
s=0; l1=4; av2 = avma;
for (i=n; i>=2; i--)
{
setlg(p4,l1); p3 = divrs(p4,i);
s -= expo(p3); p1 = mulrr(p3,p2); setlg(p1,l1);
l1 += s>>TWOPOTBITS_IN_LONG; if (l1>l2) l1=l2;
s %= BITS_IN_LONG;
setlg(unr,l1); p1 = addrr(unr,p1);
setlg(p2,l1); affrr(p1,p2); avma = av2;
}
setlg(p2,l2); setlg(p4,l2); p2 = mulrr(p4,p2);
for (i=1; i<=m; i++)
{
setlg(p2,l2);
p2 = mulrr(addsr(2,p2),p2);
}
if (sx == -1)
{
setlg(unr,l2); p2 = addrr(unr,p2);
setlg(p2,l2); p2 = ginv(p2);
p2 = subrr(p2,unr);
}
affrr(p2,y); avma = av; return y;
}
#undef EXMAX
/********************************************************************/
/** **/
/** FONCTION EXPONENTIELLE **/
/** **/
/********************************************************************/
GEN
mpexp(GEN x)
{
long av, sx = signe(x);
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); }
return gerepileupto(av,y);
}
static GEN
paexp(GEN x)
{
long k,av, e = valp(x), pp = precp(x), n = e + pp;
GEN y,r,p1;
if (gcmp0(x)) return gaddgs(x,1);
if (e<=0 || (!cmpis((GEN)x[2],2) && e==1))
err(talker,"p-adic argument out of range in gexp");
av = avma;
if (egalii(gdeux,(GEN)x[2]))
{
n--; e--; k = n/e;
if (n%e==0) k--;
}
else
{
p1 = subis((GEN)x[2], 1);
k = itos(dvmdii(subis(mulis(p1,n), 1), subis(mulis(p1,e), 1), &r));
if (!signe(r)) k--;
}
for (y=gun; k; k--) y = gaddsg(1, gdivgs(gmul(y,x), k));
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:
return mpexp(x);
case t_COMPLEX:
y=cgetg(3,t_COMPLEX); av=avma;
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);
}
void
gexpz(GEN x, GEN y)
{
long av=avma, prec = precision(y);
if (!prec) err(infprecer,"gexpz");
gaffect(gexp(x,prec),y); avma=av;
}
/********************************************************************/
/** **/
/** FONCTION LOGARITHME **/
/** **/
/********************************************************************/
GEN
mplog(GEN x)
{
long l,l1,l2,m,m1,n,i,ex,s,sgn,ltop,av;
double alpha,beta,a,b;
GEN y,p1,p2,p3,p4,p5,unr;
if (typ(x)!=t_REAL) err(typeer,"mplog");
if (signe(x)<=0) err(talker,"non positive argument in mplog");
l=lg(x); sgn = cmpsr(1,x); if (!sgn) return realzero(l);
y=cgetr(l); ltop=avma;
l2 = l+1; p2=p1=cgetr(l2); affrr(x,p1);
av=avma;
if (sgn > 0) p2 = divsr(1,p1); /* x<1 changer x en 1/x */
for (m1=1; expo(p2)>0; m1++) p2 = mpsqrt(p2);
if (m1 > 1 || sgn > 0) { affrr(p2,p1); avma=av; }
alpha = 1+p1[2]/C31; if (!alpha) alpha = 0.00000001;
l -= 2; alpha = -log(alpha);
beta = (BITS_IN_LONG/2)*l*LOG2;
a = alpha/LOG2;
b = sqrt((BITS_IN_LONG/2)*l/3.0);
if (a<=b)
{
n=(long)(1+sqrt((BITS_IN_LONG/2)*3.0*l));
m=(long)(1+b-a);
l2 += m>>TWOPOTBITS_IN_LONG;
p4=cgetr(l2); affrr(p1,p4);
p1 = p4; av=avma;
for (i=1; i<=m; i++) p1 = mpsqrt(p1);
affrr(p1,p4); avma=av;
}
else
{
n=(long)(1+beta/alpha);
m=0; p4 = p1;
}
unr=realun(l2);
p2=cgetr(l2); p3=cgetr(l2); av=avma;
/* affrr needed here instead of setlg since prec may DECREASE */
p1 = cgetr(l2); affrr(subrr(p4,unr), p1);
p5 = addrr(p4,unr); setlg(p5,l2);
affrr(divrr(p1,p5), p2);
affrr(mulrr(p2,p2), p3);
affrr(divrs(unr,2*n+1), p4); setlg(p4,4); avma=av;
s=0; ex=expo(p3); l1 = 4;
for (i=n; i>=1; i--)
{
setlg(p3,l1); p5 = mulrr(p4,p3);
setlg(unr,l1); p1 = divrs(unr,2*i-1);
s -= ex;
l1 += s>>TWOPOTBITS_IN_LONG; if (l1>l2) l1=l2;
s %= BITS_IN_LONG;
setlg(p1,l1); setlg(p4,l1); setlg(p5,l1);
p1 = addrr(p1,p5); affrr(p1,p4); avma=av;
}
setlg(p4,l2); affrr(mulrr(p2,p4), y);
setexpo(y, expo(y)+m+m1);
if (sgn > 0) setsigne(y, -1);
avma=ltop; return y;
}
GEN
teich(GEN x)
{
GEN aux,y,z,p1;
long av,n,k;
if (typ(x)!=t_PADIC) err(talker,"not a p-adic argument in teichmuller");
if (!signe(x[4])) return gcopy(x);
y=cgetp(x); z=(GEN)x[4];
if (!cmpis((GEN)x[2],2))
{
n = mod4(z);
if (n & 2)
addsiz(-1,(GEN)x[3],(GEN)y[4]);
else
affsi(1,(GEN)y[4]);
return y;
}
av = avma; p1 = addsi(-1,(GEN)x[2]);
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;
}
static GEN
palogaux(GEN x)
{
long av1=avma,tetpil,k,e,pp;
GEN y,s,y2;
if (egalii(gun,(GEN)x[4]))
{
y=gaddgs(x,-1);
if (egalii(gdeux,(GEN)x[2]))
{
setvalp(y,valp(y)-1);
if (!gcmp1((GEN)y[3])) y[3]=lshifti((GEN)y[3],-1);
}
tetpil=avma; return gerepile(av1,tetpil,gcopy(y));
}
y=gdiv(gaddgs(x,-1),gaddgs(x,1)); e=valp(y); pp=e+precp(y);
if (egalii(gdeux,(GEN)x[2])) pp--;
else
{
long av=avma;
GEN p1;
for (p1=stoi(e); cmpsi(pp,p1)>0; pp++)
p1 = mulii(p1,(GEN)x[2]);
avma=av; pp-=2;
}
k=pp/e; if (!odd(k)) k--;
y2=gsqr(y); s=gdivgs(gun,k);
while (k>=3)
{
k-=2; s=gadd(gmul(y2,s),gdivgs(gun,k));
}
tetpil=avma; return gerepile(av1,tetpil,gmul(s,y));
}
GEN
palog(GEN x)
{
long av=avma,tetpil;
GEN p1,y;
if (!signe(x[4])) err(talker,"zero argument in palog");
if (cmpis((GEN)x[2],2))
{
y=cgetp(x); p1=gsubgs((GEN)x[2],1);
affii(powmodulo((GEN)x[4],p1,(GEN)x[3]),(GEN)y[4]);
y=gmulgs(palogaux(y),2); tetpil=avma;
return gerepile(av,tetpil,gdiv(y,p1));
}
y=gsqr(x); setvalp(y,0); tetpil=avma;
return gerepile(av,tetpil,palogaux(y));
}
GEN
log0(GEN x, long flag,long prec)
{
switch(flag)
{
case 0: return glog(x,prec);
case 1: return glogagm(x,prec);
default: err(flagerr,"log");
}
return NULL; /* not reached */
}
GEN
glog(GEN x, long prec)
{
long av,tetpil;
GEN y,p1,p2;
switch(typ(x))
{
case t_REAL:
if (signe(x)>=0) return mplog(x);
y=cgetg(3,t_COMPLEX); y[2]=lmppi(lg(x));
setsigne(x,1); y[1]=lmplog(x);
setsigne(x,-1); return y;
case t_COMPLEX:
y=cgetg(3,t_COMPLEX); y[2]=larg(x,prec);
av=avma; p1=glog(gnorm(x),prec); tetpil=avma;
y[1]=lpile(av,tetpil,gmul2n(p1,-1));
return y;
case t_PADIC:
return palog(x);
case t_SER:
av=avma; if (valp(x)) err(negexper,"glog");
p1=gdiv(derivser(x),x);
tetpil=avma; p1=integ(p1,varn(x));
if (gcmp1((GEN)x[2])) return gerepile(av,tetpil,p1);
p2=glog((GEN)x[2],prec); tetpil=avma;
return gerepile(av,tetpil,gadd(p1,p2));
case t_INTMOD:
err(typeer,"glog");
}
return transc(glog,x,prec);
}
void
glogz(GEN x, GEN y)
{
long av=avma, prec = precision(y);
if (!prec) err(infprecer,"glogz");
gaffect(glog(x,prec),y); avma=av;
}
/********************************************************************/
/** **/
/** FONCTION SINCOS-1 **/
/** **/
/********************************************************************/
static GEN
mpsc1(GEN x, long *ptmod8)
{
const long mmax = 23169;
/* on a 64-bit machine with true 128 bit/64 bit division, one could
* take mmax=1518500248; on the alpha it does not seem worthwhile
*/
long l,l0,l1,l2,l4,ee,i,n,m,s,t,av;
double alpha,beta,a,b,c,d;
GEN y,p1,p2,p3,p4,pitemp;
if (typ(x)!=t_REAL) err(typeer,"mpsc1");
if (!signe(x))
{
y=cgetr(3);
y[1] = evalexpo((expo(x)<<1) - 1);
y[2] = 0; *ptmod8=0; return y;
}
l=lg(x); y=cgetr(l); av=avma;
l++; pitemp = mppi(l+1); setexpo(pitemp,-1);
p1 = addrr(x,pitemp); setexpo(pitemp,0);
if (expo(p1) >= bit_accuracy(min(l,lg(p1))) + 3)
err(talker,"loss of precision in mpsc1");
p3 = divrr(p1,pitemp); p2 = mpent(p3);
if (signe(p2)) x = subrr(x, mulir(p2,pitemp));
p1 = cgetr(l+1); affrr(x, p1);
*ptmod8 = (signe(p1) < 0)? 4: 0;
if (signe(p2))
{
long mod4 = mod4(p2);
if (signe(p2) < 0 && mod4) mod4 = 4-mod4;
*ptmod8 += mod4;
}
if (gcmp0(p1)) alpha=1000000.0;
else
{
m=expo(p1); alpha=(m<-1023)? -1-m*LOG2: -1-log(fabs(rtodbl(p1)));
}
beta = 5 + bit_accuracy(l)*LOG2;
a=0.5/LOG2; b=0.5*a;
c = a+sqrt((beta+b)/LOG2);
d = ((beta/c)-alpha-log(c))/LOG2;
if (d>=0)
{
m=(long)(1+d);
n=(long)((1+c)/2.0);
setexpo(p1,expo(p1)-m);
}
else { m=0; n=(long)((1+beta/alpha)/2.0); }
l2=l+1+(m>>TWOPOTBITS_IN_LONG);
p2=realun(l2); setlg(p2,4);
p4=cgetr(l2); av = avma;
affrr(gsqr(p1),p4); setlg(p4,4);
if (n>mmax)
p3 = divrs(divrs(p4,2*n+2),2*n+1);
else
p3 = divrs(p4, (2*n+2)*(2*n+1));
ee = -expo(p3); s=0;
l4 = l1 = 3 + (ee>>TWOPOTBITS_IN_LONG);
if (l4<=l2) { setlg(p2,l4); setlg(p4,l4); }
for (i=n; i>mmax; i--)
{
p3 = divrs(divrs(p4,2*i),2*i-1); s -= expo(p3);
t=s&(BITS_IN_LONG-1); l0=(s>>TWOPOTBITS_IN_LONG);
if (t) l0++;
l1 += l0; if (l1>l2) { l0 += l2-l1; l1=l2; }
l4 += l0;
p3 = mulrr(p3,p2);
if (l4<=l2) { setlg(p2,l4); setlg(p4,l4); }
subsrz(1,p3,p2); avma=av;
}
for ( ; i>=2; i--)
{
p3 = divrs(p4, 2*i*(2*i-1)); s -= expo(p3);
t=s&(BITS_IN_LONG-1); l0=(s>>TWOPOTBITS_IN_LONG);
if (t) l0++;
l1 += l0; if (l1>l2) { l0 += l2-l1; l1=l2; }
l4 += l0;
p3 = mulrr(p3,p2);
if (l4<=l2) { setlg(p2,l4); setlg(p4,l4); }
subsrz(1,p3,p2); avma=av;
}
if (l4<=l2) { setlg(p2,l4); setlg(p4,l4); }
setexpo(p4,expo(p4)-1); setsigne(p4, -signe(p4));
p2 = mulrr(p4,p2);
for (i=1; i<=m; i++)
{
p2 = mulrr(p2,addsr(2,p2)); setexpo(p2,expo(p2)+1);
}
affrr(p2,y); avma=av; return y;
}
/********************************************************************/
/** **/
/** FONCTION sqrt(-x*(x+2)) **/
/** **/
/********************************************************************/
static GEN
mpaut(GEN x)
{
long av = avma;
GEN p1;
cgetr(2*lg(x)+3); /* bound for 2*lg(p1)+1 */
p1=mulrr(x,addsr(2,x)); setsigne(p1,-signe(p1));
avma = av; return mpsqrt(p1); /* safe ! */
}
/********************************************************************/
/** **/
/** FONCTION COSINUS **/
/** **/
/********************************************************************/
static GEN
mpcos(GEN x)
{
long mod8,av,tetpil;
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;
switch(mod8)
{
case 0: case 4:
y=addsr(1,p1); break;
case 1: case 7:
y=mpaut(p1); setsigne(y,-signe(y)); break;
case 2: case 6:
y=subsr(-1,p1); break;
case 3: case 5:
y=mpaut(p1); break;
default:;
}
return gerepile(av,tetpil,y);
}
GEN
gcos(GEN x, long prec)
{
long av,tetpil;
GEN r,u,v,y,p1,p2;
switch(typ(x))
{
case t_REAL:
return mpcos(x);
case t_COMPLEX:
y=cgetg(3,t_COMPLEX); av=avma;
r=gexp((GEN)x[2],prec); p1=ginv(r);
p2=gmul2n(mpadd(p1,r),-1);
p1=mpsub(p2,r);
gsincos((GEN)x[1],&u,&v,prec);
tetpil=avma;
y[1]=lmul(p2,v); y[2]=lmul(p1,u);
gerepilemanyvec(av,tetpil,y+1,2);
return y;
case t_SER:
if (gcmp0(x)) return gaddsg(1,x);
if (valp(x)<0) err(negexper,"gcos");
av=avma; gsincos(x,&u,&v,prec); tetpil=avma;
return gerepile(av,tetpil,gcopy(v));
case t_INTMOD: case t_PADIC:
err(typeer,"gcos");
}
return transc(gcos,x,prec);
}
void
gcosz(GEN x, GEN y)
{
long av = avma, prec = precision(y);
if (!prec) err(infprecer,"gcosz");
gaffect(gcos(x,prec),y); avma=av;
}
/********************************************************************/
/** **/
/** FONCTION SINUS **/
/** **/
/********************************************************************/
GEN
mpsin(GEN x)
{
long mod8,av,tetpil;
GEN y,p1;
if (typ(x)!=t_REAL) err(typeer,"mpsin");
if (!signe(x))
{
y=cgetr(3); y[1]=x[1]; y[2]=0;
return y;
}
av=avma; p1=mpsc1(x,&mod8); tetpil=avma;
switch(mod8)
{
case 0: case 6:
y=mpaut(p1); break;
case 1: case 5:
y=addsr(1,p1); break;
case 2: case 4:
y=mpaut(p1); setsigne(y,-signe(y)); break;
case 3: case 7:
y=subsr(-1,p1); break;
default:;
}
return gerepile(av,tetpil,y);
}
GEN
gsin(GEN x, long prec)
{
long av,tetpil;
GEN r,u,v,y,p1,p2;
switch(typ(x))
{
case t_REAL:
return mpsin(x);
case t_COMPLEX:
y=cgetg(3,t_COMPLEX); av=avma;
r=gexp((GEN)x[2],prec); p1=ginv(r);
p2=gmul2n(mpadd(p1,r),-1);
p1=mpsub(p2,p1);
gsincos((GEN)x[1],&u,&v,prec);
tetpil=avma;
y[1]=lmul(p2,u); y[2]=lmul(p1,v);
gerepilemanyvec(av,tetpil,y+1,2);
return y;
case t_SER:
if (gcmp0(x)) return gcopy(x);
if (valp(x)<0) err(negexper,"gsin");
av=avma; gsincos(x,&u,&v,prec); tetpil=avma;
return gerepile(av,tetpil,gcopy(u));
case t_INTMOD: case t_PADIC:
err(typeer,"gsin");
}
return transc(gsin,x,prec);
}
void
gsinz(GEN x, GEN y)
{
long av=avma, prec = precision(y);
if (!prec) err(infprecer,"gsinz");
gaffect(gsin(x,prec),y); avma=av;
}
/********************************************************************/
/** **/
/** PROCEDURE SINUS,COSINUS **/
/** **/
/********************************************************************/
void
mpsincos(GEN x, GEN *s, GEN *c)
{
long av,tetpil,mod8;
GEN p1, *gptr[2];
if (typ(x)!=t_REAL) err(typeer,"mpsincos");
if (!signe(x))
{
p1=cgetr(3); *s=p1; p1[1]=x[1];
p1[2]=0; *c=addsr(1,x); return;
}
av=avma; p1=mpsc1(x,&mod8); tetpil=avma;
switch(mod8)
{
case 0: *c=addsr( 1,p1); *s=mpaut(p1); break;
case 1: *s=addsr( 1,p1); *c=mpaut(p1); setsigne(*c,-signe(*c)); break;
case 2: *c=subsr(-1,p1); *s=mpaut(p1); setsigne(*s,-signe(*s)); break;
case 3: *s=subsr(-1,p1); *c=mpaut(p1); break;
case 4: *c=addsr( 1,p1); *s=mpaut(p1); setsigne(*s,-signe(*s)); break;
case 5: *s=addsr( 1,p1); *c=mpaut(p1); break;
case 6: *c=subsr(-1,p1); *s=mpaut(p1); break;
case 7: *s=subsr(-1,p1); *c=mpaut(p1); setsigne(*c,-signe(*c)); break;
}
gptr[0]=s; gptr[1]=c;
gerepilemanysp(av,tetpil,gptr,2);
}
void
gsincos(GEN x, GEN *s, GEN *c, long prec)
{
long av,tetpil,ii,i,j,ex,ex2,lx,ly;
GEN r,u,v,u1,v1,p1,p2,p3,p4,ps,pc;
GEN *gptr[4];
switch(typ(x))
{
case t_INT: case t_FRAC: case t_FRACN:
av=avma; p1=cgetr(prec); gaffect(x,p1); tetpil=avma;
mpsincos(p1,s,c); gptr[0]=s; gptr[1]=c;
gerepilemanysp(av,tetpil,gptr,2);
return;
case t_REAL:
mpsincos(x,s,c); return;
case t_COMPLEX:
ps=cgetg(3,t_COMPLEX); pc=cgetg(3,t_COMPLEX);
*s=ps; *c=pc; av=avma;
r=gexp((GEN)x[2],prec); p1=ginv(r);
p2=gmul2n(mpadd(p1,r),-1);
p1=mpsub(p2,p1); r=mpsub(p2,r);
gsincos((GEN)x[1],&u,&v,prec);
tetpil=avma;
p1=gmul(p2,u); p2=gmul(p1,v);
p3=gmul(p2,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;
return;
case t_SER:
if (gcmp0(x)) { *c=gaddsg(1,x); *s=gcopy(x); return; }
ex=valp(x); lx=lg(x); ex2=2*ex+2;
if (ex < 0) err(talker,"non zero exponent in gsincos");
if (ex2 > lx)
{
*s=gcopy(x); av=avma; p1=gdivgs(gsqr(x),2);
tetpil=avma; *c=gerepile(av,tetpil,gsubsg(1,p1));
return;
}
if (!ex)
{
av=avma; p1=gcopy(x); p1[2]=zero;
gsincos(normalize(p1),&u,&v,prec);
gsincos((GEN)x[2],&u1,&v1,prec);
p1=gmul(v1,v); p2=gmul(u1,u); p3=gmul(v1,u);
p4=gmul(u1,v); tetpil=avma;
*c=gsub(p1,p2); *s=gadd(p3,p4);
gptr[0]=s; gptr[1]=c;
gerepilemanysp(av,tetpil,gptr,2);
return;
}
ly=lx+2*ex;
pc=cgetg(ly,t_SER); *c=pc;
ps=cgetg(lx,t_SER); *s=ps;
pc[1] = evalsigne(1) | evalvalp(0) | evalvarn(varn(x));
pc[2]=un; ps[1]=x[1];
for (i=2; i<ex+2; i++) ps[i]=lcopy((GEN)x[i]);
for (i=3; i<ex2; i++) pc[i]=zero;
for (i=ex2; i<ly; i++)
{
ii=i-ex; av=avma; p1=gzero;
for (j=ex; j<ii-1; 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++)
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]));
}
}
return;
}
err(typeer,"gsincos");
}
/********************************************************************/
/** **/
/** FONCTIONS TANGENTE ET COTANGENTE **/
/** **/
/********************************************************************/
static GEN
mptan(GEN x)
{
long av=avma,tetpil;
GEN s,c;
mpsincos(x,&s,&c); tetpil=avma;
return gerepile(av,tetpil,divrr(s,c));
}
GEN
gtan(GEN x, long prec)
{
long av,tetpil;
GEN s,c;
switch(typ(x))
{
case t_REAL:
return mptan(x);
case t_SER:
if (gcmp0(x)) return gcopy(x);
if (valp(x)<0) err(negexper,"gtan"); /* fall through */
case t_COMPLEX:
av=avma; gsincos(x,&s,&c,prec); tetpil=avma;
return gerepile(av,tetpil,gdiv(s,c));
case t_INTMOD: case t_PADIC:
err(typeer,"gtan");
}
return transc(gtan,x,prec);
}
void
gtanz(GEN x, GEN y)
{
long av=avma, prec = precision(y);
if (!prec) err(infprecer,"gtanz");
gaffect(gtan(x,prec),y); avma=av;
}
static GEN
mpcotan(GEN x)
{
long av=avma,tetpil;
GEN s,c;
mpsincos(x,&s,&c); tetpil=avma;
return gerepile(av,tetpil,divrr(c,s));
}
GEN
gcotan(GEN x, long prec)
{
long av,tetpil;
GEN s,c;
switch(typ(x))
{
case t_REAL:
return mpcotan(x);
case t_SER: err(negexper,"gcotan"); /* fall through */
case t_COMPLEX:
av=avma; gsincos(x,&s,&c,prec); tetpil=avma;
return gerepile(av,tetpil,gdiv(c,s));
case t_INTMOD: case t_PADIC:
err(typeer,"gcotan");
}
return transc(gcotan,x,prec);
}
void
gcotanz(GEN x, GEN y)
{
long av=avma, prec = precision(y);
if (!prec) err(infprecer,"gcotanz");
gaffect(gcotan(x,prec),y); avma=av;
}
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to trans1.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari / src / basemath