/* $Id: gen3.c,v 1.239.2.3 2007/02/13 22:30:25 kb Exp $
Copyright (C) 2000 The PARI group.
This file is part of the PARI/GP package.
PARI/GP is free software; you can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software
Foundation. It is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY WHATSOEVER.
Check the License for details. You should have received a copy of it, along
with the package; see the file 'COPYING'. If not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
/********************************************************************/
/** **/
/** GENERIC OPERATIONS **/
/** (third part) **/
/** **/
/********************************************************************/
#include "pari.h"
#include "paripriv.h"
/********************************************************************/
/** **/
/** PRINCIPAL VARIABLE NUMBER **/
/** **/
/********************************************************************/
long
gvar(GEN x)
{
long i, v, w;
switch(typ(x))
{
case t_POL: case t_SER: return varn(x);
case t_POLMOD: return varn(gel(x,1));
case t_RFRAC: return varn(gel(x,2));
case t_VEC: case t_COL: case t_MAT:
v = BIGINT;
for (i=1; i < lg(x); i++) { w=gvar(gel(x,i)); if (w<v) v=w; }
return v;
case t_VECSMALL:
case t_STR:
case t_LIST:
pari_err(typeer, "gvar");
}
return BIGINT;
}
/* T main polynomial in R[X], A auxiliary in R[X] (possibly degree 0).
* Guess and return the main variable of R */
static long
var2_aux(GEN T, GEN A)
{
long a = gvar2(T);
long b = (typ(A) == t_POL && varn(A) == varn(T))? gvar2(A): gvar(A);
if (a < b) a = b;
return a;
}
static long
var2_rfrac(GEN x) { return var2_aux(gel(x,2), gel(x,1)); }
static long
var2_polmod(GEN x) { return var2_aux(gel(x,1), gel(x,2)); }
/* main variable of x, with the convention that the "natural" main
* variable of a POLMOD is mute, so we want the next one. */
static long
gvar9(GEN x)
{
return (typ(x) == t_POLMOD)? var2_polmod(x): gvar(x);
}
/* main variable of the ring over wich x is defined */
long
gvar2(GEN x)
{
long i, v, w;
switch(typ(x))
{
case t_POLMOD:
return var2_polmod(x);
case t_POL: case t_SER:
v = BIGINT;
for (i=2; i < lg(x); i++) { w = gvar9(gel(x,i)); if (w<v) v=w; }
return v;
case t_RFRAC:
return var2_rfrac(x);
case t_VEC: case t_COL: case t_MAT:
v = BIGINT;
for (i=1; i < lg(x); i++) { w = gvar2(gel(x,i)); if (w<v) v=w; }
return v;
}
return BIGINT;
}
GEN
gpolvar(GEN x)
{
long v;
if (typ(x)==t_PADIC) return gcopy( gel(x,2) );
v = gvar(x);
if (v==BIGINT) pari_err(typeer,"gpolvar");
return pol_x[v];
}
/*******************************************************************/
/* */
/* PRECISION OF SCALAR OBJECTS */
/* */
/*******************************************************************/
static long
prec0(long e) { return (e < 0)? 2 - (e >> TWOPOTBITS_IN_LONG): 2; }
static long
precREAL(GEN x) { return signe(x) ? lg(x): prec0(expo(x)); }
/* t t_REAL, s an exact non-complex type. Return precision(|t| + |s|) */
static long
precrealexact(GEN t, GEN s) {
long l, e = gexpo(s);
if (e == -(long)HIGHEXPOBIT) return precREAL(t);
if (e < 0) e = 0;
e -= expo(t);
if (!signe(t)) return prec0(-e);
l = lg(t);
return (e > 0)? l + (e >> TWOPOTBITS_IN_LONG): l;
}
long
precision(GEN z)
{
long tx = typ(z), e, ex, ey, lz, lx, ly;
if (tx == t_REAL) return precREAL(z);
if (tx == t_COMPLEX)
{ /* ~ precision(|x| + |y|) */
GEN x = gel(z,1), y = gel(z,2);
if (typ(x) != t_REAL) {
if (typ(y) != t_REAL) return 0;
return precrealexact(y, x);
}
if (typ(y) != t_REAL) return precrealexact(x, y);
/* x, y are t_REALs, cf addrr_sign */
ex = expo(x);
ey = expo(y);
e = ey - ex;
if (!signe(x)) {
if (!signe(y)) return prec0( min(ex,ey) );
if (e < 0) return prec0(ex);
lz = 3 + (e >> TWOPOTBITS_IN_LONG);
ly = lg(y); if (lz > ly) lz = ly;
return lz;
}
if (!signe(y)) {
if (e > 0) return prec0(ey);
lz = 3 + ((-e)>>TWOPOTBITS_IN_LONG);
lx = lg(x); if (lz > lx) lz = lx;
return lz;
}
if (e < 0) { swap(x, y); e = -e; }
lx = lg(x);
ly = lg(y);
if (e) {
long d = e >> TWOPOTBITS_IN_LONG, l = ly-d;
return (l > lx)? lx + d: ly;
}
return min(lx, ly);
}
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: case t_VEC: case t_COL: case t_MAT:
k=VERYBIGINT;
for (i=lontyp[tx]; i<lx; i++)
{
l = gprecision(gel(x,i));
if (l && l<k) k = l;
}
return (k==VERYBIGINT)? 0: k;
case t_RFRAC:
{
k=gprecision(gel(x,1));
l=gprecision(gel(x,2)); if (l && (!k || l<k)) k=l;
return k;
}
case t_QFR:
return gprecision(gel(x,4));
}
return 0;
}
GEN
ggprecision(GEN x)
{
long a = gprecision(x);
return utoipos(a ? prec2ndec(a): VERYBIGINT);
}
GEN
precision0(GEN x, long n) { return n? gprec(x,n): 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:
return VERYBIGINT;
case t_INTMOD:
return Z_pval(gel(x,1),p);
case t_PADIC:
if (!gequal(gel(x,2),p))
pari_err(talker,"not the same prime in padicprec");
return precp(x)+valp(x);
case t_POL:
case t_COMPLEX: case t_QUAD: case t_POLMOD: case t_SER: case t_RFRAC:
case t_VEC: case t_COL: case t_MAT:
for (s=VERYBIGINT, i=lontyp[tx]; i<lx; i++)
{
t = padicprec(gel(x,i),p); if (t<s) s = t;
}
return s;
}
pari_err(typeer,"padicprec");
return 0; /* not reached */
}
#define DEGREE0 -VERYBIGINT
/* Degree of x (scalar, t_POL, t_RFRAC) wrt variable v if v >= 0,
* wrt to main variable if v < 0.
*/
long
poldegree(GEN x, long v)
{
long tx = typ(x), lx,w,i,d;
if (is_scalar_t(tx)) return gcmp0(x)? DEGREE0: 0;
switch(tx)
{
case t_POL:
if (!signe(x)) return DEGREE0;
w = varn(x);
if (v < 0 || v == w) return degpol(x);
if (v < w) return 0;
lx = lg(x); d = DEGREE0;
for (i=2; i<lx; i++)
{
long e = poldegree(gel(x,i), v);
if (e > d) d = e;
}
return d;
case t_RFRAC:
if (gcmp0(gel(x,1))) return DEGREE0;
return poldegree(gel(x,1),v) - poldegree(gel(x,2),v);
}
pari_err(typeer,"degree");
return 0; /* not reached */
}
long
degree(GEN x)
{
return poldegree(x,-1);
}
/* If v<0, leading coeff with respect to the main variable, otherwise wrt v. */
GEN
pollead(GEN x, long v)
{
long l, tx = typ(x), w;
pari_sp av;
GEN xinit;
if (is_scalar_t(tx)) return gcopy(x);
w = varn(x);
switch(tx)
{
case t_POL:
if (v < 0 || v == w)
{
l=lg(x);
return (l==2)? gen_0: gcopy(gel(x,l-1));
}
break;
case t_SER:
if (v < 0 || v == w) return signe(x)? gcopy(gel(x,2)): gen_0;
break;
default:
pari_err(typeer,"pollead");
return NULL; /* not reached */
}
if (v < w) return gcopy(x);
av = avma; xinit = x;
x = gsubst(gsubst(x,w,pol_x[MAXVARN]),v,pol_x[0]);
if (gvar(x)) { avma = av; return gcopy(xinit);}
tx = typ(x);
if (tx == t_POL) {
l = lg(x); if (l == 2) { avma = av; return gen_0; }
x = gel(x,l-1);
}
else if (tx == t_SER) {
if (!signe(x)) { avma = av; return gen_0;}
x = gel(x,2);
} else pari_err(typeer,"pollead");
return gerepileupto(av, gsubst(x,MAXVARN,pol_x[w]));
}
/* 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_const_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_POLMOD:
return isinexactreal(gel(x,1)) || isinexactreal(gel(x,2));
}
if (is_noncalc_t(tx)) return 0;
lx = lg(x);
for (i=lontyp[tx]; i<lx; i++)
if (isinexactreal(gel(x,i))) return 1;
return 0;
}
/* returns 1 if there's an inexact component in the structure, and
* 0 otherwise. */
int
isinexact(GEN x)
{
long tx = typ(x), lx, i;
switch(tx)
{
case t_REAL: case t_PADIC: case t_SER:
return 1;
case t_INT: case t_INTMOD: case t_FRAC: case t_QFR: case t_QFI:
return 0;
case t_COMPLEX: case t_QUAD: case t_RFRAC: case t_POLMOD:
return isinexact(gel(x,1)) || isinexact(gel(x,2));
case t_POL: case t_VEC: case t_COL: case t_MAT:
lx = lg(x);
for (i=lontyp[tx]; i<lx; i++)
if (isinexact(gel(x,i))) return 1;
return 0;
case t_LIST:
lx = lgeflist(x);
for (i=lontyp[tx]; i<lx; i++)
if (isinexact(gel(x,i))) return 1;
return 0;
}
return 0;
}
int
isexactzeroscalar(GEN g)
{
switch (typ(g))
{
case t_INT:
return !signe(g);
case t_INTMOD: case t_POLMOD:
return isexactzeroscalar(gel(g,2));
case t_FRAC: case t_RFRAC:
return isexactzeroscalar(gel(g,1));
case t_COMPLEX:
return isexactzeroscalar(gel(g,1)) && isexactzeroscalar(gel(g,2));
case t_QUAD:
return isexactzeroscalar(gel(g,2)) && isexactzeroscalar(gel(g,3));
case t_POL: return lg(g) == 2;
}
return 0;
}
int
isexactzero(GEN g)
{
long i;
switch (typ(g))
{
case t_INT:
return !signe(g);
case t_INTMOD: case t_POLMOD:
return isexactzero(gel(g,2));
case t_COMPLEX:
return isexactzero(gel(g,1)) && isexactzero(gel(g,2));
case t_QUAD:
return isexactzero(gel(g,2)) && isexactzero(gel(g,3));
case t_POL: return lg(g) == 2;
case t_VEC: case t_COL: case t_MAT:
for (i=lg(g)-1; i; i--)
if (!isexactzero(gel(g,i))) return 0;
return 1;
}
return 0;
}
int
iscomplex(GEN x)
{
switch(typ(x))
{
case t_INT: case t_REAL: case t_FRAC:
return 0;
case t_COMPLEX:
return !gcmp0(gel(x,2));
case t_QUAD:
return signe(gmael(x,1,2)) > 0;
}
pari_err(typeer,"iscomplex");
return 0; /* not reached */
}
int
ismonome(GEN x)
{
long i;
if (typ(x)!=t_POL || !signe(x)) return 0;
for (i=lg(x)-2; i>1; i--)
if (!isexactzero(gel(x,i))) return 0;
return 1;
}
/*******************************************************************/
/* */
/* GENERIC REMAINDER */
/* */
/*******************************************************************/
/* euclidean quotient for scalars of admissible types */
static GEN
_quot(GEN x, GEN y)
{
GEN q = gdiv(x,y), f = gfloor(q);
if (gsigne(y) < 0 && !gequal(f,q)) f = gadd(f,gen_1);
return f;
}
static GEN
_quotgs(GEN x, long y)
{
GEN q = gdivgs(x,y), f = gfloor(q);
if (y < 0 && !gequal(f,q)) f = gadd(f,gen_1);
return f;
}
static GEN
quot(GEN x, GEN y)
{
pari_sp av = avma;
return gerepileupto(av, _quot(x, y));
}
static GEN
quotgs(GEN x, long y)
{
pari_sp av = avma;
return gerepileupto(av, _quotgs(x, y));
}
GEN
gmod(GEN x, GEN y)
{
pari_sp av, tetpil;
long i,lx,ty, tx = typ(x);
GEN z,p1;
if (is_matvec_t(tx))
{
lx=lg(x); z=cgetg(lx,tx);
for (i=1; i<lx; i++) gel(z,i) = gmod(gel(x,i),y);
return z;
}
ty = typ(y);
switch(ty)
{
case t_INT:
switch(tx)
{
case t_INT:
return modii(x,y);
case t_INTMOD: z=cgetg(3,tx);
gel(z,1) = gcdii(gel(x,1),y);
gel(z,2) = modii(gel(x,2),gel(z,1)); return z;
case t_FRAC:
av=avma;
p1=mulii(gel(x,1),Fp_inv(gel(x,2),y));
tetpil=avma; return gerepile(av,tetpil,modii(p1,y));
case t_QUAD: z=cgetg(4,tx);
gel(z,1) = gcopy(gel(x,1));
gel(z,2) = gmod(gel(x,2),y);
gel(z,3) = gmod(gel(x,3),y); return z;
case t_PADIC: return padic_to_Fp(x, y);
case t_POLMOD: case t_POL:
return gen_0;
case t_REAL: /* NB: conflicting semantic with lift(x * Mod(1,y)). */
av = avma;
return gerepileupto(av, gadd(x, gneg(gmul(_quot(x,y),y))));
default: pari_err(operf,"%",x,y);
}
case t_REAL: case t_FRAC:
switch(tx)
{
case t_INT: case t_REAL: case t_FRAC:
av = avma;
return gerepileupto(av, gsub(x, gmul(_quot(x,y),y)));
case t_POLMOD: case t_POL:
return gen_0;
default: pari_err(operf,"%",x,y);
}
case t_POL:
if (is_scalar_t(tx))
{
if (tx!=t_POLMOD || varncmp(varn(x[1]), varn(y)) > 0)
return degpol(y)? gcopy(x): gen_0;
if (varn(x[1])!=varn(y)) return gen_0;
z=cgetg(3,t_POLMOD);
gel(z,1) = ggcd(gel(x,1),y);
gel(z,2) = grem(gel(x,2),gel(z,1)); return z;
}
switch(tx)
{
case t_POL:
return grem(x,y);
case t_RFRAC:
av=avma;
p1=gmul(gel(x,1),ginvmod(gel(x,2),y)); tetpil=avma;
return gerepile(av,tetpil,grem(p1,y));
case t_SER:
if (ismonome(y) && varn(x) == varn(y))
{
long d = degpol(y);
if (lg(x)-2 + valp(x) < d) pari_err(operi,"%",x,y);
av = avma;
return gerepileupto(av, gmod(ser2rfrac_i(x), y));
}
default: pari_err(operf,"%",x,y);
}
}
pari_err(operf,"%",x,y);
return NULL; /* not reached */
}
GEN
gmodgs(GEN x, long y)
{
ulong u;
long i, lx, tx = typ(x);
GEN z;
if (is_matvec_t(tx))
{
lx=lg(x); z=cgetg(lx,tx);
for (i=1; i<lx; i++) gel(z,i) = gmodgs(gel(x,i),y);
return z;
}
switch(tx)
{
case t_INT:
return modis(x,y);
case t_INTMOD: z=cgetg(3,tx);
i = cgcd(smodis(gel(x,1), y), y);
gel(z,1) = utoi(i);
gel(z,2) = modis(gel(x,2), i); return z;
case t_FRAC:
u = (ulong)labs(y);
return utoi( Fl_div(umodiu(gel(x,1), u),
umodiu(gel(x,2), u), u) );
case t_QUAD: z=cgetg(4,tx);
gel(z,1) = gcopy(gel(x,1));
gel(z,2) = gmodgs(gel(x,2),y);
gel(z,3) = gmodgs(gel(x,3),y); return z;
case t_PADIC: return padic_to_Fp(x, stoi(y));
case t_POLMOD: case t_POL:
return gen_0;
}
pari_err(operf,"%",x,stoi(y));
return NULL; /* not reached */
}
GEN
gmodulsg(long x, GEN y)
{
GEN z;
switch(typ(y))
{
case t_INT: z = cgetg(3,t_INTMOD);
gel(z,1) = absi(y);
gel(z,2) = modsi(x,y); return z;
case t_POL: z = cgetg(3,t_POLMOD);
gel(z,1) = gcopy(y);
gel(z,2) = stoi(x); return z;
}
pari_err(operf,"%",stoi(x),y); return NULL; /* not reached */
}
GEN
gmodulss(long x, long y)
{
GEN z = cgetg(3,t_INTMOD);
y = labs(y);
gel(z,1) = stoi(y);
gel(z,2) = modss(x, y); return z;
}
static GEN
specialmod(GEN x, GEN y)
{
GEN z = gmod(x,y);
if (varncmp(gvar(z), varn(y)) < 0) z = gen_0;
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++) gel(z,i) = gmodulo(gel(x,i),y);
return z;
}
switch(typ(y))
{
case t_INT: z = cgetg(3,t_INTMOD);
gel(z,1) = absi(y);
gel(z,2) = Rg_to_Fp(x,y); return z;
case t_POL: z = cgetg(3,t_POLMOD);
gel(z,1) = gcopy(y);
if (is_scalar_t(tx))
{
gel(z,2) = (lg(y) > 3)? gcopy(x): gmod(x,y);
return z;
}
if (tx!=t_POL && tx != t_RFRAC && tx!=t_SER) break;
gel(z,2) = specialmod(x,y); return z;
}
pari_err(operf,"%",x,y); return NULL; /* not reached */
}
GEN
Mod0(GEN x,GEN y,long flag)
{
switch(flag)
{
case 0:
case 1: return gmodulo(x,y);
default: pari_err(flagerr,"Mod");
}
return NULL; /* not reached */
}
/*******************************************************************/
/* */
/* GENERIC EUCLIDEAN DIVISION */
/* */
/*******************************************************************/
GEN
gdivent(GEN x, GEN y)
{
long tx = typ(x);
if (is_matvec_t(tx))
{
long i, lx = lg(x);
GEN z = cgetg(lx,tx);
for (i=1; i<lx; i++) gel(z,i) = gdivent(gel(x,i),y);
return z;
}
switch(typ(y))
{
case t_INT:
switch(tx)
{ /* equal to, but more efficient than, quot(x,y) */
case t_INT: return truedivii(x,y);
case t_REAL: case t_FRAC: return quot(x,y);
case t_POL: return gdiv(x,y);
}
break;
case t_REAL: case t_FRAC: return quot(x,y);
case t_POL:
if (is_scalar_t(tx))
{
if (tx == t_POLMOD) break;
return degpol(y)? gen_0: gdiv(x,y);
}
if (tx == t_POL) return gdeuc(x,y);
}
pari_err(operf,"\\",x,y);
return NULL; /* not reached */
}
GEN
gdiventgs(GEN x, long y)
{
long tx = typ(x);
if (is_matvec_t(tx))
{
long i, lx = lg(x);
GEN z = cgetg(lx,tx);
for (i=1; i<lx; i++) gel(z,i) = gdiventgs(gel(x,i),y);
return z;
}
switch(tx)
{ /* equal to, but more efficient than, quotgs(x,y) */
case t_INT: return truedivis(x,y);
case t_REAL: case t_FRAC: return quotgs(x,y);
case t_POL: return gdivgs(x,y);
}
pari_err(operf,"\\",x,stoi(y));
return NULL; /* not reached */
}
/* with remainder */
static GEN
quotrem(GEN x, GEN y, GEN *r)
{
pari_sp av;
GEN q = quot(x,y);
av = avma;
*r = gerepileupto(av, gsub(x, gmul(q,y)));
return q;
}
GEN
gdiventres(GEN x, GEN y)
{
long tx = typ(x);
GEN z,q,r;
if (is_matvec_t(tx))
{
long i, lx = lg(x);
z = cgetg(lx,tx);
for (i=1; i<lx; i++) gel(z,i) = gdiventres(gel(x,i),y);
return z;
}
z = cgetg(3,t_COL);
switch(typ(y))
{
case t_INT:
switch(tx)
{ /* equal to, but more efficient than next case */
case t_INT:
gel(z,1) = truedvmdii(x,y,(GEN*)(z+2));
return z;
case t_REAL: case t_FRAC:
q = quotrem(x,y,&r);
gel(z,1) = q;
gel(z,2) = r; return z;
case t_POL:
gel(z,1) = gdiv(x,y);
gel(z,2) = gen_0; return z;
}
break;
case t_REAL: case t_FRAC:
q = quotrem(x,y,&r);
gel(z,1) = q;
gel(z,2) = r; return z;
case t_POL:
if (is_scalar_t(tx))
{
if (tx == t_POLMOD) break;
if (degpol(y))
{
q = gen_0;
r = gcopy(x);
}
else
{
q = gdiv(x,y);
r = gen_0;
}
gel(z,1) = q;
gel(z,2) = r; return z;
}
if (tx == t_POL)
{
gel(z,1) = poldivrem(x,y,(GEN *)(z+2));
return z;
}
}
pari_err(operf,"\\",x,y);
return NULL; /* not reached */
}
GEN
divrem(GEN x, GEN y, long v)
{
pari_sp av = avma;
long vx, vy;
GEN q, r;
if (v < 0 || typ(y) != t_POL || typ(x) != t_POL) return gdiventres(x,y);
vx = varn(x); if (vx != v) x = swap_vars(x,v);
vy = varn(y); if (vy != v) y = swap_vars(y,v);
q = poldivrem(x,y, &r);
if (v && (vx != v || vy != v))
{
q = gsubst(q, v, pol_x[v]); /* poleval broken for t_RFRAC, subst is safe */
r = gsubst(r, v, pol_x[v]);
}
return gerepilecopy(av, mkcol2(q, r));
}
static int
is_scal(long t) { return t==t_INT || t==t_FRAC; }
GEN
diviiround(GEN x, GEN y)
{
pari_sp av1, av = avma;
GEN q,r;
int fl;
q = dvmdii(x,y,&r); /* q = x/y rounded towards 0, sgn(r)=sgn(x) */
if (r==gen_0) 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) q = gerepileuptoint(av, addis(q,sz));
}
return q;
}
/* If x and y are not both scalars, same as gdivent.
* Otherwise, compute the quotient x/y, rounded to the nearest integer
* (towards +oo in case of tie). */
GEN
gdivround(GEN x, GEN y)
{
pari_sp av1, av;
long tx=typ(x),ty=typ(y);
GEN q,r;
int fl;
if (tx==t_INT && ty==t_INT) return diviiround(x,y);
av = avma;
if (is_scal(tx) && is_scal(ty))
{ /* same as diviiround but less efficient */
q = quotrem(x,y,&r);
av1 = avma;
fl = gcmp(gmul2n(gabs(r,0),1), gabs(y,0));
avma = av1; cgiv(r);
if (fl >= 0) /* If 2*|r| >= |y| */
{
long sz = gsigne(y);
if (fl || sz > 0) q = gerepileupto(av, gaddgs(q, sz));
}
return q;
}
if (is_matvec_t(tx))
{
long i,lx = lg(x);
GEN z = cgetg(lx,tx);
for (i=1; i<lx; i++) gel(z,i) = gdivround(gel(x,i),y);
return z;
}
return gdivent(x,y);
}
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 gen_0; }
pari_err(typeer,"gdivmod");
}
if (tx!=t_POL) pari_err(typeer,"gdivmod");
return poldivrem(x,y,pr);
}
/*******************************************************************/
/* */
/* SHIFT */
/* */
/*******************************************************************/
/* Shift tronque si n<0 (multiplication tronquee par 2^n) */
GEN
gshift(GEN x, long n)
{
long i, 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);
for (i=1; i<lx; i++) gel(y,i) = gshift(gel(x,i),n);
return y;
}
return gmul2n(x,n);
}
/*******************************************************************/
/* */
/* INVERSE */
/* */
/*******************************************************************/
GEN
mpinv(GEN b)
{
long i, l1, l = lg(b), e = expo(b), s = signe(b);
GEN x = cgetr(l), a = mpcopy(b);
double t;
a[1] = evalexpo(0) | evalsigne(1);
for (i = 3; i < l; i++) x[i] = 0;
t = (((double)HIGHBIT) * HIGHBIT) / (double)(ulong)a[2];
if (((ulong)t) & HIGHBIT)
x[1] = evalexpo(0) | evalsigne(1);
else {
t *= 2;
x[1] = evalexpo(-1) | evalsigne(1);
}
x[2] = (ulong)t;
l1 = 1; l -= 2;
while (l1 < l)
{
l1 <<= 1; if (l1 > l) l1 = l;
setlg(a, l1 + 2);
setlg(x, l1 + 2);
/* TODO: mulrr(a,x) should be a half product (the higher half is known).
* mulrr(x, ) already is */
affrr(addrr(x, mulrr(x, subsr(1, mulrr(a,x)))), x);
avma = (pari_sp)a;
}
x[1] = evalexpo(expo(x)-e) | evalsigne(s);
avma = (pari_sp)x; return x;
}
GEN
inv_ser(GEN b)
{
pari_sp av = avma, av2, lim;
long i, j, le, l = lg(b), e = valp(b), v = varn(b);
GEN E, y, x = cgetg(l, t_SER), a = shallowcopy(b);
if (!signe(b)) pari_err(gdiver);
for (j = 3; j < l; j++) gel(x,j) = gen_0;
gel(x,2) = ginv(gel(b,2));
a[1] = x[1] = evalvalp(0) | evalvarn(v) | evalsigne(1);
E = Newton_exponents(l - 2);
av2 = avma; lim = stack_lim(av2, 2);
le = lg(E)-1;
for (i = le; i > 1; i--) {
long l1 = E[i-1], l2 = E[i];
setlg(a, l1 + 2);
setlg(x, l1 + 2);
/* TODO: gmul(a,x) should be a half product (the higher half is known) */
y = gmul(x, gsubsg(1, gmul(a,x))) - l2;
for (j = l2+2; j < l1+2; j++) x[j] = y[j];
if (low_stack(lim, stack_lim(av2,2)))
{
if(DEBUGMEM>1) pari_warn(warnmem,"inv_ser");
y = gerepilecopy(av2, x);
for (j = 2; j < l1+2; j++) x[j] = y[j];
}
}
x[1] = evalvalp(valp(x)-e) | evalvarn(v) | evalsigne(1);
return gerepilecopy(av, x);
}
GEN
ginv(GEN x)
{
long s;
pari_sp av, tetpil;
GEN X, z, y, p1, p2;
switch(typ(x))
{
case t_INT:
if (is_pm1(x)) return icopy(x);
s = signe(x); if (!s) pari_err(gdiver);
z = cgetg(3,t_FRAC);
gel(z,1) = s<0? gen_m1: gen_1;
gel(z,2) = absi(x); return z;
case t_REAL:
return divsr(1,x);
case t_INTMOD: z=cgetg(3,t_INTMOD);
gel(z,1) = icopy(gel(x,1));
gel(z,2) = Fp_inv(gel(x,2),gel(x,1)); return z;
case t_FRAC:
s = signe(x[1]); if (!s) pari_err(gdiver);
if (is_pm1(x[1]))
return s>0? icopy(gel(x,2)): negi(gel(x,2));
z = cgetg(3,t_FRAC);
gel(z,1) = icopy(gel(x,2));
gel(z,2) = icopy(gel(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])) pari_err(gdiver);
z[1] = evalprecp(precp(x)) | evalvalp(-valp(x));
gel(z,2) = icopy(gel(x,2));
gel(z,3) = icopy(gel(x,3));
gel(z,4) = Fp_inv(gel(x,4),gel(z,3)); return z;
case t_POLMOD: z = cgetg(3,t_POLMOD);
X = gel(x,1); gel(z,1) = gcopy(X);
if (degpol(X) == 2) { /* optimized for speed */
av = avma;
gel(z,2) = gerepileupto(av, gdiv(quad_polmod_conj(gel(x,2), X),
quad_polmod_norm(gel(x,2), X)) );
}
else gel(z,2) = ginvmod(gel(x,2), X);
return z;
case t_POL: return gred_rfrac_simple(gen_1,x);
case t_SER: return gdiv(gen_1,x);
case t_RFRAC:
{
GEN n = gel(x,1), d = gel(x,2);
pari_sp av = avma, ltop;
if (gcmp0(n)) pari_err(gdiver);
n = simplify_i(n);
if (typ(n) != t_POL || varn(n) != varn(d))
{
if (gcmp1(n)) { avma = av; return gcopy(d); }
ltop = avma;
z = RgX_Rg_div(d,n);
} else {
ltop = avma;
z = cgetg(3,t_RFRAC);
gel(z,1) = gcopy(d);
gel(z,2) = gcopy(n);
}
stackdummy(av, ltop);
return z;
}
case t_QFR:
av = avma; z = cgetg(5, t_QFR);
gel(z,1) = gel(x,1);
gel(z,2) = negi( gel(x,2) );
gel(z,3) = gel(x,3);
gel(z,4) = negr( gel(x,4) );
return gerepileupto(av, redreal(z));
case t_QFI:
y = gcopy(x);
if (!equalii(gel(x,1),gel(x,2)) && !equalii(gel(x,1),gel(x,3)))
setsigne(y[2],-signe(y[2]));
return y;
case t_MAT:
return (lg(x)==1)? cgetg(1,t_MAT): invmat(x);
case t_VECSMALL:
{
long i,lx = lg(x);
y = cgetg(lx,t_VECSMALL);
for (i=1; i<lx; i++)
{
long xi=x[i];
if (xi<1 || xi>=lx) pari_err(talker,"incorrect permtuation to inverse");
y[xi] = i;
}
return y;
}
}
pari_err(typeer,"inverse");
return NULL; /* not reached */
}
/*******************************************************************/
/* */
/* SUBSTITUTION DANS UN POLYNOME OU UNE SERIE */
/* */
/*******************************************************************/
/* Convert t_SER --> t_POL, ignoring valp. INTERNAL ! */
GEN
ser2pol_i(GEN x, long lx)
{
long i = lx-1;
GEN y;
while (i > 1 && isexactzero(gel(x,i))) i--;
y = cgetg(i+1, t_POL); y[1] = x[1] & ~VALPBITS;
for ( ; i > 1; i--) y[i] = x[i];
return y;
}
/*
subst_poly(pol, from, to) =
{ local(t='subst_poly_t, M);
\\ if fraction
M = numerator(from) - t * denominator(from);
\\ else
M = from - t;
subst(pol % M, t, to)
}
*/
GEN
gsubst_expr(GEN pol, GEN from, GEN to)
{
pari_sp av = avma;
long v = fetch_var(); /* XXX Need fetch_var_low_priority() */
GEN tmp;
from = simplify_i(from);
switch (typ(from)) {
case t_RFRAC: /* M= numerator(from) - t * denominator(from) */
tmp = gsub(gel(from,1), gmul(pol_x[v], gel(from,2)));
break;
default:
tmp = gsub(from, pol_x[v]); /* M = from - t */
}
if (v <= gvar(from)) pari_err(talker, "subst: unexpected variable precedence");
tmp = gmul(pol, mkpolmod(gen_1, tmp));
if (typ(tmp) == t_POLMOD)
tmp = gel(tmp,2); /* optimize lift */
else /* Vector? */
tmp = lift0(tmp, gvar(from));
tmp = gsubst(tmp, v, to);
(void)delete_var();
return gerepilecopy(av, tmp);
}
GEN
gsubstpol(GEN x, GEN T, GEN y)
{
pari_sp av;
long d, v;
GEN deflated;
if (typ(T) != t_POL || !ismonome(T) || !gcmp1(leading_term(T)))
return gsubst_expr(x,T,y);
d = degpol(T); v = varn(T);
if (d == 1) return gsubst(x, v, y);
av = avma;
CATCH(cant_deflate) {
avma = av;
return gsubst_expr(x,T,y);
} TRY {
deflated = gdeflate(x, v, d);
} ENDCATCH
return gerepilecopy(av, gsubst(deflated, v, y));
}
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, i, j, k, jb;
pari_sp av, lim;
GEN t,p1,p2,z;
if (ty==t_MAT)
{
if (ly==1) return cgetg(1,t_MAT);
if (ly != lg(y[1]))
pari_err(talker,"forbidden substitution by a non square matrix");
} else if (is_graphicvec_t(ty))
pari_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(gel(x,1),v,y); vx=varn(p1);
p2=gsubst(gel(x,2),v,y); vy=gvar(p2);
if (typ(p1)!=t_POL)
pari_err(talker,"forbidden substitution in a scalar type");
if (varncmp(vy, vx) >= 0) return gerepileupto(av, gmodulo(p2,p1));
lx = lg(p2);
z = cgetg(lx,t_POL); z[1] = p2[1];
for (i=2; i<lx; i++) gel(z,i) = gmodulo(gel(p2,i),p1);
return gerepileupto(av, normalizepol_i(z,lx));
}
switch(tx)
{
case t_POL:
if (lx==2)
return (ty==t_MAT)? gscalmat(gen_0,ly-1): gen_0;
vx = varn(x);
if (varncmp(vx, v) < 0)
{
if (varncmp(gvar(y), vx) > 0)
{ /* easy special case */
z = cgetg(lx, t_POL); z[1] = x[1];
for (i=2; i<lx; i++) gel(z,i) = gsubst(gel(x,i),v,y);
return normalizepol_i(z,lx);
}
/* general case */
av=avma; p1=pol_x[vx]; z= gsubst(gel(x,lx-1),v,y);
for (i=lx-1; i>2; i--) z=gadd(gmul(z,p1),gsubst(gel(x,i-1),v,y));
return gerepileupto(av,z);
}
/* v <= vx */
if (ty!=t_MAT)
return varncmp(vx,v) > 0 ? gcopy(x): poleval(x,y);
if (varncmp(vx, v) > 0) return gscalmat(x,ly-1);
if (lx==3) return gscalmat(gel(x,2),ly-1);
av=avma; z=gel(x,lx-1);
for (i=lx-1; i>2; i--) z=gaddmat(gel(x,i-1),gmul(z,y));
return gerepileupto(av,z);
case t_SER:
vx = varn(x);
if (varncmp(vx, v) > 0) return (ty==t_MAT)? gscalmat(x,ly-1): gcopy(x);
ex = valp(x);
if (varncmp(vx, v) < 0)
{ /* FIXME: improve this */
av = avma; p1 = ser2pol_i(x, lx);
z = tayl(gsubst(p1,v,y), vx, lx-2);
if (ex) z = gmul(z, monomial(gen_1,ex,vx));
return gerepileupto(av, 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));
if (lg(x) == 2) return zeroser(vy, ey*ex);
if (vy != vx)
{
av = avma; z = zeroser(vy,0);
for (i=lx-1; i>=2; i--) z = gadd(gel(x,i), gmul(y,z));
if (ex) z = gmul(z, gpowgs(y,ex));
return gerepileupto(av,z);
}
l = (lx-2)*ey+2;
if (ex) { if (l>ly) l = ly; }
else if (lx != 3)
{
long l2;
for (i = 3; i < lx; i++)
if (!isexactzero(gel(x,i))) break;
l2 = (i-2)*ey + (gcmp0(y)? 2 : ly);
if (l > l2) l = l2;
}
p2 = ex? gpowgs(y, ex): NULL;
av = avma; lim=stack_lim(av,1);
t = shallowcopy(y);
if (l < ly) setlg(t, l);
z = scalarser(gel(x,2),varn(y),l-2);
for (i=3,jb=ey; jb<=l-2; i++,jb+=ey)
{
if (i < lx) {
for (j=jb+2; j<min(l, jb+ly); j++)
gel(z,j) = gadd(gel(z,j), gmul(gel(x,i),gel(t,j-jb)));
}
for (j=l-1-jb-ey; j>1; j--)
{
p1 = gen_0;
for (k=2; k<j; k++)
p1 = gadd(p1, gmul(gel(t,j-k+2),gel(y,k)));
gel(t,j) = gadd(p1, gmul(gel(t,2),gel(y,j)));
}
if (low_stack(lim, stack_lim(av,1)))
{
if(DEBUGMEM>1) pari_warn(warnmem,"gsubst");
gerepileall(av,2, &z,&t);
}
}
if (!p2) return gerepilecopy(av,z);
return gerepileupto(av, gmul(z,p2));
case t_POL: case t_RFRAC:
if (isexactzero(y)) return scalarser(gel(x,2),v,lx-2);
vy = gvar(y); e = gval(y,vy);
if (e <= 0)
pari_err(talker,"non positive valuation in a series substitution");
av = avma; p1 = gsubst(ser2pol_i(x, lg(x)), v, y);
z = gmul(gpowgs(y, ex), tayl(p1, vy, e*(lx-2)));
return gerepileupto(av, z);
default:
pari_err(talker,"non polynomial or series type substituted in a series");
}
break;
case t_RFRAC: av=avma;
p1=gsubst(gel(x,1),v,y);
p2=gsubst(gel(x,2),v,y); return gerepileupto(av, gdiv(p1,p2));
case t_VEC: case t_COL: case t_MAT: z=cgetg(lx,tx);
for (i=1; i<lx; i++) gel(z,i) = gsubst(gel(x,i),v,y);
return z;
}
return gcopy(x);
}
GEN
gsubstvec(GEN e, GEN v, GEN r)
{
pari_sp ltop=avma;
long i,l=lg(v);
GEN w,z;
if ( !is_vec_t(typ(v)) || !is_vec_t(typ(r)) )
pari_err(typeer,"substvec");
if (lg(r)!=l)
pari_err(talker,"different number of variables and values in substvec");
w=cgetg(l,t_VECSMALL);
z=cgetg(l,t_VECSMALL);
for(i=1;i<l;i++)
{
GEN T=(GEN)v[i];
if (typ(T) != t_POL || !ismonome(T) || !gcmp1(leading_term(T)))
pari_err(talker,"not a variable in substvec");
w[i]=varn(T);
z[i]=fetch_var();
}
for(i=1;i<l;i++) e = gsubst(e,w[i],pol_x[z[i]]);
for(i=1;i<l;i++) e = gsubst(e,z[i],(GEN)r[i]);
for(i=1;i<l;i++) (void)delete_var();
return gerepileupto(ltop,e);
}
/*******************************************************************/
/* */
/* SERIE RECIPROQUE D'UNE SERIE */
/* */
/*******************************************************************/
GEN
recip(GEN x)
{
long v=varn(x), lx = lg(x);
pari_sp tetpil, av=avma;
GEN p1,p2,a,y,u;
if (typ(x)!=t_SER) pari_err(talker,"not a series in serreverse");
if (valp(x)!=1 || lx < 3)
pari_err(talker,"valuation not equal to 1 in serreverse");
a=gel(x,2);
if (gcmp1(a))
{
long i, j, k, mi;
pari_sp lim=stack_lim(av, 2);
mi = lx-1; while (mi>=3 && gcmp0(gel(x,mi))) mi--;
u = cgetg(lx,t_SER);
y = cgetg(lx,t_SER);
u[1] = y[1] = evalsigne(1) | evalvalp(1) | evalvarn(v);
gel(u,2) = gel(y,2) = gen_1;
if (lx > 3)
{
gel(u,3) = gmulsg(-2,gel(x,3));
gel(y,3) = gneg(gel(x,3));
}
for (i=3; i<lx-1; )
{
pari_sp av2;
for (j=3; j<i+1; j++)
{
av2 = avma; p1 = gel(x,j);
for (k = max(3,j+2-mi); k < j; k++)
p1 = gadd(p1, gmul(gel(u,k),gel(x,j-k+2)));
p1 = gneg(p1);
gel(u,j) = gerepileupto(av2, gadd(gel(u,j), p1));
}
av2 = avma;
p1 = gmulsg(i,gel(x,i+1));
for (k = 2; k < min(i,mi); k++)
{
p2 = gmul(gel(x,k+1),gel(u,i-k+2));
p1 = gadd(p1, gmulsg(k,p2));
}
i++;
gel(u,i) = gerepileupto(av2, gneg(p1));
gel(y,i) = gdivgs(gel(u,i), i-1);
if (low_stack(lim, stack_lim(av,2)))
{
if(DEBUGMEM>1) pari_warn(warnmem,"recip");
for(k=i+1; k<lx; k++) gel(u,k) = gel(y,k) = gen_0; /* dummy */
gerepileall(av,2, &u,&y);
}
}
return gerepilecopy(av,y);
}
y = gdiv(x,a); gel(y,2) = gen_1; y = recip(y);
a = gdiv(pol_x[v],a); tetpil = avma;
return gerepile(av,tetpil, gsubst(y,v,a));
}
/*******************************************************************/
/* */
/* DERIVATION ET INTEGRATION */
/* */
/*******************************************************************/
GEN
derivpol(GEN x)
{
long i,lx = lg(x)-1;
GEN y;
if (lx<3) return zeropol(varn(x));
y = cgetg(lx,t_POL);
for (i=2; i<lx ; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
y[1] = x[1]; return normalizepol_i(y,i);
}
GEN
derivser(GEN x)
{
long i, vx = varn(x), e = valp(x), lx = lg(x);
GEN y;
if (lx == 2) return zeroser(vx,e? e-1: 0);
if (e)
{
y = cgetg(lx,t_SER); y[1] = evalvalp(e-1) | evalvarn(vx);
for (i=2; i<lx; i++) gel(y,i) = gmulsg(i+e-2,gel(x,i));
} else {
if (lx == 3) return zeroser(vx, 0);
lx--;
y = cgetg(lx,t_SER); y[1] = evalvalp(0) | evalvarn(vx);
for (i=2; i<lx; i++) gel(y,i) = gmulsg(i-1,gel(x,i+1));
}
return normalize(y);
}
GEN
deriv(GEN x, long v)
{
long lx, vx, tx, i, j;
pari_sp av;
GEN y;
tx = typ(x); if (is_const_t(tx)) return gen_0;
if (v < 0) v = gvar9(x);
switch(tx)
{
case t_POLMOD:
if (v<=varn(x[1])) return gen_0;
y = cgetg(3,t_POLMOD);
gel(y,1) = gcopy(gel(x,1));
gel(y,2) = deriv(gel(x,2),v); return y;
case t_POL:
vx = varn(x);
if (varncmp(vx, v) > 0) return gen_0;
if (varncmp(vx, v) == 0) return derivpol(x);
lx = lg(x); y = cgetg(lx,t_POL);
y[1] = x[1];
for (i=2; i<lx; i++) gel(y,i) = deriv(gel(x,i),v);
return normalizepol_i(y,i);
case t_SER:
vx = varn(x);
if (varncmp(vx, v) > 0) return gen_0;
if (varncmp(vx, v) == 0) return derivser(x);
lx = lg(x); y = cgetg(lx, t_SER);
y[1] = x[1];
for (j=2; j<lx; j++) gel(y,j) = deriv(gel(x,j),v);
return normalize(y);
case t_RFRAC: {
GEN a = gel(x,1), b = gel(x,2), bp, b0, d, t;
y = cgetg(3,t_RFRAC); av = avma;
bp = deriv(b, v);
d = ggcd(bp, b);
if (gcmp1(d)) {
d = gadd(gmul(b, deriv(a,v)), gmul(gneg_i(a), bp));
gel(y,1) = gerepileupto(av, d);
gel(y,2) = gsqr(b); return y;
}
b0 = gdivexact(b, d);
bp = gdivexact(bp,d);
a = gadd(gmul(b0, deriv(a,v)), gmul(gneg_i(a), bp));
t = ggcd(a, d);
if (!gcmp1(t)) { a = gdivexact(a, t); d = gdivexact(d, t); }
gel(y,1) = a;
gel(y,2) = gmul(d, gsqr(b0));
return gerepilecopy((pari_sp)(y+3), y);
}
case t_VEC: case t_COL: case t_MAT: lx=lg(x); y=cgetg(lx,tx);
for (i=1; i<lx; i++) gel(y,i) = deriv(gel(x,i),v);
return y;
}
pari_err(typeer,"deriv");
return NULL; /* not reached */
}
/********************************************************************/
/** **/
/** TAYLOR SERIES **/
/** **/
/********************************************************************/
static GEN
tayl_vec(long v, long vx) {
GEN y = cgetg(v+2,t_VEC);
long i;
for (i=0; i<v; i++) gel(y,i+1) = pol_x[i];
gel(y,vx+1) = pol_x[v];
gel(y,v+1) = pol_x[vx]; return y;
}
GEN
tayl(GEN x, long v, long precS)
{
long vx = gvar9(x);
pari_sp av = avma;
GEN y, t;
if (v <= vx) return gadd(zeroser(v,precS),x);
y = tayl_vec(v, vx);
t = tayl(changevar(x,y), vx,precS);
return gerepileupto(av, changevar(t,y));
}
GEN
ggrando(GEN x, long n)
{
long m, v;
switch(typ(x))
{
case t_INT:/* bug 3 + O(1). We suppose x is a truc() */
if (signe(x) <= 0) pari_err(talker,"non-positive argument in O()");
if (!is_pm1(x)) return zeropadic(x,n);
/* +/-1 = x^0 */
v = m = 0; break;
case t_POL:
if (!signe(x)) pari_err(talker,"zero argument in O()");
v = varn(x); if ((ulong)v > MAXVARN) pari_err(talker,"incorrect object in O()");
m = n * polvaluation(x, NULL); break;
case t_RFRAC:
if (!gcmp0((GEN)x[1])) pari_err(talker,"zero argument in O()");
v = gvar(x); if ((ulong)v > MAXVARN) pari_err(talker,"incorrect object in O()");
m = n * gval(x,v); break;
default: pari_err(talker,"incorrect argument in O()");
v = m = 0; /* not reached */
}
return zeroser(v,m);
}
/*******************************************************************/
/* */
/* FORMAL INTEGRATION */
/* */
/*******************************************************************/
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++) gel(y,i) = integ(gel(x,i),v);
return y;
}
GEN
integ(GEN x, long v)
{
long lx, tx, e, i, vx, n;
pari_sp av = avma;
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);
gel(y,1) = gcopy(gel(x,1));
gel(y,2) = integ(gel(x,2),v); return y;
}
if (gcmp0(x)) return gen_0;
y = cgetg(4,t_POL);
y[1] = evalsigne(1) | evalvarn(v);
gel(y,2) = gen_0;
gel(y,3) = gcopy(x); return y;
}
switch(tx)
{
case t_POL:
vx = varn(x); lx = lg(x);
if (varncmp(vx, v) < 0) v = vx;
if (lx == 2) return zeropol(v);
if (varncmp(vx, v) > 0)
{
y = cgetg(4,t_POL);
y[1] = x[1];
gel(y,2) = gen_0;
gel(y,3) = gcopy(x); return y;
}
if (varncmp(vx, v) < 0) return triv_integ(x,v,tx,lx);
y = cgetg(lx+1,tx); y[1] = x[1]; gel(y,2) = gen_0;
for (i=3; i<=lx; i++) gel(y,i) = gdivgs(gel(x,i-1),i-2);
return y;
case t_SER:
lx = lg(x); vx = varn(x); e = valp(x);
if (lx == 2)
{
if (vx == v) e++; else if (varncmp(vx, v) < 0) v = vx;
return zeroser(v, e);
}
if (varncmp(vx, v) > 0)
{
y = cgetg(4,t_POL);
y[1] = evalvarn(v) | evalsigne(1);
gel(y,2) = gen_0;
gel(y,3) = gcopy(x); return y;
}
if (varncmp(vx, v) < 0) return triv_integ(x,v,tx,lx);
y = cgetg(lx,tx);
for (i=2; i<lx; i++)
{
long j = i+e-1;
if (!j)
{ /* should be isexactzero, but try to avoid error */
if (gcmp0(gel(x,i))) { gel(y,i) = gen_0; continue; }
pari_err(talker, "a log appears in intformal");
}
else gel(y,i) = gdivgs(gel(x,i),j);
}
y[1] = evalsigne(1) | evalvarn(vx) | evalvalp(e+1); return y;
case t_RFRAC:
vx = gvar(x);
if (varncmp(vx, v) > 0)
{
y=cgetg(4,t_POL);
y[1] = signe(x[1])? evalvarn(v) | evalsigne(1)
: evalvarn(v);
gel(y,2) = gen_0;
gel(y,3) = gcopy(x); return y;
}
if (varncmp(vx, v) < 0)
{
p1 = tayl_vec(v, vx);
y = integ(changevar(x,p1),vx);
return gerepileupto(av, changevar(y,p1));
}
tx = typ(x[1]); i = is_scalar_t(tx)? 0: degpol(x[1]);
tx = typ(x[2]); n = is_scalar_t(tx)? 0: degpol(x[2]);
n = i+n + 2;
y = gdiv(gtrunc(gmul(gel(x,2), integ(tayl(x,v,n),v))), gel(x,2));
if (!gequal(deriv(y,v),x)) pari_err(talker,"a log/atan appears in intformal");
if (typ(y)==t_RFRAC && lg(y[1]) == lg(y[2]))
{
GEN p2;
tx=typ(y[1]); p1=is_scalar_t(tx)? gel(y,1): leading_term(gel(y,1));
tx=typ(y[2]); p2=is_scalar_t(tx)? gel(y,2): leading_term(gel(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++) gel(y,i) = integ(gel(x,i),v);
return y;
}
pari_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 floorr(x);
case t_FRAC: return truedivii(gel(x,1),gel(x,2));
case t_RFRAC: return gdeuc(gel(x,1),gel(x,2));
case t_VEC: case t_COL: case t_MAT:
lx = lg(x); y = cgetg(lx,tx);
for (i=1; i<lx; i++) gel(y,i) = gfloor(gel(x,i));
return y;
}
pari_err(typeer,"gfloor");
return NULL; /* not reached */
}
GEN
gfrac(GEN x)
{
pari_sp av = avma, tetpil;
GEN p1 = gneg_i(gfloor(x));
tetpil = avma; return gerepile(av,tetpil,gadd(x,p1));
}
/* assume x t_REAL */
GEN
ceilr(GEN x) {
pari_sp av = avma;
GEN y = floorr(x);
if (cmpri(x, y)) return gerepileuptoint(av, addsi(1,y));
return y;
}
GEN
gceil(GEN x)
{
GEN y, p1;
long i, lx, tx = typ(x);
pari_sp av;
switch(tx)
{
case t_INT: case t_POL: return gcopy(x);
case t_REAL: return ceilr(x);
case t_FRAC:
av = avma; y = dvmdii(gel(x,1),gel(x,2),&p1);
if (p1 != gen_0 && gsigne(x) > 0)
{
cgiv(p1);
return gerepileuptoint(av, addsi(1,y));
}
return y;
case t_RFRAC:
return gdeuc(gel(x,1),gel(x,2));
case t_VEC: case t_COL: case t_MAT:
lx = lg(x); y = cgetg(lx,tx);
for (i=1; i<lx; i++) gel(y,i) = gceil(gel(x,i));
return y;
}
pari_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);
}
/* assume x a t_REAL */
GEN
roundr(GEN x)
{
long ex, s = signe(x);
pari_sp av;
GEN t;
if (!s || (ex=expo(x)) < -1) return gen_0;
if (ex == -1) return s>0? gen_1:
absrnz_egal2n(x)? gen_0: gen_m1;
av = avma;
t = real2n(-1, nbits2prec(ex+1)); /* = 0.5 */
return gerepileuptoint(av, floorr( addrr(x,t) ));
}
GEN
ground(GEN x)
{
GEN y;
long i, lx, tx=typ(x);
pari_sp av;
switch(tx)
{
case t_INT: case t_INTMOD: case t_QUAD: return gcopy(x);
case t_REAL: return roundr(x);
case t_FRAC: return diviiround(gel(x,1), gel(x,2));
case t_POLMOD: y=cgetg(3,t_POLMOD);
gel(y,1) = gcopy(gel(x,1));
gel(y,2) = ground(gel(x,2)); return y;
case t_COMPLEX:
av = avma; y = cgetg(3, t_COMPLEX);
gel(y,2) = ground(gel(x,2));
if (!signe(y[2])) { avma = av; return ground(gel(x,1)); }
gel(y,1) = ground(gel(x,1)); return y;
case t_POL:
y = init_gen_op(x, tx, &lx, &i);
for (; i<lx; i++) gel(y,i) = ground(gel(x,i));
return normalizepol_i(y, lx);
case t_SER:
y = init_gen_op(x, tx, &lx, &i);
for (; i<lx; i++) gel(y,i) = ground(gel(x,i));
return normalize(y);
case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
y = init_gen_op(x, tx, &lx, &i);
for (; i<lx; i++) gel(y,i) = ground(gel(x,i));
return y;
}
pari_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, ex, e1;
pari_sp av;
*e = -(long)HIGHEXPOBIT;
switch(tx)
{
case t_INT: case t_INTMOD: case t_QUAD: return gcopy(x);
case t_FRAC: return diviiround(gel(x,1), gel(x,2));
case t_REAL:
ex = expo(x);
if (!signe(x) || ex < -1) { *e = ex; return gen_0; }
av = avma;
p1 = addrr(real2n(-1,nbits2prec(ex+2)), x); e1 = expo(p1);
if (e1 < 0)
{
if (signe(p1) >= 0) { *e = ex; avma = av; return gen_0; }
*e = expo(addsr(1,x)); avma = av; return gen_m1;
}
lx = lg(x);
e1 = e1 - bit_accuracy(lx) + 1;
y = ishiftr_lg(p1, lx, e1);
if (signe(x) < 0) y = addsi(-1,y);
y = gerepileuptoint(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);
gel(y,1) = gcopy(gel(x,1));
gel(y,2) = grndtoi(gel(x,2), e); return y;
case t_COMPLEX:
av = avma; y = cgetg(3, t_COMPLEX);
gel(y,2) = grndtoi(gel(x,2), e);
if (!signe(y[2])) {
avma = av;
y = grndtoi(gel(x,1), &e1);
}
else
gel(y,1) = grndtoi(gel(x,1), &e1);
if (e1 > *e) *e = e1;
return y;
case t_POL:
y = init_gen_op(x, tx, &lx, &i);
for (; i<lx; i++)
{
gel(y,i) = grndtoi(gel(x,i),&e1);
if (e1 > *e) *e = e1;
}
return normalizepol_i(y, lx);
case t_SER:
y = init_gen_op(x, tx, &lx, &i);
for (; i<lx; i++)
{
gel(y,i) = grndtoi(gel(x,i),&e1);
if (e1 > *e) *e = e1;
}
return normalize(y);
case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
y = init_gen_op(x, tx, &lx, &i);
for (; i<lx; i++)
{
gel(y,i) = grndtoi(gel(x,i),&e1);
if (e1 > *e) *e = e1;
}
return y;
}
pari_err(typeer,"grndtoi");
return NULL; /* not reached */
}
/* floor(x * 2^s) */
GEN
gfloor2n(GEN x, long s)
{
GEN z;
switch(typ(x))
{
case t_INT:
return shifti(x, s);
case t_REAL:
return ishiftr(x, s);
case t_COMPLEX:
z = cgetg(3, t_COMPLEX);
gel(z,1) = gfloor2n(gel(x,1), s);
gel(z,2) = gfloor2n(gel(x,2), s);
return z;
default: pari_err(typeer,"gfloor2n");
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, e1;
GEN y;
if (tx == t_REAL)
{
ex = expo(x); if (ex < 0) { *e = ex; return gen_0; }
lx = lg(x); e1 = ex - bit_accuracy(lx) + 1;
y = ishiftr_lg(x, lx, e1);
if (e1 <= 0) { pari_sp av = avma; e1 = expo(subri(x,y)); avma = av; }
*e = e1; return y;
}
*e = -(long)HIGHEXPOBIT;
if (is_matvec_t(tx))
{
lx = lg(x); y = cgetg(lx,tx);
for (i=1; i<lx; i++)
{
gel(y,i) = gcvtoi(gel(x,i),&e1);
if (e1 > *e) *e = e1;
}
return y;
}
return gtrunc(x);
}
int
isint(GEN n, GEN *ptk)
{
switch(typ(n))
{
case t_INT: *ptk = n; return 1;
case t_REAL: {
pari_sp av0 = avma;
GEN z = floorr(n);
pari_sp av = avma;
long s = signe(subri(n, z));
if (s) { avma = av0; return 0; }
*ptk = z; avma = av; return 1;
}
case t_FRAC: return 0;
case t_COMPLEX: return gcmp0(gel(n,2)) && isint(gel(n,1),ptk);
case t_QUAD: return gcmp0(gel(n,3)) && isint(gel(n,2),ptk);
default: pari_err(typeer,"isint"); return 0; /* not reached */
}
}
int
issmall(GEN n, long *ptk)
{
pari_sp av = avma;
GEN z;
long k;
if (!isint(n, &z)) return 0;
k = itos_or_0(z); avma = av;
if (k || lgefint(z) == 2) { *ptk = k; return k; }
return 0;
}
/* smallest integer greater than any incarnations of the real x
* [avoid mpfloor() and "precision loss in truncation"] */
GEN
ceil_safe(GEN x)
{
pari_sp av = avma;
long e;
GEN y = gcvtoi(x,&e);
if (e < 0) e = 0;
y = addii(y, int2n(e));
return gerepileuptoint(av, y);
}
GEN
ser2rfrac_i(GEN x)
{
long e = valp(x);
GEN a = ser2pol_i(x, lg(x));
if (e) {
if (e > 0) a = RgX_shift_shallow(a, e);
else a = gred_rfrac_simple(a, monomial(gen_1, -e, varn(a)));
}
return a;
}
static GEN
ser2rfrac(GEN x)
{
pari_sp av = avma;
return gerepilecopy(av, ser2rfrac_i(x));
}
GEN
gtrunc(GEN x)
{
long tx=typ(x), i, v;
pari_sp av;
GEN y;
switch(tx)
{
case t_INT: case t_POL:
return gcopy(x);
case t_REAL:
return mptrunc(x);
case t_FRAC:
return divii(gel(x,1),gel(x,2));
case t_PADIC:
if (!signe(x[4])) return gen_0;
v = valp(x);
if (!v) return gcopy(gel(x,4));
if (v>0)
{ /* here p^v is an integer */
av = avma; y = powiu(gel(x,2),v);
return gerepileuptoint(av, mulii(y,gel(x,4)));
}
y=cgetg(3,t_FRAC);
gel(y,1) = icopy(gel(x,4));
gel(y,2) = gpowgs(gel(x,2),-v);
return y;
case t_RFRAC:
return gdeuc(gel(x,1),gel(x,2));
case t_SER:
return ser2rfrac(x);
case t_VEC: case t_COL: case t_MAT:
{
long lx = lg(x); y = cgetg(lx,tx);
for (i=1; i<lx; i++) gel(y,i) = gtrunc(gel(x,i));
return y;
}
}
pari_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 --> INT, POL & SER */
/* */
/*******************************************************************/
/* return a_(n-1) B^(n-1) + ... + a_0, where B = 2^32.
* The a_i are 32bits integers */
GEN
mkintn(long n, ...)
{
va_list ap;
GEN x, y;
long i;
#ifdef LONG_IS_64BIT
long e = (n&1);
n = (n+1) >> 1;
#endif
va_start(ap,n);
x = cgeti(n+2);
x[1] = evallgefint(n+2) | evalsigne(1);
y = int_MSW(x);
for (i=0; i <n; i++)
{
#ifdef LONG_IS_64BIT
ulong a = (e && !i)? 0: va_arg(ap, long);
ulong b = va_arg(ap, long);
*y = (a << 32) | b;
#else
*y = va_arg(ap, long);
#endif
y = int_precW(y);
}
va_end(ap);
return int_normalize(x, 0);
}
/* 2^32 a + b */
GEN
u2toi(ulong a, ulong b)
{
GEN x;
if (!a && !b) return gen_0;
#ifdef LONG_IS_64BIT
x = cgeti(3);
x[1] = evallgefint(3)|evalsigne(1);
x[2] = ((a << 32) | b);
#else
if (a) {
x = cgeti(4);
x[1] = evallgefint(4)|evalsigne(1);
*(int_MSW(x)) = (long)a;
*(int_LSW(x)) = (long)b;
} else {
x = cgeti(3);
x[1] = evallgefint(3)|evalsigne(1);
x[2] = (long)b;
}
#endif
return x;
}
/* return a_(n-1) x^(n-1) + ... + a_0 */
GEN
mkpoln(long n, ...)
{
va_list ap;
GEN x, y;
long i;
va_start(ap,n);
x = cgetg(n+2, t_POL); y = x + 2;
x[1] = evalvarn(0);
for (i=n-1; i >= 0; i--) gel(y,i) = va_arg(ap, GEN);
va_end(ap); return normalizepol(x);
}
/* return [a_1, ..., a_n] */
GEN
mkvecn(long n, ...)
{
va_list ap;
GEN x;
long i;
va_start(ap,n);
x = cgetg(n+1, t_VEC);
for (i=1; i <= n; i++) gel(x,i) = va_arg(ap, GEN);
va_end(ap); return x;
}
GEN
mkcoln(long n, ...)
{
va_list ap;
GEN x;
long i;
va_start(ap,n);
x = cgetg(n+1, t_COL);
for (i=1; i <= n; i++) gel(x,i) = va_arg(ap, GEN);
va_end(ap); return x;
}
GEN
scalarpol(GEN x, long v)
{
GEN y;
if (isexactzero(x)) return zeropol(v);
y = cgetg(3,t_POL);
y[1] = gcmp0(x)? evalvarn(v)
: evalvarn(v) | evalsigne(1);
gel(y,2) = gcopy(x); return y;
}
/* deg1pol(a,b,x)=a*x+b, assumes a != 0 */
GEN
deg1pol(GEN x1, GEN x0,long v)
{
GEN x = cgetg(4,t_POL);
x[1] = evalsigne(1) | evalvarn(v);
gel(x,2) = gcopy(x0);
gel(x,3) = gcopy(x1); return normalizepol_i(x,4);
}
/* same, no copy */
GEN
deg1pol_i(GEN x1, GEN x0,long v)
{
GEN x = cgetg(4,t_POL);
x[1] = evalsigne(1) | evalvarn(v);
gel(x,2) = x0;
gel(x,3) = x1; return normalizepol_i(x,4);
}
static GEN
_gtopoly(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:
if (varncmp(varn(x), v) < 0)
pari_err(talker,"variable must have higher priority in gtopoly");
y=gcopy(x); break;
case t_SER:
if (varncmp(varn(x), v) < 0)
pari_err(talker,"variable must have higher priority in gtopoly");
y = ser2rfrac(x);
if (typ(y) != t_POL)
pari_err(talker,"t_SER with negative valuation in gtopoly");
break;
case t_RFRAC:
if (varncmp(varn(gel(x,2)), v) < 0)
pari_err(talker,"variable must have higher priority in gtopoly");
y=gdeuc(gel(x,1),gel(x,2)); break;
case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
lx = lg(x); if (tx == t_QFR) lx--;
if (varncmp(gvar(x), v) <= 0)
pari_err(talker,"variable must have higher priority in gtopoly");
if (reverse)
{
while (lx-- && isexactzero(gel(x,lx)));
i = lx+2; y = cgetg(i,t_POL);
y[1] = gcmp0(x)? 0: evalsigne(1);
for (j=2; j<i; j++) gel(y,j) = gcopy(gel(x,j-1));
}
else
{
i=1; j=lx; while (lx-- && isexactzero(gel(x,i++)));
i = lx+2; y = cgetg(i,t_POL);
y[1] = gcmp0(x)? 0: evalsigne(1);
lx = j-1;
for (j=2; j<i; j++) gel(y,j) = gcopy(gel(x,lx--));
}
break;
default: pari_err(typeer,"gtopoly");
return NULL; /* not reached */
}
setvarn(y,v); return y;
}
GEN
gtopolyrev(GEN x, long v) { return _gtopoly(x,v,1); }
GEN
gtopoly(GEN x, long v) { return _gtopoly(x,v,0); }
GEN
scalarser(GEN x, long v, long prec)
{
long i, l;
GEN y;
if (isexactzero(x)) return zeroser(v,0);
l = prec + 2; y = cgetg(l, t_SER);
y[1] = evalsigne(1) | evalvalp(0) | evalvarn(v);
gel(y,2) = gcopy(x); for (i=3; i<l; i++) gel(y,i) = gen_0;
return y;
}
static GEN _gtoser(GEN x, long v, long prec);
/* assume x a t_[SER|POL], apply gtoser to all coeffs */
static GEN
coefstoser(GEN x, long v, long prec)
{
long i, tx = typ(x), lx = lg(x);
GEN y = cgetg(lx, tx);
y[1] = x[1];
for (i=2; i<lx; i++) gel(y,i) = _gtoser(gel(x,i), v, prec);
return y;
}
/* assume x a scalar or t_POL. Not stack-clean */
GEN
poltoser(GEN x, long v, long prec)
{
long tx = typ(x), vx = varn(x);
GEN y;
if (is_scalar_t(tx) || varncmp(vx, v) > 0) return scalarser(x, v, prec);
if (varncmp(vx, v) < 0) return coefstoser(x, v, prec);
if (!lgpol(x)) return zeroser(v, prec);
y = greffe(x, prec+2, 1);
setvarn(y, v); return y;
}
/* x a t_RFRAC[N]. Not stack-clean */
GEN
rfractoser(GEN x, long v, long prec)
{
return gdiv(poltoser(gel(x,1), v, prec),
poltoser(gel(x,2), v, prec));
}
GEN
toser_i(GEN x)
{
switch(typ(x))
{
case t_SER: return x;
case t_POL: return poltoser(x, varn(x), precdl);
case t_RFRAC: return rfractoser(x, gvar(x), precdl);
}
return NULL;
}
static GEN
_gtoser(GEN x, long v, long prec)
{
long tx=typ(x), lx, i, j, l;
pari_sp av;
GEN y;
if (v < 0) v = 0;
if (tx == t_SER)
{
long vx = varn(x);
if (varncmp(vx, v) < 0) y = coefstoser(x, v, prec);
else if (varncmp(vx, v) > 0) y = scalarser(x, v, prec);
else y = gcopy(x);
return y;
}
if (is_scalar_t(tx)) return scalarser(x,v,prec);
switch(tx)
{
case t_POL:
if (varncmp(varn(x), v) < 0)
pari_err(talker,"main variable must have higher priority in gtoser");
y = poltoser(x, v, prec); l = lg(y);
for (i=2; i<l; i++)
if (gel(y,i) != gen_0) gel(y,i) = gcopy(gel(y,i));
break;
case t_RFRAC:
if (varncmp(varn(gel(x,2)), v) < 0)
pari_err(talker,"main variable must have higher priority in gtoser");
av = avma;
return gerepileupto(av, rfractoser(x, v, prec));
case t_QFR: case t_QFI: case t_VEC: case t_COL:
if (varncmp(gvar(x), v) < 0)
pari_err(talker,"main variable must have higher priority in gtoser");
lx = lg(x); if (tx == t_QFR) lx--;
i = 1; while (i<lx && isexactzero(gel(x,i))) i++;
if (i == lx) return zeroser(v, lx-1);
lx -= i-2; x += i-2;
y = cgetg(lx,t_SER);
y[1] = evalsigne(1) | evalvalp(i-1) | evalvarn(v);
for (j=2; j<lx; j++) gel(y,j) = gcopy(gel(x,j));
break;
default: pari_err(typeer,"gtoser");
return NULL; /* not reached */
}
return y;
}
GEN
gtoser(GEN x, long v) { return _gtoser(x,v,precdl); }
/* assume typ(x) = t_STR */
static GEN
str_to_vecsmall(GEN x)
{
char *s = GSTR(x);
long i, l = strlen(s);
GEN y = cgetg(l+1, t_VECSMALL);
s--;
for (i=1; i<=l; i++) y[i] = (long)s[i];
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) || tx == t_RFRAC) return mkveccopy(x);
if (tx == t_STR)
{
char t[2] = {0,0};
y = str_to_vecsmall(x);
lx = lg(y);
for (i=1; i<lx; i++)
{
t[0] = y[i];
gel(y,i) = strtoGENstr(t);
}
settyp(y,t_VEC); return y;
}
if (is_graphicvec_t(tx))
{
lx=lg(x); y=cgetg(lx,t_VEC);
for (i=1; i<lx; i++) gel(y,i) = gcopy(gel(x,i));
return y;
}
if (tx==t_POL)
{
lx=lg(x); y=cgetg(lx-1,t_VEC);
for (i=1; i<=lx-2; i++) gel(y,i) = gcopy(gel(x,lx-i));
return y;
}
if (tx==t_LIST)
{
lx=lgeflist(x); y=cgetg(lx-1,t_VEC); x++;
for (i=1; i<=lx-2; i++) gel(y,i) = gcopy(gel(x,i));
return y;
}
if (tx == t_VECSMALL) return vecsmall_to_vec(x);
if (!signe(x)) return mkvec(gen_0);
lx=lg(x); y=cgetg(lx-1,t_VEC); x++;
for (i=1; i<=lx-2; i++) gel(y,i) = gcopy(gel(x,i));
return y;
}
GEN
gtocol(GEN x)
{
long lx, tx, i, j, h;
GEN y;
if (!x) return cgetg(1,t_COL);
tx = typ(x);
if (tx != t_MAT) { y = gtovec(x); settyp(y, t_COL); return y; }
lx = lg(x); if (lx == 1) return cgetg(1, t_COL);
h = lg(x[1]); y = cgetg(h, t_COL);
for (i = 1 ; i < h; i++) {
gel(y,i) = cgetg(lx, t_VEC);
for (j = 1; j < lx; j++) gmael(y,i,j) = gcopy(gcoeff(x,i,j));
}
return y;
}
GEN
gtovecsmall(GEN x)
{
GEN V;
long tx, l,i;
if (!x) return cgetg(1,t_VECSMALL);
tx = typ(x);
if (tx == t_VECSMALL) return gcopy(x);
if (tx == t_INT) return mkvecsmall(itos(x));
if (tx == t_STR) return str_to_vecsmall(x);
if (!is_vec_t(tx)) pari_err(typeer,"vectosmall");
l = lg(x);
V = cgetg(l,t_VECSMALL);
for(i=1;i<l;i++) V[i] = itos(gel(x,i));
return V;
}
GEN
compo(GEN x, long n)
{
long tx = typ(x);
ulong l, lx = (ulong)lg(x);
if (!is_recursive_t(tx))
pari_err(talker, "this object is a leaf. It has no components");
if (n < 1) pari_err(talker,"nonexistent component");
if (tx == t_POL && (ulong)n+1 >= lx) return gen_0;
if (tx == t_LIST) lx = (ulong)lgeflist(x);
l = (ulong)lontyp[tx] + (ulong)n-1; /* beware overflow */
if (l >= lx) pari_err(talker,"nonexistent component");
return gcopy(gel(x,l));
}
/* assume v > varn(x), extract coeff of pol_x[v]^n */
static GEN
multi_coeff(GEN x, long n, long v, long dx)
{
long i, lx = dx+3;
GEN z = cgetg(lx, t_POL); z[1] = x[1];
for (i = 2; i < lx; i++) gel(z,i) = polcoeff_i(gel(x,i), n, v);
return normalizepol_i(z, lx);
}
/* assume x a t_POL */
static GEN
_polcoeff(GEN x, long n, long v)
{
long w, dx;
dx = degpol(x);
if (dx < 0) return gen_0;
if (v < 0 || v == (w=varn(x)))
return (n < 0 || n > dx)? gen_0: gel(x,n+2);
if (w > v) return n? gen_0: x;
/* w < v */
return multi_coeff(x, n, v, dx);
}
/* assume x a t_SER */
static GEN
_sercoeff(GEN x, long n, long v)
{
long w, dx = degpol(x), ex = valp(x), N = n - ex;
GEN z;
if (dx < 0)
{
if (N >= 0) pari_err(talker,"non existent component in truecoeff");
return gen_0;
}
if (v < 0 || v == (w=varn(x)))
{
if (N > dx) pari_err(talker,"non existent component in truecoeff");
return (N < 0)? gen_0: gel(x,N+2);
}
if (w > v) return N? gen_0: x;
/* w < v */
z = multi_coeff(x, n, v, dx);
if (ex) z = gmul(z, monomial(gen_1,ex, w));
return z;
}
/* assume x a t_RFRAC(n) */
static GEN
_rfraccoeff(GEN x, long n, long v)
{
GEN P,Q, p = gel(x,1), q = gel(x,2);
long vp = gvar(p), vq = gvar(q);
if (v < 0) v = min(vp, vq);
P = (vp == v)? p: swap_vars(p, v);
Q = (vq == v)? q: swap_vars(q, v);
if (!ismonome(Q)) pari_err(typeer, "polcoeff");
n += degpol(Q);
return gdiv(_polcoeff(P, n, v), leading_term(Q));
}
GEN
polcoeff_i(GEN x, long n, long v)
{
switch(typ(x))
{
case t_POL: return _polcoeff(x,n,v);
case t_SER: return _sercoeff(x,n,v);
case t_RFRAC: return _rfraccoeff(x,n,v);
default: return n? gen_0: x;
}
}
/* 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);
pari_sp av;
if (is_scalar_t(tx)) return n? gen_0: gcopy(x);
av = avma;
switch(tx)
{
case t_POL: x = _polcoeff(x,n,v); break;
case t_SER: x = _sercoeff(x,n,v); break;
case t_RFRAC: x = _rfraccoeff(x,n,v); break;
case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
if (n>=1 && n<lg(x)) return gcopy(gel(x,n));
/* fall through */
default: pari_err(talker,"nonexistent component in truecoeff");
}
if (x == gen_0) return x;
if (avma == av) return gcopy(x);
return gerepilecopy(av, x);
}
GEN
truecoeff(GEN x, long n)
{
return polcoeff0(x,n,-1);
}
GEN
denom(GEN x)
{
long lx, i;
pari_sp av, tetpil;
GEN s,t;
switch(typ(x))
{
case t_INT: case t_REAL: case t_INTMOD: case t_PADIC: case t_SER:
return gen_1;
case t_FRAC:
return icopy(gel(x,2));
case t_COMPLEX:
av=avma; t=denom(gel(x,1)); s=denom(gel(x,2)); tetpil=avma;
return gerepile(av,tetpil,glcm(s,t));
case t_QUAD:
av=avma; t=denom(gel(x,2)); s=denom(gel(x,3)); tetpil=avma;
return gerepile(av,tetpil,glcm(s,t));
case t_POLMOD:
return denom(gel(x,2));
case t_RFRAC:
return gcopy(gel(x,2));
case t_POL:
return pol_1[varn(x)];
case t_VEC: case t_COL: case t_MAT:
lx=lg(x); if (lx==1) return gen_1;
av = tetpil = avma; s = denom(gel(x,1));
for (i=2; i<lx; i++)
{
t = denom(gel(x,i));
if (t != gen_1) { tetpil=avma; s=glcm(s,t); }
}
return gerepile(av,tetpil,s);
}
pari_err(typeer,"denom");
return NULL; /* not reached */
}
GEN
numer(GEN x)
{
pari_sp 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:
return (signe(x[2])>0)? icopy(gel(x,1)): negi(gel(x,1));
case t_POLMOD:
av=avma; s=numer(gel(x,2)); tetpil=avma;
return gerepile(av,tetpil,gmodulo(s,gel(x,1)));
case t_RFRAC:
return gcopy(gel(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));
}
pari_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(gel(x,2));
case t_POLMOD:
if (v < 0 || v == varn(gel(x,1))) return gcopy(gel(x,2));
y = cgetg(3,tx);
gel(y,1) = lift0(gel(x,1),v);
gel(y,2) = lift0(gel(x,2),v); return y;
case t_PADIC:
return gtrunc(x);
case t_FRAC: case t_COMPLEX: case t_RFRAC:
case t_POL: case t_SER: case t_VEC: case t_COL: case t_MAT:
y = init_gen_op(x, tx, &lx, &i);
for (; i<lx; i++) gel(y,i) = lift0(gel(x,i), v);
return y;
case t_QUAD:
y=cgetg(4,t_QUAD); gel(y,1) = gcopy(gel(x,1));
for (i=2; i<4; i++) gel(y,i) = lift0(gel(x,i), v);
return y;
}
pari_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 gel(x,2);
case t_POLMOD:
if (v < 0 || v == varn(gel(x,1))) return gel(x,2);
gel(x,1) = lift_intern0(gel(x,1),v);
gel(x,2) = lift_intern0(gel(x,2),v);
return x;
case t_SER:
case t_FRAC: case t_COMPLEX: case t_QUAD: case t_POL:
case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
lx = lg(x);
for (i = lx-1; i>=lontyp[tx]; i--)
gel(x,i) = lift_intern0(gel(x,i),v);
return x;
}
pari_err(typeer,"lift_intern");
return NULL; /* not reached */
}
/* memes conventions pour v que lift */
GEN
centerlift0(GEN x, long v)
{
long i, lx, tx = typ(x);
pari_sp av;
GEN y;
switch(tx)
{
case t_INT:
return icopy(x);
case t_INTMOD:
av = avma; i = cmpii(shifti(gel(x,2),1), gel(x,1)); avma = av;
return (i > 0)? subii(gel(x,2),gel(x,1)): icopy(gel(x,2));
case t_POLMOD:
if (v < 0 || v == varn(gel(x,1))) return gcopy(gel(x,2));
y = cgetg(3, t_POLMOD);
gel(y,1) = centerlift0(gel(x,1),v);
gel(y,2) = centerlift0(gel(x,2),v); return y;
case t_POL: case t_SER:
case t_FRAC: case t_COMPLEX: case t_RFRAC:
case t_VEC: case t_COL: case t_MAT:
y = init_gen_op(x, tx, &lx, &i);
for (; i<lx; i++) gel(y,i) = centerlift0(gel(x,i),v);
return y;
case t_QUAD:
y=cgetg(4,t_QUAD); gel(y,1) = gcopy(gel(x,1));
gel(y,2) = centerlift0(gel(x,2),v);
gel(y,3) = centerlift0(gel(x,3),v); return y;
}
pari_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, tx = typ(x);
pari_sp av;
GEN z;
switch(tx)
{
case t_POL:
lx = lg(x); z = cgetg(lx,t_POL); z[1] = x[1];
for (j=2; j<lx; j++) gel(z,j) = f(gel(x,j));
return normalizepol_i(z, lx);
case t_SER:
lx = lg(x); z = cgetg(lx,t_SER); z[1] = x[1];
for (j=2; j<lx; j++) gel(z,j) = f(gel(x,j));
return normalize(z);
case t_RFRAC: {
GEN dxb, n, d;
av = avma; dxb = gconj(gel(x,2));
n = gmul(gel(x,1), dxb);
d = gmul(gel(x,2), dxb);
return gerepileupto(av, gdiv(f(n), d));
}
case t_VEC: case t_COL: case t_MAT:
lx=lg(x); z=cgetg(lx,tx);
for (i=1; i<lx; i++) gel(z,i) = f(gel(x,i));
return z;
}
pari_err(typeer,"greal/gimag");
return NULL; /* not reached */
}
GEN
real_i(GEN x)
{
switch(typ(x))
{
case t_INT: case t_REAL: case t_FRAC:
return x;
case t_COMPLEX:
return gel(x,1);
case t_QUAD:
return gel(x,2);
}
return op_ReIm(real_i,x);
}
GEN
imag_i(GEN x)
{
switch(typ(x))
{
case t_INT: case t_REAL: case t_FRAC:
return gen_0;
case t_COMPLEX:
return gel(x,2);
case t_QUAD:
return gel(x,3);
}
return op_ReIm(imag_i,x);
}
GEN
greal(GEN x)
{
switch(typ(x))
{
case t_INT: case t_REAL: case t_FRAC:
return gcopy(x);
case t_COMPLEX:
return gcopy(gel(x,1));
case t_QUAD:
return gcopy(gel(x,2));
}
return op_ReIm(greal,x);
}
GEN
gimag(GEN x)
{
switch(typ(x))
{
case t_INT: case t_REAL: case t_FRAC:
return gen_0;
case t_COMPLEX:
return gcopy(gel(x,2));
case t_QUAD:
return gcopy(gel(x,3));
}
return op_ReIm(gimag,x);
}
/*******************************************************************/
/* */
/* LOGICAL OPERATIONS */
/* */
/*******************************************************************/
static long
_egal(GEN x, GEN y)
{
pari_sp av = avma;
long r = gequal(simplify_i(x), simplify_i(y));
avma = av; return r;
}
GEN
glt(GEN x, GEN y) { return gcmp(x,y)<0? gen_1: gen_0; }
GEN
gle(GEN x, GEN y) { return gcmp(x,y)<=0? gen_1: gen_0; }
GEN
gge(GEN x, GEN y) { return gcmp(x,y)>=0? gen_1: gen_0; }
GEN
ggt(GEN x, GEN y) { return gcmp(x,y)>0? gen_1: gen_0; }
GEN
geq(GEN x, GEN y) { return _egal(x,y)? gen_1: gen_0; }
GEN
gne(GEN x, GEN y) { return _egal(x,y)? gen_0: gen_1; }
GEN
gand(GEN x, GEN y) { return gcmp0(x)? gen_0: (gcmp0(y)? gen_0: gen_1); }
GEN
gor(GEN x, GEN y) { return gcmp0(x)? (gcmp0(y)? gen_0: gen_1): gen_1; }
GEN
gnot(GEN x) { return gcmp0(x)? gen_1: gen_0; }
/*******************************************************************/
/* */
/* FORMAL SIMPLIFICATIONS */
/* */
/*******************************************************************/
GEN
geval(GEN x)
{
long lx, i, tx = typ(x);
pari_sp av, tetpil;
GEN y,z;
if (is_const_t(tx)) return gcopy(x);
if (is_graphicvec_t(tx))
{
lx=lg(x); y=cgetg(lx, tx);
for (i=1; i<lx; i++) gel(y,i) = geval(gel(x,i));
return y;
}
switch(tx)
{
case t_STR:
return gp_read_str(GSTR(x));
case t_POLMOD: y=cgetg(3,tx);
gel(y,1) = geval(gel(x,1));
av=avma; z=geval(gel(x,2)); tetpil=avma;
gel(y,2) = gerepile(av,tetpil,gmod(z,gel(y,1)));
return y;
case t_POL:
lx=lg(x); if (lx==2) return gen_0;
{
long vx = varn(x);
entree *ep = varentries[vx];
if (!ep) return gcopy(x);
z = (GEN)ep->value;
if (gequal(x, pol_x[vx])) return gcopy(z);
}
y=gen_0; av=avma;
for (i=lx-1; i>1; i--)
y = gadd(geval(gel(x,i)), gmul(z,y));
return gerepileupto(av, y);
case t_SER:
pari_err(impl, "evaluation of a power series");
case t_RFRAC:
return gdiv(geval(gel(x,1)),geval(gel(x,2)));
}
pari_err(typeer,"geval");
return NULL; /* not reached */
}
GEN
simplify_i(GEN x)
{
long tx=typ(x),i,lx;
GEN y;
switch(tx)
{
case t_INT: case t_REAL: case t_FRAC:
case t_INTMOD: case t_PADIC: case t_QFR: case t_QFI:
case t_LIST: case t_STR: case t_VECSMALL:
return x;
case t_COMPLEX:
if (isexactzero(gel(x,2))) return simplify_i(gel(x,1));
y=cgetg(3,t_COMPLEX);
gel(y,1) = simplify_i(gel(x,1));
gel(y,2) = simplify_i(gel(x,2)); return y;
case t_QUAD:
if (isexactzero(gel(x,3))) return simplify_i(gel(x,2));
y=cgetg(4,t_QUAD);
y[1]=x[1];
gel(y,2) = simplify_i(gel(x,2));
gel(y,3) = simplify_i(gel(x,3)); return y;
case t_POLMOD: y=cgetg(3,t_POLMOD);
gel(y,1) = simplify_i(gel(x,1));
if (typ(y[1]) != t_POL) y[1] = x[1]; /* invalid object otherwise */
gel(y,2) = simplify_i(gel(x,2)); return y;
case t_POL:
lx = lg(x); if (lx==2) return gen_0;
if (lx==3) return simplify_i(gel(x,2));
y = cgetg(lx,t_POL); y[1] = x[1];
for (i=2; i<lx; i++) gel(y,i) = simplify_i(gel(x,i));
return y;
case t_SER:
lx = lg(x); y = cgetg(lx,t_SER); y[1] = x[1];
for (i=2; i<lx; i++) gel(y,i) = simplify_i(gel(x,i));
return y;
case t_RFRAC: y=cgetg(3,t_RFRAC);
gel(y,1) = simplify_i(gel(x,1));
gel(y,2) = simplify_i(gel(x,2));
if (typ(gel(y,2)) != t_POL) return gdiv(gel(y,1), gel(y,2));
return y;
case t_VEC: case t_COL: case t_MAT:
lx=lg(x); y=cgetg(lx,tx);
for (i=1; i<lx; i++) gel(y,i) = simplify_i(gel(x,i));
return y;
}
pari_err(typeer,"simplify_i");
return NULL; /* not reached */
}
GEN
simplify(GEN x)
{
pari_sp av = avma;
return gerepilecopy(av, simplify_i(x));
}
/*******************************************************************/
/* */
/* EVALUATION OF SOME SIMPLE OBJECTS */
/* */
/*******************************************************************/
static GEN
qfeval0_i(GEN q, GEN x, long n)
{
long i, j;
pari_sp av=avma;
GEN res=gen_0;
for (i=2;i<n;i++)
for (j=1;j<i;j++)
res = gadd(res, gmul(gcoeff(q,i,j), mulii(gel(x,i),gel(x,j))) );
res=gshift(res,1);
for (i=1;i<n;i++)
res = gadd(res, gmul(gcoeff(q,i,i), sqri(gel(x,i))) );
return gerepileupto(av,res);
}
#if 0
static GEN
qfeval0(GEN q, GEN x, long n)
{
long i, j;
pari_sp av=avma;
GEN res=gen_0;
for (i=2;i<n;i++)
for (j=1;j<i;j++)
res = gadd(res, gmul(gcoeff(q,i,j), gmul(gel(x,i),gel(x,j))) );
res=gshift(res,1);
for (i=1;i<n;i++)
res = gadd(res, gmul(gcoeff(q,i,i), gsqr(gel(x,i))) );
return gerepileupto(av,res);
}
#else
static GEN
qfeval0(GEN q, GEN x, long n)
{
long i, j;
pari_sp av=avma;
GEN res = gmul(gcoeff(q,1,1), gsqr(gel(x,1)));
for (i=2; i<n; i++)
{
GEN l = gel(q,i);
GEN sx = gmul(gel(l,1), gel(x,1));
for (j=2; j<i; j++)
sx = gadd(sx, gmul(gel(l,j),gel(x,j)));
sx = gadd(gshift(sx,1), gmul(gel(l,i),gel(x,i)));
res = gadd(res, gmul(gel(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)
pari_err(talker,"invalid data in qfeval");
return gen_0;
}
if (typ(q) != t_MAT || lg(q[1]) != n)
pari_err(talker,"invalid quadratic form in qfeval");
if (typ(x) != t_COL || lg(x) != n)
pari_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;
pari_sp av=avma;
GEN res = gmul(gcoeff(q,1,1), mulii(gel(x,1),gel(y,1)));
for (i=2;i<n;i++)
{
if (!signe(x[i]))
{
if (!signe(y[i])) continue;
for (j=1;j<i;j++)
res = gadd(res, gmul(gcoeff(q,i,j), mulii(gel(x,j),gel(y,i))));
}
else if (!signe(y[i]))
{
for (j=1;j<i;j++)
res = gadd(res, gmul(gcoeff(q,i,j), mulii(gel(x,i),gel(y,j))));
}
else
{
for (j=1;j<i;j++)
{
GEN p1 = addii(mulii(gel(x,i),gel(y,j)), mulii(gel(x,j),gel(y,i)));
res = gadd(res, gmul(gcoeff(q,i,j),p1));
}
res = gadd(res, gmul(gcoeff(q,i,i), mulii(gel(x,i),gel(y,i))));
}
}
return gerepileupto(av,res);
}
#if 0
static GEN
qfbeval0(GEN q, GEN x, GEN y, long n)
{
long i, j;
pari_sp av=avma;
GEN res = gmul(gcoeff(q,1,1), gmul(gel(x,1),gel(y,1)));
for (i=2;i<n;i++)
{
for (j=1;j<i;j++)
{
GEN p1 = gadd(gmul(gel(x,i),gel(y,j)), gmul(gel(x,j),gel(y,i)));
res = gadd(res, gmul(gcoeff(q,i,j),p1));
}
res = gadd(res, gmul(gcoeff(q,i,i), gmul(gel(x,i),gel(y,i))));
}
return gerepileupto(av,res);
}
#else
static GEN
qfbeval0(GEN q, GEN x, GEN y, long n)
{
long i, j;
pari_sp av=avma;
GEN res = gmul(gcoeff(q,1,1), gmul(gel(x,1), gel(y,1)));
for (i=2; i<n; i++)
{
GEN l = gel(q,i);
GEN sx = gmul(gel(l,1), gel(y,1));
GEN sy = gmul(gel(l,1), gel(x,1));
for (j=2; j<i; j++)
{
sx = gadd(sx, gmul(gel(l,j),gel(y,j)));
sy = gadd(sy, gmul(gel(l,j),gel(x,j)));
}
sx = gadd(sx, gmul(gel(l,i),gel(y,i)));
sx = gmul(gel(x,i), sx);
sy = gmul(gel(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)
pari_err(talker,"invalid data in qfbeval");
return gen_0;
}
if (typ(q) != t_MAT || lg(q[1]) != n)
pari_err(talker,"invalid quadratic form in qfbeval");
if (typ(x) != t_COL || lg(x) != n || typ(y) != t_COL || lg(y) != n)
pari_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)
pari_err(talker,"invalid data in qf_base_change");
return res;
}
if (typ(M) != t_MAT || k == 1 || lg(M[1]) != n)
pari_err(talker,"invalid base change matrix in qf_base_change");
for (i=1;i<k;i++)
{
gel(res,i) = cgetg(k,t_COL);
gcoeff(res,i,i) = qf(q,gel(M,i),n);
}
for (i=2;i<k;i++)
for (j=1;j<i;j++)
gcoeff(res,i,j)=gcoeff(res,j,i) = qfb(q,gel(M,i),gel(M,j),n);
return res;
}
/* return Re(x * y), x and y scalars */
GEN
mul_real(GEN x, GEN y)
{
if (typ(x) == t_COMPLEX)
{
if (typ(y) == t_COMPLEX)
{
pari_sp av=avma, tetpil;
GEN p1 = gmul(gel(x,1), gel(y,1));
GEN p2 = gneg(gmul(gel(x,2), gel(y,2)));
tetpil=avma; return gerepile(av,tetpil,gadd(p1,p2));
}
x = gel(x,1);
}
else if (typ(y) == t_COMPLEX) y = gel(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 i, j, k, lx = lg(x), ly = lg(y), l = lg(x[1]);
pari_sp av;
GEN p1, z = cgetg(ly,t_MAT);
for (j=1; j<ly; j++)
{
gel(z,j) = cgetg(l,t_COL);
for (i=1; i<l; i++)
{
p1 = gen_0; av=avma;
for (k=1; k<lx; k++)
p1 = gadd(p1, mul_real(gcoeff(x,i,k),gcoeff(y,k,j)));
gcoeff(z,i,j) = gerepileupto(av, p1);
}
}
return z;
}
static GEN
hqfeval0(GEN q, GEN x, long n)
{
long i, j;
pari_sp av=avma;
GEN res=gen_0;
for (i=2;i<n;i++)
for (j=1;j<i;j++)
{
GEN p1 = gmul(gel(x,i),gconj(gel(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(gel(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)
pari_err(talker,"invalid data in hqfeval");
return gen_0;
}
if (typ(q) != t_MAT || lg(q[1]) != n)
pari_err(talker,"invalid quadratic form in hqfeval");
if (typ(x) != t_COL || lg(x) != n)
pari_err(talker,"invalid vector in hqfeval");
return hqfeval0(q,x,n);
}
GEN
poleval(GEN x, GEN y)
{
long i, j, imin, tx = typ(x);
pari_sp av0 = avma, av, lim;
GEN p1, p2, r, s;
if (is_scalar_t(tx)) return gcopy(x);
switch(tx)
{
case t_POL:
i = lg(x)-1; imin = 2; break;
case t_RFRAC:
p1 = poleval(gel(x,1),y);
p2 = poleval(gel(x,2),y);
return gerepileupto(av0, gdiv(p1,p2));
case t_VEC: case t_COL:
i = lg(x)-1; imin = 1; break;
default: pari_err(typeer,"poleval");
return NULL; /* not reached */
}
if (i<=imin)
return (i==imin)? gcopy(gel(x,imin)): gen_0;
lim = stack_lim(av0,2);
p1 = gel(x,i); i--;
if (typ(y)!=t_COMPLEX)
{
#if 0 /* standard Horner's rule */
for ( ; i>=imin; i--)
p1 = gadd(gmul(p1,y),gel(x,i));
#endif
/* specific attention to sparse polynomials */
for ( ; i>=imin; i=j-1)
{
for (j=i; isexactzero(gel(x,j)); j--)
if (j==imin)
{
if (i!=j) y = gpowgs(y, i-j+1);
return gerepileupto(av0, gmul(p1,y));
}
r = (i==j)? y: gpowgs(y, i-j+1);
p1 = gadd(gmul(p1,r), gel(x,j));
if (low_stack(lim, stack_lim(av0,2)))
{
if (DEBUGMEM>1) pari_warn(warnmem,"poleval: i = %ld",i);
p1 = gerepileupto(av0, p1);
}
}
return gerepileupto(av0,p1);
}
p2 = gel(x,i); i--; r = gtrace(y); s = gneg_i(gnorm(y));
av = avma;
for ( ; i>=imin; i--)
{
GEN p3 = gadd(p2, gmul(r, p1));
p2 = gadd(gel(x,i), gmul(s, p1)); p1 = p3;
if (low_stack(lim, stack_lim(av0,2)))
{
if (DEBUGMEM>1) pari_warn(warnmem,"poleval: i = %ld",i);
gerepileall(av, 2, &p1, &p2);
}
}
return gerepileupto(av0, gadd(p2, gmul(y,p1)));
}
syntax highlighted by Code2HTML, v. 0.9.1