/********************************************************************/
/** **/
/** GENERIC OPERATIONS **/
/** (first part) **/
/** **/
/********************************************************************/
/* $Id: gen1.c,v 1.2 1999/09/21 10:57:55 karim Exp $ */
#include "pari.h"
#define swapspec(x,y, nx,ny) {long _a=nx;GEN _z=x; nx=ny; ny=_a; x=y; y=_z;}
GEN quickmul(GEN a, GEN b, long na, long nb);
#define cpifstack(x) isonstack(x)?gcopy(x):x
/* y is a polmod, f is gadd or gmul */
static GEN
op_polmod(GEN f(GEN,GEN), GEN x, GEN y, long tx)
{
GEN mod,k,l, z=cgetg(3,t_POLMOD);
long av,tetpil;
l=(GEN)y[1];
if (tx==t_POLMOD)
{
k=(GEN)x[1];
if (gegal(k,l))
{ mod=cpifstack(k); x=(GEN)x[2]; y=(GEN)y[2]; }
else
{
long vx=varn(k), vy=varn(l);
if (vx==vy) { mod=srgcd(k,l); x=(GEN)x[2]; y=(GEN)y[2]; }
else
if (vx<vy) { mod=cpifstack(k); x=(GEN)x[2]; }
else { mod=cpifstack(l); y=(GEN)y[2]; }
}
}
else
{
mod=cpifstack(l); y=(GEN)y[2];
if (is_scalar_t(tx))
{
z[2] = (long)f(x,y);
z[1] = (long)mod; return z;
}
}
av=avma; x = f(x,y); tetpil=avma;
z[2] = lpile(av,tetpil,gmod(x,mod));
z[1] = (long)mod; return z;
}
#undef cpifstack
/********************************************************************/
/** **/
/** SUBTRACTION **/
/** **/
/********************************************************************/
GEN
gsub(GEN x, GEN y)
{
long tetpil, av = avma;
y=gneg_i(y); tetpil=avma;
return gerepile(av,tetpil,gadd(x,y));
}
/********************************************************************/
/** **/
/** ADDITION **/
/** **/
/********************************************************************/
static GEN
addpadic(GEN x, GEN y)
{
long c,e,r,d,r1,r2,av,tetpil;
GEN z,p1,p2, p = (GEN)x[2];
z=cgetg(5,t_PADIC); icopyifstack(p, z[2]); av=avma;
e=valp(x); r=valp(y); d = r-e;
if (d<0) { p1=x; x=y; y=p1; e=r; d = -d; }
r1=precp(x); r2=precp(y);
if (d)
{
r = d+r2;
p1 = (d==1)? p: gclone(gpuigs(p,d));
avma=av;
if (r<r1) z[3]=lmulii(p1,(GEN)y[3]);
else
{
r=r1; z[3]=licopy((GEN)x[3]);
}
av=avma; p2=mulii(p1,(GEN)y[4]);
if (d!=1) gunclone(p1);
p1=addii(p2,(GEN)x[4]); tetpil=avma;
z[4]=lpile(av,tetpil, modii(p1,(GEN)z[3]));
z[1]=evalprecp(r) | evalvalp(e); return z;
}
if (r2<r1) { r=r2; p1=x; x=y; y=p1; } else r=r1;
p1 = addii((GEN)x[4],(GEN)y[4]);
if (!signe(p1) || (c = pvaluation(p1,p,&p2)) >=r)
{
avma=av; z[4]=zero; z[3]=un;
z[1]=evalvalp(e+r); return z;
}
if (c)
{
p2=gclone(p2); avma=av;
if (c==1)
z[3] = ldivii((GEN)x[3], p);
else
{
p1 = gpuigs(p,c); tetpil=avma;
z[3] = lpile(av,tetpil, divii((GEN)x[3], p1));
}
z[4]=lmodii(p2,(GEN)z[3]); gunclone(p2);
z[1]=evalprecp(r-c) | evalvalp(e+c); return z;
}
tetpil=avma;
z[4]=lpile(av,tetpil,modii(p1,(GEN)x[3]));
z[3]=licopy((GEN)x[3]);
z[1]=evalprecp(r) | evalvalp(e); return z;
}
/* return x + y, where x is t_INT or t_FRAC(N), y t_PADIC */
static GEN
gaddpex(GEN x, GEN y)
{
long tx,e1,e2,e3,av,tetpil;
GEN z,p,p1,p2;
if (gcmp0(x)) return gcopy(y);
av=avma; p=(GEN)y[2]; tx=typ(x);
z=cgetg(5,t_PADIC); z[2]=(long)p;
e3 = (tx == t_INT)? pvaluation(x,p,&p1)
: pvaluation((GEN)x[1],p,&p1) -
pvaluation((GEN)x[2],p,&p2);
e1 = valp(y)-e3; e2 = signe(y[4])? e1+precp(y): e1;
if (e2<=0)
{
z[1] = evalprecp(0) | evalvalp(e3);
z[3] = un;
z[4] = zero;
}
else
{
if (tx != t_INT && !is_pm1(p2)) p1 = gdiv(p1,p2);
z[1] = evalprecp(e2) | evalvalp(e3);
z[3] = e1? lmul((GEN)y[3], gpuigs(p,e1)): y[3];
z[4] = lmod(p1,(GEN)z[3]);
}
tetpil=avma; return gerepile(av,tetpil,addpadic(z,y));
}
static long
kro_quad(GEN x, GEN y)
{
long k, av=avma;
x = subii(sqri((GEN)x[3]), shifti((GEN)x[2],2));
k = kronecker(x,y); avma=av; return k;
}
static GEN
addfrac(GEN x, GEN y)
{
GEN x1 = (GEN)x[1], x2 = (GEN)x[2];
GEN y1 = (GEN)y[1], y2 = (GEN)y[2], p1,p2,n,d,delta,z;
z = cgetg(3,t_FRAC);
(void)new_chunk((lgefint(x1) + lgefint(x2) +
lgefint(y1) + lgefint(y2)) << 1);
delta = mppgcd(x2,y2);
if (is_pm1(delta))
{
p1 = mulii(x1,y2);
p2 = mulii(y1,x2); avma = (long)z;
z[1] = laddii(p1,p2);
z[2] = lmulii(x2,y2); return z;
}
x2 = divii(x2,delta);
y2 = divii(y2,delta);
n = addii(mulii(x1,y2), mulii(y1,x2));
if (!signe(n)) { avma = (long)(z+3); return gzero; }
d = mulii(x2, y2);
p1 = dvmdii(n, delta, &p2);
if (p2 == gzero)
{
if (is_pm1(d)) { avma = (long)(z+3); return icopy(p1); }
avma = (long)z;
z[1] = licopy(p1);
z[2] = licopy(d); return z;
}
p1 = mppgcd(delta, p2);
if (is_pm1(p1))
{
avma = (long)z;
z[1] = licopy(n);
}
else
{
delta = divii(delta, p1);
avma = (long)z;
z[1] = ldivii(n,p1);
}
z[2] = lmulii(d,delta); return z;
}
static GEN
addrfrac(GEN x, GEN y)
{
GEN z = cgetg(3,t_RFRAC);
GEN x1 = (GEN)x[1], x2 = (GEN)x[2];
GEN y1 = (GEN)y[1], y2 = (GEN)y[2], p1,p2,n,d,delta;
long tetpil;
delta = ggcd(x2,y2);
if (gcmp1(delta))
{
p1 = gmul(x1,y2);
p2 = gmul(y1,x2);
tetpil = avma; /* numerator is non-zero */
z[1] = lpile((long)z,tetpil, gadd(p1,p2));
z[2] = lmul(x2, y2); return z;
}
x2 = gdeuc(x2,delta);
y2 = gdeuc(y2,delta);
n = gadd(gmul(x1,y2), gmul(y1,x2));
if (!signe(n)) return gerepileupto((long)(z+3), n);
tetpil = avma; d = gmul(x2, y2);
p1 = poldivres(n, delta, &p2);
if (!signe(p2))
{
if (gcmp1(d)) return gerepileupto((long)(z+3), p1);
z[1]=(long)p1; z[2]=(long)d;
gerepilemanyvec((long)z,tetpil,z+1,2); return z;
}
p1 = ggcd(delta, p2);
if (gcmp1(p1))
{
tetpil = avma;
z[1] = lcopy(n);
}
else
{
delta = gdeuc(delta, p1);
tetpil = avma;
z[1] = ldeuc(n,p1);
}
z[2] = lmul(d,delta);
gerepilemanyvec((long)z,tetpil,z+1,2); return z;
}
static GEN
addscalrfrac(GEN x, GEN y)
{
GEN p1,num, z = cgetg(3,t_RFRAC);
long tetpil, av;
p1 = gmul(x,(GEN)y[2]); tetpil = avma;
num = gadd(p1,(GEN)y[1]);
av = avma;
p1 = content((GEN)y[2]);
if (!gcmp1(p1))
{
p1 = ggcd(p1, content(num));
if (!gcmp1(p1))
{
tetpil = avma;
z[1] = ldiv(num, p1);
z[2] = ldiv((GEN)y[2], p1);
gerepilemanyvec((long)z,tetpil,z+1,2); return z;
}
}
avma = av;
z[1]=lpile((long)z,tetpil, num);
z[2]=lcopy((GEN)y[2]); return z;
}
/* assume gvar(x) = varn(mod) */
static GEN
to_polmod(GEN x, GEN mod)
{
long tx = typ(x);
GEN z = cgetg(3, t_POLMOD);
if (tx == t_RFRACN) { x = gred_rfrac(x); tx = t_RFRAC; }
if (tx == t_RFRAC) x = gmul((GEN)x[1], ginvmod((GEN)x[2],mod));
z[1] = (long)mod;
z[2] = (long)x;
return z;
}
GEN
gadd(GEN x, GEN y)
{
long vx,vy,lx,ly,tx,ty,i,j,k,l,av,tetpil;
GEN z,p1,p2;
tx=typ(x); ty=typ(y);
if (is_const_t(tx) && is_const_t(ty))
{
if (tx>ty) { p1=x; x=y; y=p1; i=tx; tx=ty; ty=i; }
}
else
{
vx=gvar(x); vy=gvar(y);
if (vx<vy || (vx==vy && tx>ty))
{
p1=x; x=y; y=p1;
i=tx; tx=ty; ty=i;
i=vx; vx=vy; vy=i;
}
if (ty==t_POLMOD) return op_polmod(gadd,x,y,tx);
}
/* here vx >= vy */
if (is_scalar_t(ty))
{
switch(tx)
{
case t_INT:
switch(ty)
{
case t_INT: return addii(x,y);
case t_REAL: return addir(x,y);
case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)y[1];
(void)new_chunk(lgefint(p2)+1);
p1 = addii(modii(x,p2),(GEN)y[2]); avma = (long)z;
z[2] = (cmpii(p1,p2) >=0)? lsubii(p1,p2): licopy(p1);
icopyifstack(p2,z[1]); return z;
case t_FRAC: case t_FRACN: z=cgetg(3,ty);
(void)new_chunk(lgefint(x) + lgefint(y[1]) + lgefint(y[2]) + 1);
p1 = mulii((GEN)y[2],x); avma = (long)z;
z[1] = laddii((GEN)y[1], p1);
z[2] = licopy((GEN)y[2]); return z;
case t_COMPLEX: z=cgetg(3,t_COMPLEX);
z[1]=ladd(x,(GEN)y[1]);
z[2]=lcopy((GEN)y[2]); return z;
case t_PADIC:
return gaddpex(x,y);
case t_QUAD: z=cgetg(4,t_QUAD);
copyifstack(y[1], z[1]);
z[2]=ladd(x,(GEN)y[2]);
z[3]=lcopy((GEN)y[3]); return z;
}
case t_REAL:
switch(ty)
{
case t_REAL: return addrr(x,y);
case t_FRAC: case t_FRACN:
if (!signe(y[1])) return rcopy(x);
if (!signe(x))
{
lx = expi((GEN)y[1]) - expi((GEN)y[2]) - expo(x);
if (lx < 0) return rcopy(x);
lx >>= TWOPOTBITS_IN_LONG;
z=cgetr(lx+3); diviiz((GEN)y[1],(GEN)y[2],z);
return z;
}
av=avma; z=addir((GEN)y[1],mulir((GEN)y[2],x)); tetpil=avma;
return gerepile(av,tetpil,divri(z,(GEN)y[2]));
case t_COMPLEX: z=cgetg(3,t_COMPLEX);
z[1]=ladd(x,(GEN)y[1]);
z[2]=lcopy((GEN)y[2]); return z;
case t_QUAD:
if (gcmp0(y)) return rcopy(x);
av=avma; i=gexpo(y)-expo(x);
if (i<=0) i=0; else i >>= TWOPOTBITS_IN_LONG;
p1=co8(y,lg(x)+i); tetpil=avma;
return gerepile(av,tetpil,gadd(p1,x));
case t_INTMOD: case t_PADIC: err(gadderf,tx,ty);
}
case t_INTMOD:
switch(ty)
{
case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)x[1]; p1=(GEN)y[1];
if (p1==p2 || egalii(p1,p2))
{
icopyifstack(p2,z[1]);
if (!is_bigint(p2))
{
z[2] = lstoi(addssmod(itos((GEN)x[2]),itos((GEN)y[2]), p2[2]));
return z;
}
}
else
{ p2 = mppgcd(p1,p2); z[1] = (long)p2; }
av=avma; (void)new_chunk((lgefint(p1)<<1) + lgefint(x[1]));
p1=addii((GEN)x[2],(GEN)y[2]); avma=av;
z[2]=lmodii(p1,p2); return z;
case t_FRAC: case t_FRACN: z=cgetg(3,t_INTMOD); p2=(GEN)x[1];
(void)new_chunk(lgefint(p2)<<2);
p1 = mulii((GEN)y[1], mpinvmod((GEN)y[2],p2));
p1 = addii(modii(p1,p2), (GEN)x[2]); avma=(long)z;
z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
case t_COMPLEX: z=cgetg(3,t_COMPLEX);
z[1]=ladd(x,(GEN)y[1]);
z[2]=lcopy((GEN)y[2]); return z;
case t_PADIC:
l=avma; p1=cgetg(3,t_INTMOD);
p1[1]=x[1]; p1[2]=lgeti(lgefint(x[1]));
gaffect(y,p1); tetpil=avma;
return gerepile(l,tetpil,gadd(p1,x));
case t_QUAD: z=cgetg(4,t_QUAD);
copyifstack(y[1], z[1]);
z[2]=ladd(x,(GEN)y[2]);
z[3]=lcopy((GEN)y[3]); return z;
}
case t_FRAC: case t_FRACN:
switch (ty)
{
case t_FRAC: return addfrac(x,y);
case t_FRACN: z=cgetg(3,t_FRACN); l=avma;
p1=mulii((GEN)x[1],(GEN)y[2]);
p2=mulii((GEN)x[2],(GEN)y[1]);
tetpil=avma; z[1]=lpile(l,tetpil,addii(p1,p2));
z[2]=lmulii((GEN)x[2],(GEN)y[2]);
return z;
case t_COMPLEX: z=cgetg(3,t_COMPLEX);
z[1]=ladd((GEN)y[1],x);
z[2]=lcopy((GEN)y[2]); return z;
case t_PADIC:
return gaddpex(x,y);
case t_QUAD: z=cgetg(4,t_QUAD);
z[1]=lcopy((GEN)y[1]);
z[2]=ladd((GEN)y[2],x);
z[3]=lcopy((GEN)y[3]); return z;
}
case t_COMPLEX:
switch(ty)
{
case t_COMPLEX: z=cgetg(3,t_COMPLEX);
z[1]=ladd((GEN)x[1],(GEN)y[1]);
z[2]=ladd((GEN)x[2],(GEN)y[2]); return z;
case t_PADIC:
if (krosg(-1,(GEN)y[2])== -1)
{
z=cgetg(3,t_COMPLEX);
z[1]=ladd((GEN)x[1],y);
z[2]=lcopy((GEN)x[2]); return z;
}
av=avma; l = signe(y[4])? precp(y): 1;
p1=cvtop(x,(GEN)y[2], l + valp(y)); tetpil=avma;
return gerepile(av,tetpil,gadd(p1,y));
case t_QUAD:
lx=precision(x); if (!lx) err(gadderi,tx,ty);
if (gcmp0(y)) return gcopy(x);
av=avma; i=gexpo(y)-gexpo(x);
if (i<=0) i=0; else i >>= TWOPOTBITS_IN_LONG;
p1=co8(y,lx+i); tetpil=avma;
return gerepile(av,tetpil,gadd(p1,x));
}
case t_PADIC:
switch(ty)
{
case t_PADIC:
if (!egalii((GEN)x[2],(GEN)y[2])) err(gadderi,tx,ty);
return addpadic(x,y);
case t_QUAD:
if (kro_quad((GEN)y[1],(GEN)x[2]) == -1)
{
z=cgetg(4,t_QUAD);
copyifstack(y[1], z[1]);
z[2]=ladd((GEN)y[2],x);
z[3]=lcopy((GEN)y[3]); return z;
}
av=avma; l = signe(x[4])? precp(x): 1;
p1=cvtop(y,(GEN)x[2],valp(x)+l); tetpil=avma;
return gerepile(av,tetpil,gadd(p1,x));
}
case t_QUAD: z=cgetg(4,t_QUAD); k=x[1]; l=y[1];
if (!gegal((GEN)k,(GEN)l)) err(gadderi,tx,ty);
copyifstack(l, z[1]);
z[2]=ladd((GEN)x[2],(GEN)y[2]);
z[3]=ladd((GEN)x[3],(GEN)y[3]); return z;
}
err(bugparier,"gadd");
}
/* here !isscalar(y) */
if ( (vx>vy && (!is_matvec_t(tx) || !is_matvec_t(ty)))
|| (vx==vy && is_scalar_t(tx)) )
{
if (tx == t_POLMOD && vx == vy && ty != t_SER)
{
av = avma;
return gerepileupto(av, op_polmod(gadd, x, to_polmod(y,(GEN)x[1]), tx));
}
switch(ty)
{
case t_POL: ly=lgef(y);
if (ly==2) return isexactzero(x)? zeropol(vy): scalarpol(x,vy);
z = cgetg(ly,t_POL); z[1] = y[1];
z[2] = ladd(x,(GEN)y[2]);
for (i=3; i<ly; i++) z[i]=lcopy((GEN)y[i]);
return normalizepol_i(z, ly);
case t_SER: l=valp(y); ly=lg(y);
if (l<3-ly) return gcopy(y);
if (l<0)
{
z=cgetg(ly,t_SER); z[1]=y[1];
for (i=2; i<=1-l; i++) z[i]=lcopy((GEN)y[i]);
for (i=3-l; i<ly; i++) z[i]=lcopy((GEN)y[i]);
z[2-l]=ladd(x,(GEN)y[2-l]); return z;
}
if (l>0)
{
if (gcmp0(x)) return gcopy(y);
if (gcmp0(y)) ly=2;
ly += l; z=cgetg(ly,t_SER);
z[1]=evalsigne(1) | evalvalp(0) | evalvarn(vy);
for (i=3; i<=l+1; i++) z[i]=zero;
for ( ; i<ly; i++) z[i]=lcopy((GEN)y[i-l]);
z[2]=lcopy(x); return z;
}
av=avma; z=cgetg(ly,t_SER);
p1=signe(y)? gadd(x,(GEN)y[2]): x;
if (!isexactzero(p1))
{
z[1] = evalvalp(0) | evalvarn(vy);
if (signe(y))
{
z[1] |= evalsigne(1); z[2]=(long)p1;
for (i=3; i<ly; i++) z[i]=lcopy((GEN)y[i]);
}
return z;
}
avma=av; /* first coeff is 0 */
i=3; while (i<ly && gcmp0((GEN)y[i])) i++;
if (i==ly) return zeroser(vy,i-2);
z=cgetg(ly-i+2,t_SER); z[1]=evalvalp(i-2)|evalvarn(vy)|evalsigne(1);
for (j=2; j<=ly-i+1; j++) z[j]=lcopy((GEN)y[j+i-2]);
return z;
case t_RFRAC: return addscalrfrac(x,y);
case t_RFRACN: z=cgetg(3,t_RFRACN);
av=avma; p1=gmul(x,(GEN)y[2]); tetpil=avma;
z[1]=lpile(av,tetpil, gadd(p1,(GEN)y[1]));
z[2]=lcopy((GEN)y[2]); return z;
case t_VEC: case t_COL: case t_MAT:
if (isexactzero(x)) return gcopy(y);
if (ty == t_MAT) return gaddmat(x,y);
/* fall through */
case t_QFR: case t_QFI: err(gadderf,tx,ty);
}
err(typeer,"addition");
}
/* here !isscalar(x) && isscalar(y) && (vx=vy || ismatvec(x and y)) */
if (tx>ty) { p1=x; x=y; y=p1; i=tx; tx=ty; ty=i; }
switch(tx)
{
case t_POL:
switch (ty)
{
case t_POL:
lx = lgef(x); ly = lgef(y); if (lx < ly) swapspec(x,y, lx,ly);
z = cgetg(lx,t_POL); z[1] = x[1];
for (i=2; i<ly; i++) z[i]=ladd((GEN)x[i],(GEN)y[i]);
for ( ; i<lx; i++) z[i]=lcopy((GEN)x[i]);
(void)normalizepol_i(z, lx);
if (lgef(z) == 2) { avma = (long)(z + lx); z = zeropol(vx); }
return z;
case t_SER:
if (gcmp0(x)) return gcopy(y);
ly = signe(y)? lg(y): 3;
i = ly+valp(y)-gval(x,vx);
if (i<3) return gcopy(y);
p1=greffe(x,i,0); y=gadd(p1,y);
free(p1); return y;
case t_RFRAC: return addscalrfrac(x,y);
case t_RFRACN: z=cgetg(3,t_RFRACN);
av=avma; p1=gmul(x,(GEN)y[2]); tetpil=avma;
z[1]=lpile(av,tetpil, gadd(p1,(GEN)y[1]));
z[2]=lcopy((GEN)y[2]); return z;
case t_QFR: case t_QFI:
case t_VEC: case t_COL: case t_MAT: err(gadderf,tx,ty);
default: err(typeer,"addition");
}
case t_SER:
switch(ty)
{
case t_SER:
l=valp(y)-valp(x);
if (l<0) { l= -l; p1=x; x=y; y=p1; }
if (gcmp0(x)) return gcopy(x);
lx = lg(x);
ly = signe(y)? lg(y): 2;
ly += l; if (lx<ly) ly=lx;
z=cgetg(ly,t_SER);
if (l)
{
if (l>=ly-2)
for (i=2; i<ly; i++) z[i]=lcopy((GEN)x[i]);
else
{
for (i=2; i<=l+1; i++) z[i]=lcopy((GEN)x[i]);
for ( ; i<ly; i++) z[i]=ladd((GEN)x[i],(GEN)y[i-l]);
}
z[1]=x[1]; return z;
}
if (ly>2)
{
l=avma;
for (i=2; i<ly; i++)
{
p1 = gadd((GEN)x[i],(GEN)y[i]);
if (!isexactzero(p1))
{
l = i-2; stackdummy(z,l); z += l;
z[0]=evaltyp(t_SER) | evallg(ly-l);
z[1]=evalvalp(valp(x)+i-2) | evalsigne(1) | evalvarn(vx);
for (j=i+1; j<ly; j++)
z[j-l]=ladd((GEN)x[j],(GEN)y[j]);
z[2]=(long)p1; return z;
}
avma=l;
}
}
return zeroser(vx,ly-2+valp(y));
case t_RFRAC: case t_RFRACN:
if (gcmp0(y)) return gcopy(x);
l = valp(x)-gval(y,vy); l += gcmp0(x)? 3: lg(x);
if (l<3) return gcopy(x);
av=avma; ty=typ(y[2]);
p1 = is_scalar_t(ty)? (GEN)y[2]: greffe((GEN)y[2],l,1);
p1 = gdiv((GEN)y[1], p1); tetpil=avma;
return gerepile(av,tetpil,gadd(p1,x));
case t_QFR: case t_QFI: case t_VEC: case t_COL: case t_MAT:
err(gadderf,tx,ty);
default: err(typeer,"addition");
}
case t_RFRAC:
if (!is_rfrac_t(ty)) err(gadderi,tx,ty);
return addrfrac(x,y);
case t_RFRACN:
if (!is_rfrac_t(ty)) err(gadderi,tx,ty);
z=cgetg(3,t_RFRACN); av=avma;
p1=gmul((GEN)x[1],(GEN)y[2]);
p2=gmul((GEN)x[2],(GEN)y[1]); tetpil=avma;
z[1]=lpile(av,tetpil, gadd(p1,p2));
z[2]=lmul((GEN)x[2],(GEN)y[2]); return z;
case t_VEC: case t_COL: case t_MAT:
lx = lg(x); ly = lg(y);
if (lx!=ly || tx!=ty) err(gadderi,tx,ty);
z=cgetg(ly,ty);
for (i=1; i<ly; i++)
z[i]=ladd((GEN)x[i],(GEN)y[i]);
return z;
case t_QFR: case t_QFI: err(gadderf,tx,ty);
}
err(typeer,"addition");
return NULL; /* not reached */
}
/********************************************************************/
/** **/
/** MULTIPLICATION **/
/** **/
/********************************************************************/
#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]);
/* assume z[1] was created last */
#define fix_frac_if_int_GC(z,tetpil) { if (is_pm1(z[2]))\
z = gerepileupto((long)(z+3), (GEN)z[1]);\
else\
gerepilemanyvec((long)z, tetpil, z+1, 2); }
GEN
fix_rfrac_if_pol(GEN x, GEN y)
{
if (gcmp1(y)) return x;
if (typ(y) != t_POL)
{
if (typ(x) != t_POL || gvar2(y) > varn(x))
return gdiv(x,y);
}
else if (varn(y) > varn(x)) return gdiv(x,y);
return NULL;
}
static long
mingvar(GEN x, GEN y)
{
long i = gvar(x);
long j = gvar(y);
return min(i,j);
}
GEN
mulscalrfrac(GEN x, GEN y)
{
GEN p1,z,y1,y2,cx,cy1,cy2;
long tetpil,tx;
if (gcmp0(x)) return gcopy(x);
y1=(GEN)y[1]; if (gcmp0(y1)) return gcopy(y1);
y2=(GEN)y[2]; tx = typ(x);
z = cgetg(3, t_RFRAC);
if (is_const_t(tx) || varn(x) > mingvar(y1,y2)) { cx = x; x = gun; }
else
{
p1 = ggcd(x,y2); if (isnonscalar(p1)) { x=gdeuc(x,p1); y2=gdeuc(y2,p1); }
cx = content(x); if (!gcmp1(cx)) x = gdiv(x,cx);
}
cy1 = content(y1); if (!gcmp1(cy1)) y1 = gdiv(y1,cy1);
cy2 = content(y2); if (!gcmp1(cy2)) y2 = gdiv(y2,cy2);
if (x != gun) y1 = gmul(y1,x);
x = gdiv(gmul(cx,cy1), cy2);
cy1 = numer(x);
cy2 = denom(x); tetpil = avma;
z[2] = lmul(y2, cy2);
z[1] = lmul(y1, cy1);
p1 = fix_rfrac_if_pol((GEN)z[1],(GEN)z[2]);
if (p1) return gerepileupto((long)(z+3), p1);
gerepilemanyvec((long)z,tetpil,z+1,2); return z;
}
static GEN
mulrfrac(GEN x, GEN y)
{
GEN z = cgetg(3,t_RFRAC), p1;
GEN x1 = (GEN)x[1], x2 = (GEN)x[2];
GEN y1 = (GEN)y[1], y2 = (GEN)y[2];
long tetpil;
p1 = ggcd(x1, y2); if (!gcmp1(p1)) { x1 = gdiv(x1,p1); y2 = gdiv(y2,p1); }
p1 = ggcd(x2, y1); if (!gcmp1(p1)) { x2 = gdiv(x2,p1); y1 = gdiv(y1,p1); }
tetpil = avma;
z[2] = lmul(x2,y2);
z[1] = lmul(x1,y1);
p1 = fix_rfrac_if_pol((GEN)z[1],(GEN)z[2]);
if (p1) return gerepileupto((long)(z+3), p1);
gerepilemanyvec((long)z,tetpil,z+1,2); return z;
}
GEN
to_Kronecker(GEN P, GEN Q)
{
/* P(X) = sum Mod(.,Q(Y)) * X^i, lift then set X := Y^(2n-1) */
long i,j,k,l, lx = lgef(P), N = ((lgef(Q)-3)<<1) + 1, vQ = varn(Q);
GEN p1, y = cgetg((N-2)*(lx-2) + 2, t_POL);
for (k=i=2; i<lx; i++)
{
p1 = (GEN)P[i];
if ((l=typ(p1)) == t_POLMOD) { p1 = (GEN)p1[2]; l = typ(p1); }
if (is_scalar_t(l) || varn(p1)<vQ) { y[k++] = (long)p1; j = 3; }
else
{
l = lgef(p1);
for (j=2; j < l; j++) y[k++] = p1[j];
}
if (i == lx-1) break;
for ( ; j < N; j++) y[k++] = zero;
}
y[1] = evalsigne(1)|evalvarn(vQ)|evallgef(k);
return y;
}
int
ff_poltype(GEN *x, GEN *p, GEN *pol)
{
GEN Q, P = *x, pr,p1,p2,y;
long i, lx;
if (!signe(P)) return 0;
lx = lgef(P); Q = *pol;
for (i=2; i<lx; i++)
{
p1 = (GEN)P[i]; if (typ(p1) != t_POLMOD) { Q = NULL; break; }
p2 = (GEN)p1[1];
if (Q==NULL) Q = p2;
else if (p2 != Q)
{
if (!gegal(p2, Q))
{
if (DEBUGMEM) err(warner,"different modulus in ff_poltype");
return 0;
}
if (DEBUGMEM > 2) err(warner,"different pointers in ff_poltype");
}
}
if (Q) {
*x = P = to_Kronecker(P, Q);
*pol = Q; lx = lgef(P);
}
pr = *p; y = cgetg(lx, t_POL);
for (i=lx-1; i>1; i--)
{
p1 = (GEN)P[i];
switch(typ(p1))
{
case t_INTMOD: break;
case t_INT:
if (*p) p1 = modii(p1, *p);
y[i] = (long)p1; continue;
default:
return (Q && !pr)? 1: 0;
}
p2 = (GEN)p1[1];
if (pr==NULL) pr = p2;
else if (p2 != pr)
{
if (!egalii(p2, pr))
{
if (DEBUGMEM) err(warner,"different modulus in ff_poltype");
return 0;
}
if (DEBUGMEM > 2) err(warner,"different pointers in ff_poltype");
}
y[i] = p1[2];
}
y[1] = evalsigne(1)|evalvarn(varn(P))|evallgef(lx);
*x = y; *p = pr; return (Q || pr);
}
GEN
gmul(GEN x, GEN y)
{
long tx,ty,lx,ly,vx,vy,i,j,k,l,av,tetpil;
GEN z,p1,p2,p3,p4;
if (x == y) return gsqr(x);
tx=typ(x); ty=typ(y);
if (is_const_t(tx) && is_const_t(ty))
{
if (tx>ty) { p1=x; x=y; y=p1; i=tx; tx=ty; ty=i; }
}
else
{
vx=gvar(x); vy=gvar(y);
if (!is_matvec_t(ty))
if (is_matvec_t(tx) || vx<vy || (vx==vy && tx>ty))
{
p1=x; x=y; y=p1;
i=tx; tx=ty; ty=i;
i=vx; vx=vy; vy=i;
}
if (ty==t_POLMOD) return op_polmod(gmul,x,y,tx);
}
if (is_scalar_t(ty))
{
switch(tx)
{
case t_INT:
switch(ty)
{
case t_INT: return mulii(x,y);
case t_REAL: return mulir(x,y);
case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)y[1];
(void)new_chunk(lgefint(p2)<<2);
p1=mulii(modii(x,p2),(GEN)y[2]); avma=(long)z;
z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
case t_FRAC:
if (!signe(x)) return gzero;
z=cgetg(3,t_FRAC);
p1 = mppgcd(x,(GEN)y[2]);
if (is_pm1(p1))
{
avma = (long)z;
z[2] = licopy((GEN)y[2]);
z[1] = lmulii((GEN)y[1], x);
}
else
{
x = divii(x,p1); tetpil = avma;
z[2] = ldivii((GEN)y[2], p1);
z[1] = lmulii((GEN)y[1], x);
fix_frac_if_int_GC(z,tetpil);
}
return z;
case t_FRACN: z=cgetg(3,t_FRACN);
z[1]=lmulii(x,(GEN)y[1]);
z[2]=licopy((GEN)y[2]); return z;
case t_COMPLEX: z=cgetg(3,t_COMPLEX);
z[1]=lmul(x,(GEN)y[1]);
z[2]=lmul(x,(GEN)y[2]); return z;
case t_PADIC:
if (!signe(x)) return gzero;
l=avma; p1=cgetp(y); gaffect(x,p1); tetpil=avma;
return gerepile(l,tetpil,gmul(p1,y));
case t_QUAD: z=cgetg(4,t_QUAD);
copyifstack(y[1], z[1]);
z[2]=lmul(x,(GEN)y[2]);
z[3]=lmul(x,(GEN)y[3]); return z;
}
case t_REAL:
switch(ty)
{
case t_REAL: return mulrr(x,y);
case t_FRAC: case t_FRACN:
l=avma; p1=cgetr(lg(x)); tetpil=avma; gaffect(y,p1);
p2=mulrr(p1,x); return gerepile(l,tetpil,p2);
case t_COMPLEX: z=cgetg(3,t_COMPLEX);
z[1]=lmul(x,(GEN)y[1]);
z[2]=lmul(x,(GEN)y[2]); return z;
case t_QUAD:
l=avma; p1=co8(y,lg(x)); tetpil=avma;
return gerepile(l,tetpil,gmul(p1,x));
case t_INTMOD: err(gmulerf,tx,ty);
}
case t_INTMOD:
switch(ty)
{
case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)x[1]; p1=(GEN)y[1];
if (p1==p2 || egalii(p1,p2))
{
icopyifstack(p2,z[1]);
if (!is_bigint(p2))
{
z[2] = lstoi(mulssmod(itos((GEN)x[2]),itos((GEN)y[2]), p2[2]));
return z;
}
}
else
{ p2 = mppgcd(p1,p2); z[1] = (long)p2; }
av=avma;
(void)new_chunk(lgefint(x[1]) + (lgefint(p1)<<1));
p1=mulii((GEN)x[2],(GEN)y[2]); avma=av;
z[2]=lmodii(p1,p2); return z;
case t_FRAC: case t_FRACN: z=cgetg(3,t_INTMOD); p2=(GEN)x[1];
(void)new_chunk(lgefint(p2)<<2);
p1 = mulii((GEN)y[1], mpinvmod((GEN)y[2],p2));
p1 = mulii(modii(p1,p2),(GEN)x[2]); avma=(long)z;
z[2] = lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
case t_COMPLEX: z=cgetg(3,t_COMPLEX);
z[1]=lmul(x,(GEN)y[1]);
z[2]=lmul(x,(GEN)y[2]); return z;
case t_PADIC:
l=avma; p1=cgetg(3,t_INTMOD);
p1[1]=x[1]; p1[2]=lgeti(lg(x[1]));
gaffect(y,p1); tetpil=avma;
return gerepile(l,tetpil,gmul(x,p1));
case t_QUAD: z=cgetg(4,t_QUAD);
copyifstack(y[1], z[1]);
z[2]=lmul(x,(GEN)y[2]);
z[3]=lmul(x,(GEN)y[3]); return z;
}
case t_FRAC: case t_FRACN:
switch(ty)
{
case t_FRAC:
{
GEN x1 = (GEN)x[1], x2 = (GEN)x[2];
GEN y1 = (GEN)y[1], y2 = (GEN)y[2];
z=cgetg(3,t_FRAC);
p1 = mppgcd(x1, y2);
if (!is_pm1(p1)) { x1 = divii(x1,p1); y2 = divii(y2,p1); }
p1 = mppgcd(x2, y1);
if (!is_pm1(p1)) { x2 = divii(x2,p1); y1 = divii(y1,p1); }
tetpil = avma;
z[2] = lmulii(x2,y2);
z[1] = lmulii(x1,y1);
fix_frac_if_int_GC(z,tetpil); return z;
}
case t_FRACN: z=cgetg(3,t_FRACN);
z[1]=lmulii((GEN)x[1],(GEN)y[1]);
z[2]=lmulii((GEN)x[2],(GEN)y[2]); return z;
case t_COMPLEX: z=cgetg(3,t_COMPLEX);
z[1]=lmul((GEN)y[1],x);
z[2]=lmul((GEN)y[2],x); return z;
case t_PADIC:
if (!signe(x[1])) return gzero;
l=avma; p1=cgetp(y); gaffect(x,p1); tetpil=avma;
return gerepile(l,tetpil,gmul(p1,y));
case t_QUAD: z=cgetg(4,t_QUAD);
copyifstack(y[1], z[1]);
z[2]=lmul((GEN)y[2],x);
z[3]=lmul((GEN)y[3],x); return z;
}
case t_COMPLEX:
switch(ty)
{
case t_COMPLEX: z=cgetg(3,t_COMPLEX); l=avma;
p1=gmul((GEN)x[1],(GEN)y[1]);
p2=gmul((GEN)x[2],(GEN)y[2]);
x=gadd((GEN)x[1],(GEN)x[2]);
y=gadd((GEN)y[1],(GEN)y[2]);
y=gmul(x,y); x=gadd(p1,p2);
tetpil=avma; z[1]=lsub(p1,p2); z[2]=lsub(y,x);
gerepilemanyvec(l,tetpil,z+1,2); return z;
case t_PADIC:
if (krosg(-1,(GEN)y[2]))
{
z=cgetg(3,t_COMPLEX);
z[1]=lmul((GEN)x[1],y);
z[2]=lmul((GEN)x[2],y); return z;
}
av=avma;
if (signe(y[4])) l=precp(y);
else
{
l=valp(y)+1; if (l<=0) l=1;
}
p1=cvtop(x,(GEN)y[2],l); tetpil=avma;
return gerepile(av,tetpil,gmul(p1,y));
case t_QUAD:
lx=precision(x); if (!lx) err(gmuleri,tx,ty);
l=avma; p1=co8(y,lx); tetpil=avma;
return gerepile(l,tetpil,gmul(p1,x));
}
case t_PADIC:
switch(ty)
{
case t_PADIC:
if (!egalii((GEN)x[2],(GEN)y[2])) err(gmuleri,tx,ty);
l = valp(x)+valp(y);
if (!signe(x[4])) { z=gcopy(x); setvalp(z,l); return z; }
if (!signe(y[4])) { z=gcopy(y); setvalp(z,l); return z; }
p1 = (precp(x) > precp(y))? y: x;
z=cgetp(p1); setvalp(z,l); av=avma;
modiiz(mulii((GEN)x[4],(GEN)y[4]),(GEN)p1[3],(GEN)z[4]);
avma=av; return z;
case t_QUAD:
if (kro_quad((GEN)y[1],(GEN)x[2])== -1)
{
z=cgetg(4,t_QUAD);
copyifstack(y[1], z[1]);
z[2]=lmul((GEN)y[2],x);
z[3]=lmul((GEN)y[3],x); return z;
}
l = signe(x[4])? precp(x): valp(x)+1;
av=avma; p1=cvtop(y,(GEN)x[2],l); tetpil=avma;
return gerepile(av,tetpil,gmul(p1,x));
}
case t_QUAD: z=cgetg(4,t_QUAD);
p1=(GEN)x[1]; p2=(GEN)y[1];
if (!gegal(p1,p2)) err(gmuleri,tx,ty);
copyifstack(p2, z[1]); l=avma;
p2=gmul((GEN)x[2],(GEN)y[2]);
p3=gmul((GEN)x[3],(GEN)y[3]);
p4=gmul(gneg_i((GEN)p1[2]),p3);
if (gcmp0((GEN)p1[3]))
{
tetpil=avma;
z[2]=lpile(l,tetpil,gadd(p4,p2)); l=avma;
p2=gmul((GEN)x[2],(GEN)y[3]);
p3=gmul((GEN)x[3],(GEN)y[2]); tetpil=avma;
z[3]=lpile(l,tetpil,gadd(p2,p3)); return z;
}
p1 = gadd(gmul((GEN)x[2],(GEN)y[3]), gmul((GEN)x[3],(GEN)y[2]));
tetpil=avma;
z[2]=ladd(p2,p4);
z[3]=ladd(p1,p3);
gerepilemanyvec(l,tetpil,z+2,2); return z;
}
err(bugparier,"multiplication");
}
/* here !isscalar(y) */
if (is_matvec_t(ty))
{
ly=lg(y);
if (!is_matvec_t(tx))
{
z=cgetg(ly,ty);
for (i=1; i<ly; i++) z[i]=lmul(x,(GEN)y[i]);
return z;
}
lx=lg(x);
switch(tx)
{
case t_VEC:
switch(ty)
{
case t_COL:
if (lx!=ly) err(gmuleri,tx,ty);
z=gzero; l=avma;
for (i=1; i<lx; i++)
{
p1=gmul((GEN)x[i],(GEN)y[i]);
tetpil=avma; z=gadd(z,p1);
}
return gerepile(l,tetpil,z);
case t_MAT:
if (ly==1) return cgetg(1,t_VEC);
l=lg(y[1]); if (lx!=l) err(gmuleri,tx,ty);
z=cgetg(ly,tx);
for (i=1; i<ly; i++)
{
p1=gzero; av=avma;
for (j=1; j<lx; j++)
{
p2=gmul((GEN)x[j],gcoeff(y,j,i));
tetpil=avma; p1=gadd(p1,p2);
}
z[i]=lpile(av,tetpil,p1);
}
return z;
case t_VEC: err(gmulerf,tx,ty);
default: err(typeer,"multiplication");
}
case t_COL:
switch(ty)
{
case t_VEC:
z=cgetg(ly,t_MAT);
for (i=1; i<ly; i++)
{
p1 = gmul((GEN)y[i],x);
if (typ(p1) != t_COL) err(gmuleri,tx,ty);
z[i]=(long)p1;
}
return z;
case t_MAT:
if (ly!=1 && lg(y[1])!=2) err(gmuleri,tx,ty);
z=cgetg(ly,t_MAT);
for (i=1; i<ly; i++) z[i]=lmul(gcoeff(y,1,i),x);
return z;
case t_COL: err(gmulerf,tx,ty);
default: err(typeer,"multiplication");
}
case t_MAT:
switch(ty)
{
case t_VEC:
if (lx!=2) err(gmuleri,tx,ty);
z=cgetg(ly,t_MAT);
for (i=1; i<ly; i++) z[i]=lmul((GEN)y[i],(GEN)x[1]);
return z;
case t_COL:
if (lx!=ly) err(gmuleri,tx,ty);
if (lx==1) return gcopy(y);
lx=lg(x[1]); z=cgetg(lx,t_COL);
for (i=1; i<lx; i++)
{
p1=gzero; l=avma;
for (j=1; j<ly; j++)
{
p2=gmul(gcoeff(x,i,j),(GEN)y[j]);
tetpil=avma; p1=gadd(p1,p2);
}
z[i]=lpile(l,tetpil,p1);
}
return z;
case t_MAT:
if (ly==1) return cgetg(ly,t_MAT);
if (lx != lg(y[1])) err(gmuleri,tx,ty);
z=cgetg(ly,t_MAT);
if (lx==1)
{
for (i=1; i<ly; i++) z[i]=lgetg(1,t_COL);
return z;
}
l=lg(x[1]);
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=gmul(gcoeff(x,i,k),gcoeff(y,k,j));
tetpil=avma; p1=gadd(p1,p2);
}
coeff(z,i,j)=lpile(av,tetpil,p1);
}
}
return z;
}
}
err(bugparier,"multiplication");
}
/* now !ismatvec(x and y) */
if (vx>vy || (vx==vy && is_scalar_t(tx)))
{
if (isexactzero(x)) return zeropol(vy);
if (tx == t_INT && is_pm1(x))
return (signe(x)>0) ? gcopy(y): gneg(y);
if (tx == t_POLMOD && vx == vy && ty != t_SER)
{
av = avma;
return gerepileupto(av, op_polmod(gmul, x, to_polmod(y,(GEN)x[1]), tx));
}
switch(ty)
{
case t_POL:
if (isexactzero(y)) return zeropol(vy);
ly = lgef(y); z = cgetg(ly,t_POL); z[1]=y[1];
for (i=2; i<ly; i++) z[i]=lmul(x,(GEN)y[i]);
return normalizepol_i(z,ly);
case t_SER:
if (!signe(y)) return gcopy(y);
ly=lg(y); z=cgetg(ly,t_SER);
for (i=2; i<ly; i++) z[i]=lmul(x,(GEN)y[i]);
z[1]=y[1]; return normalize(z);
case t_RFRAC: return mulscalrfrac(x,y);
case t_RFRACN: av=avma; z=cgetg(3,t_RFRACN);
z[1]=lmul(x,(GEN)y[1]);
z[2]=lcopy((GEN)y[2]); return z;
}
err(typeer,"multiplication");
}
if (tx>ty) { p1=x; x=y; y=p1; i=tx; tx=ty; ty=i; }
switch(tx)
{
case t_POL:
switch (ty)
{
case t_POL:
{
GEN a = x,b = y, p = NULL, pol = NULL;
long av = avma;
if (ff_poltype(&x,&p,&pol) && ff_poltype(&y,&p,&pol))
{
/* fprintferr("HUM"); */
if (pol && varn(x) != varn(y))
x = to_Kronecker(x,pol);
z = quickmul(x+2, y+2, lgef(x)-2, lgef(y)-2);
if (p) z = Fp_pol(z,p);
if (pol) z = from_Kronecker(z,pol);
z = gerepileupto(av, z);
}
else
{
avma = av;
z = quickmul(a+2, b+2, lgef(a)-2, lgef(b)-2);
}
setvarn(z,vx); return z;
}
case t_SER:
if (gcmp0(x)) return zeropol(vx);
if (gcmp0(y)) return zeroser(vx, valp(y)+gval(x,vx));
p1=greffe(x,lg(y),0); p2=gmul(p1,y);
free(p1); return p2;
case t_RFRAC: return mulscalrfrac(x,y);
case t_RFRACN: av=avma; z=cgetg(3,t_RFRACN);
z[1]=lmul(x,(GEN)y[1]);
z[2]=lcopy((GEN)y[2]); return z;
default: err(typeer,"multiplication");
}
case t_SER:
switch (ty)
{
case t_SER:
if (gcmp0(x) || gcmp0(y)) return zeroser(vx, valp(x)+valp(y));
lx=lg(x); ly=lg(y);
if (lx>ly) { k=ly; ly=lx; lx=k; p1=y; y=x; x=p1; }
z = cgetg(lx,t_SER);
z[1] = evalvalp(valp(x)+valp(y)) | evalvarn(vx) | evalsigne(1);
x += 2; y += 2; z += 2; lx -= 3;
p2 = (GEN)gpmalloc((lx+1)*sizeof(long));
for (i=0; i<=lx; i++)
{
p2[i] = !isexactzero((GEN)y[i]);
p1 = gzero; av = avma;
for (j=0; j<=i; j++)
if (p2[j])
p1 = gadd(p1, gmul((GEN)y[j],(GEN)x[i-j]));
z[i] = lpileupto(av,p1);
}
z -= 2; /* back to normalcy */
free(p2); return normalize(z);
case t_RFRAC: case t_RFRACN:
if (gcmp0(y)) return zeropol(vx);
if (gcmp0(x)) return zeroser(vx, valp(x)+gval(y,vx));
l=avma; p1=gmul((GEN)y[1],x); tetpil=avma;
return gerepile(l,tetpil,gdiv(p1,(GEN)y[2]));
default: err(typeer,"multiplication");
}
/* (tx,ty) == t_RFRAC <==> ty == t_RFRAC */
case t_RFRAC: return mulrfrac(x,y);
case t_RFRACN:
if (!is_rfrac_t(ty)) err(typeer,"multiplication");
av=avma; z=cgetg(3,ty);
z[1]=lmul((GEN)x[1],(GEN)y[1]);
z[2]=lmul((GEN)x[2],(GEN)y[2]); return z;
}
if (tx==ty)
switch(tx)
{
case t_QFI: return compimag(x,y);
case t_QFR: return compreal(x,y);
}
err(typeer,"multiplication");
return NULL; /* not reached */
}
/********************************************************************/
/** **/
/** DIVISION **/
/** **/
/********************************************************************/
static
GEN divrfracscal(GEN x, GEN y)
{
long Y[3]; Y[1]=un; Y[2]=(long)y;
return mulrfrac(x,Y);
}
static
GEN divscalrfrac(GEN x, GEN y)
{
long Y[3]; Y[1]=y[2]; Y[2]=y[1];
return mulscalrfrac(x,Y);
}
static
GEN divrfrac(GEN x, GEN y)
{
long Y[3]; Y[1]=y[2]; Y[2]=y[1];
return mulrfrac(x,Y);
}
GEN
gdiv(GEN x, GEN y)
{
long tx = typ(x), ty = typ(y), lx,ly,vx,vy,i,j,k,l,av,tetpil;
GEN z,p1,p2,p3;
if (tx==t_INT && is_const_t(ty))
{
switch (signe(x))
{
case 0:
if (gcmp0(y)) err(gdiver2);
if (ty != t_INTMOD) return gzero;
z = cgetg(3,t_INTMOD); icopyifstack(y[1],z[1]); z[2]=zero;
return z;
case 1:
if (is_pm1(x)) return ginv(y);
break;
case -1:
if (is_pm1(x)) { av = avma; return gerepileupto(av, ginv(gneg(y))); }
}
switch(ty)
{
case t_INT:
av=avma; z=dvmdii(x,y,&p1);
if (p1==gzero) return z;
(void)new_chunk((lgefint(x) + lgefint(y)) << 2);
p1 = mppgcd(y,p1);
avma=av; z=cgetg(3,t_FRAC);
if (is_pm1(p1)) {
z[1]=licopy(x);
z[2]=licopy(y);
}
else {
z[1]=ldivii(x,p1);
z[2]=ldivii(y,p1);
}
fix_frac(z); return z;
case t_REAL:
return divir(x,y);
case t_INTMOD:
z=cgetg(3,t_INTMOD); p2=(GEN)y[1];
(void)new_chunk(lgefint(p2)<<2);
p1=mulii(modii(x,p2), mpinvmod((GEN)y[2],p2)); avma=(long)z;
z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
case t_FRAC:
z=cgetg(3,t_FRAC);
p1 = mppgcd(x,(GEN)y[1]);
if (is_pm1(p1))
{
avma = (long)z; tetpil = 0;
z[2] = licopy((GEN)y[1]);
}
else
{
x = divii(x,p1); tetpil = avma;
z[2] = ldivii((GEN)y[1], p1);
}
z[1] = lmulii((GEN)y[2], x);
fix_frac(z);
if (tetpil)
{ fix_frac_if_int_GC(z,tetpil); }
else
fix_frac_if_int(z);
return z;
case t_FRACN:
z=cgetg(3,t_FRACN);
z[1]=lmulii((GEN)y[2], x);
z[2]=licopy((GEN)y[1]);
fix_frac(z); return z;
case t_PADIC:
l=avma; p1=cgetp(y); gaffect(x,p1); tetpil=avma;
return gerepile(l,tetpil,gdiv(p1,y));
case t_COMPLEX: case t_QUAD:
l=avma; p1=gnorm(y); p2=gmul(x,gconj(y)); tetpil=avma;
return gerepile(l,tetpil,gdiv(p2,p1));
}
}
if (gcmp0(y)) err(gdiver2);
if (is_const_t(tx) && is_const_t(ty))
{
switch(tx)
{
case t_REAL:
switch(ty)
{
case t_INT:
return divri(x,y);
case t_REAL:
return divrr(x,y);
case t_FRAC: case t_FRACN:
l=avma; p1=cgetg(lg(x),t_REAL); gaffect(y,p1);
return gerepile(l,(long)p1,divrr(x,p1));
case t_COMPLEX: z=cgetg(3,t_COMPLEX);
l=avma; p1=gnorm(y);
p2=gmul(x,(GEN)y[1]);
p3=gmul(x,(GEN)y[2]);
if (!gcmp0(p3)) p3 = gneg_i(p3);
tetpil=avma;
z[1]=ldiv(p2,p1);
z[2]=ldiv(p3,p1);
gerepilemanyvec(l,tetpil,z+1,2); return z;
case t_QUAD:
l=avma; p1=co8(y,lg(x)); tetpil=avma;
return gerepile(l,tetpil,gdiv(x,p1));
case t_INTMOD: case t_PADIC: err(gdiverf,tx,ty);
}
case t_INTMOD:
switch(ty)
{
case t_INT: z=cgetg(3,t_INTMOD); p2=(GEN)x[1];
(void)new_chunk(lgefint(p2)<<2);
p1=mulii((GEN)x[2], mpinvmod(y,p2)); avma=(long)z;
z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)x[1]; p1=(GEN)y[1];
if (p1==p2 || egalii(p1,p2))
{ icopyifstack(p2,z[1]); }
else
{ p2 = mppgcd(p1,p2); z[1] = (long)p2; }
av=avma; (void)new_chunk(lgefint(x[1]) + (lgefint(p1) << 1));
p1=mulii((GEN)x[2], mpinvmod((GEN)y[2],p2)); avma=av;
z[2]=lmodii(p1,p2); return z;
case t_FRAC: z=cgetg(3,t_INTMOD); p2=(GEN)x[1];
(void)new_chunk(lgefint(p2)<<2);
p1=mulii((GEN)y[2], mpinvmod((GEN)y[1],p2));
p1=mulii(modii(p1,p2),(GEN)x[2]); avma=(long)z;
z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
case t_FRACN:
l=avma; p1=gred(y); tetpil=avma;
return gerepile(l,tetpil,gdiv(x,p1));
case t_COMPLEX: case t_QUAD:
l=avma; p1=gnorm(y); p2=gmul(x,gconj(y)); tetpil=avma;
return gerepile(l,tetpil,gdiv(p2,p1));
case t_PADIC:
l=avma; p1=cgetg(3,t_INTMOD); p1[1]=x[1]; p1[2]=lgeti(lg(x[1]));
gaffect(y,p1); tetpil=avma; return gerepile(l,tetpil,gdiv(x,p1));
case t_REAL: err(gdiverf,tx,ty);
}
case t_FRAC: case t_FRACN:
switch(ty)
{
case t_INT:
z = cgetg(3, tx);
if (tx == t_FRAC)
{
p1 = mppgcd(y,(GEN)x[1]);
if (is_pm1(p1))
{
avma = (long)z; tetpil = 0;
z[1] = licopy((GEN)x[1]);
}
else
{
y = divii(y,p1); tetpil = avma;
z[1] = ldivii((GEN)x[1], p1);
}
}
else
{
tetpil = 0;
z[1] = licopy((GEN)x[1]);
}
z[2] = lmulii((GEN)x[2],y);
fix_frac(z);
if (tetpil) fix_frac_if_int_GC(z,tetpil);
return z;
case t_REAL:
l=avma; p1=cgetg(lg(y),t_REAL); gaffect(x,p1);
p2=divrr(p1,y); return gerepile(l,(long)p1,p2);
case t_INTMOD: z=cgetg(3,t_INTMOD); p2=(GEN)y[1];
(void)new_chunk(lgefint(p2)<<2);
p1=mulii((GEN)y[2],(GEN)x[2]);
p1=mulii(mpinvmod(p1,p2), modii((GEN)x[1],p2)); avma=(long)z;
z[2]=lmodii(p1,p2); icopyifstack(p2,z[1]); return z;
case t_FRAC: if (tx == t_FRACN) ty=t_FRACN;
case t_FRACN:
z = cgetg(3,ty);
if (ty == t_FRAC)
{
GEN x1 = (GEN)x[1], x2 = (GEN)x[2];
GEN y1 = (GEN)y[1], y2 = (GEN)y[2];
p1 = mppgcd(x1, y1);
if (!is_pm1(p1)) { x1 = divii(x1,p1); y1 = divii(y1,p1); }
p1 = mppgcd(x2, y2);
if (!is_pm1(p1)) { x2 = divii(x2,p1); y2 = divii(y2,p1); }
tetpil = avma;
z[2] = lmulii(x2,y1);
z[1] = lmulii(x1,y2);
fix_frac(z);
fix_frac_if_int_GC(z,tetpil);
}
else
{
z[1]=lmulii((GEN)x[1],(GEN)y[2]);
z[2]=lmulii((GEN)x[2],(GEN)y[1]);
fix_frac(z);
}
return z;
case t_COMPLEX: z=cgetg(3,t_COMPLEX);
l=avma; p1=gnorm(y);
p2=gmul(x,(GEN)y[1]);
p3=gmul(x,(GEN)y[2]);
if(!gcmp0(p3)) p3 = gneg_i(p3);
tetpil=avma;
z[1]=ldiv(p2,p1); z[2]=ldiv(p3,p1);
gerepilemanyvec(l,tetpil,z+1,2); return z;
case t_PADIC:
if (!signe(x[1])) return gzero;
l=avma; p1=cgetp(y); gaffect(x,p1);
tetpil=avma; return gerepile(l,tetpil,gdiv(p1,y));
case t_QUAD:
l=avma; p1=gnorm(y); p2=gmul(x,gconj(y)); tetpil=avma;
return gerepile(l,tetpil,gdiv(p2,p1));
}
case t_COMPLEX:
switch(ty)
{
case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_FRACN:
z=cgetg(3,t_COMPLEX);
z[1]=ldiv((GEN)x[1],y);
z[2]=ldiv((GEN)x[2],y); return z;
case t_COMPLEX:
l=avma; p1=gnorm(y); p2=gconj(y); p2=gmul(x,p2); tetpil=avma;
return gerepile(l,tetpil, gdiv(p2,p1));
case t_PADIC:
if (krosg(-1,(GEN)y[2])== -1)
{
z=cgetg(3,t_COMPLEX);
z[1]=ldiv((GEN)x[1],y);
z[2]=ldiv((GEN)x[2],y); return z;
}
av=avma; p1=cvtop(x,(GEN)y[2],precp(y)); tetpil=avma;
return gerepile(av,tetpil,gdiv(p1,y));
case t_QUAD:
lx=precision(x); if (!lx) err(gdiveri,tx,ty);
l=avma; p1=co8(y,lx); tetpil=avma;
return gerepile(l,tetpil,gdiv(x,p1));
}
case t_PADIC:
switch(ty)
{
case t_INT: case t_FRAC: case t_FRACN:
l=avma;
if (signe(x[4])) { p1=cgetp(x); gaffect(y,p1); }
else p1=cvtop(y,(GEN)x[2],(valp(x)>0)?valp(x):1);
tetpil=avma; return gerepile(l,tetpil,gdiv(x,p1));
case t_INTMOD:
l=avma; p1=cgetg(3,t_INTMOD);
p1[1]=y[1]; p1[2]=lgeti(lg(y[1]));
gaffect(x,p1); tetpil=avma;
return gerepile(l,tetpil,gdiv(p1,y));
case t_PADIC:
if (!egalii((GEN)x[2],(GEN)y[2])) err(gdiveri,tx,ty);
if (!signe(x[4]))
{
z=gcopy(x); setvalp(z,valp(x)-valp(y));
return z;
}
p1=(precp(x)>precp(y)) ? y : x;
z=cgetp(p1); l=avma;
setvalp(z,valp(x)-valp(y));
p2=mpinvmod((GEN)y[4],(GEN)p1[3]);
modiiz(mulii((GEN)x[4],p2),(GEN)p1[3],(GEN)z[4]);
avma=l; return z;
case t_COMPLEX: case t_QUAD:
l=avma; p1=gmul(x,gconj(y)); p2=gnorm(y); tetpil=avma;
return gerepile(l,tetpil,gdiv(p1,p2));
case t_REAL:
err(talker,"forbidden division p-adic/R");
}
case t_QUAD:
switch (ty)
{
case t_INT: case t_INTMOD: case t_FRAC: case t_FRACN:
z=cgetg(4,t_QUAD);
copyifstack(x[1], z[1]);
for (i=2; i<4; i++) z[i]=ldiv((GEN)x[i],y);
return z;
case t_REAL:
l=avma; p1=co8(x,lg(y)); tetpil=avma;
return gerepile(l,tetpil,gdiv(p1,y));
case t_PADIC:
l=avma; p1=cvtop(x,(GEN)y[2],precp(y));
tetpil=avma; return gerepile(l,tetpil,gdiv(p1,y));
case t_COMPLEX:
ly=precision(y); if (!ly) err(gdiveri,tx,ty);
l=avma; p1=co8(x,ly); tetpil=avma;
return gerepile(l,tetpil,gdiv(p1,y));
case t_QUAD:
k=x[1]; l=y[1];
if (!gegal((GEN)k,(GEN)l)) err(gdiveri,tx,ty);
l=avma; p1=gnorm(y); p2=gmul(x,gconj(y)); tetpil=avma;
return gerepile(l,tetpil,gdiv(p2,p1));
}
}
err(bugparier,"division");
}
vx=gvar(x); vy=gvar(y);
if (ty==t_POLMOD && (tx==t_POLMOD || vy<vx))
{
z=cgetg(3,t_POLMOD);
if (tx==t_POLMOD)
{
k=x[1]; l=y[1];
if (gegal((GEN)k,(GEN)l))
{
copyifstack(k, z[1]); av=avma;
p1 = ginvmod((GEN)y[2],(GEN)z[1]);
p2 = gmul((GEN)x[2],p1);
}
else
{
vx=varn(x[1]); vy=varn(y[1]);
if (vx==vy)
{
z[1]=lgcd((GEN)k,(GEN)l); av=avma;
p1=ginvmod((GEN)y[2],(GEN)z[1]);
p2=gmul((GEN)x[2],p1);
}
else
{
if (vx<vy)
{ copyifstack(k,z[1]); av=avma; p2=gdiv((GEN)x[2],y); }
else
{
copyifstack(l,z[1]); av=avma;
p1 = ginvmod((GEN)y[2],(GEN)z[1]);
p2 = gmul(x, p1);
}
}
}
p2 = gmod(p2,(GEN)z[1]);
}
else
{
copyifstack(y[1],z[1]); av=avma;
p1 = ginvmod((GEN)y[2],(GEN)y[1]);
p2 = gmul(x,p1);
}
z[2]=lpileupto(av, p2); return z;
}
if (tx == t_POLMOD && vx<vy)
{
z=cgetg(3,t_POLMOD);
copyifstack(x[1],z[1]);
z[2]=ldiv((GEN)x[2],y); return z;
}
if (vx == vy)
{
av = avma;
if (tx == t_POLMOD)
return gerepileupto(av, gdiv(x, to_polmod(y,(GEN)x[1])));
if (ty == t_POLMOD)
return gerepileupto(av, gdiv(to_polmod(x,(GEN)y[1]), y));
}
/* now x and y are not both is_scalar_t */
lx = lg(x);
if ((vx<vy && (!is_matvec_t(tx) || !is_matvec_t(ty)))
|| (vx==vy && is_scalar_t(ty)) || (is_matvec_t(tx) && !is_matvec_t(ty)))
{
if (tx == t_RFRAC) return divrfracscal(x,y);
z = cgetg(lx,tx);
if (tx == t_RFRACN)
{
z[2]=lmul((GEN)x[2],y);
z[1]=lcopy((GEN)x[1]); return z;
}
switch(tx)
{
case t_POL: lx = lgef(x);
case t_SER: z[1] = x[1];
case t_VEC: case t_COL: case t_MAT:
if (ty == t_POLMOD || ty == t_INTMOD)
{
if (!gcmp1(y)) y = ginv(y); /* garbage, left alone */
for (i=lontyp[tx]; i<lx; i++) z[i]=lmul((GEN)x[i],y);
return z;
}
else
for (i=lontyp[tx]; i<lx; i++) z[i]=ldiv((GEN)x[i],y);
return z;
}
err(typeer,"division");
}
ly=lg(y);
if (vy<vx || (vy==vx && is_scalar_t(tx)))
{
switch(ty)
{
case t_POL:
if (lgef(y)==3) return gdiv(x,(GEN)y[2]);
if (isexactzero(x)) return zeropol(vy);
av=avma; z=cgetg(3,t_RFRAC); z[1]=(long)x; z[2]=(long)y;
return gerepileupto(av,gred_rfrac(z));
case t_SER:
if (gcmp0(x))
{
l=avma; p1=ginv(y); tetpil=avma; /* a ameliorer !!!! */
return gerepile(l,tetpil,gmul(x,p1));
}
p1 = (GEN)gpmalloc(ly*sizeof(long));
p1[0] = evaltyp(t_SER) | evallg(ly);
p1[1] = evalsigne(1) | evalvalp(0) | evalvarn(vy);
p1[2] = (long)x; for (i=3; i<ly; i++) p1[i]=zero;
y = gdiv(p1,y); free(p1); return y;
case t_RFRAC: return divscalrfrac(x,y);
case t_RFRACN: z=cgetg(ly,t_RFRACN);
z[1]=lmul(x,(GEN)y[2]);
z[2]=lcopy((GEN)y[1]); return z;
case t_MAT:
if (ly==1 || lg(y[1])!=ly) err(gdiveri,tx,ty);
l=avma; p1=invmat(y); tetpil=avma;
return gerepile(l,tetpil,gmul(x,p1));
case t_VEC: case t_COL: err(gdiverf,tx,ty);
}
err(typeer,"division");
}
/* ici vx=vy et tx>=10 et ty>=10*/
switch(tx)
{
case t_POL:
switch(ty)
{
case t_POL:
if (lgef(y)==3) return gdiv(x,(GEN)y[2]);
if (isexactzero(x)) return zeropol(vy);
av=avma; z=cgetg(3,t_RFRAC); z[1]=(long)x; z[2]=(long)y;
return gerepileupto(av,gred_rfrac(z));
case t_SER:
if (gcmp0(x)) return zeropol(vx);
p1=greffe(x,ly,0); p2=gdiv(p1,y);
free(p1); return p2;
case t_RFRAC: return divscalrfrac(x,y);
case t_RFRACN: z=cgetg(ly,t_RFRACN);
z[1]=lmul(x,(GEN)y[2]);
z[2]=lcopy((GEN)y[1]); return z;
case t_VEC: case t_COL: case t_MAT: err(gdiverf,tx,ty);
default: err(typeer,"division");
}
case t_SER:
switch(ty)
{
case t_POL:
p1=greffe(y,lx,0); p2=gdiv(x,p1);
free(p1); return p2;
case t_SER:
{
GEN y_lead;
l = valp(x) - valp(y);
if (gcmp0(x)) return zeroser(vx,l);
y_lead = (GEN)y[2];
if (gcmp0(y_lead)) /* normalize denominator if leading term is 0 */
{
err(warner,"normalizing a series with 0 leading term");
for (i=3,y++; i<ly; i++,y++)
{
y_lead = (GEN)y[2]; ly--; l--;
if (!gcmp0(y_lead)) break;
}
if (i==ly) err(gdiver2);
}
if (ly < lx) lx = ly;
p2 = (GEN)gpmalloc(lx*sizeof(long));
for (i=3; i<lx; i++)
{
p1 = (GEN)y[i];
if (isexactzero(p1)) p2[i] = 0;
else
{
av = avma; p2[i] = lclone(gneg_i(p1));
avma = av;
}
}
z = cgetg(lx,t_SER);
z[1] = evalvalp(l) | evalvarn(vx) | evalsigne(1);
z[2] = ldiv((GEN)x[2], y_lead);
for (i=3; i<lx; i++)
{
av=avma; p1 = (GEN)x[i];
for (j=2; j<i; j++)
{
l = i-j+2;
if (p2[l])
p1 = gadd(p1, gmul((GEN)z[j], (GEN)p2[l]));
}
tetpil=avma; z[i]=lpile(av,tetpil, gdiv(p1,y_lead));
}
for (i=3; i<lx; i++)
if (p2[i]) gunclone((GEN)p2[i]);
free(p2); return z;
}
case t_RFRAC: case t_RFRACN:
l=avma; p2=gmul(x,(GEN)y[2]); tetpil=avma;
return gerepile(l,tetpil,gdiv(p2,(GEN)y[1]));
case t_VEC: case t_COL: case t_MAT: err(gdiverf,tx,ty);
default: err(typeer,"division");
}
case t_RFRAC: case t_RFRACN:
switch(ty)
{
case t_POL:
if (tx==t_RFRAC) return divrfracscal(x,y);
z=cgetg(3,t_RFRACN);
z[2]=lmul((GEN)x[2],y);
z[1]=lcopy((GEN)x[1]); return z;
case t_SER:
l=avma; p2=gmul((GEN)x[2],y); tetpil=avma;
return gerepile(l,tetpil, gdiv((GEN)x[1],p2));
case t_RFRAC: case t_RFRACN:
if (tx == t_RFRACN) ty=t_RFRACN;
if (ty != t_RFRACN) return divrfrac(x,y);
z=cgetg(3,t_RFRACN);
z[1]=lmul((GEN)x[1],(GEN)y[2]);
z[2]=lmul((GEN)x[2],(GEN)y[1]); return z;
case t_VEC: case t_COL: case t_MAT: err(gdiverf,tx,ty);
default: err(typeer,"division");
}
case t_VEC: case t_COL: case t_MAT:
if (!is_matvec_t(ty))
{
z=cgetg(lx,tx);
for (i=1; i<lx; i++) z[i]=ldiv((GEN)x[i],y);
return z;
}
if (ty!=t_MAT || ly==1 || lg(y[1])!=ly) err(gdiveri,tx,ty);
l=avma; p1=invmat(y); tetpil=avma;
return gerepile(l,tetpil,gmul(x,p1));
}
if (tx==ty)
{
l=signe(y[2]); setsigne(y[2],-l);
switch(tx)
{
case t_QFI: z = compimag(x,y);
setsigne(y[2],l); return z;
case t_QFR:
k=signe(y[4]); setsigne(y[4],-k); z=compreal(x,y);
setsigne(y[2],l); setsigne(y[4],k); return z;
}
}
err(typeer,"division");
return NULL; /* not reached */
}
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gen1.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