File: [local] / OpenXM_contrib / pari / src / basemath / Attic / gen3.c (download)
Revision 1.1.1.1 (vendor branch), Sun Jan 9 17:35:31 2000 UTC (24 years, 5 months ago) by maekawa
Branch: PARI_GP
CVS Tags: maekawa-ipv6, VERSION_2_0_17_BETA, RELEASE_20000124, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, RELEASE_1_1_3, RELEASE_1_1_2
Changes since 1.1: +0 -0
lines
Import PARI/GP 2.0.17 beta.
|
/********************************************************************/
/** **/
/** GENERIC OPERATIONS **/
/** (third part) **/
/** **/
/********************************************************************/
/* $Id: gen3.c,v 1.2 1999/09/23 17:50:56 karim Exp $ */
#include "pari.h"
/********************************************************************/
/** **/
/** PRINCIPAL VARIABLE NUMBER **/
/** **/
/********************************************************************/
int
gvar(GEN x)
{
long tx=typ(x),i,v,w;
if (is_polser_t(tx)) return varn(x);
if (tx==t_POLMOD) return varn(x[1]);
if (is_const_t(tx) || is_qf_t(tx) || tx == t_STR || tx == t_LIST)
return BIGINT;
v=BIGINT;
for (i=1; i < lg(x); i++) { w=gvar((GEN)x[i]); if (w<v) v=w; }
return v;
}
int
gvar2(GEN x)
{
long tx=typ(x),i,v,w;
if (is_const_t(tx) || is_qf_t(tx)) return BIGINT;
v = BIGINT;
switch(tx)
{
case t_POLMOD:
v=gvar2((GEN)x[1]);
w=gvar2((GEN)x[2]); if (w<v) v=w;
return v;
case t_SER:
for (i=2; i < lg(x); i++) { w=gvar((GEN)x[i]); if (w<v) v=w; }
return v;
case t_POL:
for (i=2; i<lgef(x); i++) { w=gvar((GEN)x[i]); if (w<v) v=w; }
return v;
}
for (i=1; i<lg(x); i++) { w=gvar2((GEN)x[i]); if (w<v) v=w; }
return v;
}
GEN
gpolvar(GEN x)
{
if (typ(x)==t_PADIC) x = (GEN)x[2];
else
{
long v=gvar(x);
if (v==BIGINT) err(typeer,"polvar");
x = polx[v];
}
return gcopy(x);
}
/*******************************************************************/
/* */
/* PRECISION OF SCALAR OBJECTS */
/* */
/*******************************************************************/
long
precision(GEN x)
{
long tx=typ(x),k,l;
if (tx==t_REAL)
{
k=2-(expo(x)>>TWOPOTBITS_IN_LONG);
l=lg(x); if (l>k) k=l;
return k;
}
if (tx==t_COMPLEX)
{
k=precision((GEN)x[1]);
l=precision((GEN)x[2]); if (l && l<k) k=l;
return k;
}
return 0;
}
long
gprecision(GEN x)
{
long tx=typ(x),lx=lg(x),i,k,l;
if (is_scalar_t(tx)) return precision(x);
switch(tx)
{
case t_POL: lx=lgef(x);
case t_VEC: case t_COL: case t_MAT:
k=VERYBIGINT;
for (i=lontyp[tx]; i<lx; i++)
{
l = gprecision((GEN)x[i]);
if (l && l<k) k = l;
}
return (k==VERYBIGINT)? 0: k;
case t_RFRAC: case t_RFRACN:
{
k=gprecision((GEN)x[1]);
l=gprecision((GEN)x[2]); if (l && l<k) k=l;
return k;
}
case t_QFR:
return gprecision((GEN)x[4]);
}
return 0;
}
GEN
ggprecision(GEN x)
{
long a=gprecision(x);
return stoi(a ? (long) ((a-2)*pariK): VERYBIGINT);
}
GEN
precision0(GEN x, long n)
{
if (n) return gprec(x,n);
return ggprecision(x);
}
/* attention: precision p-adique absolue */
long
padicprec(GEN x, GEN p)
{
long i,s,t,lx=lg(x),tx=typ(x);
switch(tx)
{
case t_INT: case t_FRAC: case t_FRACN:
return VERYBIGINT;
case t_INTMOD:
return ggval((GEN)x[1],p);
case t_PADIC:
if (!gegal((GEN)x[2],p))
err(talker,"not the same prime in padicprec");
return precp(x)+valp(x);
case t_POL:
lx=lgef(x);
case t_COMPLEX: case t_QUAD: case t_POLMOD: case t_SER: case t_RFRAC:
case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
for (s=VERYBIGINT, i=lontyp[tx]; i<lx; i++)
{
t = padicprec((GEN)x[i],p); if (t<s) s = t;
}
return s;
}
err(typeer,"padicprec");
return 0; /* not reached */
}
/* Degre de x par rapport a la variable v si v>=0, par rapport a la variable
* principale si v<0. On suppose x fraction rationnelle ou polynome.
* Convention deg(0)=-1.
*/
long
poldegree(GEN x, long v)
{
long tx=typ(x), av, w, d;
if (is_scalar_t(tx)) return gcmp0(x)? -1: 0;
switch(tx)
{
case t_POL:
w = varn(x);
if (v < 0 || v == w) return lgef(x)-3;
if (v < w) return signe(x)? 0: -1;
av = avma; x = gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
if (gvar(x)) { d = gcmp0(x)? -1: 0; } else d = lgef(x)-3;
avma = av; return d;
case t_RFRAC: case t_RFRACN:
if (gcmp0((GEN)x[1])) return -1;
return poldegree((GEN)x[1],v) - poldegree((GEN)x[2],v);
}
err(typeer,"degree");
return 0; /* not reached */
}
long
degree(GEN x)
{
return poldegree(x,-1);
}
/* si v<0, par rapport a la variable principale, sinon par rapport a v.
* On suppose que x est un polynome ou une serie.
*/
GEN
pollead(GEN x, long v)
{
long l,tx = typ(x),av,tetpil,w;
GEN xinit;
if (is_scalar_t(tx)) return gcopy(x);
w = varn(x);
switch(tx)
{
case t_POL:
if (v < 0 || v == w)
{
l=lgef(x);
return (l==2)? gzero: gcopy((GEN)x[l-1]);
}
if (v < w) return gcopy(x);
av = avma; xinit = x;
x = gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
if (gvar(x)) { avma = av; return gcopy(xinit); }
l = lgef(x); if (l == 2) { avma = av; return gzero; }
tetpil = avma; x = gsubst((GEN)x[l-1],MAXVARN,polx[w]);
return gerepile(av,tetpil,x);
case t_SER:
if (v < 0 || v == w) return (!signe(x))? gzero: gcopy((GEN)x[2]);
if (v < w) return gcopy(x);
av = avma; xinit = x;
x = gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
if (gvar(x)) { avma = av; return gcopy(xinit);}
if (!signe(x)) { avma = av; return gzero;}
tetpil = avma; x = gsubst((GEN)x[2],MAXVARN,polx[w]);
return gerepile(av,tetpil,x);
}
err(typeer,"pollead");
return NULL; /* not reached */
}
/* returns 1 if there's a real component in the structure, 0 otherwise */
int
isinexactreal(GEN x)
{
long tx=typ(x),lx,i;
if (is_scalar_t(tx))
{
switch(tx)
{
case t_REAL:
return 1;
case t_COMPLEX: case t_QUAD:
return (typ(x[1])==t_REAL || typ(x[2])==t_REAL);
}
return 0;
}
switch(tx)
{
case t_QFR: case t_QFI:
return 0;
case t_RFRAC: case t_RFRACN:
return isinexactreal((GEN)x[1]) || isinexactreal((GEN)x[2]);
}
lx = (tx==t_POL)? lgef(x): lg(x);
for (i=lontyp[tx]; i<lx; i++)
if (isinexactreal((GEN)x[i])) return 1;
return 0;
}
int
isexactzero(GEN g)
{
long i;
switch (typ(g))
{
case t_INT:
return !signe(g);
case t_REAL: case t_PADIC: case t_SER:
return 0;
case t_INTMOD: case t_POLMOD:
return isexactzero((GEN)g[2]);
case t_FRAC: case t_FRACN: case t_RFRAC: case t_RFRACN:
return isexactzero((GEN)g[1]);
case t_COMPLEX:
return isexactzero((GEN)g[1]) && isexactzero((GEN)g[2]);
case t_QUAD:
return isexactzero((GEN)g[2]) && isexactzero((GEN)g[3]);
case t_POL:
for (i=lgef(g)-1; i>1; i--)
if (!isexactzero((GEN)g[i])) return 0;
return 1;
case t_VEC: case t_COL: case t_MAT:
for (i=lg(g)-1; i; i--)
if (!isexactzero((GEN)g[i])) return 0;
return 1;
}
return 0;
}
int
iscomplex(GEN x)
{
switch(typ(x))
{
case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
return 0;
case t_COMPLEX:
return !gcmp0((GEN)x[2]);
case t_QUAD:
return signe(mael(x,1,2)) > 0;
}
err(typeer,"iscomplex");
return 0; /* not reached */
}
int
ismonome(GEN x)
{
long i;
if (typ(x)!=t_POL || !signe(x)) return 0;
for (i=lgef(x)-2; i>1; i--)
if (!isexactzero((GEN)x[i])) return 0;
return 1;
}
/********************************************************************/
/** **/
/** MULTIPLICATION SIMPLE **/
/** **/
/********************************************************************/
#define fix_frac(z) if (signe(z[2])<0)\
{\
setsigne(z[1],-signe(z[1]));\
setsigne(z[2],1);\
}
/* assume z[1] was created last */
#define fix_frac_if_int(z) if (is_pm1(z[2]))\
z = gerepileupto((long)(z+3), (GEN)z[1]);
GEN
gmulsg(long s, GEN y)
{
long ty=typ(y),ly=lg(y),i,av,tetpil;
GEN z,p1,p2;
switch(ty)
{
case t_INT:
return mulsi(s,y);
case t_REAL:
return mulsr(s,y);
case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)y[1];
(void)new_chunk(lgefint(p2)<<2);
p1=mulsi(s,(GEN)y[2]); avma=(long)z;
z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
case t_FRAC:
if (!s) return gzero;
z = cgetg(3,t_FRAC);
i = cgcd(s, smodis((GEN)y[2], s));
if (i == 1)
{
z[2] = licopy((GEN)y[2]);
z[1] = lmulis((GEN)y[1], s);
}
else
{
z[2] = ldivis((GEN)y[2], i);
z[1] = lmulis((GEN)y[1], s/i);
fix_frac_if_int(z);
}
return z;
case t_FRACN: z=cgetg(3,ty);
z[1]=lmulsi(s,(GEN)y[1]);
z[2]=licopy((GEN)y[2]);
return z;
case t_COMPLEX: z=cgetg(ly,ty);
z[1]=lmulsg(s,(GEN)y[1]);
z[2]=lmulsg(s,(GEN)y[2]); return z;
case t_PADIC:
if (!s) return gzero;
av=avma; p1=cgetp(y); gaffsg(s,p1); tetpil=avma;
return gerepile(av,tetpil,gmul(p1,y));
case t_QUAD: z=cgetg(ly,ty);
copyifstack(y[1],z[1]);
z[2]=lmulsg(s,(GEN)y[2]);
z[3]=lmulsg(s,(GEN)y[3]); return z;
case t_POLMOD: z=cgetg(ly,ty);
z[2]=lmulsg(s,(GEN)y[2]);
copyifstack(y[1],z[1]); return z;
case t_POL:
if (!s || !signe(y)) return zeropol(varn(y));
ly=lgef(y); z=cgetg(ly,t_POL); z[1]=y[1];
for (i=2; i<ly; i++) z[i]=lmulsg(s,(GEN)y[i]);
return normalizepol_i(z, ly);
case t_SER:
if (!s) return zeropol(varn(y));
if (gcmp0(y)) return gcopy(y);
z=cgetg(ly,ty);
for (i=2; i<ly; i++) z[i]=lmulsg(s,(GEN)y[i]);
z[1]=y[1]; return normalize(z);
case t_RFRAC:
if (!s) return zeropol(gvar(y));
z = cgetg(3, t_RFRAC);
i = ggcd(stoi(s),(GEN)y[2])[2];
avma = (long)z;
if (i == 1)
{
z[1]=lmulgs((GEN)y[1], s);
z[2]= lcopy((GEN)y[2]);
}
else
{
z[1] = lmulgs((GEN)y[1], s/i);
z[2] = ldivgs((GEN)y[2], i);
}
return z;
case t_RFRACN:
if (!s) return zeropol(gvar(y));
z=cgetg(ly,t_RFRACN);
z[1]=lmulsg(s,(GEN)y[1]);
z[2]=lcopy((GEN)y[2]); return z;
case t_VEC: case t_COL: case t_MAT:
z=cgetg(ly,ty);
for (i=1; i<ly; i++) z[i]=lmulsg(s,(GEN)y[i]);
return z;
}
err(typeer,"gmulsg");
return NULL; /* not reached */
}
/********************************************************************/
/** **/
/** DIVISION SIMPLE **/
/** **/
/********************************************************************/
GEN
gdivgs(GEN x, long s)
{
static long court[] = { evaltyp(t_INT) | m_evallg(3),0,0 };
long tx=typ(x),lx=lg(x),i,av;
GEN z,p1;
if (!s) err(gdiver2);
switch(tx)
{
case t_INT:
av=avma; z=dvmdis(x,s,&p1);
if (p1==gzero) return z;
i = cgcd(s, smodis(x,s));
avma=av; z=cgetg(3,t_FRAC);
if (i == 1)
x = icopy(x);
else
{
s /= i;
x = divis(x, i);
}
z[1]=(long)x; z[2]=lstoi(s);
fix_frac(z); return z;
case t_REAL:
return divrs(x,s);
case t_FRAC: z=cgetg(3,tx);
i = cgcd(s, smodis((GEN)x[1], s));
if (i == 1)
{
z[2] = lmulsi(s, (GEN)x[2]);
z[1] = licopy((GEN)x[1]);
}
else
{
z[2] = lmulsi(s/i, (GEN)x[2]);
z[1] = ldivis((GEN)x[1], i);
}
fix_frac(z);
fix_frac_if_int(z); return z;
case t_FRACN: z = cgetg(3,tx);
z[1]=licopy((GEN)x[1]);
z[2]=lmulsi(s,(GEN)x[2]);
fix_frac(z); return z;
case t_COMPLEX: z=cgetg(lx,tx);
z[1]=ldivgs((GEN)x[1],s);
z[2]=ldivgs((GEN)x[2],s);
return z;
case t_QUAD: z=cgetg(lx,tx);
copyifstack(x[1],z[1]);
for (i=2; i<4; i++) z[i]=ldivgs((GEN)x[i],s);
return z;
case t_POLMOD: z=cgetg(lx,tx);
copyifstack(x[1],z[1]);
z[2]=ldivgs((GEN)x[2],s);
return z;
case t_POL: lx=lgef(x); z=cgetg(lx,tx);
for (i=2; i<lx; i++) z[i]=ldivgs((GEN)x[i],s);
z[1]=x[1]; return z;
case t_SER: z=cgetg(lx,tx);
for (i=2; i<lx; i++) z[i]=ldivgs((GEN)x[i],s);
z[1]=x[1]; return z;
case t_RFRAC:
z = cgetg(3, t_RFRAC);
i = ggcd(stoi(s),(GEN)x[1])[2];
avma = (long)z;
if (i == 1)
{
z[2]=lmulsg(s,(GEN)x[2]);
z[1]=lcopy((GEN)x[1]); return z;
}
z[1] = ldivgs((GEN)x[1], i);
z[2] = lmulgs((GEN)x[2], s/i); return z;
case t_RFRACN: z=cgetg(3,t_RFRACN);
z[2]=lmulsg(s,(GEN)x[2]);
z[1]=lcopy((GEN)x[1]); return z;
case t_VEC: case t_COL: case t_MAT: z=cgetg(lx,tx);
for (i=1; i<lx; i++) z[i]=ldivgs((GEN)x[i],s);
return z;
}
affsi(s,court); return gdiv(x,court);
}
/*******************************************************************/
/* */
/* MODULO GENERAL */
/* */
/*******************************************************************/
GEN
gmod(GEN x, GEN y)
{
long av,tetpil,i, tx=typ(x);
GEN z,p1;
if (is_matvec_t(tx))
{
i=lg(x); z=cgetg(i,tx);
for (i--; i; i--) z[i]=lmod((GEN)x[i],y);
return z;
}
switch(typ(y))
{
case t_INT:
switch(tx)
{
case t_INT:
return modii(x,y);
case t_INTMOD: z=cgetg(3,tx);
z[1]=lmppgcd((GEN)x[1],y);
z[2]=lmodii((GEN)x[2],(GEN)z[1]); return z;
case t_FRAC: case t_FRACN:
av=avma; if (tx==t_FRACN) x=gred(x);
p1=mulii((GEN)x[1],mpinvmod((GEN)x[2],y));
tetpil=avma; return gerepile(av,tetpil,modii(p1,y));
case t_QUAD: z=cgetg(4,tx);
copyifstack(x[1],z[1]);
z[2]=lmod((GEN)x[2],y);
z[3]=lmod((GEN)x[3],y); return z;
case t_PADIC:
{
long p1[] = {evaltyp(t_INTMOD)|m_evallg(3),0,0};
p1[1] = (long)y;
p1[2] = lgeti(lgefint(y));
gaffect(x,p1); return (GEN)p1[2];
}
case t_POLMOD: case t_POL:
return gzero;
default: err(gmoder1);
}
case t_REAL: case t_FRAC: case t_FRACN:
switch(tx)
{
case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
av=avma; p1 = gfloor(gdiv(x,y)); p1 = gneg_i(gmul(p1,y));
tetpil=avma; return gerepile(av,tetpil,gadd(x,p1));
case t_POLMOD: case t_POL:
return gzero;
default: err(gmoder1);
}
case t_POL:
if (is_scalar_t(tx))
{
if (tx!=t_POLMOD || varn(x[1]) > varn(y))
return (lgef(y) < 4)? gzero: gcopy(x);
if (varn(x[1])!=varn(y)) return gzero;
z=cgetg(3,t_POLMOD);
z[1]=lgcd((GEN)x[1],y);
z[2]=lres((GEN)x[2],(GEN)z[1]); return z;
}
switch(tx)
{
case t_POL:
return gres(x,y);
case t_RFRAC: case t_RFRACN:
av=avma; if (tx==t_RFRACN) x=gred_rfrac(x);
p1=gmul((GEN)x[1],ginvmod((GEN)x[2],y)); tetpil=avma;
return gerepile(av,tetpil,gres(p1,y));
default: err(talker,"type mod polynomial forbidden in gmod");
}
}
err(talker,"modulus type forbidden in gmod");
return NULL; /* not reached */
}
GEN
gmodulsg(long x, GEN y)
{
GEN z;
switch(typ(y))
{
case t_INT: z=cgetg(3,t_INTMOD);
z[1]=labsi(y); z[2]=lmodsi(x,y); return z;
case t_POL: z=cgetg(3,t_POLMOD);
z[1]=lcopy(y); z[2]=lstoi(x); return z;
}
err(gmoder1); return NULL; /* not reached */
}
GEN
gmodulss(long x, long y)
{
GEN z=cgetg(3,t_INTMOD);
y=labs(y); z[1]=lstoi(y); z[2]=lstoi(x % y); return z;
}
GEN
gmodulo(GEN x,GEN y)
{
long tx=typ(x),l,i;
GEN z;
if (is_matvec_t(tx))
{
l=lg(x); z=cgetg(l,tx);
for (i=1; i<l; i++) z[i] = (long) gmodulo((GEN)x[i],y);
return z;
}
switch(typ(y))
{
case t_INT:
if (tx!=t_INT && !is_frac_t(tx) && tx!=t_PADIC) break;
z=cgetg(3,t_INTMOD);
if (!signe(y)) err(talker,"zero modulus in gmodulo");
y = gclone(y); setsigne(y,1);
z[1]=(long)y;
z[2]=lmod(x,y); return z;
case t_POL: z=cgetg(3,t_POLMOD);
z[1] = lclone(y);
if (is_scalar_t(tx)) { z[2]=lcopy(x); return z; }
if (tx==t_POL || is_rfrac_t(tx)) { z[2]=lmod(x,y); return z; }
}
err(gmoder1); return NULL; /* not reached */
}
GEN
gmodulcp(GEN x,GEN y)
{
long tx=typ(x),l,i;
GEN z;
if (is_matvec_t(tx))
{
l=lg(x); z=cgetg(l,tx);
for (i=1; i<l; i++) z[i] = lmodulcp((GEN)x[i],y);
return z;
}
switch(typ(y))
{
case t_INT:
if (tx!=t_INT && !is_frac_t(tx) && tx!=t_PADIC) break;
z=cgetg(3,t_INTMOD);
z[1]=labsi(y);
z[2]=lmod(x,y); return z;
case t_POL: z=cgetg(3,t_POLMOD);
z[1]=lcopy(y);
if (is_scalar_t(tx)) { z[2]=lcopy(x); return z; }
if (tx==t_POL || is_rfrac_t(tx)) { z[2]=lmod(x,y); return z; }
}
err(gmoder1); return NULL; /* not reached */
}
GEN
Mod0(GEN x,GEN y,long flag)
{
switch(flag)
{
case 0: return gmodulcp(x,y);
case 1: return gmodulo(x,y);
default: err(flagerr,"Mod");
}
return NULL; /* not reached */
}
/*******************************************************************/
/* */
/* DIVISION ENTIERE GENERALE */
/* DIVISION ENTIERE AVEC RESTE GENERALE */
/* */
/*******************************************************************/
GEN
gdivent(GEN x, GEN y)
{
long tx=typ(x),ty=typ(y);
if (tx == t_INT)
{
if (ty == t_INT)
{
GEN z,p1;
long av,tetpil;
if (signe(x)>=0) return divii(x,y);
av=avma; z=dvmdii(x,y,&p1);
if (p1==gzero) return z;
tetpil=avma; return gerepile(av,tetpil,addsi(-signe(y),z));
}
if (ty!=t_POL) err(typeer,"gdivent");
return gzero;
}
if (ty != t_POL) err(typeer,"gdivent");
if (tx == t_POL) return gdeuc(x,y);
if (! is_scalar_t(tx)) err(typeer,"gdivent");
return gzero;
}
GEN
gdiventres(GEN x, GEN y)
{
long tx=typ(x),ty=typ(y);
GEN z = cgetg(3,t_COL);
if (tx==t_INT)
{
if (ty==t_INT) z[1]=(long)truedvmdii(x,y,(GEN*)(z+2));
else
{
if (ty!=t_POL) err(typeer,"gdiventres");
z[1]=zero; z[2]=licopy(x);
}
return z;
}
if (ty != t_POL) err(typeer,"gdiventres");
if (tx == t_POL)
{
z[1]=ldivres(x,y,(GEN *)(z+2));
return z;
}
if (! is_scalar_t(tx)) err(typeer,"gdiventres");
z[1]=zero; z[2]=lcopy(x); return z;
}
GEN
gdivmod(GEN x, GEN y, GEN *pr)
{
long ty,tx=typ(x);
if (tx==t_INT)
{
ty=typ(y);
if (ty==t_INT) return dvmdii(x,y,pr);
if (ty==t_POL) { *pr=gcopy(x); return gzero; }
err(typeer,"gdivmod");
}
if (tx!=t_POL) err(typeer,"gdivmod");
return poldivres(x,y,pr);
}
/* When x and y are integers, compute the quotient x/y, rounded to the
* nearest integer. If there is a tie, the quotient is rounded towards zero.
* If x and y are not both integers, same as gdivent.
*/
GEN
gdivround(GEN x, GEN y)
{
long tx=typ(x),ty=typ(y);
if (tx==t_INT)
{
if (ty==t_INT)
{
long av = avma, av1,fl;
GEN r, q=dvmdii(x,y,&r); /* q = x/y rounded towards 0, sqn(r)=sgn(x) */
if (r==gzero) return q;
av1 = avma;
fl = absi_cmp(shifti(r,1),y);
avma = av1; cgiv(r);
if (fl >= 0) /* If 2*|r| >= |y| */
{
long sz = signe(x)*signe(y);
if (fl || sz > 0)
{ av1=avma; q = gerepile(av,av1,addis(q,sz)); }
}
return q;
}
if (ty!=t_POL) err(typeer,"gdivround");
return gzero;
}
if (ty != t_POL) err(typeer,"gdivround");
if (tx == t_POL) return gdeuc(x,y);
if (! is_scalar_t(tx)) err(typeer,"gdivround");
return gzero;
}
/*******************************************************************/
/* */
/* SHIFT D'UN GEN */
/* */
/*******************************************************************/
/* Shift tronque si n<0 (multiplication tronquee par 2^n) */
GEN
gshift(GEN x, long n)
{
long i,l,lx, tx = typ(x);
GEN y;
switch(tx)
{
case t_INT:
return shifti(x,n);
case t_REAL:
return shiftr(x,n);
case t_VEC: case t_COL: case t_MAT:
lx=lg(x); y=cgetg(lx,tx); l=lontyp[tx];
for (i=1; i<l ; i++) y[i]=x[i];
for ( ; i<lx; i++) y[i]=lshift((GEN)x[i],n);
return y;
}
return gmul2n(x,n);
}
GEN mulscalrfrac(GEN x, GEN y);
/* Shift vrai (multiplication exacte par 2^n) */
GEN
gmul2n(GEN x, long n)
{
long tx,lx,i,k,l,av,tetpil;
GEN p2,p1,y;
tx=typ(x);
switch(tx)
{
case t_INT:
if (n>=0) return shifti(x,n);
if (!signe(x)) return gzero;
l = vali(x); n = -n;
if (n<=l) return shifti(x,-n);
y=cgetg(3,t_FRAC);
y[1]=lshifti(x,-l);
y[2]=lshifti(gun,n-l); return y;
case t_REAL:
return shiftr(x,n);
case t_INTMOD:
if (n > 0)
{
y=cgetg(3,t_INTMOD); p2=(GEN)x[1];
av=avma; new_chunk(lgefint(p2) * (3 + (n>>TWOPOTBITS_IN_LONG)));
p1 = shifti((GEN)x[2],n); avma=av;
y[2]=lmodii(p1,p2); icopyifstack(p2,y[1]); return y;
}
l=avma; y=gmul2n(gun,n); tetpil=avma;
return gerepile(l,tetpil,gmul(y,x));
case t_FRAC: case t_FRACN:
l = vali((GEN)x[1]);
k = vali((GEN)x[2]);
if (n+l-k>=0)
{
if (expi((GEN)x[2]) == k) /* x[2] power of 2 */
return shifti((GEN)x[1],n-k);
l = n-k; k = -k;
}
else
{
k = -l-n; l = -l;
}
y=cgetg(3,t_FRAC);
y[1]=lshifti((GEN)x[1],l);
y[2]=lshifti((GEN)x[2],k); return y;
case t_QUAD: y=cgetg(4,t_QUAD);
copyifstack(x[1],y[1]);
for (i=2; i<4; i++) y[i]=lmul2n((GEN)x[i],n);
return y;
case t_POLMOD: y=cgetg(3,t_POLMOD);
copyifstack(x[1],y[1]);
y[2]=lmul2n((GEN)x[2],n); return y;
case t_POL: case t_COMPLEX: case t_SER:
case t_VEC: case t_COL: case t_MAT:
lx = (tx==t_POL)? lgef(x): lg(x);
y=cgetg(lx,tx); l=lontyp[tx];
for (i=1; i<l ; i++) y[i]=x[i];
for ( ; i<lx; i++) y[i]=lmul2n((GEN)x[i],n);
return y;
case t_RFRAC: av=avma; p1 = gmul2n(gun,n); tetpil = avma;
return gerepile(av,tetpil, mulscalrfrac(p1,x));
case t_RFRACN: y=cgetg(3,tx);
if (n>=0)
{
y[1]=lmul2n((GEN)x[1],n); y[2]=lcopy((GEN)x[2]);
}
else
{
y[2]=lmul2n((GEN)x[2],-n); y[1]=lcopy((GEN)x[1]);
}
return y;
case t_PADIC:
l=avma; y=gmul2n(gun,n); tetpil=avma;
return gerepile(l,tetpil,gmul(y,x));
}
err(typeer,"gmul2n");
return NULL; /* not reached */
}
/*******************************************************************/
/* */
/* INVERSE D'UN GEN */
/* */
/*******************************************************************/
GEN fix_rfrac_if_pol(GEN x, GEN y);
GEN
ginv(GEN x)
{
long tx=typ(x),av,tetpil,s;
GEN z,y,p1,p2;
switch(tx)
{
case t_INT:
if (is_pm1(x)) return icopy(x);
if (!signe(x)) err(gdiver2);
z=cgetg(3,t_FRAC);
z[1] = (signe(x)<0)? lnegi(gun): un;
z[2]=labsi(x); return z;
case t_REAL:
return divsr(1,x);
case t_INTMOD: z=cgetg(3,t_INTMOD);
icopyifstack(x[1],z[1]);
z[2]=lmpinvmod((GEN)x[2],(GEN)x[1]); return z;
case t_FRAC: case t_FRACN: z=cgetg(3,tx);
s = signe(x[1]);
if (!s) err(gdiver2);
if (is_pm1(x[1]))
return s>0? icopy((GEN)x[2]): negi((GEN)x[2]);
z[1]=licopy((GEN)x[2]);
z[2]=licopy((GEN)x[1]);
if (s < 0)
{
setsigne(z[1],-signe(z[1]));
setsigne(z[2],1);
}
return z;
case t_COMPLEX: case t_QUAD:
av=avma; p1=gnorm(x); p2=gconj(x); tetpil=avma;
return gerepile(av,tetpil,gdiv(p2,p1));
case t_PADIC: z = cgetg(5,t_PADIC);
if (!signe(x[4])) err(gdiver2);
z[1] = evalprecp(precp(x)) | evalvalp(-valp(x));
icopyifstack(x[2], z[2]);
z[3] = licopy((GEN)x[3]);
z[4] = lmpinvmod((GEN)x[4],(GEN)z[3]); return z;
case t_POLMOD: z=cgetg(3,t_POLMOD);
copyifstack(x[1],z[1]);
z[2]=linvmod((GEN)x[2],(GEN)x[1]); return z;
case t_POL: case t_SER:
return gdiv(gun,x);
case t_RFRAC: case t_RFRACN:
if (gcmp0((GEN) x[1])) err(gdiver2);
p1 = fix_rfrac_if_pol((GEN)x[2],(GEN)x[1]);
if (p1) return p1;
z=cgetg(3,tx);
z[1]=lcopy((GEN)x[2]);
z[2]=lcopy((GEN)x[1]); return z;
case t_QFR:
{
long k,l;
l=signe(x[2]); setsigne(x[2],-l);
k=signe(x[4]); setsigne(x[4],-k); z=redreal(x);
setsigne(x[2],l); setsigne(x[4],k); return z;
}
case t_QFI:
y=gcopy(x);
if (!egalii((GEN)x[1],(GEN)x[2]) && !egalii((GEN)x[1],(GEN)x[3]))
setsigne(y[2],-signe(y[2]));
return y;
case t_MAT:
return (lg(x)==1)? cgetg(1,t_MAT): invmat(x);
}
err(typeer,"inverse");
return NULL; /* not reached */
}
/*******************************************************************/
/* */
/* SUBSTITUTION DANS UN POLYNOME OU UNE SERIE */
/* */
/*******************************************************************/
/* Convert t_SER --> t_POL / t_RFRAC */
static GEN
gconvsp(GEN x, int flpile)
{
long v = varn(x), av,tetpil,i;
GEN p1,y;
if (gcmp0(x)) return zeropol(v);
av=avma; y=dummycopy(x); settyp(y,t_POL);
i=lg(x)-1; while (i>1 && gcmp0((GEN)y[i])) i--;
setlgef(y,i+1);
p1=gpuigs(polx[v],valp(x));
tetpil=avma; p1=gmul(p1,y);
return flpile? gerepile(av,tetpil,p1): p1;
}
GEN
gsubst(GEN x, long v, GEN y)
{
long tx = typ(x), ty = typ(y), lx = lg(x), ly = lg(y);
long l,vx,vy,e,ex,ey,tetpil,av,i,j,k,jb,limite;
GEN t,p1,p2,z;
if (ty==t_MAT)
{
if (ly==1) return cgetg(1,t_MAT);
if (ly != lg(y[1]))
err(talker,"forbidden substitution by a non square matrix");
} else if (is_graphicvec_t(ty))
err(talker,"forbidden substitution by a vector");
if (is_scalar_t(tx))
{
if (tx!=t_POLMOD || v <= varn(x[1]))
{
if (ty==t_MAT) return gscalmat(x,ly-1);
return gcopy(x);
}
av=avma;
p1=gsubst((GEN)x[1],v,y); vx=varn(p1);
p2=gsubst((GEN)x[2],v,y); vy=gvar(p2);
if (typ(p1)!=t_POL)
err(talker,"forbidden substitution in a scalar type");
if (vy>=vx)
{
tetpil=avma;
return gerepile(av,tetpil,gmodulcp(p2,p1));
}
lx=lgef(p2); tetpil=avma; z=cgetg(lx,t_POL); z[1]=p2[1];
for (i=2; i<lx; i++) z[i]=lmodulcp((GEN)p2[i],p1);
return gerepile(av,tetpil, normalizepol_i(z,lx));
}
switch(tx)
{
case t_POL:
l=lgef(x);
if (l==2)
return (ty==t_MAT)? gscalmat(gzero,ly-1): gzero;
vx=varn(x);
if (vx<v)
{
av=avma; p1=polx[vx]; z= gsubst((GEN)x[l-1],v,y);
for (i=l-1; i>2; i--) z=gadd(gmul(z,p1),gsubst((GEN)x[i-1],v,y));
return gerepileupto(av,z);
}
if (ty!=t_MAT)
return (vx>v)? gcopy(x): poleval(x,y);
if (vx>v) return gscalmat(x,ly-1);
if (l==3) return gscalmat((GEN)x[2],ly-1);
av=avma; z=(GEN)x[l-1];
for (i=l-1; i>2; i--) z=gaddmat((GEN)x[i-1],gmul(z,y));
return gerepileupto(av,z);
case t_SER:
vx=varn(x);
if (vx > v)
return (ty==t_MAT)? gscalmat(x,ly-1): gcopy(x);
ex = valp(x);
if (vx < v)
{
if (!signe(x)) return gcopy(x);
/* a ameliorer */
av=avma; setvalp(x,0); p1=gconvsp(x,0); setvalp(x,ex);
p2=gsubst(p1,v,y); tetpil=avma; z=tayl(p2,vx,lx-2);
if (ex)
{
p1=gpuigs(polx[vx],ex); tetpil=avma; z=gmul(z,p1);
}
return gerepile(av,tetpil,z);
}
switch(ty) /* here vx == v */
{
case t_SER:
ey=valp(y); vy=varn(y);
if (ey<1) return zeroser(vy,ey*(ex+lx-2));
l=(lx-2)*ey+2;
if (ex) { if (l>ly) l=ly; }
else
{
if (gcmp0(y)) l=ey+2;
else { if (l>ly) l=ey+ly; }
}
if (vy!=vx)
{
av=avma; z = zeroser(vy,0);
for (i=lx-1; i>=2; i--)
z = gadd((GEN)x[i], gmul(y,z));
if (ex) z = gmul(z, gpuigs(y,ex));
return gerepileupto(av,z);
}
av=avma; limite=stack_lim(av,1);
t=cgetg(ly,t_SER);
z = scalarser((GEN)x[2],varn(y),l-2);
for (i=2; i<ly; i++) t[i]=y[i];
for (i=3,jb=ey; jb<=l-2; i++,jb+=ey)
{
for (j=jb+2; j<l; j++)
{
p1=gmul((GEN)x[i],(GEN)t[j-jb]);
z[j]=ladd((GEN)z[j],p1);
}
for (j=l-1-jb-ey; j>1; j--)
{
p1=gzero;
for (k=2; k<j; k++)
p1=gadd(p1,gmul((GEN)t[j-k+2],(GEN)y[k]));
p2=gmul((GEN)t[2],(GEN)y[j]);
t[j]=ladd(p1,p2);
}
if (low_stack(limite, stack_lim(av,1)))
{
GEN *gptr[2];
if(DEBUGMEM>1) err(warnmem,"gsubst");
gptr[0]=&z; gptr[1]=&t; gerepilemany(av,gptr,2);
}
}
if (!ex) { tetpil=avma; return gerepile(av,tetpil,gcopy(z)); }
if (l<ly) { setlg(y,l); p1=gpuigs(y,ex); setlg(y,ly); }
else p1=gpuigs(y,ex);
tetpil=avma; return gerepile(av,tetpil,gmul(z,p1));
case t_POL: case t_RFRAC: case t_RFRACN:
if (isexactzero(y)) return scalarser((GEN)x[2],v,lx-2);
vy=gvar(y); e=gval(y,vy);
if (e<=0)
err(talker,"non positive valuation in a series substitution");
av=avma; p1=gconvsp(x,0); p2=gsubst(p1,v,y); tetpil=avma;
return gerepile(av,tetpil,tayl(p2,vy,e*(lx-2+ex)));
default:
err(talker,"non polynomial or series type substituted in a series");
}
break;
case t_RFRAC: case t_RFRACN: av=avma;
p1=gsubst((GEN)x[1],v,y);
p2=gsubst((GEN)x[2],v,y); tetpil=avma;
return gerepile(av,tetpil,gdiv(p1,p2));
case t_VEC: case t_COL: case t_MAT: z=cgetg(lx,tx);
for (i=1; i<lx; i++) z[i]=lsubst((GEN)x[i],v,y);
}
return z;
}
/*******************************************************************/
/* */
/* SERIE RECIPROQUE D'UNE SERIE */
/* */
/*******************************************************************/
GEN
recip(GEN x)
{
long tetpil, av=avma, v=varn(x);
GEN p1,p2,a,y,u;
if (typ(x)!=t_SER) err(talker,"not a series in serreverse");
if (valp(x)!=1) err(talker,"valuation not equal to 1 in serreverse");
a=(GEN)x[2];
if (gcmp1(a))
{
long i,j,k, lx=lg(x), lim=stack_lim(av,2);
u=cgetg(lx,t_SER); y=cgetg(lx,t_SER);
u[1]=y[1]=evalsigne(1) | evalvalp(1) | evalvarn(v);
u[2]=un; u[3]=lmulsg(-2,(GEN)x[3]);
y[2]=un; y[3]=lneg((GEN)x[3]);
for (i=3; i<lx-1; )
{
for (j=3; j<i+1; j++)
{
p1=(GEN)u[j];
for (k=j-1; k>2; k--)
p1=gsub(p1,gmul((GEN)u[k],(GEN)x[j-k+2]));
u[j]=lsub(p1,(GEN)x[j]);
}
p1=gmulsg(i,(GEN)x[i+1]);
for (k=2; k<i; k++)
{
p2=gmul((GEN)x[k+1],(GEN)u[i-k+2]);
p1=gadd(p1,gmulsg(k,p2));
}
i++; u[i]=lneg(p1); y[i]=ldivgs((GEN)u[i],i-1);
if (low_stack(lim, stack_lim(av,2)))
{
GEN *gptr[2];
if(DEBUGMEM>1) err(warnmem,"recip");
for(k=i+1; k<lx; k++) u[k]=y[k]=zero; /* dummy */
gptr[0]=&u; gptr[1]=&y; gerepilemany(av,gptr,2);
}
}
tetpil=avma; return gerepile(av,tetpil,gcopy(y));
}
y=gdiv(x,a); y[2]=un; y=recip(y);
a=gdiv(polx[v],a); tetpil=avma;
return gerepile(av,tetpil,gsubst(y,v,a));
}
/*******************************************************************/
/* */
/* DERIVATION ET INTEGRATION */
/* */
/*******************************************************************/
GEN
derivpol(GEN x)
{
long i,lx = lgef(x)-1;
GEN y;
if (lx<3) return gzero;
y = cgetg(lx,t_POL);
for (i=2; i<lx ; i++) y[i] = lmulsg(i-1,(GEN)x[i+1]);
y[1] = x[1]; return normalizepol_i(y,i);
}
GEN
derivser(GEN x)
{
long i,j,vx = varn(x), e = valp(x), lx = lg(x);
GEN y;
if (gcmp0(x)) return zeroser(vx,e-1);
if (e)
{
y=cgetg(lx,t_SER); y[1] = evalsigne(1) | evalvalp(e-1) | evalvarn(vx);
for (i=2; i<lx; i++) y[i]=lmulsg(i+e-2,(GEN)x[i]);
return y;
}
i=3; while (i<lx && gcmp0((GEN)x[i])) i++;
if (i==lx) return zeroser(vx,lx-3);
lx--; if (lx<3) lx=3;
lx = lx-i+3; y=cgetg(lx,t_SER);
y[1]=evalsigne(1) | evalvalp(i-3) | evalvarn(vx);
for (j=2; j<lx; j++) y[j]=lmulsg(j+i-4,(GEN)x[i+j-2]);
return y;
}
GEN
deriv(GEN x, long v)
{
long lx,vx,tx,e,i,j,l,av,tetpil;
GEN y,p1,p2;
tx=typ(x); if (is_const_t(tx)) return gzero;
if (v < 0) v = gvar(x);
switch(tx)
{
case t_POLMOD:
if (v<=varn(x[1])) return gzero;
y=cgetg(3,t_POLMOD); copyifstack(x[1],y[1]);
y[2]=lderiv((GEN)x[2],v); return y;
case t_POL:
vx=varn(x); if (vx>v) return gzero;
if (vx<v)
{
lx = lgef(x);
y = cgetg(lx,t_POL);
for (i=2; i<lx; i++) y[i] = lderiv((GEN)x[i],v);
y[1] = evalvarn(vx);
return normalizepol_i(y,i);
}
return derivpol(x);
case t_SER:
vx=varn(x); if (vx>v) return gzero;
if (vx<v)
{
if (!signe(x)) return gcopy(x);
lx=lg(x); e=valp(x);
l=avma;
for (i=2; i<lx; i++)
{
if (!gcmp0(deriv((GEN)x[i],v))) break;
avma=l;
}
if (i==lx) return ggrando(polx[vx],e+lx-2);
y=cgetg(lx-i+2,t_SER);
y[1] = evalsigne(1) | evalvalp(e+i-2) | evalvarn(vx);
for (j=2; i<lx; j++,i++) y[j]=lderiv((GEN)x[i],v);
return y;
}
return derivser(x);
case t_RFRAC: case t_RFRACN: av=avma; y=cgetg(3,tx);
y[2]=lsqr((GEN)x[2]); l=avma;
p1=gmul((GEN)x[2],deriv((GEN)x[1],v));
p2=gmul(gneg_i((GEN)x[1]),deriv((GEN)x[2],v));
tetpil=avma; p1=gadd(p1,p2);
if (tx==t_RFRACN) { y[1]=lpile(l,tetpil,p1); return y; }
y[1]=(long)p1; return gerepileupto(av,gred_rfrac(y));
case t_VEC: case t_COL: case t_MAT: lx=lg(x); y=cgetg(lx,tx);
for (i=1; i<lx; i++) y[i]=lderiv((GEN)x[i],v);
return y;
}
err(typeer,"deriv");
return NULL; /* not reached */
}
/*******************************************************************/
/* */
/* INTEGRATION FORMELLE */
/* */
/*******************************************************************/
static GEN
triv_integ(GEN x, long v, long tx, long lx)
{
GEN y=cgetg(lx,tx);
long i;
y[1]=x[1];
for (i=2; i<lx; i++) y[i]=linteg((GEN)x[i],v);
return y;
}
GEN
integ(GEN x, long v)
{
long lx,tx,e,i,j,vx,n,av=avma,tetpil;
GEN y,p1;
tx = typ(x);
if (v < 0) v = gvar(x);
if (is_scalar_t(tx))
{
if (tx == t_POLMOD && v>varn(x[1]))
{
y=cgetg(3,t_POLMOD); copyifstack(x[1],y[1]);
y[2]=linteg((GEN)x[2],v); return y;
}
if (gcmp0(x)) return gzero;
y=cgetg(4,t_POL); y[1] = evalsigne(1) | evallgef(4) | evalvarn(v);
y[2]=zero; y[3]=lcopy(x); return y;
}
switch(tx)
{
case t_POL:
vx=varn(x); lx=lgef(x);
if (lx==2) return zeropol(min(v,vx));
if (vx>v)
{
y=cgetg(4,t_POL);
y[1] = signe(x)? evallgef(4) | evalvarn(v) | evalsigne(1)
: evallgef(4) | evalvarn(v);
y[2]=zero; y[3]=lcopy(x); return y;
}
if (vx<v) return triv_integ(x,v,tx,lx);
y=cgetg(lx+1,tx); y[2]=zero;
for (i=3; i<=lx; i++)
y[i]=ldivgs((GEN)x[i-1],i-2);
y[1] = signe(x)? evallgef(lx+1) | evalvarn(v) | evalsigne(1)
: evallgef(lx+1) | evalvarn(v);
return y;
case t_SER:
lx=lg(x); e=valp(x); vx=varn(x);
if (!signe(x))
{
if (vx == v) e++; else if (vx < v) v = vx;
return zeroser(v,e);
}
if (vx>v)
{
y=cgetg(4,t_POL);
y[1] = evallgef(4) | evalvarn(v) | evalsigne(1);
y[2]=zero; y[3]=lcopy(x); return y;
}
if (vx<v) return triv_integ(x,v,tx,lx);
y=cgetg(lx,tx);
for (i=2; i<lx; i++)
{
j=i+e-1;
if (!j)
{
if (!gcmp0((GEN)x[i])) err(inter2);
y[i]=zero;
}
else y[i] = ldivgs((GEN)x[i],j);
}
y[1]=x[1]+1; return y;
case t_RFRAC: case t_RFRACN:
vx = gvar(x);
if (vx>v)
{
y=cgetg(4,t_POL);
y[1] = signe(x[1])? evallgef(4) | evalvarn(v) | evalsigne(1)
: evallgef(4) | evalvarn(v);
y[2]=zero; y[3]=lcopy(x); return y;
}
if (vx<v)
{
p1=cgetg(v+2,t_VEC);
for (i=0; i<vx; i++) p1[i+1] = lpolx[i];
for (i=vx+1; i<v; i++) p1[i+1] = lpolx[i];
p1[v+1] = lpolx[vx]; p1[vx+1] = lpolx[v];
y=integ(changevar(x,p1),vx); tetpil=avma;
return gerepile(av,tetpil,changevar(y,p1));
}
n=lgef(x[1])+lgef(x[2])-4;
y=gdiv(gtrunc(gmul((GEN)x[2], integ(tayl(x,v,n),v))), (GEN)x[2]);
if (!gegal(deriv(y,v),x)) err(inter2);
if (typ(y)==t_RFRAC && lgef(y[1]) == lgef(y[2]))
{
GEN p2;
tx=typ(y[1]); p1=is_scalar_t(tx)? (GEN)y[1]: leading_term(y[1]);
tx=typ(y[2]); p2=is_scalar_t(tx)? (GEN)y[2]: leading_term(y[2]);
y=gsub(y, gdiv(p1,p2));
}
return gerepileupto(av,y);
case t_VEC: case t_COL: case t_MAT:
lx=lg(x); y=cgetg(lx,tx);
for (i=1; i<lg(x); i++) y[i]=linteg((GEN)x[i],v);
return y;
}
err(typeer,"integ");
return NULL; /* not reached */
}
/*******************************************************************/
/* */
/* PARTIES ENTIERES */
/* */
/*******************************************************************/
GEN
gfloor(GEN x)
{
GEN y;
long i,lx, tx = typ(x);
switch(tx)
{
case t_INT: case t_POL:
return gcopy(x);
case t_REAL:
return mpent(x);
case t_FRAC: case t_FRACN:
return truedvmdii((GEN)x[1],(GEN)x[2],NULL);
case t_RFRAC: case t_RFRACN:
return gdeuc((GEN)x[1],(GEN)x[2]);
case t_VEC: case t_COL: case t_MAT:
lx=lg(x); y=cgetg(lx,tx);
for (i=1; i<lx; i++) y[i]=lfloor((GEN)x[i]);
return y;
}
err(typeer,"gfloor");
return NULL; /* not reached */
}
GEN
gfrac(GEN x)
{
long av = avma,tetpil;
GEN p1 = gneg_i(gfloor(x));
tetpil=avma; return gerepile(av,tetpil,gadd(x,p1));
}
GEN
gceil(GEN x)
{
GEN y,p1;
long i,lx,tx=typ(x),av,tetpil;
switch(tx)
{
case t_INT: case t_POL:
return gcopy(x);
case t_REAL:
av=avma; y=mpent(x);
if (!gegal(x,y))
{
tetpil=avma; return gerepile(av,tetpil,addsi(1,y));
}
return y;
case t_FRAC: case t_FRACN:
av=avma; y=dvmdii((GEN)x[1],(GEN)x[2],&p1);
if (p1 != gzero && gsigne(x)>0)
{
cgiv(p1); tetpil=avma;
return gerepile(av,tetpil,addsi(1,y));
}
return y;
case t_RFRAC: case t_RFRACN:
return gdeuc((GEN)x[1],(GEN)x[2]);
case t_VEC: case t_COL: case t_MAT:
lx=lg(x); y=cgetg(lx,tx);
for (i=1; i<lx; i++) y[i]=lceil((GEN)x[i]);
return y;
}
err(typeer,"gceil");
return NULL; /* not reached */
}
GEN
round0(GEN x, GEN *pte)
{
if (pte) { long e; x = grndtoi(x,&e); *pte = stoi(e); }
return ground(x);
}
GEN
ground(GEN x)
{
GEN y,p1;
long i,lx,tx=typ(x),av,tetpil;
switch(tx)
{
case t_INT: case t_INTMOD: case t_QUAD:
return gcopy(x);
case t_REAL:
{
long ex, s = signe(x);
if (s==0 || (ex=expo(x)) < -1) return gzero;
if (ex < 0) return s>0? gun: negi(gun);
av=avma; p1 = realun(3 + (ex>>TWOPOTBITS_IN_LONG));
setexpo(p1,-1); /* p1 = 0.5 */
p1 = addrr(x,p1); tetpil = avma;
return gerepile(av,tetpil,mpent(p1));
}
case t_FRAC: case t_FRACN:
av=avma; p1 = addii(shifti((GEN)x[2], -1), (GEN)x[1]);
return gerepileuptoint(av, truedvmdii(p1, (GEN)x[2], NULL));
case t_POLMOD: y=cgetg(3,t_POLMOD);
copyifstack(x[1],y[1]);
y[2]=lround((GEN)x[2]); return y;
case t_COMPLEX: case t_POL: case t_SER: case t_RFRAC: case t_RFRACN:
case t_VEC: case t_COL: case t_MAT:
lx = (tx==t_POL)? lgef(x): lg(x);
y=cgetg(lx,tx);
for (i=1; i<lontyp[tx]; i++) y[i]=x[i];
for ( ; i<lx; i++) y[i]=lround((GEN)x[i]);
if (tx==t_POL) return normalizepol_i(y, lx);
if (tx==t_SER) return normalize(y);
return y;
}
err(typeer,"ground");
return NULL; /* not reached */
}
/* e = number of error bits on integral part */
GEN
grndtoi(GEN x, long *e)
{
GEN y,p1;
long i,tx=typ(x), lx,av,ex,e1;
*e = -HIGHEXPOBIT;
switch(tx)
{
case t_INT: case t_INTMOD: case t_QUAD:
case t_FRAC: case t_FRACN:
return ground(x);
case t_REAL:
av=avma; p1=gadd(ghalf,x); ex=expo(p1);
if (ex<0)
{
if (signe(p1)>=0) { *e=expo(x); avma=av; return gzero; }
*e=expo(addsr(1,x)); avma=av; return negi(gun);
}
lx=lg(x); e1 = ex - bit_accuracy(lx) + 1;
settyp(p1,t_INT); setlgefint(p1,lx);
y=shifti(p1,e1); if (signe(x)<0) y=addsi(-1,y);
y = gerepileupto(av,y);
if (e1<=0) { av=avma; e1=expo(subri(x,y)); avma=av; }
*e=e1; return y;
case t_POLMOD: y=cgetg(3,t_POLMOD);
copyifstack(x[1],y[1]);
y[2]=lrndtoi((GEN)x[2],e); return y;
case t_COMPLEX: case t_POL: case t_SER: case t_RFRAC: case t_RFRACN:
case t_VEC: case t_COL: case t_MAT:
lx=(tx==t_POL)? lgef(x): lg(x);
y=cgetg(lx,tx);
for (i=1; i<lontyp[tx]; i++) y[i]=x[i];
for ( ; i<lx; i++)
{
y[i]=lrndtoi((GEN)x[i],&e1);
if (e1>*e) *e=e1;
}
return y;
}
err(typeer,"grndtoi");
return NULL; /* not reached */
}
/* e = number of error bits on integral part */
GEN
gcvtoi(GEN x, long *e)
{
long tx=typ(x), lx,i,ex,av,e1;
GEN y;
*e = -HIGHEXPOBIT;
if (tx == t_REAL)
{
long x0,x1;
ex=expo(x); if (ex<0) { *e=ex; return gzero; }
lx=lg(x); e1 = ex - bit_accuracy(lx) + 1;
x0=x[0]; x1=x[1]; settyp(x,t_INT); setlgefint(x,lx);
y=shifti(x,e1); x[0]=x0; x[1]=x1;
if (e1<=0) { av=avma; e1=expo(subri(x,y)); avma=av; }
*e=e1; return y;
}
if (is_matvec_t(tx))
{
lx=lg(x); y=cgetg(lx,tx);
for (i=1; i<lx; i++)
{
y[i]=lcvtoi((GEN)x[i],&e1);
if (e1>*e) *e=e1;
}
return y;
}
return gtrunc(x);
}
GEN
gtrunc(GEN x)
{
long tx=typ(x),av,tetpil,i,v;
GEN y;
switch(tx)
{
case t_INT: case t_POL:
return gcopy(x);
case t_REAL:
return mptrunc(x);
case t_FRAC: case t_FRACN:
return divii((GEN)x[1],(GEN)x[2]);
case t_PADIC:
if (!signe(x[4])) return gzero;
v = valp(x);
if (!v) return gcopy((GEN)x[4]);
if (v>0)
{ /* here p^v is an integer */
av=avma; y=gpuigs((GEN)x[2],v); tetpil=avma;
return gerepile(av,tetpil, mulii(y,(GEN)x[4]));
}
y=cgetg(3,t_FRAC);
y[1]=licopy((GEN)x[4]);
y[2]=lpuigs((GEN)x[2],-v);
return y;
case t_RFRAC: case t_RFRACN:
return gdeuc((GEN)x[1],(GEN)x[2]);
case t_SER:
return gconvsp(x,1);
case t_VEC: case t_COL: case t_MAT:
{
long lx=lg(x);
y=cgetg(lx,tx);
for (i=1; i<lx; i++) y[i]=ltrunc((GEN)x[i]);
return y;
}
}
err(typeer,"gtrunc");
return NULL; /* not reached */
}
GEN
trunc0(GEN x, GEN *pte)
{
if (pte) { long e; x = gcvtoi(x,&e); *pte = stoi(e); }
return gtrunc(x);
}
/*******************************************************************/
/* */
/* CONVERSIONS GEN --> POLYNOMIALS & SERIES */
/* */
/*******************************************************************/
GEN
zeropol(long v)
{
GEN x = cgetg(2,t_POL);
x[1] = evallgef(2) | evalvarn(v); return x;
}
GEN
scalarpol(GEN x, long v)
{
GEN y=cgetg(3,t_POL);
y[1] = gcmp0(x)? evallgef(3) | evalvarn(v)
: evallgef(3) | evalvarn(v) | evalsigne(1);
y[2]=lcopy(x); return y;
}
static GEN
gtopoly0(GEN x, long v, int reverse)
{
long tx=typ(x),lx,i,j;
GEN y;
if (v<0) v = 0;
if (isexactzero(x)) return zeropol(v);
if (is_scalar_t(tx)) return scalarpol(x,v);
switch(tx)
{
case t_POL:
y=gcopy(x); break;
case t_SER:
y=gconvsp(x,1); break;
case t_RFRAC: case t_RFRACN:
y=gdeuc((GEN)x[1],(GEN)x[2]); break;
case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
lx=lg(x);
if (reverse)
{
while (lx-- && isexactzero((GEN)x[lx]));
i=lx+2; y=cgetg(i,t_POL);
y[1]=evallgef(i); if (!gcmp0(x)) y[1] |= evalsigne(1);
for (j=2; j<i; j++) y[j]=lcopy((GEN)x[j-1]);
}
else
{
i=1; j=lx; while (lx-- && isexactzero((GEN)x[i++]));
i=lx+2; y=cgetg(i,t_POL);
y[1]=evallgef(i); if (!gcmp0(x)) y[1] |= evalsigne(1);
lx = j-1;
for (j=2; j<i; j++) y[j]=lcopy((GEN)x[lx--]);
}
break;
default: err(typeer,"gtopoly");
}
setvarn(y,v); return y;
}
GEN
gtopolyrev(GEN x, long v) { return gtopoly0(x,v,1); }
GEN
gtopoly(GEN x, long v) { return gtopoly0(x,v,0); }
GEN
zeroser(long v, long val)
{
GEN x = cgetg(2,t_SER);
x[1]=evalvalp(val) | evalvarn(v); return x;
}
GEN
scalarser(GEN x, long v, long prec)
{
GEN y=cgetg(prec+2,t_SER);
long i;
y[1]=evalsigne(1) | evalvalp(0) | evalvarn(v);
y[2]=lcopy(x); for (i=3; i<=prec+1; i++) y[i]=zero;
return y;
}
GEN
gtoser(GEN x, long v)
{
long tx=typ(x),lx,i,j,l,av,tetpil;
GEN y,p1,p2;
if (v<0) v = 0;
if (tx==t_SER) { y=gcopy(x); setvarn(y,v); return y; }
if (isexactzero(x)) return zeroser(v,precdl);
if (is_scalar_t(tx)) return scalarser(x,v,precdl);
switch(tx)
{
case t_POL:
lx=lgef(x); i=2; while (i<lx && gcmp0((GEN)x[i])) i++;
l=lx-i; if (precdl>l) l=precdl;
y=cgetg(l+2,t_SER);
y[1] = evalsigne(1) | evalvalp(i-2) | evalvarn(v);
for (j=2; j<=lx-i+1; j++) y[j]=lcopy((GEN)x[j+i-2]);
for ( ; j<=l+1; j++) y[j]=zero;
break;
case t_RFRAC: case t_RFRACN:
av=avma; p1=gtoser((GEN)x[1],v); p2=gtoser((GEN)x[2],v);
tetpil=avma; return gerepile(av,tetpil,gdiv(p1,p2));
case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
lx=lg(x); i=1; while (i<lx && isexactzero((GEN)x[i])) i++;
y = cgetg(lx-i+2,t_SER);
y[1] = evalsigne(1) | evalvalp(i-1) | evalvarn(v);
for (j=2; j<=lx-i+1; j++) y[j]=lcopy((GEN)x[j+i-2]);
break;
default: err(typeer,"gtoser");
}
return y;
}
GEN
gtovec(GEN x)
{
long tx,lx,i;
GEN y;
if (!x) return cgetg(1,t_VEC);
tx = typ(x);
if (is_scalar_t(tx) || is_rfrac_t(tx))
{
y=cgetg(2,t_VEC); y[1]=lcopy(x);
return y;
}
if (is_graphicvec_t(tx))
{
lx=lg(x); y=cgetg(lx,t_VEC);
for (i=1; i<lx; i++) y[i]=lcopy((GEN)x[i]);
return y;
}
if (tx==t_POL)
{
lx=lgef(x); y=cgetg(lx-1,t_VEC);
for (i=1; i<=lx-2; i++) y[i]=lcopy((GEN)x[lx-i]);
return y;
}
if (tx==t_LIST)
{
lx=lgef(x); y=cgetg(lx-1,t_VEC); x++;
for (i=1; i<=lx-2; i++) y[i]=lcopy((GEN)x[i]);
return y;
}
if (!signe(x)) { y=cgetg(2,t_VEC); y[1]=zero; return y; }
lx=lg(x); y=cgetg(lx-1,t_VEC); x++;
for (i=1; i<=lx-2; i++) y[i]=lcopy((GEN)x[i]);
return y;
}
GEN
compo(GEN x, long n)
{
long l,tx=typ(x);
if (tx==t_POL && n+1 >= lgef(x)) return gzero;
if (tx==t_SER && !signe(x)) return gzero;
if (!is_recursive_t(tx))
err(talker, "this object doesn't have components (is a leaf)");
l=lontyp[tx]+n-1;
if (n<1 || l>=lg(x))
err(talker,"nonexistent component");
return gcopy((GEN)x[l]);
}
/* with respect to the main variable if v<0, with respect to the variable v
otherwise. v ignored if x is not a polynomial/series. */
GEN
polcoeff0(GEN x, long n, long v)
{
long tx=typ(x),lx,ex,w,av,tetpil;
GEN xinit;
if (is_scalar_t(tx)) return n? gzero: gcopy(x);
switch(tx)
{
case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
if (n<1 || n>=lg(x)) break;
return gcopy((GEN)x[n]);
case t_POL:
if (n<0) return gzero;
w=varn(x);
if (v < 0 || v == w)
return (n>=lgef(x)-2)? gzero: gcopy((GEN)x[n+2]);
if (v < w) return n? gzero: gcopy(x);
av=avma; xinit=x;
x=gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
if (gvar(x)) { avma=av; return n? gzero: gcopy(xinit); }
if (typ(x) == t_POL)
{
if (n>=lgef(x)-2) { avma=av; return gzero; }
x = (GEN)x[n+2];
}
else
x = polcoeff0(x, n, 0);
tetpil=avma; return gerepile(av,tetpil,gsubst(x,MAXVARN,polx[w]));
case t_SER:
w=varn(x);
if (v < 0 || v == w)
{
if (!signe(x)) return gzero;
lx=lg(x); ex=valp(x); if (n<ex) return gzero;
if (n>=ex+lx-2) break;
return gcopy((GEN)x[n-ex+2]);
}
if (v < w) return n? gzero: gcopy(x);
av=avma; xinit=x;
x=gsubst(gsubst(x,w,polx[MAXVARN]),v,polx[0]);
if (gvar(x)) { avma=av; return n? gzero: gcopy(xinit); }
if (gcmp0(x)) { avma=av; return gzero; }
if (typ(x) == t_SER)
{
lx=lg(x); ex=valp(x); if (n<ex) return gzero;
if (n>=ex+lx-2) break;
x = (GEN)x[n-ex+2];
}
else
x = polcoeff0(x, n, 0);
tetpil=avma; return gerepile(av,tetpil,gsubst(x,MAXVARN,polx[w]));
case t_RFRAC: case t_RFRACN:
w = precdl; av = avma;
if (v<0) v = gvar(x);
ex = ggval((GEN)x[2], polx[v]);
precdl = n + ex + 1; x = gtoser(x, v); precdl = w;
return gerepileupto(av, polcoeff0(x, n, v));
}
err(talker,"nonexistent component in truecoeff");
return NULL; /* not reached */
}
GEN
truecoeff(GEN x, long n)
{
return polcoeff0(x,n,-1);
}
GEN
denom(GEN x)
{
long lx,i,av,tetpil;
GEN s,t;
switch(typ(x))
{
case t_INT: case t_REAL: case t_INTMOD: case t_PADIC: case t_SER:
return gun;
case t_FRAC: case t_FRACN:
return absi((GEN)x[2]);
case t_COMPLEX:
av=avma; t=denom((GEN)x[1]); s=denom((GEN)x[2]); tetpil=avma;
return gerepile(av,tetpil,glcm(s,t));
case t_QUAD:
av=avma; t=denom((GEN)x[2]); s=denom((GEN)x[3]); tetpil=avma;
return gerepile(av,tetpil,glcm(s,t));
case t_POLMOD:
return denom((GEN)x[2]);
case t_RFRAC: case t_RFRACN:
return gcopy((GEN)x[2]);
case t_POL:
return polun[varn(x)];
case t_VEC: case t_COL: case t_MAT:
lx=lg(x); if (lx==1) return gun;
av = tetpil = avma; s = denom((GEN)x[1]);
for (i=2; i<lx; i++)
{
t = denom((GEN)x[i]);
/* t != gun est volontaire */
if (t != gun) { tetpil=avma; s=glcm(s,t); }
}
return gerepile(av,tetpil,s);
}
err(typeer,"denom");
return NULL; /* not reached */
}
GEN
numer(GEN x)
{
long av,tetpil;
GEN s;
switch(typ(x))
{
case t_INT: case t_REAL: case t_INTMOD:
case t_PADIC: case t_POL: case t_SER:
return gcopy(x);
case t_FRAC: case t_FRACN:
return (signe(x[2])>0)? gcopy((GEN)x[1]): gneg((GEN)x[1]);
case t_POLMOD:
av=avma; s=numer((GEN)x[2]); tetpil=avma;
return gerepile(av,tetpil,gmodulcp(s,(GEN)x[1]));
case t_RFRAC: case t_RFRACN:
return gcopy((GEN)x[1]);
case t_COMPLEX: case t_QUAD: case t_VEC: case t_COL: case t_MAT:
av=avma; s=denom(x); tetpil=avma;
return gerepile(av,tetpil,gmul(s,x));
}
err(typeer,"numer");
return NULL; /* not reached */
}
/* Lift only intmods if v does not occur in x, lift with respect to main
* variable of x if v < 0, with respect to variable v otherwise.
*/
GEN
lift0(GEN x, long v)
{
long lx,tx=typ(x),i;
GEN y;
switch(tx)
{
case t_INT: case t_REAL:
return gcopy(x);
case t_INTMOD:
return gcopy((GEN)x[2]);
case t_POLMOD:
if (v < 0 || v == varn((GEN)x[1])) return gcopy((GEN)x[2]);
y = cgetg(3,tx);
y[1] = (long)lift0((GEN)x[1],v);
y[2] = (long)lift0((GEN)x[2],v); return y;
case t_SER:
if (!signe(x)) return gcopy(x);
lx=lg(x); y=cgetg(lx,tx); y[1]=x[1];
for (i=2; i<lx; i++) y[i]=(long)lift0((GEN)x[i],v);
return y;
case t_FRAC: case t_FRACN: case t_COMPLEX: case t_RFRAC:
case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
lx=lg(x); y=cgetg(lx,tx);
for (i=1; i<lx; i++) y[i]=(long)lift0((GEN)x[i],v);
return y;
case t_POL:
y=cgetg(lx=lgef(x),tx); y[1]=x[1];
for (i=2; i<lx; i++) y[i]=(long)lift0((GEN)x[i],v);
return y;
case t_QUAD:
y=cgetg(4,tx); copyifstack(x[1],y[1]);
for (i=2; i<4; i++) y[i]=(long)lift0((GEN)x[i],v);
return y;
}
err(typeer,"lift");
return NULL; /* not reached */
}
GEN
lift(GEN x)
{
return lift0(x,-1);
}
/* same as lift, without copy. May DESTROY x. For internal use only.
Conventions on v as for lift. */
GEN
lift_intern0(GEN x, long v)
{
long i,lx,tx=typ(x);
switch(tx)
{
case t_INT: case t_REAL:
return x;
case t_INTMOD:
return (GEN)x[2];
case t_POLMOD:
if (v < 0 || v == varn((GEN)x[1])) return (GEN)x[2];
x[1]=(long)lift_intern0((GEN)x[1],v);
x[2]=(long)lift_intern0((GEN)x[2],v);
return x;
case t_SER: if (!signe(x)) return x; /* fall through */
case t_FRAC: case t_FRACN: case t_COMPLEX: case t_QUAD: case t_POL:
case t_RFRAC: case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
lx = (tx==t_POL)? lgef(x): lg(x);
for (i = lx-1; i>=lontyp[tx]; i--)
x[i] = (long) lift_intern0((GEN)x[i],v);
return x;
}
err(typeer,"lift_intern");
return NULL; /* not reached */
}
/* memes conventions pour v que lift */
GEN
centerlift0(GEN x, long v)
{
long lx,tx=typ(x),i,av;
GEN y;
switch(tx)
{
case t_INT:
return icopy(x);
case t_INTMOD:
av=avma; i=cmpii(shifti((GEN)x[2],1),(GEN)x[1]); avma=av;
return (i>0)? subii((GEN)x[2],(GEN)x[1]): icopy((GEN)x[2]);
case t_POLMOD:
if (v < 0 || v == varn((GEN)x[1])) return gcopy((GEN)x[2]);
y=cgetg(3,tx);
y[1]=(long)centerlift0((GEN)x[1],v); y[2]=(long)centerlift0((GEN)x[2],v);
return y;
case t_SER:
if (!signe(x)) return gcopy(x);
lx=lg(x); y=cgetg(lx,tx); y[1]=x[1];
for (i=2; i<lx; i++) y[i]=(long)centerlift0((GEN)x[i],v);
return y;
case t_POL:
case t_FRAC: case t_FRACN: case t_COMPLEX: case t_RFRAC:
case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
lx = (tx==t_POL)? lgef(x): lg(x);
y=cgetg(lx,tx); y[1]=x[1];
for (i=lontyp[tx]; i<lx; i++) y[i]=(long)centerlift0((GEN)x[i],v);
return y;
case t_QUAD:
y=cgetg(4,tx); copyifstack(x[1],y[1]);
for (i=2; i<4; i++) y[i]=(long)centerlift0((GEN)x[i],v);
return y;
}
err(typeer,"centerlift");
return NULL; /* not reached */
}
GEN
centerlift(GEN x)
{
return centerlift0(x,-1);
}
/*******************************************************************/
/* */
/* PARTIES REELLE ET IMAGINAIRES */
/* */
/*******************************************************************/
static GEN
op_ReIm(GEN f(GEN), GEN x)
{
long lx,i,j,av,tetpil, tx = typ(x);
GEN p1,p2,r2,i2,z;
switch(tx)
{
case t_POL:
lx=lgef(x); av=avma;
for (i=lx-1; i>=2; i--)
if (!gcmp0(f((GEN)x[i]))) break;
avma=av; if (i==1) return zeropol(varn(x));
z=cgetg(i+1,tx); z[1]=evalsigne(1)|evallgef(1+i)|evalvarn(varn(x));
for (j=2; j<=i; j++) z[j] = (long)f((GEN)x[j]);
return z;
case t_SER:
if (gcmp0(x)) { z=cgetg(2,tx); z[1]=x[1]; return z; }
lx=lg(x); av=avma;
for (i=2; i<lx; i++)
if (!gcmp0(f((GEN)x[i]))) break;
avma=av; if (i==lx) return zeroser(varn(x),lx-2+valp(x));
z=cgetg(lx-i+2,tx); z[1]=x[1]; setvalp(z, valp(x)+i-2);
for (j=2; i<lx; j++,i++) z[j] = (long) f((GEN)x[i]);
return z;
case t_RFRAC: case t_RFRACN:
{
av=avma; r2=greal((GEN)x[2]); i2=gimag((GEN)x[2]);
if (f==greal)
p1 = gadd(gmul(greal((GEN)x[1]),r2),
gmul(gimag((GEN)x[1]),i2));
else
p1 = gsub(gmul(gimag((GEN)x[1]),r2),
gmul(greal((GEN)x[1]),i2));
p2=gadd(gsqr(r2), gsqr(i2));
tetpil=avma; return gerepile(av,tetpil,gdiv(p1,p2));
}
case t_VEC: case t_COL: case t_MAT:
lx=lg(x); z=cgetg(lx,tx);
for (i=1; i<lx; i++) z[i] = (long) f((GEN)x[i]);
return z;
}
err(typeer,"greal/gimag");
return NULL; /* not reached */
}
GEN
greal(GEN x)
{
switch(typ(x))
{
case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
return gcopy(x);
case t_COMPLEX:
return gcopy((GEN)x[1]);
case t_QUAD:
return gcopy((GEN)x[2]);
}
return op_ReIm(greal,x);
}
GEN
gimag(GEN x)
{
switch(typ(x))
{
case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
return gzero;
case t_COMPLEX:
return gcopy((GEN)x[2]);
case t_QUAD:
return gcopy((GEN)x[3]);
}
return op_ReIm(gimag,x);
}
/*******************************************************************/
/* */
/* LOGICAL OPERATIONS */
/* */
/*******************************************************************/
GEN
glt(GEN x, GEN y) { return gcmp(x,y)<0? gun: gzero; }
GEN
gle(GEN x, GEN y) { return gcmp(x,y)<=0? gun: gzero; }
GEN
gge(GEN x, GEN y) { return gcmp(x,y)>=0? gun: gzero; }
GEN
ggt(GEN x, GEN y) { return gcmp(x,y)>0? gun: gzero; }
GEN
geq(GEN x, GEN y) { return gegal(x,y)? gun: gzero; }
GEN
gne(GEN x, GEN y) { return gegal(x,y)? gzero: gun; }
GEN
gand(GEN x, GEN y) { return gcmp0(x)? gzero: (gcmp0(y)? gzero: gun); }
GEN
gor(GEN x, GEN y) { return gcmp0(x)? (gcmp0(y)? gzero: gun): gun; }
GEN
gnot(GEN x) { return gcmp0(x)? gun: gzero; }
/*******************************************************************/
/* */
/* FORMAL SIMPLIFICATIONS */
/* */
/*******************************************************************/
GEN
geval(GEN x)
{
long av,tetpil,lx,i, tx = typ(x);
GEN y,z;
if (is_const_t(tx)) return gcopy(x);
if (is_graphicvec_t(tx) || tx == t_RFRACN)
{
lx=lg(x); y=cgetg(lx, tx);
for (i=1; i<lx; i++) y[i] = (long)geval((GEN)x[i]);
return y;
}
switch(tx)
{
case t_STR:
return flisexpr(GSTR(x));
case t_POLMOD: y=cgetg(3,tx);
y[1]=(long)geval((GEN)x[1]);
av=avma; z=geval((GEN)x[2]); tetpil=avma;
y[2]=lpile(av,tetpil,gmod(z,(GEN)y[1]));
return y;
case t_POL:
lx=lgef(x); if (lx==2) return gzero;
{
entree *ep = varentries[varn(x)];
if (!ep) return gcopy(x);
z = (GEN)ep->value;
}
y=gzero; av=avma;
for (i=lx-1; i>1; i--)
{
tetpil=avma;
y = gadd(geval((GEN)x[i]), gmul(z,y));
}
return gerepile(av,tetpil,y);
case t_SER:
err(impl, "evaluation of a power series");
case t_RFRAC:
return gdiv(geval((GEN)x[1]),geval((GEN)x[2]));
}
err(typeer,"geval");
return NULL; /* not reached */
}
GEN
simplify(GEN x)
{
long tx=typ(x),av,tetpil,i,lx;
GEN p1,y;
switch(tx)
{
case t_INT: case t_REAL: case t_FRAC:
return is_universal_constant(x)? x: gcopy(x);
case t_INTMOD: case t_PADIC: case t_QFR: case t_QFI:
case t_LIST: case t_STR: case t_VECSMALL:
return gcopy(x);
case t_FRACN:
return gred(x);
case t_COMPLEX:
if (is_universal_constant(x)) return x;
if (isexactzero((GEN)x[2])) return simplify((GEN)x[1]);
y=cgetg(3,tx);
y[1]=(long)simplify((GEN)x[1]);
y[2]=(long)simplify((GEN)x[2]); return y;
case t_QUAD:
if (isexactzero((GEN)x[3])) return simplify((GEN)x[2]);
y=cgetg(4,tx);
y[1]=lcopy((GEN)x[1]);
y[2]=(long)simplify((GEN)x[2]);
y[3]=(long)simplify((GEN)x[3]); return y;
case t_POLMOD: y=cgetg(3,tx);
y[1]=(long)simplify((GEN)x[1]);
av=avma; p1=gmod((GEN)x[2],(GEN)y[1]); tetpil=avma;
y[2]=lpile(av,tetpil,simplify(p1)); return y;
case t_POL:
lx=lgef(x); if (lx==2) return gzero;
if (lx==3) return simplify((GEN)x[2]);
y=cgetg(lx,tx); y[1]=x[1];
for (i=2; i<lx; i++) y[i]=(long)simplify((GEN)x[i]);
return y;
case t_SER:
if (!signe(x)) return gcopy(x);
lx=lg(x);
y=cgetg(lx,tx); y[1]=x[1];
for (i=2; i<lx; i++) y[i]=(long)simplify((GEN)x[i]);
return y;
case t_RFRAC: y=cgetg(3,tx);
y[1]=(long)simplify((GEN)x[1]);
y[2]=(long)simplify((GEN)x[2]); return y;
case t_RFRACN:
av=avma; y=gred_rfrac(x); tetpil=avma;
return gerepile(av,tetpil,simplify(y));
case t_VEC: case t_COL: case t_MAT:
lx=lg(x); y=cgetg(lx,tx);
for (i=1; i<lx; i++) y[i]=(long)simplify((GEN)x[i]);
return y;
}
err(typeer,"simplify");
return NULL; /* not reached */
}
/*******************************************************************/
/* */
/* EVALUATION OF SOME SIMPLE FORMAL OBJECTS */
/* */
/*******************************************************************/
static GEN
qfeval0_i(GEN q, GEN x, long n)
{
long i,j,av=avma;
GEN res=gzero;
for (i=2;i<n;i++)
for (j=1;j<i;j++)
res = gadd(res, gmul(gcoeff(q,i,j), mulii((GEN)x[i],(GEN)x[j])) );
res=gshift(res,1);
for (i=1;i<n;i++)
res = gadd(res, gmul(gcoeff(q,i,i), sqri((GEN)x[i])) );
return gerepileupto(av,res);
}
#if 0
static GEN
qfeval0(GEN q, GEN x, long n)
{
long i,j,av=avma;
GEN res=gzero;
for (i=2;i<n;i++)
for (j=1;j<i;j++)
res = gadd(res, gmul(gcoeff(q,i,j), gmul((GEN)x[i],(GEN)x[j])) );
res=gshift(res,1);
for (i=1;i<n;i++)
res = gadd(res, gmul(gcoeff(q,i,i), gsqr((GEN)x[i])) );
return gerepileupto(av,res);
}
#else
static GEN
qfeval0(GEN q, GEN x, long n)
{
long i,j,av=avma;
GEN res = gmul(gcoeff(q,1,1), gsqr((GEN)x[1]));
for (i=2; i<n; i++)
{
GEN l = (GEN)q[i];
GEN sx = gmul((GEN)l[1], (GEN)x[1]);
for (j=2; j<i; j++)
sx = gadd(sx, gmul((GEN)l[j],(GEN)x[j]));
sx = gadd(gshift(sx,1), gmul((GEN)l[i],(GEN)x[i]));
res = gadd(res, gmul((GEN)x[i], sx));
}
return gerepileupto(av,res);
}
#endif
/* We assume q is a real symetric matrix */
GEN
qfeval(GEN q, GEN x)
{
long n=lg(q);
if (n==1)
{
if (typ(q) != t_MAT || lg(x) != 1)
err(talker,"invalid data in qfeval");
return gzero;
}
if (typ(q) != t_MAT || lg(q[1]) != n)
err(talker,"invalid quadratic form in qfeval");
if (typ(x) != t_COL || lg(x) != n)
err(talker,"invalid vector in qfeval");
return qfeval0(q,x,n);
}
/* the Horner-type evaluation (mul x 2/3) would force us to use gmul and not
* mulii (more than 4 x slower for small entries). Not worth it.
*/
static GEN
qfbeval0_i(GEN q, GEN x, GEN y, long n)
{
long i,j,av=avma;
GEN res = gmul(gcoeff(q,1,1), mulii((GEN)x[1],(GEN)y[1]));
for (i=2;i<n;i++)
{
for (j=1;j<i;j++)
{
GEN p1 = addii(mulii((GEN)x[i],(GEN)y[j]), mulii((GEN)x[j],(GEN)y[i]));
res = gadd(res, gmul(gcoeff(q,i,j),p1));
}
res = gadd(res, gmul(gcoeff(q,i,i), mulii((GEN)x[i],(GEN)y[i])));
}
return gerepileupto(av,res);
}
#if 0
static GEN
qfbeval0(GEN q, GEN x, GEN y, long n)
{
long i,j,av=avma;
GEN res = gmul(gcoeff(q,1,1), gmul((GEN)x[1],(GEN)y[1]));
for (i=2;i<n;i++)
{
for (j=1;j<i;j++)
{
GEN p1 = gadd(gmul((GEN)x[i],(GEN)y[j]), gmul((GEN)x[j],(GEN)y[i]));
res = gadd(res, gmul(gcoeff(q,i,j),p1));
}
res = gadd(res, gmul(gcoeff(q,i,i), gmul((GEN)x[i],(GEN)y[i])));
}
return gerepileupto(av,res);
}
#else
static GEN
qfbeval0(GEN q, GEN x, GEN y, long n)
{
long i,j,av=avma;
GEN res = gmul(gcoeff(q,1,1), gmul((GEN)x[1], (GEN)y[1]));
for (i=2; i<n; i++)
{
GEN l = (GEN)q[i];
GEN sx = gmul((GEN)l[1], (GEN)y[1]);
GEN sy = gmul((GEN)l[1], (GEN)x[1]);
for (j=2; j<i; j++)
{
sx = gadd(sx, gmul((GEN)l[j],(GEN)y[j]));
sy = gadd(sy, gmul((GEN)l[j],(GEN)x[j]));
}
sx = gadd(sx, gmul((GEN)l[i],(GEN)y[i]));
sx = gmul((GEN)x[i], sx);
sy = gmul((GEN)y[i], sy);
res = gadd(res, gadd(sx,sy));
}
return gerepileupto(av,res);
}
#endif
/* We assume q is a real symetric matrix */
GEN
qfbeval(GEN q, GEN x, GEN y)
{
long n=lg(q);
if (n==1)
{
if (typ(q) != t_MAT || lg(x) != 1 || lg(y) != 1)
err(talker,"invalid data in qfbeval");
return gzero;
}
if (typ(q) != t_MAT || lg(q[1]) != n)
err(talker,"invalid quadratic form in qfbeval");
if (typ(x) != t_COL || lg(x) != n || typ(y) != t_COL || lg(y) != n)
err(talker,"invalid vector in qfbeval");
return qfbeval0(q,x,y,n);
}
/* yield X = M'.q.M, assuming q is symetric.
* X_ij are X_ji identical, not copies
* if flag is set, M has integer entries
*/
GEN
qf_base_change(GEN q, GEN M, int flag)
{
long i,j, k = lg(M), n=lg(q);
GEN res = cgetg(k,t_MAT);
GEN (*qf)(GEN,GEN,long) = flag? qfeval0_i: qfeval0;
GEN (*qfb)(GEN,GEN,GEN,long) = flag? qfbeval0_i: qfbeval0;
if (n==1)
{
if (typ(q) != t_MAT || k != 1)
err(talker,"invalid data in qf_base_change");
return res;
}
if (typ(M) != t_MAT || k == 1 || lg(M[1]) != n)
err(talker,"invalid base change matrix in qf_base_change");
for (i=1;i<k;i++)
{
res[i] = lgetg(k,t_COL);
coeff(res,i,i) = (long) qf(q,(GEN)M[i],n);
}
for (i=2;i<k;i++)
for (j=1;j<i;j++)
coeff(res,i,j)=coeff(res,j,i) = (long) qfb(q,(GEN)M[i],(GEN)M[j],n);
return res;
}
/* compute M'.M */
GEN
gram_matrix(GEN M)
{
long n=lg(M),i,j,k, av;
GEN res = cgetg(n,t_MAT),p1;
if (n==1)
{
if (typ(M) != t_MAT)
err(talker,"invalid data in gram_matrix");
return res;
}
if (typ(M) != t_MAT || lg(M[1]) != n)
err(talker,"not a square matrix in gram_matrix");
for (i=1;i<n;i++) res[i]=lgetg(n,t_COL);
av=avma;
for (i=1;i<n;i++)
{
p1 = gzero;
for (k=1;k<n;k++)
p1 = gadd(p1, gsqr(gcoeff(M,k,i)));
coeff(res,i,i) = (long) gerepileupto(av,p1);
av=avma;
}
for (i=2;i<n;i++)
for (j=1;j<i;j++)
{
p1=gzero;
for (k=1;k<n;k++)
p1 = gadd(p1, gmul(gcoeff(M,k,i),gcoeff(M,k,j)));
coeff(res,i,j)=coeff(res,j,i) = lpileupto(av,p1);
av=avma;
}
return res;
}
/* return Re(x * y), x and y scalars */
static GEN
mul_real(GEN x, GEN y)
{
if (typ(x) == t_COMPLEX)
{
if (typ(y) == t_COMPLEX)
{
long av=avma, tetpil;
GEN p1 = gmul((GEN)x[1], (GEN) y[1]);
GEN p2 = gneg(gmul((GEN)x[2], (GEN) y[2]));
tetpil=avma; return gerepile(av,tetpil,gadd(p1,p2));
}
x = (GEN)x[1];
}
else if (typ(y) == t_COMPLEX) y = (GEN)y[1];
return gmul(x,y);
}
/* Compute Re(x * y), x and y matrices of compatible dimensions
* assume lx, ly > 1, and scalar entries */
GEN
mulmat_real(GEN x, GEN y)
{
long av,tetpil,i,j,k, lx = lg(x), ly = lg(y), l = lg(x[1]);
GEN p1,p2, z = cgetg(ly,t_MAT);
for (j=1; j<ly; j++)
{
z[j] = lgetg(l,t_COL);
for (i=1; i<l; i++)
{
p1=gzero; av=avma;
for (k=1; k<lx; k++)
{
p2=mul_real(gcoeff(x,i,k),gcoeff(y,k,j));
tetpil=avma; p1=gadd(p1,p2);
}
coeff(z,i,j)=lpile(av,tetpil,p1);
}
}
return z;
}
static GEN
hqfeval0(GEN q, GEN x, long n)
{
long i,j, av=avma;
GEN res=gzero;
for (i=2;i<n;i++)
for (j=1;j<i;j++)
{
GEN p1 = gmul((GEN)x[i],gconj((GEN)x[j]));
res = gadd(res, mul_real(gcoeff(q,i,j),p1));
}
res=gshift(res,1);
for (i=1;i<n;i++)
res = gadd(res, gmul(gcoeff(q,i,i), gnorm((GEN)x[i])) );
return gerepileupto(av,res);
}
/* We assume q is a hermitian complex matrix */
GEN
hqfeval(GEN q, GEN x)
{
long n=lg(q);
if (n==1)
{
if (typ(q) != t_MAT || lg(x) != 1)
err(talker,"invalid data in hqfeval");
return gzero;
}
if (typ(q) != t_MAT || lg(q[1]) != n)
err(talker,"invalid quadratic form in hqfeval");
if (typ(x) != t_COL || lg(x) != n)
err(talker,"invalid vector in hqfeval");
return hqfeval0(q,x,n);
}
GEN
poleval(GEN x, GEN y)
{
long av,tetpil,i,j,imin,tx=typ(x);
GEN p1,p2,p3,r,s;
if (is_scalar_t(tx)) return gcopy(x);
switch(tx)
{
case t_POL:
i=lgef(x)-1; imin=2; break;
case t_RFRAC: case t_RFRACN: av=avma;
p1=poleval((GEN)x[1],y);
p2=poleval((GEN)x[2],y); tetpil=avma;
return gerepile(av,tetpil,gdiv(p1,p2));
case t_VEC: case t_COL:
i=lg(x)-1; imin=1; break;
default: err(typeer,"poleval");
}
if (i<=imin)
return (i==imin)? gcopy((GEN)x[imin]): gzero;
av=avma; p1=(GEN)x[i]; i--;
if (typ(y)!=t_COMPLEX)
{
#if 0 /* standard Horner's rule */
for ( ; i>=imin; i--)
p1 = gadd(gmul(p1,y),(GEN)x[i]);
#endif
/* specific attention to sparse polynomials */
for ( ; i>=imin; i=j-1)
{
for (j=i; gcmp0((GEN)x[j]); j--)
if (j==imin)
{
if (i!=j) y = gpuigs(y,i-j+1);
tetpil=avma; return gerepile(av,tetpil,gmul(p1,y));
}
r = (i==j)? y: gpuigs(y,i-j+1);
p1 = gadd(gmul(p1,r), (GEN)x[j]);
}
return gerepileupto(av,p1);
}
p2=(GEN)x[i]; i--; r=gtrace(y); s=gneg_i(gnorm(y));
for ( ; i>=imin; i--)
{
p3=gadd(p2,gmul(r,p1));
p2=gadd((GEN)x[i],gmul(s,p1)); p1=p3;
}
p1=gmul(y,p1); tetpil=avma;
return gerepile(av,tetpil,gadd(p1,p2));
}
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gen3.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