OpenXM_contrib/pari-2.2/src/basemath/polarit2.c - diff

Return to polarit2.c CVS log

Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath

Diff for /OpenXM_contrib/pari-2.2/src/basemath/Attic/polarit2.c between version 1.1 and 1.2

version 1.1, 2001/10/02 11:17:05

version 1.2, 2002/09/11 07:26:52

Line 24 Foundation, Inc., 59 Temple Place - Suite 330, Boston,

#define swap(a,b) { GEN _x = a; a = b; b = _x; }

#define lswap(a,b) { long _x = a; a = b; b = _x; }

#define pswap(a,b) { GEN *_x = a; a = b; b = _x; }

#define both_even(a,b) ((((a)|(b))&1) == 0)

#define both_odd(a,b) ((a)&(b)&1)

#define addshift(x,y) addshiftpol((x),(y),1)

GEN gassoc_proto(GEN f(GEN,GEN),GEN,GEN);

extern GEN FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p);

extern GEN ZY_ZXY_resultant(GEN A, GEN B0, long *lambda);

extern GEN addshiftpol(GEN x, GEN y, long d);

extern GEN cauchy_bound(GEN p);

extern GEN gassoc_proto(GEN f(GEN,GEN),GEN,GEN);

extern GEN gauss_intern(GEN a, GEN b);

extern GEN indexpartial(GEN P, GEN DP);

extern GEN qf_disc(GEN x, GEN y, GEN z);

extern GEN to_polmod(GEN x, GEN mod);

extern GEN vconcat(GEN Q1, GEN Q2);

extern int approx_0(GEN x, GEN y);

extern long FpX_split_berlekamp(GEN *t, GEN pp);

extern long u_center(ulong u, ulong p, ulong ps2);

extern void gerepilemanycoeffs2(gpmem_t av, GEN x, int n, GEN y, int o);

GEN matratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN denom);

GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den);

/* compute Newton sums S_1(P), ... , S_n(P). S_k(P) = sum a_j^k, a_j root of P

* If N != NULL, assume p-adic roots and compute mod N [assume integer coeffs]

* If T != NULL, compute mod (T,N) [assume integer coeffs if N != NULL]

* If y0!= NULL, precomputed i-th powers, i=1..m, m = length(y0) */

GEN

polsym(GEN x, long n)

polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N)

{

long av1,av2,dx=degpol(x),i,k;

long dP=degpol(P), i, k, m;

GEN s,y,x_lead;

gpmem_t av1, av2;

GEN s,y,P_lead;

if (n<0) err(impl,"polsym of a negative n");

if (typ(x) != t_POL) err(typeer,"polsym");

if (typ(P) != t_POL) err(typeer,"polsym");

if (!signe(x)) err(zeropoler,"polsym");

if (!signe(P)) err(zeropoler,"polsym");

y=cgetg(n+2,t_COL); y[1]=lstoi(dx);

y = cgetg(n+2,t_COL);

x_lead=(GEN)x[dx+2]; if (gcmp1(x_lead)) x_lead=NULL;

if (y0)

for (k=1; k<=n; k++)

{

av1=avma; s = (dx>=k)? gmulsg(k,(GEN)x[dx+2-k]): gzero;

if (typ(y0) != t_COL) err(typeer,"polsym_gen");

for (i=1; i<k && i<=dx; i++)

m = lg(y0)-1;

s = gadd(s,gmul((GEN)y[k-i+1],(GEN)x[dx+2-i]));

for (i=1; i<=m; i++) y[i] = lcopy((GEN)y0[i]);

if (x_lead) s = gdiv(s,x_lead);

av2=avma; y[k+1]=lpile(av1,av2,gneg(s));

}

else

{

m = 1;

y[1] = lstoi(dP);

}

P += 2; /* strip codewords */

P_lead = (GEN)P[dP]; if (gcmp1(P_lead)) P_lead = NULL;

if (P_lead)

{

if (N) P_lead = FpXQ_inv(P_lead,T,N);

else if (T) P_lead = QX_invmod(P_lead,T);

}

for (k=m; k<=n; k++)

{

av1 = avma; s = (dP>=k)? gmulsg(k,(GEN)P[dP-k]): gzero;

for (i=1; i<k && i<=dP; i++)

s = gadd(s, gmul((GEN)y[k-i+1],(GEN)P[dP-i]));

if (N)

{

if (T) s = FpX_res(FpX_red(s,N), T, N);

else s = modii(s, N);

if (P_lead) s = Fq_mul(s, P_lead, T, N);

}

else if (T)

{

s = gres(s, T);

if (P_lead) s = gres(gmul(s, P_lead), T);

}

else

if (P_lead) s = gdiv(s, P_lead);

av2 = avma; y[k+1] = lpile(av1,av2, gneg(s));

}

return y;

}

GEN

polsym(GEN x, long n)

{

return polsym_gen(x, NULL, n, NULL,NULL);

}

static int (*polcmp_coeff_cmp)(GEN,GEN);

/* assume x and y are polynomials in the same variable whose coeffs can be

Line 60 polcmp(GEN x, GEN y)

Line 118 polcmp(GEN x, GEN y)

{

long i, lx = lgef(x), ly = lgef(y);

int fl;

#if 0 /* for efficiency */

if (typ(x) != t_POL || typ(y) != t_POL)

err(talker,"not a polynomial in polcmp");

#endif

if (lx > ly) return 1;

if (lx < ly) return -1;

for (i=lx-1; i>1; i--)

Line 78 sort_factor(GEN y, int (*cmp)(GEN,GEN))

Line 132 sort_factor(GEN y, int (*cmp)(GEN,GEN))

{

int (*old)(GEN,GEN) = polcmp_coeff_cmp;

polcmp_coeff_cmp = cmp;

(void)sort_factor_gen(y,polcmp);

(void)sort_factor_gen(y,&polcmp);

polcmp_coeff_cmp = old; return y;

}

Line 86 sort_factor(GEN y, int (*cmp)(GEN,GEN))

Line 140 sort_factor(GEN y, int (*cmp)(GEN,GEN))

GEN

sort_factor_gen(GEN y, int (*cmp)(GEN,GEN))

{

long n,i, av = avma;

long n, i;

gpmem_t av = avma;

GEN a,b,A,B,w;

a = (GEN)y[1]; n = lg(a); A = new_chunk(n);

b = (GEN)y[2]; B = new_chunk(n);

Line 96 sort_factor_gen(GEN y, int (*cmp)(GEN,GEN))

Line 151 sort_factor_gen(GEN y, int (*cmp)(GEN,GEN))

avma = av; return y;

}

GEN

sort_vecpol(GEN a)

{

long n, i;

gpmem_t av = avma;

GEN A,w;

n = lg(a); A = new_chunk(n);

w = gen_sort(a, cmp_IND | cmp_C, cmp_pol);

for (i=1; i<n; i++) A[i] = a[w[i]];

for (i=1; i<n; i++) a[i] = A[i];

avma = av; return a;

}

GEN

centermodii(GEN x, GEN p, GEN po2)

{

GEN y = modii(x,p);

if (cmpii(y,po2) > 0) return subii(y,p);

return y;

}

long

s_centermod(long x, ulong pp, ulong pps2)

{

long y = x % (long)pp;

if (y < 0) y += pp;

return u_center(y, pp,pps2);

}

/* for internal use */

GEN

centermod_i(GEN x, GEN p, GEN ps2)

{

long av,i,lx;

long i, lx;

GEN y,p1;

gpmem_t av;

GEN y;

if (!ps2) ps2 = shifti(p,-1);

switch(typ(x))

{

case t_INT:

case t_INT: return centermodii(x,p,ps2);

y=modii(x,p);

if (cmpii(y,ps2)>0) return subii(y,p);

return y;

case t_POL: lx=lgef(x);

case t_POL: lx = lgef(x);

y=cgetg(lx,t_POL); y[1]=x[1];

y = cgetg(lx,t_POL); y[1] = x[1];

for (i=2; i<lx; i++)

{

av=avma; p1=modii((GEN)x[i],p);

av = avma;

if (cmpii(p1,ps2)>0) p1=subii(p1,p);

y[i] = lpileuptoint(av, centermodii((GEN)x[i],p,ps2));

y[i]=lpileupto(av,p1);

}

return normalizepol_i(y, lx);

case t_COL: lx = lg(x);

y = cgetg(lx,t_COL);

for (i=1; i<lx; i++) y[i] = (long)centermodii((GEN)x[i],p,ps2);

return y;

case t_COL: lx=lg(x);

case t_MAT: lx = lg(x);

y=cgetg(lx,t_COL);

y = cgetg(lx,t_MAT);

for (i=1; i<lx; i++)

for (i=1; i<lx; i++) y[i] = (long)centermod_i((GEN)x[i],p,ps2);

{

p1=modii((GEN)x[i],p);

if (cmpii(p1,ps2)>0) p1=subii(p1,p);

y[i]=(long)p1;

}

return y;

case t_VECSMALL: lx = lg(x);

{

ulong pp = itou(p), pps2 = itou(ps2);

y = cgetg(lx,t_VECSMALL);

for (i=1; i<lx; i++) y[i] = s_centermod(x[i], pp, pps2);

return y;

}

return x;

}

Line 141 centermod(GEN x, GEN p) { return centermod_i(x,p,NULL)

Line 230 centermod(GEN x, GEN p) { return centermod_i(x,p,NULL)

static GEN

polidivis(GEN x, GEN y, GEN bound)

{

long dx,dy,dz,i,j,av, vx = varn(x);

long dx, dy, dz, i, j, vx = varn(x);

gpmem_t av;

GEN z,p1,y_lead;

dy=degpol(y);

Line 153 polidivis(GEN x, GEN y, GEN bound)

Line 243 polidivis(GEN x, GEN y, GEN bound)

if (gcmp1(y_lead)) y_lead = NULL;

p1 = (GEN)x[dx];

z[dz]=y_lead? ldivii(p1,y_lead): licopy(p1);

z[dz]=y_lead? (long)diviiexact(p1,y_lead): licopy(p1);

for (i=dx-1; i>=dy; i--)

{

av=avma; p1=(GEN)x[i];

av = avma; p1 = (GEN)x[i];

for (j=i-dy+1; j<=i && j<=dz; j++)

p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-j]));

if (y_lead) { p1 = gdiv(p1,y_lead); if (typ(p1)!=t_INT) return NULL; }

if (y_lead) p1 = diviiexact(p1,y_lead);

if (absi_cmp(p1, bound) > 0) return NULL;

p1 = gerepileupto(av,p1);

z[i-dy] = (long)p1;

Line 176 polidivis(GEN x, GEN y, GEN bound)

Line 266 polidivis(GEN x, GEN y, GEN bound)

if (!i) break;

}

z -= 2;

z[1]=evalsigne(1) | evallgef(dz+3) | evalvarn(vx);

z[1] = evalsigne(1) | evallgef(dz+3) | evalvarn(vx);

return z;

}

static long

min_deg(long jmax,ulong tabbit[])

{

long j, k = 1, j1 = 2;

for (j=jmax; j>=0; j--)

{

for ( ; k<=15; k++)

{

if (tabbit[j]&j1) return ((jmax-j)<<4)+k;

j1<<=1;

}

k = 0; j1 = 1;

}

return (jmax<<4)+15;

}

/* tabkbit is a bit vector (only lowest 32 bits of each word are used

* on 64bit architecture): reading from right to left, bit i+1 is set iff

* degree i is attainable from the factorisation mod p.

* record N modular factors of degree d. */

static void

record_factors(long N, long d, long jmax, ulong *tabkbit, ulong *tmp)

{

long n,j, a = d>>4, b = d&15; /* d = 16 a + b */

ulong *tmp2 = tmp - a;

for (n = 1; n <= N; n++)

{

ulong pro, rem = 0;

for (j=jmax; j>=a; j--)

{

pro = tabkbit[j] << b;

tmp2[j] = (pro&0xffff) + rem; rem = pro >> 16;

}

for (j=jmax-a; j>=0; j--) tabkbit[j] |= tmp[j];

}

/***********************************************************************/

/** **/

/** QUADRATIC HENSEL LIFT (adapted from V. Shoup's NTL) **/

/** **/

/***********************************************************************/

/* Setup for divide/conquer quadratic Hensel lift

* a = set of k t_POL in Z[X] corresponding to factors over Fp

* a = set of k t_POL in Z[X] = factors over Fp (T=NULL) or Fp[Y]/(T)

* V = set of products of factors built as follows

* 1) V[1..k] = initial a

* 2) iterate:

Line 238 record_factors(long N, long d, long jmax, ulong *tabkb

Line 288 record_factors(long N, long d, long jmax, ulong *tabkb

* link[i] = -j if V[i] = a[j]

* j if V[i] = V[j] * V[j+1]

* Arrays (link, V, W) pre-allocated for 2k - 2 elements */

#if 0 /* restricted to Fp */

static void

BuildTree(GEN link, GEN V, GEN W, GEN a, GEN p)

{

long k = lg(a)-1;

long i, j, s, minp, mind;

for (i=1; i<=k; i++) { V[i] = a[i]; link[i] = -i; }

for (j=1; j <= 2*k-5; j+=2,i++)

{

minp = j;

mind = degpol(V[j]);

for (s=j+1; s<i; s++)

if (degpol(V[s]) < mind) { minp = s; mind = degpol(V[s]); }

lswap(V[j], V[minp]);

lswap(link[j], link[minp]);

minp = j+1;

mind = degpol(V[j+1]);

for (s=j+2; s<i; s++)

if (degpol(V[s]) < mind) { minp = s; mind = degpol(V[s]); }

lswap(V[j+1], V[minp]);

lswap(link[j+1], link[minp]);

V[i] = (long)FpX_mul((GEN)V[j], (GEN)V[j+1], p);

link[i] = j;

}

for (j=1; j <= 2*k-3; j+=2)

{

GEN d, u, v;

d = FpX_extgcd((GEN)V[j], (GEN)V[j+1], p, &u, &v);

if (degpol(d) > 0) err(talker, "relatively prime polynomials expected");

d = (GEN)d[2];

if (!gcmp1(d))

{

d = mpinvmod(d, p);

u = FpX_Fp_mul(u, d, p);

v = FpX_Fp_mul(v, d, p);

}

W[j] = (long)u;

W[j+1] = (long)v;

}

#endif

/* same over Fp[Y] / (T) [copy-paste] */

static void

BuildTree(GEN link, GEN V, GEN W, GEN a, GEN T, GEN p)

{

long k = lg(a)-1;

Line 340 BuildTree(GEN link, GEN V, GEN W, GEN a, GEN T, GEN p)

static void

HenselLift(GEN V, GEN W, long j, GEN f, GEN T, GEN pd, GEN p0, int noinv)

{

ulong av = avma;

gpmem_t av = avma;

long space = lgef(f) * (lgefint(pd) + lgefint(p0) - 2);

GEN a2,b2,g,z,s,t;

GEN a = (GEN)V[j], b = (GEN)V[j+1];

Line 351 HenselLift(GEN V, GEN W, long j, GEN f, GEN T, GEN pd,

(void)new_chunk(space); /* HACK */

g = gadd(f, gneg_i(gmul(a,b)));

if (T) g = FpXQX_red(g, T, mulii(p0,pd));

g = gdivexact(g, p0); g = FpXQX_red(g, NULL, pd);

g = gdivexact(g, p0); if (!T) g = FpXQX_red(g, NULL, pd);

z = FpXQX_mul(v,g, T,pd);

t = FpXQX_divres(z,a, T,pd, &s);

t = gadd(gmul(u,g), gmul(t,b)); t = FpXQX_red(t, T, pd);

Line 370 HenselLift(GEN V, GEN W, long j, GEN f, GEN T, GEN pd,

g = gadd(gneg_i(g), gun);

if (T) g = FpXQX_red(g, T, mulii(p0,pd));

g = gdivexact(g, p0); g = FpXQX_red(g, NULL, pd);

g = gdivexact(g, p0); if (!T) g = FpXQX_red(g, NULL, pd);

z = FpXQX_mul(v,g, T,pd);

t = FpXQX_divres(z,a, T,pd, &s);

t = gadd(gmul(u,g), gmul(t,b)); t = FpXQX_red(t, T, pd);

Line 422 MultiLift(GEN f, GEN a, GEN T, GEN p, long e0, int fla

l = 0; E[++l] = e;

while (e > 1) { e = (e+1)/2; E[++l] = e; }

if (DEBUGLEVEL > 3) timer2();

if (DEBUGLEVEL > 3) (void)timer2();

if (flag != 2)

{

Line 528 GEN

polhensellift(GEN pol, GEN fct, GEN p, long exp)

{

GEN p1, p2;

long av = avma, i, j, l;

long i, j, l;

gpmem_t av = avma;

/* we check the arguments */

if (typ(pol) != t_POL)

err(talker, "not a polynomial in polhensellift");

if ((typ(fct) != t_COL && typ(fct) != t_VEC) || (lg(fct) < 3))

err(talker, "not a factorization in polhensellift");

if (typ(p) != t_INT || !isprime(p))

if (typ(p) != t_INT || !BSW_psp(p))

err(talker, "not a prime number in polhensellift");

if (exp < 1)

err(talker, "not a positive exponent in polhensellift");

Line 564 polhensellift(GEN pol, GEN fct, GEN p, long exp)

Line 565 polhensellift(GEN pol, GEN fct, GEN p, long exp)

{

for (i = 1; i <= l; i++)

for (j = 1; j < i; j++)

if (lgef(FpX_gcd((GEN)p1[i], (GEN)p1[j], p)) != 3)

if (degpol(FpX_gcd((GEN)p1[i], (GEN)p1[j], p)))

err(talker, "polhensellift: factors %Z and %Z are not coprime",

p1[i], p1[j]);

}

Line 579 polhensellift(GEN pol, GEN fct, GEN p, long exp)

Line 580 polhensellift(GEN pol, GEN fct, GEN p, long exp)

static GEN

two_factor_bound(GEN x)

{

long av = avma, i,j, n = lgef(x) - 3;

long i, j, n = lgef(x) - 3;

gpmem_t av = avma;

GEN *invbin, c, r = cgetr(3), z;

x += 2; invbin = (GEN*)new_chunk(n+1);

Line 604 two_factor_bound(GEN x)

Line 606 two_factor_bound(GEN x)

}

#endif

/* A | S ==> |a_i| <= binom(d-1, i-1) || S ||_2 + binom(d-1, i) lc(S) */

static GEN

uniform_Mignotte_bound(GEN x)

Mignotte_bound(GEN S)

{

long e, n = degpol(x);

long i, d = degpol(S);

GEN p1, N2 = mpsqrt(QuickNormL2(x,DEFAULTPREC));

GEN lS, C, binlS, bin, N2, p1;

N2 = mpsqrt(QuickNormL2(S,DEFAULTPREC));

binlS = bin = vecbinome(d-1);

lS = leading_term(S);

if (!is_pm1(lS)) binlS = gmul(lS, bin);

p1 = grndtoi(gmul2n(N2, n), &e);

/* i = 0 */

if (e>=0) p1 = addii(p1, shifti(gun, e));

C = (GEN)binlS[1];

return p1;

/* i = d */

p1 = N2; if (gcmp(C, p1) < 0) C = p1;

for (i = 1; i < d; i++)

{

p1 = gadd(gmul((GEN)bin[i], N2), (GEN)binlS[i+1]);

if (gcmp(C, p1) < 0) C = p1;

}

return C;

}

/* A | S ==> |a_i|^2 <= 3^{3/2 + d} / (4 \pi d) [P]_2^2,

* where [P]_2 is Bombieri's 2-norm */

static GEN

Beauzamy_bound(GEN S)

{

const long prec = DEFAULTPREC;

long i, d = degpol(S);

GEN bin, lS, s, C;

bin = vecbinome(d);

s = realzero(prec);

for (i=0; i<=d; i++)

{

GEN c = (GEN)S[i+2];

if (gcmp0(c)) continue;

/* s += P_i^2 / binomial(d,i) */

s = addrr(s, gdiv(itor(sqri(c), prec), (GEN)bin[i+1]));

}

/* s = [S]_2^2 */

C = gpow(stoi(3), dbltor(1.5 + d), prec); /* 3^{3/2 + d} */

C = gdiv(gmul(C, s), gmulsg(4*d, mppi(prec)));

lS = absi(leading_term(S));

return mulir(lS, mpsqrt(C));

}

static GEN

factor_bound(GEN S)

{

gpmem_t av = avma;

GEN a = Mignotte_bound(S);

GEN b = Beauzamy_bound(S);

if (DEBUGLEVEL>2)

{

fprintferr("Mignotte bound: %Z\n",a);

fprintferr("Beauzamy bound: %Z\n",b);

}

return gerepileupto(av, ceil_safe(gmin(a, b)));

}

#if 0

/* all factors have coeffs less than the bound */

static GEN

all_factor_bound(GEN x)

{

long i, n = lgef(x) - 3;

GEN t, z = gzero;

for (i=2; i<=n+2; i++) z = addii(z,sqri((GEN)x[i]));

for (i=2; i<=n+2; i++) z = addii(z, sqri((GEN)x[i]));

t = absi((GEN)x[n+2]);

z = addii(t, addsi(1, racine(z)));

z = mulii(z, binome(stoi(n-1), n>>1));

return shifti(mulii(t,z),1);

}

#endif

/* x defined mod p^a, return mods ( (x - mods(x, p^b)) / p^b , p^(b-a) ) */

static GEN

TruncTrace(GEN x, GEN pb, GEN pa_b, GEN pa_bs2, GEN pbs2)

{

GEN r, q = dvmdii(x, pb, &r);

if (cmpii(r, pbs2) > 0) q = addis(q,1);

if (pa_bs2 && cmpii(q,pa_bs2) > 0) q = subii(q,pa_b);

return q;

}

/* Naive recombination of modular factors: combine up to maxK modular

* factors, degree <= klim and divisible by hint

* target = polynomial we want to factor

* famod = array of modular factors. Product should be congruent to

* target/lc(target) modulo p^a

* All rational factors are bounded by p^b, b <= a and p^(b-a) < 2^(BIL-1) */

* For true factors: S1,S2 <= p^b, with b <= a and p^(b-a) < 2^31 */

static GEN

cmbf(GEN target, GEN famod, GEN p, long b, long a,

cmbf(GEN pol, GEN famod, GEN bound, GEN p, long a, long b,

long maxK, long klim,long hint)

{

long K = 1, cnt = 1, i,j,k, curdeg, lfamod = lg(famod)-1;

long spa_b, spa_bs2;

long spa_b, spa_bs2, Sbound;

GEN lc, lctarget, pa = gpowgs(p,a), pas2 = shifti(pa,-1);

GEN lc, lcpol, pa = gpowgs(p,a), pas2 = shifti(pa,-1);

GEN trace = cgetg(lfamod+1, t_VECSMALL);

GEN trace1 = cgetg(lfamod+1, t_VECSMALL);

GEN trace2 = cgetg(lfamod+1, t_VECSMALL);

GEN ind = cgetg(lfamod+1, t_VECSMALL);

GEN degpol = cgetg(lfamod+1, t_VECSMALL);

GEN degsofar = cgetg(lfamod+1, t_VECSMALL);

GEN listmod = cgetg(lfamod+1, t_COL);

GEN fa = cgetg(lfamod+1, t_COL);

GEN res = cgetg(3, t_VEC);

GEN bound = all_factor_bound(target);

if (maxK < 0) maxK = lfamod-1;

lc = absi(leading_term(target));

lc = absi(leading_term(pol));

lctarget = gmul(lc,target);

if (is_pm1(lc)) lc = NULL;

lcpol = lc? gmul(lc,pol): pol;

{

GEN pa_b,pa_bs2,pb,pbs2;

GEN pa_b,pa_bs2,pb, lc2 = lc? sqri(lc): NULL;

pa_b = gpowgs(p, a-b); /* make sure p^(a-b) < 2^(BIL-1) */

while (is_bigint(pa_b)) { b++; pa_b = divii(pa_b, p); }

pa_b = gpowgs(p, a-b); /* < 2^31 */

pa_bs2 = shifti(pa_b,-1);

pb= gpowgs(p, b); pbs2 = shifti(pb,-1);

pb= gpowgs(p, b);

for (i=1; i <= lfamod; i++)

{

GEN T, p1 = (GEN)famod[i];

GEN T1,T2, P = (GEN)famod[i];

degpol[i] = degpol(p1);

long d = degpol(P);

if (!gcmp1(lc))

T = modii(mulii(lc, (GEN)p1[degpol[i]+1]), pa);

degpol[i] = d; P += 2;

else

T1 = (GEN)P[d-1];/* = - S_1 */

T = (GEN)p1[degpol[i]+1]; /* d-1 term */

T2 = sqri(T1);

trace[i] = itos( TruncTrace(T, pb,pa_b,NULL,pbs2) );

if (d > 1) T2 = subii(T2, shifti((GEN)P[d-2],1));

T2 = modii(T2, pa); /* = S_2 Newton sum */

if (lc)

{

T1 = modii(mulii(lc, T1), pa);

T2 = modii(mulii(lc2,T2), pa);

}

trace1[i] = itos(diviiround(T1, pb));

trace2[i] = itos(diviiround(T2, pb));

}

spa_b = pa_b[2]; /* < 2^(BIL-1) */

spa_b = pa_b[2]; /* < 2^31 */

spa_bs2 = pa_bs2[2]; /* < 2^(BIL-1) */

spa_bs2 = pa_bs2[2]; /* < 2^31 */

}

degsofar[0] = 0; /* sentinel */

Line 691 cmbf(GEN target, GEN famod, GEN p, long b, long a,

Line 745 cmbf(GEN target, GEN famod, GEN p, long b, long a,

* 1 <= ind[i] <= lfamod */

nextK:

if (K > maxK || 2*K > lfamod) goto END;

if (DEBUGLEVEL > 4)

if (DEBUGLEVEL > 3)

fprintferr("\n### K = %d, %Z combinations\n", K,binome(stoi(lfamod), K));

setlg(ind, K+1); ind[1] = 1;

Sbound = ((K+1)>>1);

i = 1; curdeg = degpol[ind[1]];

for(;;)

{ /* try all combinations of K factors */

Line 704 nextK:

Line 759 nextK:

}

if (curdeg <= klim && curdeg % hint == 0) /* trial divide */

{

GEN y, q, cont, list;

GEN y, q, list;

ulong av;

gpmem_t av;

long t;

/* d - 1 test, overflow is not a problem (correct mod 2^BIL) */

/* d - 1 test */

for (t=trace[ind[1]],i=2; i<=K; i++)

for (t=trace1[ind[1]],i=2; i<=K; i++)

t = addssmod(t, trace[ind[i]], spa_b);

t = addssmod(t, trace1[ind[i]], spa_b);

if (t > spa_bs2) t -= spa_b;

if (labs(t) > ((K+1)>>1))

if (labs(t) > Sbound)

{

if (DEBUGLEVEL>6) fprintferr(".");

goto NEXT;

}

av = avma;

/* d - 2 test */

for (t=trace2[ind[1]],i=2; i<=K; i++)

t = addssmod(t, trace2[ind[i]], spa_b);

if (t > spa_bs2) t -= spa_b;

if (labs(t) > Sbound)

{

if (DEBUGLEVEL>6) fprintferr("|");

goto NEXT;

}

av = avma;

/* check trailing coeff */

y = lc;

for (i=1; i<=K; i++)

y = centermod_i(mulii(y, gmael(famod,ind[i],2)), pa, pas2);

if (!signe(y) || resii((GEN)lctarget[2], y) != gzero)

{

if (DEBUGLEVEL>6) fprintferr("T");

GEN q = constant_term((GEN)famod[ind[i]]);

if (y) q = mulii(y, q);

y = centermod_i(q, pa, pas2);

}

if (!signe(y) || resii((GEN)lcpol[2], y) != gzero)

{

if (DEBUGLEVEL>3) fprintferr("T");

avma = av; goto NEXT;

}

y = lc; /* full computation */

for (i=1; i<=K; i++)

y = centermod_i(gmul(y, (GEN)famod[ind[i]]), pa, pas2);

{

GEN q = (GEN)famod[ind[i]];

if (y) q = gmul(y, q);

y = centermod_i(q, pa, pas2);

}

/* y is the candidate factor */

if (! (q = polidivis(lctarget,y,bound)) )

if (! (q = polidivis(lcpol,y,bound)) )

{

if (DEBUGLEVEL>6) fprintferr("*");

if (DEBUGLEVEL>3) fprintferr("*");

avma = av; goto NEXT;

}

/* found a factor */

cont = content(y);

if (signe(leading_term(y)) < 0) cont = negi(cont);

y = gdiv(y, cont);

list = cgetg(K+1, t_VEC);

listmod[cnt] = (long)list;

for (i=1; i<=K; i++) list[i] = famod[ind[i]];

y = primpart(y);

fa[cnt++] = (long)y;

/* fix up target */

/* fix up pol */

target = gdiv(q, leading_term(y));

pol = q;

if (lc) pol = gdivexact(pol, leading_term(y));

for (i=j=k=1; i <= lfamod; i++)

{ /* remove used factors */

if (j <= K && i == ind[j]) j++;

else

{

famod[k] = famod[i];

trace[k] = trace[i];

trace1[k] = trace1[i];

trace2[k] = trace2[i];

degpol[k] = degpol[i]; k++;

}

lfamod -= K;

if (lfamod < 2*K) goto END;

i = 1; curdeg = degpol[ind[1]];

bound = all_factor_bound(target);

bound = factor_bound(pol);

lc = absi(leading_term(target));

if (lc) lc = absi(leading_term(pol));

lctarget = gmul(lc,target);

lcpol = lc? gmul(lc,pol): pol;

if (DEBUGLEVEL > 2)

{

fprintferr("\nfound factor %Z\nremaining modular factor(s): %ld\n",

fprintferr("\n"); msgtimer("to find factor %Z",y);

y, lfamod);

fprintferr("remaining modular factor(s): %ld\n", lfamod);

}

continue;

}

Line 786 NEXT:

Line 855 NEXT:

}

END:

if (degpol(target) > 0)

if (degpol(pol) > 0)

{ /* leftover factor */

if (signe(leading_term(target)) < 0) target = gneg_i(target);

if (signe(leading_term(pol)) < 0) pol = gneg_i(pol);

setlg(famod, lfamod+1);

listmod[cnt] = (long)dummycopy(famod);

fa[cnt++] = (long)target;

fa[cnt++] = (long)pol;

}

if (DEBUGLEVEL>6) fprintferr("\n");

setlg(listmod, cnt); setlg(fa, cnt);

Line 822 factor_quad(GEN x, GEN res, long *ptcnt)

Line 891 factor_quad(GEN x, GEN res, long *ptcnt)

long

logint(GEN B, GEN y, GEN *ptq)

{

ulong av = avma;

gpmem_t av = avma;

long e,i,fl;

GEN q,pow2, r = y;

Line 858 END:

Line 927 END:

if (ptq) *ptq = gerepileuptoint(av, icopy(r)); else avma = av;

return e;

}

/* recombination of modular factors: van Hoeij's algorithm */

/* compute Newton sums of P (i-th powers of roots, i=1..n)

/* return integer y such that all |a| <= y if P(a) = 0 */

* If N != NULL, assume p-adic roots and compute mod N [assume integer coeffs]

* If y0!= NULL, precomputed i-th powers, i=1..m, m = length(y0) */

GEN

polsym_gen(GEN P, GEN y0, long n, GEN N)

{

long av1,av2,dP=degpol(P),i,k,m;

GEN s,y,P_lead;

if (n<0) err(impl,"polsym of a negative n");

if (typ(P) != t_POL) err(typeer,"polsym");

if (!signe(P)) err(zeropoler,"polsym");

y = cgetg(n+2,t_COL);

if (y0)

{

if (typ(y0) != t_COL) err(typeer,"polsym_gen");

m = lg(y0)-1;

for (i=1; i<=m; i++) y[i] = lcopy((GEN)y0[i]);

}

else

{

m = 1;

y[1] = lstoi(dP);

}

P += 2; /* strip codewords */

P_lead=(GEN)P[dP]; if (gcmp1(P_lead)) P_lead=NULL;

if (N && P_lead)

if (!invmod(P_lead,N,&P_lead))

err(talker,"polsyn: non-invertible leading coeff: %Z", P_lead);

for (k=m; k<=n; k++)

{

av1=avma; s = (dP>=k)? gmulsg(k,(GEN)P[dP-k]): gzero;

for (i=1; i<k && i<=dP; i++)

s = gadd(s, gmul((GEN)y[k-i+1],(GEN)P[dP-i]));

if (N)

{

s = modii(s, N);

if (P_lead) { s = mulii(s,P_lead); s = modii(s,N); }

}

else

if (P_lead) s = gdiv(s,P_lead);

av2=avma; y[k+1]=lpile(av1,av2,gneg(s));

}

return y;

}

/* return integer y such that all roots of P are less than y */

static GEN

root_bound(GEN P)

root_bound(GEN P0)

{

GEN P0 = dummycopy(P), lP = absi(leading_term(P)), x,y,z;

GEN Q = dummycopy(P0), lP = absi(leading_term(Q)), P,x,y,z;

long k,d = degpol(P);

long k, d = degpol(Q);

setlgef(P0, d+2); /* P = lP x^d + P0 */

setlgef(Q, d+2); /* P = lP x^d + Q */

P = P0+2; /* strip codewords */

P = Q+2; /* strip codewords */

for (k=0; k<d; k++) P[k] = labsi((GEN)P[k]);

x = y = gun;

k = gexpo(cauchy_bound(P0)) - 5;

if (k < 0) k = 0;

x = y = shifti(gun, k);

for (k=0; ; k++)

{

if (cmpii(poleval(P0,y), mulii(lP, gpowgs(y, d))) < 0) break;

gpmem_t av = avma;

if (cmpii(poleval(Q,y), mulii(lP, gpowgs(y, d))) < 0) break;

avma = av;

x = y; y = shifti(y,1);

}

for(;;)

for(k=0; ; k++)

{

z = shifti(addii(x,y), -1);

if (egalii(x,z)) break;

if (egalii(x,z) || k == 20) break;

if (cmpii(poleval(P0,z), mulii(lP, gpowgs(z, d))) < 0)

if (cmpii(poleval(Q,z), mulii(lP, gpowgs(z, d))) < 0)

y = z;

else

x = z;

Line 937 root_bound(GEN P)

Line 963 root_bound(GEN P)

return y;

}

extern GEN gscal(GEN x,GEN y);

static GEN

extern GEN sqscal(GEN x);

ZM_HNFimage(GEN x)

/* return Gram-Schmidt orthogonal basis (f) associated to (e), B is the

* vector of the (f_i . f_i)*/

GEN

gram_schmidt(GEN e, GEN *ptB)

{

long i,j,lx = lg(e);

return (lg(x) > 50)? hnflll_i(x, NULL, 1): hnfall_i(x, NULL, 1);

GEN f = dummycopy(e), B, iB;

B = cgetg(lx, t_VEC);

iB= cgetg(lx, t_VEC);

for (i=1;i<lx;i++)

{

GEN p1 = NULL;

ulong av;

B[i] = (long)sqscal((GEN)f[i]);

iB[i]= linv((GEN)B[i]); av = avma;

for (j=1; j<i; j++)

{

GEN mu = gmul(gscal((GEN)e[i],(GEN)f[j]), (GEN)iB[j]);

GEN p2 = gmul(mu, (GEN)f[j]);

p1 = p1? gadd(p1,p2): p2;

}

p1 = p1? gerepileupto(av, gsub((GEN)e[i], p1)): (GEN)e[i];

f[i] = (long)p1;

}

*ptB = B; return f;

}

extern GEN hnfall_i(GEN A, GEN *ptB, long remove);

GEN

special_pivot(GEN x)

{

GEN t, H = hnfall_i(x, NULL, 1);

GEN t, H = ZM_HNFimage(x);

long i,j, l = lg(H), h = lg(H[1]);

for (i=1; i<h; i++)

{

Line 992 special_pivot(GEN x)

Line 990 special_pivot(GEN x)

return H;

}

/* x matrix: compute a bound for | \sum e_i x_i | ^ 2, e_i = 0,1 */

/* B from lllint_i: return [ |b_i^*|^2, i ] */

GEN

GS_norms(GEN B, long prec)

{

long i, l = lg(B);

GEN v = gmul(B, realun(prec));

l--; setlg(v, l);

for (i=1; i<l; i++)

v[i] = ldivrr((GEN)v[i+1], (GEN)v[i]);

return v;

}

static GEN

my_norml2(GEN x)

check_factors(GEN P, GEN M_L, GEN bound, GEN famod, GEN pa)

{

long i,j, co = lg(x), li;

long i, j, r, n0;

GEN S = gzero;

GEN pol = P, lcpol, lc, list, piv, y, pas2;

if (co == 1) return S;

li = lg(x[1]);

piv = special_pivot(M_L);

for (i=1; i<li; i++)

if (!piv) return NULL;

if (DEBUGLEVEL>3) fprintferr("special_pivot output:\n%Z\n",piv);

pas2 = shifti(pa,-1);

r = lg(piv)-1;

n0 = lg(piv[1])-1;

list = cgetg(r+1, t_COL);

lc = absi(leading_term(pol));

if (is_pm1(lc)) lc = NULL;

lcpol = lc? gmul(lc, pol): pol;

for (i=1; i<r; i++)

{

GEN p = gzero;

GEN c = (GEN)piv[i];

GEN m = gzero;

if (DEBUGLEVEL) fprintferr("LLL_cmbf: checking factor %ld\n",i);

for (j=1; j<co; j++)

y = lc;

for (j=1; j<=n0; j++)

if (signe(c[j]))

{

GEN q = (GEN)famod[j];

if (y) q = gmul(y, q);

y = centermod_i(q, pa, pas2);

}

/* y is the candidate factor */

if (! (pol = polidivis(lcpol,y,bound)) ) return NULL;

y = primpart(y);

if (lc)

{

GEN u = gcoeff(x,i,j);

pol = gdivexact(pol, leading_term(y));

long s = gsigne(u);

lc = absi(leading_term(pol));

if (s < 0) m = gadd(m,u);

else if (s > 0) p = gadd(p,u);

}

if (m != gzero) m = gneg(m);

lcpol = lc? gmul(lc, pol): pol;

if (gcmp(p,m) > 0) m = p;

list[i] = (long)y;

}

y = primpart(pol);

list[i] = (long)y; return list;

}

S = gadd(S, gsqr(m));

GEN

LLL_check_progress(GEN Bnorm, long n0, GEN m, int final,

pari_timer *T, long *ti_LLL)

{

GEN B, norm, u;

long i, R;

if (DEBUGLEVEL>2) (void)TIMER(T);

u = lllint_i(m, final? 1000: 4, 0, NULL, NULL, &B);

if (DEBUGLEVEL>2) *ti_LLL += TIMER(T);

norm = GS_norms(B, DEFAULTPREC);

for (R=lg(m)-1; R > 0; R--)

if (cmprr((GEN)norm[R], Bnorm) < 0) break;

for (i=1; i<=R; i++) setlg(u[i], n0+1);

if (R <= 1)

{

if (!R) err(bugparier,"LLL_cmbf [no factor]");

return NULL; /* irreducible */

}

return S;

setlg(u, R+1); return u;

}

/* each entry < 2^e */

static ulong

static long

next2pow(ulong a)

init_padic_prec(long e, int BitPerFactor, long n0, double LOGp2)

{

/* long b, goodb = (long) ((e - 0.175 * n0 * n0) * LOGp2); */

ulong b = 1;

long b, goodb = (long) (((e - BitPerFactor*n0)) * LOGp2);

while (b < a) b <<= 1;

b = (long) ((e - 32)*LOGp2); if (b < goodb) goodb = b;

return b;

/* b = (long) ((e - Max*32)*LOGp2); if (b > goodb) goodb = b; */

return goodb;

}

extern GEN sindexrank(GEN x);

extern GEN vconcat(GEN Q1, GEN Q2);

extern GEN gauss_intern(GEN a, GEN b);

/* Recombination phase of Berlekamp-Zassenhaus algorithm using a variant of

* van Hoeij's knapsack

* P = monic squarefree in Z[X].

* P = squarefree in Z[X].

* famod = array of (lifted) modular factors mod p^a

* bound = Mignotte bound for the size of divisors of P (sor the sup norm)

* bound = Mignotte bound for the size of divisors of P (for the sup norm)

* previously recombined all set of factors with less than rec elts

* previously recombined all set of factors with less than rec elts */

static GEN

GEN

LLL_cmbf(GEN P, GEN famod, GEN p, GEN pa, GEN bound, long a, long rec)

{

long s = 2; /* # of traces added at each step */

const long N0 = 1; /* # of traces added at each step */

long BitPerFactor = 3; /* nb bits in p^(a-b) / modular factor */

double BitPerFactor = 0.5; /* nb bits in p^(a-b) / modular factor */

long i,j,C,r,S,tmax,n0,n,dP = degpol(P);

long i,j,tmax,n0,C, dP = degpol(P);

double logp = log(gtodouble(p)), LOGp2 = LOG2/logp;

double logp = log((double)itos(p)), LOGp2 = LOG2/logp;

double b0 = log((double)dP*2) / logp;

double b0 = log((double)dP*2) / logp, logBr;

double k = gtodouble(glog(root_bound(P), DEFAULTPREC)) / logp;

GEN lP, Br, Bnorm, Tra, T2, TT, CM_L, m, list, ZERO;

GEN y, T, T2, TT, BL, m, u, norm, target, M, piv, list;

gpmem_t av, av2, lim;

GEN run = realun(DEFAULTPREC);

long ti_LLL = 0, ti_CF = 0;

ulong av,av2, lim;

pari_timer ti, ti2, TI;

int did_recompute_famod = 0;

n0 = n = r = lg(famod) - 1;

if (DEBUGLEVEL>2) (void)TIMER(&ti);

list = cgetg(n0+1, t_COL);

lP = absi(leading_term(P));

if (is_pm1(lP)) lP = NULL;

Br = root_bound(P);

if (lP) Br = gmul(lP, Br);

logBr = gtodouble(glog(Br, DEFAULTPREC)) / logp;

n0 = lg(famod) - 1;

C = (long)ceil( sqrt(N0 * n0 / 4.) ); /* > 1 */

Bnorm = dbltor(n0 * (C*C + N0*n0/4.) * 1.00001);

ZERO = zeromat(n0, N0);

av = avma; lim = stack_lim(av, 1);

TT = cgetg(n0+1, t_VEC);

T = cgetg(n0+1, t_MAT);

Tra = cgetg(n0+1, t_MAT);

for (i=1; i<=n0; i++)

{

TT[i] = 0;

T [i] = lgetg(s+1, t_COL);

Tra[i] = lgetg(N0+1, t_COL);

}

BL = idmat(n0);

CM_L = gscalsmat(C, n0);

/* tmax = current number of traces used (and computed so far)

/* tmax = current number of traces used (and computed so far) */

* S = number of traces used at the round's end = tmax + s */

for (tmax = 0;; tmax += N0)

for(tmax = 0;; tmax = S)

{

GEN pas2, pa_b, BE;

long b, bmin, bgood, delta, tnew = tmax + N0, r = lg(CM_L)-1;

long b, goodb;

GEN M_L, q, CM_Lp, oldCM_L;

double Nx;

int first = 1;

bmin = (long)ceil(b0 + tnew*logBr);

if (DEBUGLEVEL>2)

fprintferr("\nLLL_cmbf: %ld potential factors (tmax = %ld, bmin = %ld)\n",

r, tmax, bmin);

if (DEBUGLEVEL>3)

/* compute Newton sums (possibly relifting first) */

fprintferr("LLL_cmbf: %ld potential factors (tmax = %ld)\n", r, tmax);

if (a <= bmin)

av2 = avma;

if (tmax > 0)

{ /* bound small vector in terms of a modified L2 norm of a

* left inverse of BL */

GEN z = gauss_intern(BL,NULL); /* 1/BL */

if (!z) /* not maximal rank */

{

avma = av2;

BL = hnfall_i(BL,NULL,1);

r = lg(BL)-1; z = invmat(BL);

av2 = avma;

}

Nx = gtodouble(my_norml2(gmul(run, z)));

avma = av2;

}

else

Nx = (double)n0; /* first time: BL = id */

C = (long)sqrt(s*n0*n0/4. / Nx);

if (C == 0) C = 1; /* frequent after a few iterations */

M = dbltor((Nx * C*C + s*n0*n0/4.) * 1.00001);

S = tmax + s;

b = (long)ceil(b0 + S*k);

if (a <= b)

{

a = (long)ceil(b + 3*s*k) + 1; /* enough for 3 more rounds */

a = (long)ceil(bmin + 3*N0*logBr) + 1; /* enough for 3 more rounds */

a = next2pow(a);

pa = gpowgs(p,a);

famod = hensel_lift_fact(P,famod,NULL,p,pa,a);

/* recompute old Newton sums to new precision */

for (i=1; i<=n0; i++) TT[i] = 0;

for (i=1; i<=n0; i++)

TT[i] = (long)polsym_gen((GEN)famod[i], NULL, tmax, pa);

did_recompute_famod = 1;

}

for (i=1; i<=n0; i++)

{

GEN p1 = (GEN)T[i];

GEN p1 = (GEN)Tra[i];

GEN p2 = polsym_gen((GEN)famod[i], (GEN)TT[i], S, pa);

GEN p2 = polsym_gen((GEN)famod[i], (GEN)TT[i], tnew, NULL, pa);

TT[i] = (long)p2;

p2 += 1+tmax; /* ignore traces number 0...tmax */

for (j=1; j<=s; j++) p1[j] = p2[j];

for (j=1; j<=N0; j++) p1[j] = p2[j];

if (lP)

{ /* make Newton sums integral */

GEN lPpow = gpowgs(lP, tmax);

for (j=1; j<=N0; j++)

{

lPpow = mulii(lPpow,lP);

p1[j] = lmulii((GEN)p1[j], lPpow);

}

/* compute truncation parameter */

if (DEBUGLEVEL>2) { TIMERstart(&ti2); TIMERstart(&TI); }

oldCM_L = CM_L;

av2 = avma;

T2 = gmod( gmul(T, BL), pa );

delta = b = 0; /* -Wall */

goodb = init_padic_prec(gexpo(T2), BitPerFactor, n0, LOGp2);

AGAIN:

if (goodb > b) b = goodb;

M_L = gdivexact(CM_L, stoi(C));

if (DEBUGLEVEL>2)

T2 = centermod( gmul(Tra, M_L), pa );

fprintferr("LLL_cmbf: (a, b) = (%ld, %ld), expo(T) = %ld\n",

if (first)

a,b,gexpo(T2));

{ /* initialize lattice, using few p-adic digits for traces */

pa_b = gpowgs(p, a-b);

double t = gexpo(T2) - max(32, BitPerFactor*r);

{

bgood = (long) (t * LOGp2);

GEN pb = gpowgs(p, b);

b = max(bmin, bgood);

GEN pa_bs2 = shifti(pa_b,-1);

delta = a - b;

GEN pbs2 = shifti(pb,-1);

for (i=1; i<=r; i++)

{

GEN p1 = (GEN)T2[i];

for (j=1; j<=s; j++)

p1[j] = (long)TruncTrace((GEN)p1[j],pb,pa_b,pa_bs2,pbs2);

}

if (gcmp0(T2)) { avma = av2; continue; }

else

{ /* add more p-adic digits and continue reduction */

long b0 = gexpo(T2) * LOGp2;

if (b0 < b) b = b0;

b = max(b-delta, bmin);

if (b - delta/2 < bmin) b = bmin; /* near there. Go all the way */

}

BE = cgetg(s+1, t_MAT);

q = gpowgs(p, b);

for (i=1; i<=s; i++)

m = vconcat( CM_L, gdivround(T2, q) );

if (first)

{

BE[i] = (long)zerocol(r+s);

GEN P1 = gscalmat(gpowgs(p, a-b), N0);

coeff(BE, i+r, i) = (long)pa_b;

first = 0;

m = concatsp( m, vconcat(ZERO, P1) );

/* [ C M_L 0 ]

* m = [ ] square matrix

* [ T2' p^(a-b) I_s ] T2' = Tra * M_L truncated

}

m = concatsp( vconcat( gscalmat(stoi(C), r), T2 ), BE );

/* [ C 0 ]

* m = [ ] square matrix

* [ T2 p^(a-b) I_S ] T2 = T * BL truncated

u = lllgramintern(gram_matrix(m), 4, 0, 0);

m = gmul(m,u);

(void)gram_schmidt(gmul(run,m), &norm);

for (i=r+s; i>0; i--)

if (cmprr((GEN)norm[i], M) < 0) break;

if (i > r)

{ /* no progress */

avma = av2; BitPerFactor += 2;

if (DEBUGLEVEL>2)

fprintferr("LLL_cmbf: increasing BitPerFactor = %ld\n", BitPerFactor);

#if 0

s++; for (i=1; i<=n; i++) T[i] = lgetg(s+1, t_COL);

if (DEBUGLEVEL>2) fprintferr("LLL_cmbf: increasing s = %ld\n", s);

#endif

continue;

}

n = r; r = i;

CM_L = LLL_check_progress(Bnorm, n0, m, b == bmin, /*dbg:*/ &ti, &ti_LLL);

if (r <= 1)

if (DEBUGLEVEL>2)

fprintferr("LLL_cmbf: (a,b) =%4ld,%4ld; r =%3ld -->%3ld, time = %ld\n",

a,b, lg(m)-1, CM_L? lg(CM_L)-1: 1, TIMER(&TI));

if (!CM_L) { list = _col(P); break; }

if (b > bmin)

{

if (r == 0) err(bugparier,"LLL_cmbf [no factor]");

CM_L = gerepilecopy(av2, CM_L);

if (DEBUGLEVEL>2) fprintferr("LLL_cmbf: 1 factor\n");

goto AGAIN;

list[1] = (long)P; setlg(list,2); return list;

}

if (DEBUGLEVEL>2) msgTIMER(&ti2, "for this block of traces");

setlg(u, r+1);

i = lg(CM_L) - 1;

for (i=1; i<=r; i++) setlg(u[i], n+1);

if (i == r && gegal(CM_L, oldCM_L))

BL = gerepileupto(av2, gmul(BL, u));

if (low_stack(lim, stack_lim(av,1)))

{

GEN *gptr[4]; gptr[0]=&BL; gptr[1]=&TT; gptr[2]=&T; gptr[3]=&famod;

CM_L = oldCM_L;

if(DEBUGMEM>1) err(warnmem,"LLL_cmbf");

avma = av2; continue;

gerepilemany(av, gptr, did_recompute_famod? 4: 3);

}

if (r*rec >= n0) continue;

av2 = avma;

CM_Lp = FpM_image(CM_L, stoi(27449)); /* inexpensive test */

piv = special_pivot(BL);

if (lg(CM_Lp) != lg(CM_L))

if (DEBUGLEVEL) fprintferr("special_pivot output:\n%Z\n",piv);

if (!piv) { avma = av2; continue; }

pas2 = shifti(pa,-1); target = P;

for (i=1; i<=r; i++)

{

GEN p1 = (GEN)piv[i];

if (DEBUGLEVEL>2) fprintferr("LLL_cmbf: rank decrease\n");

if (DEBUGLEVEL) fprintferr("LLL_cmbf: checking factor %ld\n",i);

CM_L = ZM_HNFimage(CM_L);

}

y = gun;

if (i <= r && i*rec < n0)

for (j=1; j<=n0; j++)

{

if (signe(p1[j]))

if (DEBUGLEVEL>2) (void)TIMER(&ti);

y = centermod_i(gmul(y, (GEN)famod[j]), pa, pas2);

list = check_factors(P, M_L, bound, famod, pa);

if (DEBUGLEVEL>2) ti_CF += TIMER(&ti);

/* y is the candidate factor */

if (list) break;

if (! (target = polidivis(target,y,bound)) ) break;

CM_L = gerepilecopy(av2, CM_L);

list[i] = (long)y;

}

if (i > r)

if (low_stack(lim, stack_lim(av,1)))

{

if (DEBUGLEVEL>2) fprintferr("LLL_cmbf: %ld factors\n", r);

if(DEBUGMEM>1) err(warnmem,"LLL_cmbf");

setlg(list,r+1); return list;

gerepileall(av, 5, &CM_L, &TT, &Tra, &famod, &pa);

}

if (DEBUGLEVEL>2)

fprintferr("* Time LLL: %ld\n* Time Check Factor: %ld\n",ti_LLL,ti_CF);

return list;

}

extern GEN primitive_pol_to_monic(GEN pol, GEN *ptlead);

/* Return P(h * x) */

GEN

/* P(hx), in place. Assume P in Z[X], h in Z */

unscale_pol(GEN P, GEN h)

void

rescale_pol_i(GEN P, GEN h)

{

GEN hi = gun;

long i, l = lgef(P);

GEN hi = gun, Q = cgetg(l, t_POL);

Q[1] = P[1];

Q[2] = lcopy((GEN)P[2]);

for (i=3; i<l; i++)

{

hi = mulii(hi,h);

hi = gmul(hi,h);

P[i] = lmulii((GEN)P[i], hi);

Q[i] = lmul((GEN)P[i], hi);

}

return Q;

}

/* P(hx) in Fp[X], in place. Assume P in Z[X], h in Z */

/* Return h^degpol(P) P(x / h) */

void

GEN

FpX_rescale_i(GEN P, GEN h, GEN p)

rescale_pol(GEN P, GEN h)

{

GEN hi = gun;

long i, l = lgef(P);

for (i=3; i<l; i++)

GEN Q = cgetg(l,t_POL), hi = gun;

Q[l-1] = P[l-1];

for (i=l-2; i>=2; i--)

{

hi = gmul(hi,h);

Q[i] = lmul((GEN)P[i], hi);

}

Q[1] = P[1]; return Q;

}

/* as above over Fp[X] */

GEN

FpX_rescale(GEN P, GEN h, GEN p)

{

long i, l = lgef(P);

GEN Q = cgetg(l,t_POL), hi = gun;

Q[l-1] = P[l-1];

for (i=l-2; i>=2; i--)

{

hi = modii(mulii(hi,h), p);

P[i] = lmodii(mulii((GEN)P[i], hi), p);

Q[i] = lmodii(mulii((GEN)P[i], hi), p);

}

Q[1] = P[1]; return Q;

}

/* Find a,b minimal such that A < q^a, B < q^b, 1 << q^(a-b) < 2^31 */

int

cmbf_precs(GEN q, GEN A, GEN B, long *pta, long *ptb, GEN *qa, GEN *qb)

{

long a,b,amin,d = (long)(31 * LOG2/gtodouble(glog(q,DEFAULTPREC)) - 1e-5);

int fl = 0;

b = logint(B, q, qb);

amin = b + d;

if (gcmp(gpowgs(q, amin), A) <= 0)

{

a = logint(A, q, qa);

b = a - d; *qb = gpowgs(q, b);

}

else

{ /* not enough room */

a = amin; *qa = gpowgs(q, a);

fl = 1;

}

if (DEBUGLEVEL > 3) {

fprintferr("S_2 bound: %Z^%ld\n", q,b);

fprintferr("coeff bound: %Z^%ld\n", q,a);

}

*pta = a;

*ptb = b; return fl;

}

/* use van Hoeij's knapsack algorithm */

static GEN

combine_factors(GEN a, GEN famod, GEN p, long klim, long hint)

combine_factors(GEN target, GEN famod, GEN p, long klim, long hint)

{

GEN B = uniform_Mignotte_bound(a), res,y,lt,L,pe,pE,listmod,p1;

GEN la, B, A, res, L, pa, pb, listmod;

long i,E,e,l, maxK = 3, nft = lg(famod)-1;

long a,b, l, maxK = 3, nft = lg(famod)-1, n = degpol(target);

pari_timer T;

e = logint(B, p, &pe);

A = factor_bound(target);

if (DEBUGLEVEL > 4) fprintferr("Mignotte bound: %Z^%ld\n", p,e);

{ /* overlift for the d-1 test */

la = absi(leading_term(target));

int over = (int) (32 * LOG2 / log((double)itos(p)));

B = mulsi(n, sqri(gmul(la, root_bound(target)))); /* = bound for S_2 */

pE = mulii(pe, gpowgs(p,over)); E = e + over;

}

famod = hensel_lift_fact(a,famod,NULL,p,pE,E);

(void)cmbf_precs(p, A, B, &a, &b, &pa, &pb);

if (DEBUGLEVEL>2) (void)TIMER(&T);

famod = hensel_lift_fact(target,famod,NULL,p,pa,a);

if (nft < 11) maxK = -1; /* few modular factors: try all posibilities */

else

if (DEBUGLEVEL>2) msgTIMER(&T, "Hensel lift (mod %Z^%ld)", p,a);

{

L = cmbf(target, famod, A, p, a, b, maxK, klim, hint);

lt = leading_term(a);

if (DEBUGLEVEL>2) msgTIMER(&T, "Naive recombination");

if (!is_pm1(lt) && nft < 13) maxK = -1;

}

L = cmbf(a, famod, p, e, E, maxK, klim, hint);

res = (GEN)L[1];

listmod = (GEN)L[2]; l = lg(listmod);

listmod = (GEN)L[2]; l = lg(listmod)-1;

famod = (GEN)listmod[l-1];

famod = (GEN)listmod[l];

if (maxK >= 0 && lg(famod)-1 > 2*maxK)

{

a = (GEN)res[l-1];

if (l!=1) A = factor_bound((GEN)res[l]);

lt = leading_term(a);

if (signe(lt) < 0) { a = gneg_i(a); lt = leading_term(a); }

if (DEBUGLEVEL > 4) fprintferr("last factor still to be checked\n");

if (gcmp1(lt))

L = LLL_cmbf((GEN)res[l], famod, p, pa, A, a, maxK);

lt = NULL;

if (DEBUGLEVEL>2) msgTIMER(&T,"Knapsack");

else

/* remove last elt, possibly unfactored. Add all new ones. */

{

setlg(res, l); res = concatsp(res, L);

GEN invlt, invLT;

}

if (DEBUGLEVEL > 4) fprintferr("making it monic\n");

return res;

a = primitive_pol_to_monic(a, &lt);

}

B = uniform_Mignotte_bound(a);

e = logint(B, p, &pe);

invlt = mpinvmod(lt,p);

#define u_FpX_div(x,y,p) u_FpX_divrem((x),(y),(p),NULL)

for (i = 1; i<lg(famod); i++)

{ /* renormalize modular factors */

extern GEN u_FpX_Kmul(GEN a, GEN b, ulong p, long na, long nb);

p1 = (GEN)famod[i]; FpX_rescale_i(p1, invlt, p);

extern GEN u_pol_to_vec(GEN x, long N);

invLT = powmodulo(lt, stoi(degpol(p1)), p);

extern GEN u_FpM_FpX_mul(GEN x, GEN y, ulong p);

famod[i] = (long)FpX_Fp_mul(p1, invLT, p);

#define u_FpX_mul(x,y, p) u_FpX_Kmul(x+2,y+2,p, lgef(x)-2,lgef(y)-2)

}

famod = hensel_lift_fact(a,famod,NULL,p,pe,e);

static GEN

}

u_FpM_Frobenius(GEN u, ulong p)

setlg(res, l-1); /* remove last elt (possibly unfactored) */

{

L = LLL_cmbf(a, famod, p, pe, B, e, maxK);

long j,N = degpol(u);

if (lt)

GEN v,w,Q;

if (DEBUGLEVEL > 7) (void)timer2();

Q = cgetg(N+1,t_MAT);

Q[1] = (long)vecsmall_const(N, 0);

coeff(Q,1,1) = 1;

w = v = u_FpXQ_pow(u_Fp_FpX(polx[varn(u)], p), utoi(p), u, p);

for (j=2; j<=N; j++)

{

Q[j] = (long)u_pol_to_vec(w, N);

if (j < N)

{

for (i=1; i<lg(L); i++)

gpmem_t av = avma;

{

w = gerepileupto(av, u_FpX_rem(u_FpX_mul(w, v, p), u, p));

y = (GEN)L[i]; rescale_pol_i(y, lt);

L[i] = (long)primitive_part(y, &p1);

}

res = concatsp(res, L);

}

return res;

if (DEBUGLEVEL > 7) msgtimer("frobenius");

return Q;

}

extern long split_berlekamp(GEN *t, GEN pp, GEN pps2);

/* Assume a squarefree, degree(a) > 0, a(0) != 0 */

#define u_FpX_div(x,y,p) u_FpX_divrem((x),(y),(p),(0),NULL)

/* assume degree(a) > 0, a(0) != 0, and a squarefree */

static GEN

squff(GEN a, long hint)

DDF(GEN a, long hint)

{

GEN PolX,lead,res,prime,primes2,famod,y,g,z,w,*tabd,*tabdnew;

GEN MP, PolX, lead, prime, famod, z, w;

long da = degpol(a), va = varn(a);

const long da = degpol(a);

long klim,chosenp,p,nfacp,lbit,j,d,e,np,nmax,nf,nft;

long chosenp, p, nfacp, d, e, np, nmax;

ulong *tabbit, *tabkbit, *tmp, av = avma;

gpmem_t av = avma, av1;

byteptr pt=diffptr;

byteptr pt = diffptr;

const int MAXNP = max(5, (long)sqrt(da));

const int MAXNP = max(5, (long)sqrt((double)da));

long ti = 0;

pari_timer T;

if (DEBUGLEVEL>2) (void)TIMER(&T);

if (hint <= 0) hint = 1;

if (DEBUGLEVEL > 2) timer2();

if (DEBUGLEVEL > 2) (void)timer2();

lbit = (da>>4)+1; nmax = da+1; klim = da>>1;

nmax = da+1;

chosenp = 0;

tabd = NULL;

tabdnew = (GEN*)new_chunk(nmax);

tabbit = (ulong*)new_chunk(lbit);

tabkbit = (ulong*)new_chunk(lbit);

tmp = (ulong*)new_chunk(lbit);

prime = icopy(gun);

lead = (GEN)a[da+2]; PolX = u_Fp_FpX(polx[0],0, 2);

lead = (GEN)a[da+2]; if (gcmp1(lead)) lead = NULL;

for (p = np = 0; np < MAXNP; )

PolX = u_Fp_FpX(polx[0], 2);

av1 = avma;

for (p = np = 0; np < MAXNP; avma = av1)

{

ulong av0 = avma;

NEXT_PRIME_VIADIFF_CHECK(p,pt);

p += *pt++; if (!*pt) err(primer1);

if (lead && !smodis(lead,p)) continue;

if (!smodis(lead,p)) continue;

z = u_Fp_FpX(a, p);

z = u_Fp_FpX(a,0, p);

if (!u_FpX_is_squarefree(z, p)) continue;

if (!u_FpX_is_squarefree(z, p)) { avma = av0; continue ; }

for (j=0; j<lbit-1; j++) tabkbit[j] = 0;

MP = u_FpM_Frobenius(z, p);

tabkbit[j] = 1;

d = 0; e = da; nfacp = 0;

w = PolX; prime[2] = p;

while (d < (e>>1))

{

long lgg;

d++; w = u_FpXQ_pow(w, prime, z, p);

GEN g;

/* here e = degpol(z) */

d++;

w = u_FpM_FpX_mul(MP, w, p); /* w^p mod (z,p) */

g = u_FpX_gcd(z, u_FpX_sub(w, PolX, p), p);

lgg = degpol(g);

if (lgg == 0) g = NULL;

if (!lgg) continue;

else

{

e -= lgg;

z = u_FpX_div(z, g, p); e = degpol(z);

nfacp += lgg/d;

w = u_FpX_rem(w, z, p);

if (DEBUGLEVEL>5)

lgg /= d; nfacp += lgg;

fprintferr(" %3ld fact. of degree %3ld\n", lgg/d, d);

if (DEBUGLEVEL>5)

if (!e) break;

fprintferr(" %3ld factor%s of degree %3ld\n", lgg, lgg==1?"":"s",d);

z = u_FpX_div(z, g, p);

record_factors(lgg, d, lbit-1, tabkbit, tmp);

w = u_FpX_rem(w, z, p);

}

tabdnew[d] = g;

}

if (e > 0)

if (e)

{

if (DEBUGLEVEL>5) fprintferr(" %3ld factor of degree %3ld\n",1,e);

tabdnew[e] = z; nfacp++;

nfacp++;

record_factors(1, e, lbit-1, tabkbit, tmp);

}

if (DEBUGLEVEL>4)

if (np)

for (j=0; j<lbit; j++) tabbit[j] &= tabkbit[j];

else

for (j=0; j<lbit; j++) tabbit[j] = tabkbit[j];

if (DEBUGLEVEL > 4)

fprintferr("...tried prime %3ld (%-3ld factor%s). Time = %ld\n",

p, nfacp, nfacp==1?"": "s", timer2());

if (min_deg(lbit-1,tabbit) > klim) { avma = av; return _col(a); }

if (nfacp < nmax)

{

nmax = nfacp; tabd = tabdnew; chosenp = p;

if (nfacp == 1) { avma = av; return _col(a); }

for (j=d+1; j<e; j++) tabd[j] = NULL;

nmax = nfacp; chosenp = p;

tabdnew = (GEN*)new_chunk(da);

if (nmax < 5) break; /* very few factors. Enough */

}

else avma = av0;

np++;

}

prime[2] = chosenp; primes2 = shifti(prime, -1);

prime[2] = chosenp;

nf = nmax; nft = 1;

famod = cgetg(nmax+1,t_COL);

y = cgetg(nf+1,t_COL); famod = cgetg(nf+1,t_COL);

famod[1] = (long)(lead? FpX_normalize(a, prime): FpX_red(a, prime));

for (d = 1; nft <= nf; d++)

if (nmax != FpX_split_berlekamp((GEN*)(famod+1),prime))

err(bugparier,"DDF: wrong numbers of factors");

if (DEBUGLEVEL>2)

{

g = tabd[d];

if (DEBUGLEVEL>4) msgtimer("splitting mod p = %ld",chosenp);

if (g)

ti = TIMER(&T);

{

fprintferr("Time setup: %ld\n", ti);

long n = degpol(g)/d;

famod[nft] = (long)small_to_pol(u_FpX_normalize(g, chosenp),va);

if (n > 1 && n != split_berlekamp((GEN*)(famod+nft),prime,primes2))

err(bugparier,"squff: wrong numbers of factors");

nft += n;

}

if (DEBUGLEVEL > 4) msgtimer("splitting mod p = %ld",chosenp);

z = combine_factors(a, famod, prime, da-1, hint);

res = combine_factors(a, famod, prime, da-1, hint);

if (DEBUGLEVEL>2)

return gerepilecopy(av, res);

fprintferr("Total Time: %ld\n===========\n", ti + TIMER(&T));

return gerepilecopy(av, z);

}

/* A(X^d) --> A(X) */

Line 1442 gdeflate(GEN x, long v, long d)

Line 1498 gdeflate(GEN x, long v, long d)

long i, lx, tx = typ(x);

GEN z;

if (is_scalar_t(tx)) return gcopy(x);

if (d <= 0) err(talker,"need positive degree in gdeflate");

if (tx == t_POL)

{

long vx = varn(x);

ulong av;

gpmem_t av;

if (vx < v)

{

lx = lgef(x);

Line 1499 polinflate(GEN x0, long d)

Line 1556 polinflate(GEN x0, long d)

return y;

}

/* Distinct Degree Factorization (deflating first)

* Assume x squarefree, degree(x) > 0, x(0) != 0 */

GEN

squff2(GEN x, long hint)

DDF2(GEN x, long hint)

{

GEN L;

long m;

x = poldeflate(x, &m);

L = squff(x, hint);

L = DDF(x, hint);

if (m > 1)

{

GEN e, v, fa = decomp(stoi(m));

Line 1525 squff2(GEN x, long hint)

Line 1584 squff2(GEN x, long hint)

{

GEN L2 = cgetg(1,t_VEC);

for (i=1; i < lg(L); i++)

L2 = concatsp(L2, squff(polinflate((GEN)L[i], v[k]), hint));

L2 = concatsp(L2, DDF(polinflate((GEN)L[i], v[k]), hint));

L = L2;

}

return L;

}

/* Factor x in Z[t]. Assume all factors have degree divisible by hint */

/* SquareFree Factorization. f = prod P^e, all e distinct, in Z[X] (char 0

* would be enough, if modulargcd --> ggcd). Return (P), set *ex = (e) */

GEN

factpol(GEN x, long hint)

ZX_squff(GEN f, GEN *ex)

{

GEN fa,p1,d,t,v,w, y = cgetg(3,t_MAT);

GEN T,V,W,P,e,cf;

long av=avma,av2,lx,vx,i,j,k,ex,nbfac,zval;

long i,k,dW,n,val;

if (typ(x)!=t_POL) err(notpoler,"factpol");

val = polvaluation(f, &f);

if (!signe(x)) err(zeropoler,"factpol");

n = 1 + degpol(f); if (val) n++;

e = cgetg(n,t_VECSMALL);

P = cgetg(n,t_COL);

ex = nbfac = 0;

cf = content(f); if (gsigne(leading_term(f)) < 0) cf = gneg_i(cf);

p1 = x+2; while (gcmp0((GEN)*p1)) p1++;

if (!gcmp1(cf)) f = gdiv(f,cf);

zval = p1 - (x + 2);

lx = lgef(x) - zval;

T = modulargcd(derivpol(f), f);

vx = varn(x);

V = gdeuc(f,T);

if (zval)

for (k=i=1;; k++)

{

x = cgetg(lx, t_POL); p1 -= 2;

W = modulargcd(T,V); T = gdeuc(T,W); dW = degpol(W);

for (i=2; i<lx; i++) x[i] = p1[i];

/* W = prod P^e, e > k; V = prod P^e, e >= k */

x[1] = evalsigne(1)|evalvarn(vx)|evallgef(lx);

if (dW != degpol(V)) { P[i] = ldeuc(V,W); e[i] = k; i++; }

nbfac++;

if (dW <= 0) break;

V = W;

}

/* now x(0) != 0 */

if (val) { P[i] = lpolx[varn(f)]; e[i] = val; i++;}

if (lx==3) { fa = NULL;/* for lint */ goto END; }

setlg(P,i);

p1 = cgetg(1,t_VEC); fa=cgetg(lx,t_VEC);

setlg(e,i); *ex=e; return P;

for (i=1; i<lx; i++) fa[i] = (long)p1;

}

d=content(x); if (gsigne(leading_term(x)) < 0) d = gneg_i(d);

if (!gcmp1(d)) x=gdiv(x,d);

if (lx==4) { nbfac++; ex++; fa[1] = (long)concatsp(p1,x); goto END; }

w=derivpol(x); t=modulargcd(x,w);

GEN

if (!gcmp1(t)) { x=gdeuc(x,t); w=gdeuc(w,t); }

fact_from_DDF(GEN fa, GEN e, long n)

k=1;

{

while (k)

GEN v,w, y = cgetg(3, t_MAT);

long i,j,k, l = lg(fa);

v = cgetg(n+1,t_COL); y[1] = (long)v;

w = cgetg(n+1,t_COL); y[2] = (long)w;

for (k=i=1; i<l; i++)

{

ex++; w=gadd(w, gneg_i(derivpol(x))); k=signe(w);

GEN L = (GEN)fa[i], ex = stoi(e[i]);

if (k) { t=modulargcd(x,w); x=gdeuc(x,t); w=gdeuc(w,t); } else t=x;

long J = lg(L);

if (degpol(t) > 0)

for (j=1; j<J; j++,k++)

{

fa[ex] = (long)squff2(t,hint);

v[k] = lcopy((GEN)L[j]);

nbfac += lg(fa[ex])-1;

w[k] = (long)ex;

}

END: av2=avma;

return y;

v=cgetg(nbfac+1,t_COL); y[1]=(long)v;

w=cgetg(nbfac+1,t_COL); y[2]=(long)w;

if (zval) { v[1]=lpolx[vx]; w[1]=lstoi(zval); k=1; } else k=0;

for (i=1; i<=ex; i++)

for (j=1; j<lg(fa[i]); j++)

{

k++; v[k]=lcopy(gmael(fa,i,j)); w[k]=lstoi(i);

}

gerepilemanyvec(av,av2,y+1,2);

return sort_factor(y, cmpii);

}

/* Factor x in Z[t]. Assume all factors have degree divisible by hint */

GEN

factpol(GEN x, long hint)

{

gpmem_t av = avma;

GEN fa,ex,y;

long n,i,l;

if (typ(x)!=t_POL) err(notpoler,"factpol");

if (!signe(x)) err(zeropoler,"factpol");

fa = ZX_squff(x, &ex);

l = lg(fa); n = 0;

for (i=1; i<l; i++)

{

fa[i] = (long)DDF2((GEN)fa[i], hint);

n += lg(fa[i])-1;

}

y = fact_from_DDF(fa,ex,n);

return gerepileupto(av, sort_factor(y, cmpii));

}

/***********************************************************************/

/** **/

/** FACTORIZATION **/

Line 1630 poltype(GEN x, GEN *ptp, GEN *ptpol, long *ptpa)

Line 1707 poltype(GEN x, GEN *ptp, GEN *ptpol, long *ptpa)

assign_or_fail((GEN)p1[1],p);

t[3]=1; break;

case t_COMPLEX:

if (!pcx)

if (!pcx) pcx = coefs_to_pol(3, gun,gzero,gun); /* x^2 + 1 */

{

pcx = cgetg(5,t_POL); /* x^2 + 1 */

pcx[1] = evalsigne(1)|evalvarn(0)|m_evallgef(5),

pcx[2]=pcx[4]=un; pcx[3]=zero;

}

for (j=1; j<=2; j++)

{

p2 = (GEN)p1[j];

Line 1756 factor0(GEN x,long flag)

Line 1828 factor0(GEN x,long flag)

return NULL; /* not reached */

}

/* assume f and g coprime integer factorizations */

GEN

merge_factor_i(GEN f, GEN g)

concat_factor(GEN f, GEN g)

{

GEN h;

if (lg(f) == 1) return g;

Line 1766 merge_factor_i(GEN f, GEN g)

Line 1837 merge_factor_i(GEN f, GEN g)

h = cgetg(3,t_MAT);

h[1] = (long)concatsp((GEN)f[1], (GEN)g[1]);

h[2] = (long)concatsp((GEN)f[2], (GEN)g[2]);

return sort_factor_gen(h, cmpii);

return h;

}

/* assume f and g coprime integer factorizations */

GEN

merge_factor_i(GEN f, GEN g)

{

if (lg(f) == 1) return g;

if (lg(g) == 1) return f;

return sort_factor_gen(concat_factor(f,g), cmpii);

}

GEN

deg1_from_roots(GEN L, long v)

{

long i, l = lg(L);

GEN z = cgetg(l,t_COL);

for (i=1; i<l; i++)

z[i] = (long)deg1pol_i(gun, gneg((GEN)z[i]), v);

return z;

}

GEN

factor(GEN x)

{

long tx=typ(x),lx,av,tetpil,i,j,pa,v,r1;

long tx=typ(x), lx, i, j, pa, v, r1;

gpmem_t av, tetpil;

GEN y,p,p1,p2,p3,p4,p5,pol;

if (is_matvec_t(tx))

Line 1794 factor(GEN x)

Line 1885 factor(GEN x)

case t_INT: return decomp(x);

case t_FRACN:

x = gred(x); /* fall through */

return gerepileupto(av, factor(gred(x)));

case t_FRAC:

p1 = decomp((GEN)x[1]);

p2 = decomp((GEN)x[2]); p2[2] = (long)gneg_i((GEN)p2[2]);

Line 1810 factor(GEN x)

Line 1901 factor(GEN x)

case t_COMPLEX: y=cgetg(3,t_MAT); lx=lgef(x)-2; v=varn(x);

av = avma; p1 = roots(x,pa); tetpil = avma;

p2=cgetg(lx,t_COL);

p2 = deg1_from_roots(p1, v);

for (i=1; i<lx; i++)

p2[i] = (long)deg1pol_i(gun, gneg((GEN)p1[i]), v);

y[1]=lpile(av,tetpil,p2);

p3=cgetg(lx,t_COL); for (i=1; i<lx; i++) p3[i] = un;

y[2]=(long)p3; return y;

Line 1857 factor(GEN x)

Line 1946 factor(GEN x)

p2[2] = (long)deg1pol_i((GEN)p1[2], (GEN)p1[1], v);

}

killv = (avma != (ulong)pol);

killv = (avma != (gpmem_t)pol);

if (killv) setvarn(pol, fetch_var());

switch (typ2(tx))

{

Line 1869 factor(GEN x)

Line 1958 factor(GEN x)

switch (typ1(tx))

{

case t_POLMOD:

if (killv) delete_var();

if (killv) (void)delete_var();

return gerepileupto(av,p1);

case t_COMPLEX: p5 = gi; break;

case t_QUAD: p5=cgetg(4,t_QUAD);

Line 1888 factor(GEN x)

Line 1977 factor(GEN x)

if(typ(p4)==t_POLMOD) p3[j]=lsubst((GEN)p4[2],v,p5);

}

if (killv) delete_var();

if (killv) (void)delete_var();

tetpil=avma; y=cgetg(3,t_MAT);

y[1]=lcopy(p2);y[2]=lcopy((GEN)p1[2]);

return gerepile(av,tetpil,y);

Line 1896 factor(GEN x)

Line 1985 factor(GEN x)

}

case t_RFRACN:

x=gred_rfrac(x); /* fall through */

return gerepileupto(av, factor(gred_rfrac(x)));

case t_RFRAC:

p1=factor((GEN)x[1]);

p2=factor((GEN)x[2]); p3=gneg_i((GEN)p2[2]);

Line 1913 factor(GEN x)

Line 2002 factor(GEN x)

#undef typs

#undef tsh

/* assume n != 0, t_INT. Compute x^|n| using left-right binary powering */

GEN

leftright_pow(GEN x, GEN n, void *data,

GEN (*sqr)(void*,GEN), GEN (*mul)(void*,GEN,GEN))

{

GEN nd = n+2, y = x;

long i, m = *nd, j = 1+bfffo((ulong)m);

gpmem_t av = avma, lim = stack_lim(av, 1);

/* normalize, i.e set highest bit to 1 (we know m != 0) */

m<<=j; j = BITS_IN_LONG-j;

/* first bit is now implicit */

for (i=lgefint(n)-2;;)

{

for (; j; m<<=1,j--)

{

y = sqr(data,y);

if (m < 0) y = mul(data,y,x); /* first bit set: multiply by base */

if (low_stack(lim, stack_lim(av,1)))

{

if (DEBUGMEM>1) err(warnmem,"leftright_pow");

y = gerepilecopy(av, y);

}

if (--i == 0) return avma==av? gcopy(y): y;

m = *++nd; j = BITS_IN_LONG;

}

GEN

divide_conquer_prod(GEN x, GEN (*mul)(GEN,GEN))

{

long i,k,lx = lg(x);

Line 1936 divide_conquer_prod(GEN x, GEN (*mul)(GEN,GEN))

Line 2054 divide_conquer_prod(GEN x, GEN (*mul)(GEN,GEN))

static GEN static_nf;

static GEN

myidealmulred(GEN x, GEN y) { return idealmulred(static_nf, x, y, 0); }

idmulred(GEN x, GEN y) { return idealmulred(static_nf, x, y, 0); }

static GEN

myidealpowred(GEN x, GEN n) { return idealpowred(static_nf, x, n, 0); }

idpowred(GEN x, GEN n) { return idealpowred(static_nf, x, n, 0); }

static GEN

myidealmul(GEN x, GEN y) { return idealmul(static_nf, x, y); }

idmul(GEN x, GEN y) { return idealmul(static_nf, x, y); }

static GEN

myidealpow(GEN x, GEN n) { return idealpow(static_nf, x, n); }

idpow(GEN x, GEN n) { return idealpow(static_nf, x, n); }

static GEN

eltmul(GEN x, GEN y) { return element_mul(static_nf, x, y); }

static GEN

eltpow(GEN x, GEN n) { return element_pow(static_nf, x, n); }

GEN

factorback_i(GEN fa, GEN nf, int red)

_factorback(GEN fa, GEN e, GEN (*_mul)(GEN,GEN), GEN (*_pow)(GEN,GEN))

{

long av=avma,k,l,lx,t=typ(fa);

gpmem_t av = avma;

GEN ex, x;

long k,l,lx,t = typ(fa);

GEN (*_mul)(GEN,GEN);

GEN p,x;

GEN (*_pow)(GEN,GEN);

if (nf)

if (e) /* supplied vector of exponents */

p = fa;

else /* genuine factorization */

{

static_nf = nf;

if (t != t_MAT || lg(fa) != 3)

if (red)

{

_mul = &myidealmulred;

if (!is_vec_t(t)) err(talker,"not a factorisation in factorback");

_pow = &myidealpowred;

return gerepileupto(av, divide_conquer_prod(fa, _mul));

}

else

p = (GEN)fa[1];

{

e = (GEN)fa[2];

_mul = &myidealmul;

_pow = &myidealpow;

}

else

lx = lg(p);

t = t_INT; /* dummy */

/* check whether e is an integral vector of correct length */

if (is_vec_t(typ(e)) && lx == lg(e))

{

_mul = &gmul;

for (k=1; k<lx; k++)

_pow = &powgi;

if (typ(e[k]) != t_INT) break;

if (k == lx) t = t_MAT;

}

if ( t == t_VEC || t == t_COL)

if (t != t_MAT) err(talker,"not a factorisation in factorback");

return gerepileupto(av, divide_conquer_prod(fa, _mul));

if (lx == 1) return gun;

if (t!=t_MAT || lg(fa)!=3)

err(talker,"not a factorisation in factorback");

ex=(GEN)fa[2]; fa=(GEN)fa[1];

lx = lg(fa); if (lx == 1) return gun;

x = cgetg(lx,t_VEC);

for (l=1,k=1; k<lx; k++)

if (signe(ex[k]))

if (signe(e[k]))

x[l++] = (long)_pow((GEN)fa[k],(GEN)ex[k]);

x[l++] = (long)_pow((GEN)p[k],(GEN)e[k]);

setlg(x,l);

return gerepileupto(av, divide_conquer_prod(x, _mul));

}

GEN

factorback_i(GEN fa, GEN e, GEN nf, int red)

{

if (!nf && e && lg(e) > 1 && typ(e[1]) != t_INT) { nf = e; e = NULL; }

if (!nf) return _factorback(fa, e, &gmul, &powgi);

static_nf = checknf(nf);

if (red)

return _factorback(fa, e, &idmulred, &idpowred);

else

return _factorback(fa, e, &idmul, &idpow);

}

GEN

factorbackelt(GEN fa, GEN e, GEN nf)

{

if (!nf && e && lg(e) > 1 && typ(e[1]) != t_INT) { nf = e; e = NULL; }

if (!nf) err(talker, "missing nf in factorbackelt");

static_nf = checknf(nf);

return _factorback(fa, e, &eltmul, &eltpow);

}

GEN

factorback0(GEN fa, GEN e, GEN nf)

{

return factorback_i(fa,e,nf,0);

}

GEN

factorback(GEN fa, GEN nf)

{

return factorback_i(fa,nf,0);

return factorback_i(fa,nf,NULL,0);

}

GEN

gisirreducible(GEN x)

{

long av=avma, tx = typ(x),l,i;

long tx = typ(x), l, i;

gpmem_t av=avma;

GEN y;

if (is_matvec_t(tx))

Line 2035 gcd0(GEN x, GEN y, long flag)

Line 2182 gcd0(GEN x, GEN y, long flag)

static GEN

triv_cont_gcd(GEN x, GEN y)

{

long av = avma, tetpil;

gpmem_t av = avma, tetpil;

GEN p1 = (typ(x)==t_COMPLEX)? ggcd((GEN)x[1],(GEN)x[2])

: ggcd((GEN)x[2],(GEN)x[3]);

tetpil=avma; return gerepile(av,tetpil,ggcd(p1,y));

Line 2044 triv_cont_gcd(GEN x, GEN y)

Line 2191 triv_cont_gcd(GEN x, GEN y)

static GEN

cont_gcd(GEN x, GEN y)

{

long av = avma, tetpil;

gpmem_t av = avma, tetpil;

GEN p1 = content(x);

tetpil=avma; return gerepile(av,tetpil,ggcd(p1,y));

Line 2056 padic_gcd(GEN x, GEN y)

Line 2203 padic_gcd(GEN x, GEN y)

{

long v = ggval(x,(GEN)y[2]), w = valp(y);

if (w < v) v = w;

return gpuigs((GEN)y[2], v);

return gpowgs((GEN)y[2], v);

}

/* x,y in Z[i], at least one of which is t_COMPLEX */

static GEN

gaussian_gcd(GEN x, GEN y)

{

gpmem_t av = avma;

GEN dx = denom(x);

GEN dy = denom(y);

GEN den = mpppcm(dx, dy);

x = gmul(x, dx);

y = gmul(y, dy);

while (!gcmp0(y))

{

GEN z = gdiv(x,y);

GEN r0 = greal(z), r = gfloor(r0);

GEN i0 = gimag(z), i = gfloor(i0);

if (gcmp(gsub(r0,r), ghalf) > 0) r = addis(r,1);

if (gcmp(gsub(i0,i), ghalf) > 0) i = addis(i,1);

if (gcmp0(i)) z = r;

else

{

z = cgetg(3, t_COMPLEX);

z[1] = (long)r;

z[2] = (long)i;

}

z = gsub(x, gmul(z,y));

x = y; y = z;

}

if (signe(greal(x)) < 0) x = gneg(x);

if (signe(gimag(x)) < 0) x = gmul(x, gi);

if (typ(x) == t_COMPLEX)

{

if (gcmp0((GEN)x[2])) x = (GEN)x[1];

else if (gcmp0((GEN)x[1])) x = (GEN)x[2];

}

return gerepileupto(av, gdiv(x, den));

}

#define fix_frac(z) if (signe(z[2])<0)\

setsigne(z[1],-signe(z[1]));\

setsigne(z[2],1);\

}

int

isrational(GEN x)

{

long i, t = typ(x);

if (t != t_POL) return is_rational_t(t);

for (i=lgef(x)-1; i>1; i--)

{

t = typ(x[i]);

if (!is_rational_t(t)) return 0;

}

return 1;

}

static int

cx_isrational(GEN x)

{

return (isrational((GEN)x[1]) && isrational((GEN)x[2]));

}

/* y == 0 */

static GEN

zero_gcd(GEN x, GEN y)

{

if (typ(y) == t_PADIC)

{

GEN p = (GEN)y[2];

long v = ggval(x,p), w = valp(y);

if (w < v) return padiczero(p, w);

if (gcmp0(x)) return padiczero(p, v);

return gpowgs(p, v);

}

switch(typ(x))

{

case t_INT: case t_FRAC: case t_FRACN:

return gabs(x,0);

case t_COMPLEX:

if (cx_isrational(x)) return gaussian_gcd(x, gzero);

/* fall through */

case t_REAL:

return gun;

default:

return gcopy(x);

}

GEN

ggcd(GEN x, GEN y)

{

long l,av,tetpil,i,vx,vy, tx = typ(x), ty = typ(y);

long l, i, vx, vy, tx = typ(x), ty = typ(y);

gpmem_t av, tetpil;

GEN p1,z;

if (tx>ty) { swap(x,y); lswap(tx,ty); }

Line 2078 ggcd(GEN x, GEN y)

Line 2311 ggcd(GEN x, GEN y)

for (i=1; i<l; i++) z[i]=lgcd(x,(GEN)y[i]);

return z;

}

if (is_noncalc_t(tx) || is_noncalc_t(ty)) err(operf,"g",tx,ty);

if (is_noncalc_t(tx) || is_noncalc_t(ty)) err(operf,"g",x,y);

if (gcmp0(x)) return gcopy(y);

if (gcmp0(x)) return zero_gcd(y, x);

if (gcmp0(y)) return gcopy(x);

if (gcmp0(y)) return zero_gcd(x, y);

if (is_const_t(tx))

{

if (ty == t_FRACN)

Line 2100 ggcd(GEN x, GEN y)

Line 2333 ggcd(GEN x, GEN y)

case t_INTMOD: z=cgetg(3,t_INTMOD);

if (egalii((GEN)x[1],(GEN)y[1]))

{ copyifstack(x[1],z[1]); }

copyifstack(x[1],z[1]);

else

z[1] = lmppgcd((GEN)x[1],(GEN)y[1]);

if (gcmp1((GEN)z[1])) z[2] = zero;

Line 2118 ggcd(GEN x, GEN y)

Line 2351 ggcd(GEN x, GEN y)

return z;

case t_COMPLEX:

av=avma; p1=gdiv(x,y);

if (cx_isrational(x) && cx_isrational(y)) return gaussian_gcd(x,y);

if (gcmp0((GEN)p1[2]))

{

p1 = (GEN)p1[1];

switch(typ(p1))

{

case t_INT:

avma=av; return gcopy(y);

case t_FRAC: case t_FRACN:

tetpil=avma; return gerepile(av,tetpil,gdiv(y,(GEN)p1[2]));

default: avma=av; return gun;

}

if (typ(p1[1])==t_INT && typ(p1[2])==t_INT) {avma=av; return gcopy(y);}

p1 = ginv(p1); avma=av;

if (typ(p1[1])==t_INT && typ(p1[2])==t_INT) return gcopy(x);

return triv_cont_gcd(y,x);

case t_PADIC:

if (!egalii((GEN)x[2],(GEN)y[2])) return gun;

vx = valp(x);

vy = valp(y); return gpuigs((GEN)y[2], min(vy,vx));

vy = valp(y); return gpowgs((GEN)y[2], min(vy,vx));

case t_QUAD:

av=avma; p1=gdiv(x,y);

Line 2171 ggcd(GEN x, GEN y)

Line 2389 ggcd(GEN x, GEN y)

z[1] = lmppgcd(x,(GEN)y[1]);

z[2] = licopy((GEN)y[2]); return z;

case t_COMPLEX: case t_QUAD:

case t_COMPLEX:

if (cx_isrational(y)) return gaussian_gcd(x,y);

case t_QUAD:

return triv_cont_gcd(y,x);

case t_PADIC:

Line 2185 ggcd(GEN x, GEN y)

Line 2405 ggcd(GEN x, GEN y)

{

case t_FRAC:

av = avma; p1=mppgcd((GEN)x[1],(GEN)y[2]); avma = av;

if (!is_pm1(p1)) err(operi,"g",tx,ty);

if (!is_pm1(p1)) err(operi,"g",x,y);

return ggcd((GEN)y[1], x);

case t_COMPLEX: case t_QUAD:

Line 2198 ggcd(GEN x, GEN y)

Line 2418 ggcd(GEN x, GEN y)

case t_FRAC:

switch(ty)

{

case t_COMPLEX: case t_QUAD:

case t_COMPLEX:

if (cx_isrational(y)) return gaussian_gcd(x,y);

case t_QUAD:

return triv_cont_gcd(y,x);

case t_PADIC:

Line 2226 ggcd(GEN x, GEN y)

Line 2448 ggcd(GEN x, GEN y)

{

case t_POLMOD: z=cgetg(3,t_POLMOD);

if (gegal((GEN)x[1],(GEN)y[1]))

{ copyifstack(x[1],z[1]); }

copyifstack(x[1],z[1]);

else

z[1] = lgcd((GEN)x[1],(GEN)y[1]);

if (lgef(z[1])<=3) z[2]=zero;

Line 2251 ggcd(GEN x, GEN y)

Line 2473 ggcd(GEN x, GEN y)

case t_RFRAC:

av = avma; p1=ggcd((GEN)x[1],(GEN)y[2]); avma = av;

if (!gcmp1(p1)) err(operi,"g",tx,ty);

if (!gcmp1(p1)) err(operi,"g",x,y);

return ggcd((GEN)y[1],x);

case t_RFRACN:

Line 2267 ggcd(GEN x, GEN y)

Line 2489 ggcd(GEN x, GEN y)

return srgcd(x,y);

case t_SER:

return gpuigs(polx[vx],min(valp(y),gval(x,vx)));

return gpowgs(polx[vx],min(valp(y),gval(x,vx)));

case t_RFRAC: case t_RFRACN: av=avma; z=cgetg(3,ty);

z[1]=lgcd(x,(GEN)y[1]);

Line 2279 ggcd(GEN x, GEN y)

Line 2501 ggcd(GEN x, GEN y)

switch(ty)

{

case t_SER:

return gpuigs(polx[vx],min(valp(x),valp(y)));

return gpowgs(polx[vx],min(valp(x),valp(y)));

case t_RFRAC: case t_RFRACN:

return gpuigs(polx[vx],min(valp(x),gval(y,vx)));

return gpowgs(polx[vx],min(valp(x),gval(y,vx)));

}

break;

case t_RFRAC: case t_RFRACN: z=cgetg(3,ty);

if (!is_rfrac_t(ty)) err(operf,"g",tx,ty);

if (!is_rfrac_t(ty)) err(operf,"g",x,y);

p1 = gdiv((GEN)y[2], ggcd((GEN)x[2], (GEN)y[2]));

tetpil = avma;

z[2] = lpile((long)z,tetpil,gmul(p1, (GEN)x[2]));

z[2] = lpile((gpmem_t)z,tetpil,gmul(p1, (GEN)x[2]));

z[1] = lgcd((GEN)x[1], (GEN)y[1]); return z;

}

err(operf,"g",tx,ty);

err(operf,"g",x,y);

return NULL; /* not reached */

}

Line 2305 GEN glcm0(GEN x, GEN y)

Line 2527 GEN glcm0(GEN x, GEN y)

GEN

glcm(GEN x, GEN y)

{

long av,tx,ty,i,l;

long tx, ty, i, l;

gpmem_t av;

GEN p1,p2,z;

ty = typ(y);

Line 2340 glcm(GEN x, GEN y)

Line 2563 glcm(GEN x, GEN y)

return gerepileupto(av,p2);

}

/* x + r ~ x ? Assume x,r are t_POL, deg(r) <= deg(x) */

static int

pol_approx0(GEN r, GEN x, int exact)

{

long i, lx,lr;

if (exact) return gcmp0(r);

lx = lgef(x);

lr = lgef(r); if (lr < lx) lx = lr;

for (i=2; i<lx; i++)

if (!approx_0((GEN)r[i], (GEN)x[i])) return 0;

return 1;

}

static GEN

polgcdnun(GEN x, GEN y)

{

long av1, av = avma, lim = stack_lim(av,1);

gpmem_t av1, av = avma, lim = stack_lim(av, 1);

GEN p1, yorig = y;

GEN r, yorig = y;

int exact = !(isinexactreal(x) || isinexactreal(y));

for(;;)

{

av1 = avma; p1 = gres(x,y);

av1 = avma; r = gres(x,y);

if (gcmp0(p1))

if (pol_approx0(r, x, exact))

{

avma = av1;

if (lgef(y) == 3 && isinexactreal(y)) { avma = av; return gun; }

if (lgef(y) == 3 && !exact) { avma = av; return gun; }

return (y==yorig)? gcopy(y): gerepileupto(av,y);

}

x = y; y = p1;

x = y; y = r;

if (low_stack(lim,stack_lim(av,1)))

{

GEN *gptr[2]; x = gcopy(x); gptr[0]=&x; gptr[1]=&y;

Line 2365 polgcdnun(GEN x, GEN y)

Line 2602 polgcdnun(GEN x, GEN y)

}

static int issimplefield(GEN x);

static int

issimplepol(GEN x)

{

long i, lx = lgef(x);

for (i=2; i<lx; i++)

if (issimplefield((GEN)x[i])) return 1;

return 0;

}

/* return 1 if coeff explosion is not possible */

static int

issimplefield(GEN x)

{

long lx,i;

switch(typ(x))

{

case t_REAL: case t_INTMOD: case t_PADIC: case t_SER:

return 1;

case t_POL:

case t_COMPLEX:

lx=lgef(x);

for (i=2; i<lx; i++)

if (issimplefield((GEN)x[i])) return 1;

return 0;

case t_COMPLEX: case t_POLMOD:

return issimplefield((GEN)x[1]) || issimplefield((GEN)x[2]);

case t_POLMOD:

return issimplepol((GEN)x[1]) || issimplepol((GEN)x[2]);

}

return 0;

}

int

isrational(GEN x)

{

long i;

for (i=lgef(x)-1; i>1; i--)

switch(typ(x[i]))

{

case t_INT: case t_FRAC: break;

default: return 0;

}

return 1;

}

static int

can_use_modular_gcd(GEN x, GEN *mod, long *v)

{

Line 2423 can_use_modular_gcd(GEN x, GEN *mod, long *v)

Line 2653 can_use_modular_gcd(GEN x, GEN *mod, long *v)

break;

case t_POL:

if (!isrational(p1)) return 0;

if (*v >= 0)

{

if (*v != varn(p1)) return 0;

} else *v = varn(p1);

Line 2456 isinexactfield(GEN x)

Line 2686 isinexactfield(GEN x)

static GEN

gcdmonome(GEN x, GEN y)

{

long tetpil,av=avma, dx=degpol(x), v=varn(x), e=gval(y,v);

long dx=degpol(x), v=varn(x), e=gval(y, v);

gpmem_t tetpil, av=avma;

GEN p1,p2;

if (dx < e) e = dx;

p1=ggcd(leading_term(x),content(y)); p2=gpuigs(polx[v],e);

p1=ggcd(leading_term(x),content(y)); p2=gpowgs(polx[v],e);

tetpil=avma; return gerepile(av,tetpil,gmul(p1,p2));

}

Line 2473 gcdmonome(GEN x, GEN y)

Line 2704 gcdmonome(GEN x, GEN y)

static GEN

polinvinexact(GEN x, GEN y)

{

ulong av = avma;

gpmem_t av = avma;

long i,dx=degpol(x),dy=degpol(y),lz=dx+dy;

GEN v,z;

Line 2490 polinvinexact(GEN x, GEN y)

Line 2721 polinvinexact(GEN x, GEN y)

static GEN

polinvmod(GEN x, GEN y)

{

long av,av1,tx,vx=varn(x),vy=varn(y);

long tx, vx=varn(x), vy=varn(y);

gpmem_t av, av1;

GEN u,v,d,p1;

while (vx!=vy)

Line 2560 GEN

Line 2792 GEN

content(GEN x)

{

long lx,i,tx=typ(x);

ulong av,tetpil;

gpmem_t av, tetpil;

GEN p1,p2;

if (is_scalar_t(tx))

return tx==t_POLMOD? content((GEN)x[2]): gcopy(x);

{

switch (tx)

{

case t_INT: return absi(x);

case t_FRAC:

case t_FRACN: return gabs(x,0);

case t_POLMOD: return content((GEN)x[2]);

}

return gcopy(x);

}

av = avma;

switch(tx)

{

Line 2614 content(GEN x)

Line 2855 content(GEN x)

}

GEN

primitive_part(GEN x, GEN *c)

primitive_part(GEN x, GEN *ptc)

{

*c = content(x);

gpmem_t av = avma;

if (gcmp1(*c)) *c = NULL; else x = gdiv(x,*c);

GEN c = content(x);

if (gcmp1(c)) { avma = av; c = NULL; }

else if (!gcmp0(c)) x = gdiv(x,c);

if (ptc) *ptc = c;

return x;

}

GEN

Q_primitive_part(GEN x, GEN *ptc)

{

gpmem_t av = avma;

GEN c = content(x);

if (gcmp1(c)) { avma = av; c = NULL; }

else if (!gcmp0(c)) x = Q_div_to_int(x,c);

if (ptc) *ptc = c;

return x;

}

GEN

primpart(GEN x) { return primitive_part(x, NULL); }

GEN

Q_primpart(GEN x) { return Q_primitive_part(x, NULL); }

/* NOT MEMORY CLEAN (because of t_FRAC).

* As denom(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead

* of Q(x2,...,xn)[x1] */

GEN

Q_denom(GEN x)

{

long i, l = lgef(x);

GEN d, D;

gpmem_t av;

switch(typ(x))

{

case t_INT: return gun;

case t_FRAC: return (GEN)x[2];

case t_VEC: case t_COL: case t_MAT:

l = lg(x); if (l == 1) return gun;

av = avma; d = Q_denom((GEN)x[1]);

for (i=2; i<l; i++)

{

D = Q_denom((GEN)x[i]);

if (D != gun) d = mpppcm(d, D);

}

return gerepileuptoint(av, d);

case t_POL:

l = lgef(x); if (l == 2) return gun;

av = avma; d = Q_denom((GEN)x[2]);

for (i=3; i<l; i++)

{

D = Q_denom((GEN)x[i]);

if (D != gun) d = mpppcm(d, D);

}

return gerepileuptoint(av, d);

}

err(typeer,"Q_denom");

return NULL; /* not reached */

}

GEN

Q_remove_denom(GEN x, GEN *ptd)

{

GEN d = Q_denom(x);

if (gcmp1(d)) d = NULL; else x = Q_muli_to_int(x,d);

if (ptd) *ptd = d;

return x;

}

/* return y = x * d, assuming x rational, and d,y integral */

GEN

Q_muli_to_int(GEN x, GEN d)

{

long i, l, t = typ(x);

GEN y, xn, xd;

gpmem_t av;

if (typ(d) != t_INT) err(typeer,"Q_muli_to_int");

switch (t)

{

case t_VEC: case t_COL: case t_MAT:

l = lg(x); y = cgetg(l, t);

for (i=1; i<l; i++) y[i] = (long)Q_muli_to_int((GEN)x[i], d);

return y;

case t_POL: l = lgef(x);

y = cgetg(l, t_POL); y[1] = x[1];

for (i=2; i<l; i++) y[i] = (long)Q_muli_to_int((GEN)x[i], d);

return y;

case t_INT:

return mulii(x,d);

case t_FRAC:

xn = (GEN)x[1];

xd = (GEN)x[2]; av = avma;

y = mulii(xn, diviiexact(d, xd));

return gerepileuptoint(av, y);

}

err(typeer,"Q_muli_to_int");

return NULL; /* not reached */

}

/* return x * n/d. x: rational; d,n,result: integral. n = NULL represents 1 */

GEN

Q_divmuli_to_int(GEN x, GEN d, GEN n)

{

long i, l, t = typ(x);

GEN y, xn, xd;

gpmem_t av;

switch(t)

{

case t_VEC: case t_COL: case t_MAT:

l = lg(x); y = cgetg(l, t);

for (i=1; i<l; i++) y[i] = (long)Q_divmuli_to_int((GEN)x[i], d,n);

return y;

case t_POL: l = lgef(x);

y = cgetg(l, t_POL); y[1] = x[1];

for (i=2; i<l; i++) y[i] = (long)Q_divmuli_to_int((GEN)x[i], d,n);

return y;

case t_INT:

av = avma; y = diviiexact(x,d);

if (n) y = gerepileuptoint(av, mulii(y,n));

return y;

case t_FRAC:

xn = (GEN)x[1];

xd = (GEN)x[2]; av = avma;

y = mulii(diviiexact(xn, d), diviiexact(n, xd));

return gerepileuptoint(av, y);

}

err(typeer,"Q_divmuli_to_int");

return NULL; /* not reached */

}

/* return y = x / c, assuming x,c rational, and y integral */

GEN

Q_div_to_int(GEN x, GEN c)

{

GEN d, n;

switch(typ(c))

{

case t_INT:

return Q_divmuli_to_int(x, c, NULL);

case t_FRAC:

n = (GEN)c[1];

d = (GEN)c[2]; if (gcmp1(n)) return Q_muli_to_int(x,d);

return Q_divmuli_to_int(x, n,d);

}

err(typeer,"Q_div_to_int");

return NULL; /* not reached */

}

/*******************************************************************/

/* */

/* SUBRESULTANT */

Line 2647 gdivexact(GEN x, GEN y)

Line 3043 gdivexact(GEN x, GEN y)

case t_INTMOD:

case t_POLMOD: return gmul(x,ginv(y));

case t_POL:

if (varn(x)==varn(y)) return poldivres(x,y, ONLY_DIVIDES_EXACT);

if (varn(x)==varn(y)) return poldivres(x,y, ONLY_DIVIDES);

}

lx = lgef(x); z = cgetg(lx,t_POL);

for (i=2; i<lx; i++) z[i]=(long)gdivexact((GEN)x[i],y);

Line 2669 init_resultant(GEN x, GEN y)

Line 3065 init_resultant(GEN x, GEN y)

tx = typ(x); ty = typ(y);

if (is_scalar_t(tx) || is_scalar_t(ty))

{

if (tx==t_POL) return gpuigs(y,degpol(x));

if (tx==t_POL) return gpowgs(y,degpol(x));

if (ty==t_POL) return gpuigs(x,degpol(y));

if (ty==t_POL) return gpowgs(x,degpol(y));

return gun;

}

if (tx!=t_POL || ty!=t_POL) err(typeer,"subresall");

if (varn(x)==varn(y)) return NULL;

return (varn(x)<varn(y))? gpuigs(y,degpol(x)): gpuigs(x,degpol(y));

return (varn(x)<varn(y))? gpowgs(y,degpol(x)): gpowgs(x,degpol(y));

}

/* return coefficients s.t x = x_0 X^n + ... + x_n */

static GEN

GEN

revpol(GEN x)

{

long i,n = degpol(x);

Line 2689 revpol(GEN x)

Line 3085 revpol(GEN x)

return y;

}

/* assume dx >= dy, y non constant */

/* assume dx >= dy, y non constant. */

GEN

static GEN

pseudorem(GEN x, GEN y)

pseudorem_i(GEN x, GEN y, GEN mod)

{

long av = avma, av2,lim, vx = varn(x),dx,dy,dz,i,lx, p;

long vx = varn(x), dx, dy, dz, i, lx, p;

gpmem_t av = avma, av2, lim;

GEN (*red)(GEN,GEN) = NULL;

if (mod) red = (typ(mod) == t_INT)? &FpX_red: &gmod;

if (!signe(y)) err(talker,"euclidean division by zero (pseudorem)");

(void)new_chunk(2);

dx=degpol(x); x = revpol(x);

Line 2704 pseudorem(GEN x, GEN y)

Line 3104 pseudorem(GEN x, GEN y)

{

x[0] = lneg((GEN)x[0]); p--;

for (i=1; i<=dy; i++)

{

x[i] = ladd(gmul((GEN)y[0], (GEN)x[i]), gmul((GEN)x[0],(GEN)y[i]));

if (red) x[i] = (long)red((GEN)x[i], mod);

}

for ( ; i<=dx; i++)

{

x[i] = lmul((GEN)y[0], (GEN)x[i]);

if (red) x[i] = (long)red((GEN)x[i], mod);

}

do { x++; dx--; } while (dx >= 0 && gcmp0((GEN)x[0]));

if (dx < dy) break;

if (low_stack(lim,stack_lim(av2,1)))

Line 2720 pseudorem(GEN x, GEN y)

Line 3126 pseudorem(GEN x, GEN y)

x[0]=evaltyp(t_POL) | evallg(lx);

x[1]=evalsigne(1) | evalvarn(vx) | evallgef(lx);

x = revpol(x) - 2;

return gerepileupto(av, gmul(x, gpowgs((GEN)y[0], p)));

if (p)

{ /* multiply by y[0]^p [beware dummy vars from FpY_FpXY_resultant] */

GEN t = (GEN)y[0];

if (mod)

{ /* assume p fairly small */

for (i=1; i<p; i++)

t = red(gmul(t, (GEN)y[0]), mod);

}

else

t = gpowgs(t, p);

for (i=2; i<lx; i++)

{

x[i] = lmul((GEN)x[i], t);

if (red) x[i] = (long)red((GEN)x[i], mod);

}

if (!red) return gerepileupto(av, x);

}

return gerepilecopy(av, x);

}

extern void gerepilemanycoeffs2(long av, GEN x, long n, GEN y, long o);

GEN

pseudorem(GEN x, GEN y) { return pseudorem_i(x,y, NULL); }

/* assume dx >= dy, y non constant

* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */

GEN

pseudodiv(GEN x, GEN y, GEN *ptr)

{

long av = avma, av2,lim, vx = varn(x),dx,dy,dz,i,iz,lx,lz,p;

long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;

gpmem_t av = avma, av2, lim;

GEN z, r, ypow;

if (!signe(y)) err(talker,"euclidean division by zero (pseudorem)");

if (!signe(y)) err(talker,"euclidean division by zero (pseudodiv)");

(void)new_chunk(2);

dx=degpol(x); x = revpol(x);

dy=degpol(y); y = revpol(y); dz=dx-dy; p = dz+1;

Line 2751 pseudodiv(GEN x, GEN y, GEN *ptr)

Line 3176 pseudodiv(GEN x, GEN y, GEN *ptr)

x[i] = ladd(gmul((GEN)y[0], (GEN)x[i]), gmul((GEN)x[0],(GEN)y[i]));

for ( ; i<=dx; i++)

x[i] = lmul((GEN)y[0], (GEN)x[i]);

x++; dx--;

while (dx >= dy && gcmp0((GEN)x[0])) { x++; dx--; z[iz++] = zero; }

if (dx < dy) break;

if (low_stack(lim,stack_lim(av2,1)))

Line 2761 pseudodiv(GEN x, GEN y, GEN *ptr)

Line 3186 pseudodiv(GEN x, GEN y, GEN *ptr)

}

while (dx >= 0 && gcmp0((GEN)x[0])) { x++; dx--; }

if (dx < 0)

x = zeropol(vx);

else

{

Line 2776 pseudodiv(GEN x, GEN y, GEN *ptr)

Line 3201 pseudodiv(GEN x, GEN y, GEN *ptr)

z[1] = evalsigne(1) | evalvarn(vx) | evallgef(lz);

z = revpol(z) - 2;

r = gmul(x, (GEN)ypow[p]);

{

gerepileall(av, 2, &z, &r);

GEN *gptr[2]; gptr[0] = &z; gptr[1] = &r;

gerepilemany(av,gptr,2);

}

*ptr = r; return z;

}

/* Return resultant(u,v). If sol != NULL: set *sol to the last non-zero

* polynomial in the prs IF the sequence was computed, and gzero otherwise

* polynomial in the prs IF the sequence was computed, and gzero otherwise */

GEN

subresall(GEN u, GEN v, GEN *sol)

{

ulong av,av2,lim;

gpmem_t av, av2, lim;

long degq,dx,dy,du,dv,dr,signh;

GEN z,g,h,r,p1,p2,cu,cv;

Line 2797 subresall(GEN u, GEN v, GEN *sol)

Line 3218 subresall(GEN u, GEN v, GEN *sol)

if ((r = init_resultant(u,v))) return r;

if (isinexactfield(u) || isinexactfield(v)) return resultant2(u,v);

dx=lgef(u); dy=lgef(v); signh=1;

dx=degpol(u); dy=degpol(v); signh=1;

if (dx<dy)

if (dx < dy)

{

swap(u,v); lswap(dx,dy);

if (both_even(dx, dy)) signh = -signh;

if (both_odd(dx, dy)) signh = -signh;

}

if (dy==3) return gpowgs((GEN)v[2],dx-3);

if (dy==0) return gpowgs((GEN)v[2],dx);

av = avma;

u = primitive_part(u, &cu);

v = primitive_part(v, &cv);

Line 2813 subresall(GEN u, GEN v, GEN *sol)

Line 3234 subresall(GEN u, GEN v, GEN *sol)

r = pseudorem(u,v); dr = lgef(r);

if (dr == 2)

{

if (sol) { avma = (long)(r+2); *sol=gerepileupto(av,v); } else avma = av;

if (sol) { avma = (gpmem_t)(r+2); *sol=gerepileupto(av,v); } else avma = av;

return gzero;

}

du = lgef(u); dv = lgef(v); degq = du-dv;

du = degpol(u); dv = degpol(v); degq = du-dv;

u = v; p1 = g; g = leading_term(u);

switch(degq)

{

Line 2827 subresall(GEN u, GEN v, GEN *sol)

Line 3248 subresall(GEN u, GEN v, GEN *sol)

p1 = gmul(gpowgs(h,degq),p1);

h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1));

}

if (both_even(du,dv)) signh = -signh;

if (both_odd(du,dv)) signh = -signh;

v = gdivexact(r,p1);

if (dr==3) break;

if (low_stack(lim,stack_lim(av2,1)))

{

GEN *gptr[4]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;

if(DEBUGMEM>1) err(warnmem,"subresall, dr = %ld",dr);

gerepilemany(av2,gptr,4);

gerepileall(av2,4, &u, &v, &g, &h);

}

z = (GEN)v[2];

if (dv > 4) z = gdivexact(gpowgs(z,dv-3), gpowgs(h,dv-4));

if (dv > 1) z = gdivexact(gpowgs(z,dv), gpowgs(h,dv-1));

if (signh < 0) z = gneg(z); /* z = resultant(ppart(x), ppart(y)) */

p2 = gun;

if (cu) p2 = gmul(p2, gpowgs(cu,dy-3));

if (cu) p2 = gmul(p2, gpowgs(cu,dy));

if (cv) p2 = gmul(p2, gpowgs(cv,dx-3));

if (cv) p2 = gmul(p2, gpowgs(cv,dx));

z = gmul(z, p2);

if (sol) u = gclone(u);

Line 2854 subresall(GEN u, GEN v, GEN *sol)

Line 3274 subresall(GEN u, GEN v, GEN *sol)

static GEN

scalar_res(GEN x, GEN y, GEN *U, GEN *V)

{

long dx=lgef(x)-4;

*V=gpowgs(y,degpol(x)-1); *U=gzero; return gmul(y,*V);

*V=gpuigs(y,dx); *U=gzero; return gmul(y,*V);

}

/* compute U, V s.t Ux + Vy = resultant(x,y) */

GEN

subresext(GEN x, GEN y, GEN *U, GEN *V)

{

ulong av,av2,tetpil,lim;

gpmem_t av, av2, tetpil, lim;

long degq,tx,ty,dx,dy,du,dv,dr,signh;

GEN q,z,g,h,r,p1,p2,cu,cv,u,v,um1,uze,vze,xprim,yprim;

GEN q,z,g,h,r,p1,p2,cu,cv,u,v,um1,uze,vze,xprim,yprim, *gptr[3];

if (gcmp0(x) || gcmp0(y)) { *U = *V = gzero; return gzero; }

tx=typ(x); ty=typ(y);

Line 2878 subresext(GEN x, GEN y, GEN *U, GEN *V)

Line 3297 subresext(GEN x, GEN y, GEN *U, GEN *V)

if (varn(x) != varn(y))

return (varn(x)<varn(y))? scalar_res(x,y,U,V)

: scalar_res(y,x,V,U);

dx = lgef(x); dy = lgef(y); signh = 1;

dx = degpol(x); dy = degpol(y); signh = 1;

if (dx < dy)

{

pswap(U,V); lswap(dx,dy); swap(x,y);

if (both_even(dx, dy)) signh = -signh;

if (both_odd(dx, dy)) signh = -signh;

}

if (dy == 3)

if (dy == 0)

{

*V = gpowgs((GEN)y[2],dx-4);

*V = gpowgs((GEN)y[2],dx-1);

*U = gzero; return gmul(*V,(GEN)y[2]);

}

av = avma;

Line 2899 subresext(GEN x, GEN y, GEN *U, GEN *V)

Line 3318 subresext(GEN x, GEN y, GEN *U, GEN *V)

q = pseudodiv(u,v, &r); dr = lgef(r);

if (dr == 2) { *U = *V = gzero; avma = av; return gzero; }

du = lgef(u); dv = lgef(v); degq = du-dv;

du = degpol(u); dv = degpol(v); degq = du-dv;

/* lead(v)^(degq + 1) * um1 - q * uze */

p1 = gsub(gmul(gpowgs((GEN)v[dv-1],degq+1),um1), gmul(q,uze));

p1 = gsub(gmul(gpowgs((GEN)v[dv+2],degq+1),um1), gmul(q,uze));

um1 = uze; uze = p1;

u = v; p1 = g; g = leading_term(u);

switch(degq)

Line 2912 subresext(GEN x, GEN y, GEN *U, GEN *V)

Line 3331 subresext(GEN x, GEN y, GEN *U, GEN *V)

p1 = gmul(gpowgs(h,degq),p1);

h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1));

}

if (both_even(du, dv)) signh = -signh;

if (both_odd(du, dv)) signh = -signh;

v = gdivexact(r,p1);

uze= gdivexact(uze,p1);

if (dr==3) break;

if (low_stack(lim,stack_lim(av2,1)))

{

GEN *gptr[6]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;

gptr[4]=&uze; gptr[5]=&um1;

if(DEBUGMEM>1) err(warnmem,"subresext, dr = %ld",dr);

gerepilemany(av2,gptr,6);

gerepileall(av2,6, &u,&v,&g,&h,&uze,&um1);

}

z = (GEN)v[2];

if (dv > 4)

if (dv > 1)

{

/* z = gdivexact(gpowgs(z,dv-3), gpowgs(h,dv-4)); */

/* z = gdivexact(gpowgs(z,dv), gpowgs(h,dv-1)); */

p1 = gpowgs(gdiv(z,h),dv-4);

p1 = gpowgs(gdiv(z,h),dv-1);

z = gmul(z,p1); uze = gmul(uze,p1);

}

if (signh < 0) { z = gneg_i(z); uze = gneg_i(uze); }

p1 = gadd(z, gneg(gmul(uze,xprim)));

if (!poldivis(p1,yprim,&vze)) err(bugparier,"subresext");

vze = poldivres(p1, yprim, &r);

if (!gcmp0(r)) err(warner,"inexact computation in subresext");

/* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */

p2 = gun;

if (cu) p2 = gmul(p2, gpowgs(cu,dy-3));

if (cu) p2 = gmul(p2, gpowgs(cu,dy));

if (cv) p2 = gmul(p2, gpowgs(cv,dx-3));

if (cv) p2 = gmul(p2, gpowgs(cv,dx));

cu = cu? gdiv(p2,cu): p2;

cv = cv? gdiv(p2,cv): p2;

tetpil = avma;

z = gmul(z,p2); uze = gmul(uze,cu); vze = gmul(vze,cv);

z = gmul(z,p2);

{

uze = gmul(uze,cu);

GEN *gptr[3]; gptr[0]=&z; gptr[1]=&uze; gptr[2]=&vze;

vze = gmul(vze,cv);

gerepilemanysp(av,tetpil,gptr,3);

gptr[0]=&z; gptr[1]=&uze; gptr[2]=&vze;

}

gerepilemanysp(av,tetpil,gptr,3);

*U = uze; *V = vze; return z;

*U = uze;

*V = vze; return z;

}

static GEN

Line 2967 zero_bezout(GEN y, GEN *U, GEN *V)

Line 3386 zero_bezout(GEN y, GEN *U, GEN *V)

GEN

bezoutpol(GEN x, GEN y, GEN *U, GEN *V)

{

ulong av,av2,tetpil,lim;

gpmem_t av, av2, tetpil, lim;

long degq,tx,ty,dx,dy,du,dv,dr;

GEN q,z,g,h,r,p1,cu,cv,u,v,um1,uze,vze,xprim,yprim, *gptr[3];

Line 2984 bezoutpol(GEN x, GEN y, GEN *U, GEN *V)

Line 3403 bezoutpol(GEN x, GEN y, GEN *U, GEN *V)

if (varn(x)!=varn(y))

return (varn(x)<varn(y))? scalar_bezout(x,y,U,V)

: scalar_bezout(y,x,V,U);

dx = lgef(x); dy = lgef(y);

dx = degpol(x); dy = degpol(y);

if (dx < dy)

{

pswap(U,V); lswap(dx,dy); swap(x,y);

}

if (dy==3) return scalar_bezout(x,y,U,V);

if (dy==0) return scalar_bezout(x,y,U,V);

u = primitive_part(x, &cu); xprim = u;

v = primitive_part(y, &cv); yprim = v;

Line 3000 bezoutpol(GEN x, GEN y, GEN *U, GEN *V)

Line 3419 bezoutpol(GEN x, GEN y, GEN *U, GEN *V)

q = pseudodiv(u,v, &r); dr = lgef(r);

if (dr == 2) break;

du = lgef(u); dv = lgef(v); degq = du-dv;

du = degpol(u); dv = degpol(v); degq = du-dv;

p1 = gsub(gmul(gpowgs((GEN)v[dv-1],degq+1),um1), gmul(q,uze));

p1 = gsub(gmul(gpowgs((GEN)v[dv+2],degq+1),um1), gmul(q,uze));

um1 = uze; uze = p1;

u = v; p1 = g; g = leading_term(u);

switch(degq)

Line 3018 bezoutpol(GEN x, GEN y, GEN *U, GEN *V)

Line 3437 bezoutpol(GEN x, GEN y, GEN *U, GEN *V)

if (dr==3) break;

if (low_stack(lim,stack_lim(av2,1)))

{

GEN *gptr[6]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;

gptr[4]=&uze; gptr[5]=&um1;

if(DEBUGMEM>1) err(warnmem,"bezoutpol, dr = %ld",dr);

gerepilemany(av2,gptr,6);

gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);

}

if (!poldivis(gsub(v,gmul(uze,xprim)),yprim, &vze))

p1 = gsub(v,gmul(uze,xprim));

err(warner,"inexact computation in bezoutpol");

vze = poldivres(p1, yprim, &r);

if (!gcmp0(r)) err(warner,"inexact computation in bezoutpol");

if (cu) uze = gdiv(uze,cu);

if (cv) vze = gdiv(vze,cv);

p1 = ginv(content(v)); tetpil = avma;

p1 = ginv(content(v));

tetpil = avma;

uze = gmul(uze,p1);

vze = gmul(vze,p1);

z = gmul(v,p1);

gptr[0]=&uze; gptr[1]=&vze; gptr[2]=&z;

gerepilemanysp(av,tetpil,gptr,3);

*U = uze; *V = vze; return z;

*U = uze;

*V = vze; return z;

}

/*******************************************************************/

Line 3075 Lazard2(GEN F, GEN x, GEN y, long n)

Line 3495 Lazard2(GEN F, GEN x, GEN y, long n)

x = Lazard(x,y,n-1); return gdivexact(gmul(x,F),y);

}

extern GEN addshiftpol(GEN x, GEN y, long d);

#define addshift(x,y) addshiftpol((x),(y),1)

static GEN

nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)

{

GEN p0, q0, z0, H, A;

long p, q, j, av, lim, v = varn(P);

long p, q, j, v = varn(P);

gpmem_t av, lim;

z0 = leading_term(Z);

p = degpol(P); p0 = leading_term(P); P = reductum(P);

Line 3099 nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)

Line 3517 nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)

A = gadd(A,gmul((GEN)P[j+2],H));

if (low_stack(lim,stack_lim(av,1)))

{

GEN *gptr[2]; gptr[0]=&A; gptr[1]=&H;

if(DEBUGMEM>1) err(warnmem,"nextSousResultant j = %ld/%ld",j,p);

gerepilemany(av,gptr,2);

gerepileall(av,2,&A,&H);

}

P = normalizepol_i(P, q+2);

Line 3115 nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)

Line 3532 nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)

GEN

resultantducos(GEN P, GEN Q)

{

ulong av = avma, av2,lim;

gpmem_t av = avma, av2, lim;

long dP,dQ,delta;

GEN cP,cQ,Z,s;

Line 3134 resultantducos(GEN P, GEN Q)

Line 3551 resultantducos(GEN P, GEN Q)

if (degpol(Q) > 0)

{

av2 = avma; lim = stack_lim(av2,1);

s = gpuigs(leading_term(Q),delta);

s = gpowgs(leading_term(Q),delta);

Z = Q;

Q = pseudorem(P, gneg(Q));

P = Z;

Line 3142 resultantducos(GEN P, GEN Q)

Line 3559 resultantducos(GEN P, GEN Q)

{

if (low_stack(lim,stack_lim(av,1)))

{

GEN *gptr[2]; gptr[0]=&P; gptr[1]=&Q;

if(DEBUGMEM>1) err(warnmem,"resultantducos, degpol Q = %ld",degpol(Q));

gerepilemany(av2,gptr,2); s = leading_term(P);

gerepileall(av2,2,&P,&Q); s = leading_term(P);

}

delta = degpol(P) - degpol(Q);

Z = Lazard2(Q, leading_term(Q), s, delta);

Line 3215 sylvestermatrix(GEN x, GEN y)

Line 3631 sylvestermatrix(GEN x, GEN y)

GEN

resultant2(GEN x, GEN y)

{

long av;

gpmem_t av;

GEN r;

if ((r = init_resultant(x,y))) return r;

av = avma; return gerepileupto(av,det(sylvestermatrix_i(x,y)));

Line 3251 fix_pol(GEN x, long v, long *mx)

Line 3667 fix_pol(GEN x, long v, long *mx)

GEN

polresultant0(GEN x, GEN y, long v, long flag)

{

long av = avma, m = 0;

long m = 0;

gpmem_t av = avma;

if (v >= 0)

{

Line 3274 polresultant0(GEN x, GEN y, long v, long flag)

Line 3691 polresultant0(GEN x, GEN y, long v, long flag)

/* GCD USING SUBRESULTANT */

/* */

/*******************************************************************/

extern GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den);

extern GEN to_polmod(GEN x, GEN mod);

GEN

srgcd(GEN x, GEN y)

{

long av,av1,tetpil,tx=typ(x),ty=typ(y),dx,dy,vx,vmod,lim;

long tx=typ(x), ty=typ(y), dx, dy, vx, vmod;

gpmem_t av, av1, tetpil, lim;

GEN d,g,h,p1,p2,u,v,mod,cx,cy;

if (!signe(y)) return gcopy(x);

Line 3309 srgcd(GEN x, GEN y)

Line 3724 srgcd(GEN x, GEN y)

if (vmod < 0) return modulargcd(x,y); /* Q[X] */

}

if (issimplefield(x) || issimplefield(y)) x = polgcdnun(x,y);

if (issimplepol(x) || issimplepol(y)) x = polgcdnun(x,y);

else

{

dx=lgef(x); dy=lgef(y);

Line 3344 srgcd(GEN x, GEN y)

Line 3759 srgcd(GEN x, GEN y)

h = g = leading_term(u);

break;

default:

v = gdiv(r,gmul(gpuigs(h,degq),g));

v = gdiv(r,gmul(gpowgs(h,degq),g));

g = leading_term(u);

h = gdiv(gpuigs(g,degq), gpuigs(h,degq-1));

h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));

}

if (low_stack(lim, stack_lim(av1,1)))

{

GEN *gptr[4]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;

if(DEBUGMEM>1) err(warnmem,"srgcd");

gerepilemany(av1,gptr,4);

gerepileall(av1,4,&u,&v,&g,&h);

}

p1 = content(v); if (!gcmp1(p1)) v = gdiv(v,p1);

Line 3368 srgcd(GEN x, GEN y)

Line 3782 srgcd(GEN x, GEN y)

return gerepileupto(av,x);

}

extern GEN qf_disc(GEN x, GEN y, GEN z);

GEN

poldisc0(GEN x, long v)

{

long av,tx=typ(x),i;

long tx=typ(x), i;

gpmem_t av;

GEN z,p1,p2;

switch(tx)

Line 3415 discsr(GEN x)

Line 3828 discsr(GEN x)

GEN

reduceddiscsmith(GEN pol)

{

long av=avma,tetpil,i,j,n;

long i, j, n;

gpmem_t av=avma, tetpil;

GEN polp,alpha,p1,m;

if (typ(pol)!=t_POL) err(typeer,"reduceddiscsmith");

Line 3446 reduceddiscsmith(GEN pol)

Line 3860 reduceddiscsmith(GEN pol)

long

sturmpart(GEN x, GEN a, GEN b)

{

long av = avma, lim = stack_lim(av,1), sl,sr,s,t,r1;

long sl, sr, s, t, r1;

gpmem_t av = avma, lim = stack_lim(av, 1);

GEN g,h,u,v;

if (typ(x) != t_POL) err(typeer,"sturm");

Line 3508 sturmpart(GEN x, GEN a, GEN b)

Line 3923 sturmpart(GEN x, GEN a, GEN b)

case 1:

p1 = gmul(h,p1); h = g; break;

default:

p1 = gmul(gpuigs(h,degq),p1);

p1 = gmul(gpowgs(h,degq),p1);

h = gdivexact(gpuigs(g,degq), gpuigs(h,degq-1));

h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1));

}

v = gdivexact(r,p1);

if (low_stack(lim,stack_lim(av,1)))

{

GEN *gptr[4]; gptr[0]=&u; gptr[1]=&v; gptr[2]=&g; gptr[3]=&h;

if(DEBUGMEM>1) err(warnmem,"polsturm, dr = %ld",dr);

gerepilemany(av,gptr,4);

gerepileall(av,4,&u,&v,&g,&h);

}

Line 3530 sturmpart(GEN x, GEN a, GEN b)

Line 3944 sturmpart(GEN x, GEN a, GEN b)

GEN

quadpoly0(GEN x, long v)

{

long res,l,tetpil,i,sx, tx = typ(x);

long res, i, l, sx, tx = typ(x);

gpmem_t av;

GEN y,p1;

if (is_matvec_t(tx))

{

l=lg(x); y=cgetg(l,tx);

l = lg(x); y = cgetg(l,tx);

for (i=1; i<l; i++) y[i]=(long)quadpoly0((GEN)x[i],v);

for (i=1; i<l; i++) y[i] = (long)quadpoly0((GEN)x[i],v);

return y;

}

if (tx!=t_INT) err(arither1);

if (tx != t_INT) err(arither1);

if (v < 0) v = 0;

sx=signe(x);

sx = signe(x);

if (!sx) err(talker,"zero discriminant in quadpoly");

y=cgetg(5,t_POL);

res = mod4(x); if (sx < 0 && res) res = 4-res;

y[1]=evalsigne(1) | evalvarn(v) | evallgef(5); y[4]=un;

if (res > 1) err(funder2,"quadpoly");

res=mod4(x); if (sx<0 && res) res=4-res;

if (res>1) err(funder2,"quadpoly");

l=avma; p1=shifti(x,-2); setsigne(p1,-signe(p1));

y = cgetg(5,t_POL);

y[1] = evalsigne(1) | evalvarn(v) | evallgef(5);

av = avma; p1 = shifti(x,-2); setsigne(p1,-signe(p1));

y[2] = (long) p1;

if (!res) y[3] = zero;

else

{

if (sx<0) { tetpil=avma; y[2] = lpile(l,tetpil,addsi(1,p1)); }

if (sx < 0) y[2] = lpileuptoint(av, addsi(1,p1));

y[3] = lnegi(gun);

}

return y;

y[4] = un; return y;

}

GEN

Line 3609 newtonpoly(GEN x, GEN p)

Line 4025 newtonpoly(GEN x, GEN p)

{

GEN y;

long n,ind,a,b,c,u1,u2,r1,r2;

long *vval, num[] = {evaltyp(t_INT) | m_evallg(3), 0, 0};

long *vval, num[] = {evaltyp(t_INT) | _evallg(3), 0, 0};

if (typ(x)!=t_POL) err(typeer,"newtonpoly");

n=degpol(x); if (n<=0) { y=cgetg(1,t_VEC); return y; }

Line 3636 newtonpoly(GEN x, GEN p)

Line 4052 newtonpoly(GEN x, GEN p)

free(vval); return y;

}

extern int cmp_pol(GEN x, GEN y);

extern GEN ZY_ZXY_resultant(GEN A, GEN B0, long *lambda);

GEN matratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN denom);

GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den);

/* Factor polynomial a on the number field defined by polynomial t */

GEN

polfnf(GEN a, GEN t)

{

ulong av = avma;

gpmem_t av = avma;

GEN x0,y,p1,p2,u,fa,n,unt,dent,alift;

GEN x0,y,p1,p2,u,G,fa,n,unt,dent,alift;

long lx,i,k,e;

int sqfree, tmonic;

Line 3655 polfnf(GEN a, GEN t)

Line 4066 polfnf(GEN a, GEN t)

a = fix_relative_pol(t,a,0);

alift = lift(a);

p1 = content(alift); if (!gcmp1(p1)) { a = gdiv(a, p1); alift = lift(a); }

p1 = content(t); if (!gcmp1(t)) t = gdiv(t, p1);

t = primpart(t);

tmonic = is_pm1(leading_term(t));

dent = ZX_disc(t); unt = gmodulsg(1,t);

dent = indexpartial(t, NULL); unt = gmodulsg(1,t);

u = nfgcd(alift,derivpol(alift), t, dent);

G = nfgcd(alift,derivpol(alift), t, dent);

sqfree = gcmp1(u);

sqfree = gcmp1(G);

u = sqfree? alift: lift_intern(gdiv(a, gmul(unt,u)));

u = sqfree? alift: lift_intern(gdiv(a, gmul(unt,G)));

k = 0; n = ZY_ZXY_resultant(t, u, &k);

if (DEBUGLEVEL > 4) fprintferr("polfnf: choosing k = %ld\n",k);

if (DEBUGLEVEL>4) fprintferr("polfnf: choosing k = %ld\n",k);

if (!sqfree)

{

G = poleval(G, gadd(polx[varn(a)], gmulsg(k, polx[varn(t)])));

G = ZY_ZXY_resultant(t, G, NULL);

}

/* n guaranteed to be squarefree */

fa = squff2(n,0); lx = lg(fa);

fa = DDF2(n,0); lx = lg(fa);

y = cgetg(3,t_MAT);

p1 = cgetg(lx,t_COL); y[1] = (long)p1;

p2 = cgetg(lx,t_COL); y[2] = (long)p2;

if (lx == 2)

{ /* P^k, k irreducible */

p1[1] = lmul(unt,u);

p2[1] = lstoi(degpol(a) / degpol(u));

return gerepilecopy(av, y);

}

x0 = gadd(polx[varn(a)], gmulsg(-k, gmodulcp(polx[varn(t)], t)));

for (i=lx-1; i>1; i--)

for (i=lx-1; i>0; i--)

{

GEN b, F = lift_intern(poleval((GEN)fa[i], x0));

GEN f = (GEN)fa[i], F = lift_intern(poleval(f, x0));

if (!tmonic) F = gdiv(F, content(F));

if (!tmonic) F = primpart(F);

F = nfgcd(u, F, t, dent);

if (typ(F) != t_POL || degpol(F) == 0)

err(talker,"reducible modulus in factornf");

F = gmul(unt, F);

e = 1;

F = gdiv(F, leading_term(F));

if (!sqfree)

u = lift_intern(gdeuc(u, F));

u = gdiv(u, content(u));

if (sqfree) e = 1;

else

{

e = 0; while (poldivis(a,F, &b)) { a = b; e++; }

while (poldivis(G,f, &G)) e++;

if (degpol(a) == degpol(u)) sqfree = 1;

if (degpol(G) == 0) sqfree = 1;

}

p1[i] = (long)F;

p1[i] = ldiv(gmul(unt,F), leading_term(F));

p2[i] = lstoi(e);

}

u = gmul(unt, u);

u = gdiv(u, leading_term(u));

p1[1] = (long)u;

e = sqfree? 1: degpol(a)/degpol(u);

p2[1] = lstoi(e);

return gerepilecopy(av, sort_factor(y, cmp_pol));

}

extern GEN FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p);

extern GEN FpM(GEN z, GEN p);

extern GEN polpol_to_mat(GEN v, long n);

extern GEN mat_to_polpol(GEN x, long v, long w);

static GEN

to_frac(GEN a, GEN b)

{

Line 3736 lift_to_frac(GEN t, GEN mod, GEN amax, GEN bmax, GEN d

Line 4144 lift_to_frac(GEN t, GEN mod, GEN amax, GEN bmax, GEN d

GEN

matratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN denom)

{

ulong ltop = avma;

gpmem_t ltop = avma;

GEN N, a;

long l2, l3;

long i, j;

Line 3759 matratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN den

Line 4167 matratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN den

GEN

polratlift(GEN P, GEN mod, GEN amax, GEN bmax, GEN denom)

{

ulong ltop = avma;

gpmem_t ltop = avma;

GEN Q,a;

long l2, j;

if (typ(P)!=t_POL) err(typeer,"polratlift");

Line 3793 polratlift(GEN P, GEN mod, GEN amax, GEN bmax, GEN den

Line 4201 polratlift(GEN P, GEN mod, GEN amax, GEN bmax, GEN den

* If not NULL, den must a a multiple of the denominator of the gcd,

* for example the discriminant of nf.

* NOTE: if nf is not irreducible, nfgcd may loop forever, especially if the

* NOTE: if nf is not irreducible, nfgcd may loop forever, esp. if gcd | nf */

* gcd divides nf !

GEN

nfgcd(GEN P, GEN Q, GEN nf, GEN den)

{

ulong ltop = avma;

gpmem_t ltop = avma;

GEN sol, mod = NULL;

long x = varn(P);

long y = varn(nf);

Line 3815 nfgcd(GEN P, GEN Q, GEN nf, GEN den)

Line 4221 nfgcd(GEN P, GEN Q, GEN nf, GEN den)

den = mulii(den, mppgcd(ZX_resultant(lP, nf), ZX_resultant(lQ, nf)));

/*Modular GCD*/

{

ulong btop = avma;

gpmem_t btop = avma, st_lim = stack_lim(btop, 1);

ulong st_lim = stack_lim(btop, 1);

long p;

long dM=0, dR;

GEN M, dsol, dens;

GEN M, dsol;

GEN R, ax, bo;

byteptr primepointer;

for (p = 27449, primepointer = diffptr + 3000; ; p += *(primepointer++))

for (p = 27449, primepointer = diffptr + 3000; ; )

{

if (!*primepointer) err(primer1);

if (!smodis(den, p))

continue;/*Discard primes dividing the disc(T) or leadingcoeff(PQ) */

goto repeat;/*Discard primes dividing disc(T) or leadingcoeff(PQ) */

if (DEBUGLEVEL>=5) fprintferr("nfgcd: p=%d\n",p);

if (DEBUGLEVEL>5) fprintferr("nfgcd: p=%d\n",p);

if ((R = FpXQX_safegcd(P, Q, nf, stoi(p))) == NULL)

continue;/*Discard primes when modular gcd does not exist*/

goto repeat;/*Discard primes when modular gcd does not exist*/

dR = degpol(R);

if (dR == 0) return scalarpol(gun, x);

if (mod && dR > dM) continue; /* p divides Res(P/gcd, Q/gcd). Discard. */

if (mod && dR > dM) goto repeat; /* p divides Res(P/gcd, Q/gcd). Discard. */

R = polpol_to_mat(R, d);

/* previous primes divided Res(P/gcd, Q/gcd)? Discard them. */

if (!mod || dR < dM) { M = R; mod = stoi(p); dM = dR; continue; }

if (!mod || dR < dM) { M = R; mod = stoi(p); dM = dR; goto repeat; }

if (low_stack(st_lim, stack_lim(btop, 1)))

{

GEN *bptr[2]; bptr[0]=&M; bptr[1]=&mod;

if (DEBUGMEM>1) err(warnmem,"nfgcd");

gerepilemany(btop, bptr, 2);

gerepileall(btop, 2, &M, &mod);

}

ax = gmulgs(mpinvmod(stoi(p), mod), p);

Line 3851 nfgcd(GEN P, GEN Q, GEN nf, GEN den)

Line 4255 nfgcd(GEN P, GEN Q, GEN nf, GEN den)

/* I suspect it must be better to take amax > bmax*/

bo = racine(shifti(mod, -1));

if ((sol = matratlift(M, mod, bo, bo, den)) == NULL)

continue;

goto repeat;

dens = denom(sol);

sol = mat_to_polpol(sol,x,y);

dsol = gmul(sol, gmodulcp(dens, nf));

dsol = primpart(sol);

if (gdivise(P, dsol) && gdivise(Q, dsol))

if (gcmp0(pseudorem_i(P, dsol, nf))

break;

&& gcmp0(pseudorem_i(Q, dsol, nf))) break;

repeat:

NEXT_PRIME_VIADIFF_CHECK(p, primepointer);

}

return gerepilecopy(ltop, sol);

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>