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, |
|
Line 24 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
#define swap(a,b) { GEN _x = a; a = b; b = _x; } |
#define swap(a,b) { GEN _x = a; a = b; b = _x; } |
#define lswap(a,b) { long _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 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 |
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 (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; |
return y; |
} |
} |
|
|
|
GEN |
|
polsym(GEN x, long n) |
|
{ |
|
return polsym_gen(x, NULL, n, NULL,NULL); |
|
} |
|
|
static int (*polcmp_coeff_cmp)(GEN,GEN); |
static int (*polcmp_coeff_cmp)(GEN,GEN); |
|
|
/* assume x and y are polynomials in the same variable whose coeffs can be |
/* 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); |
long i, lx = lgef(x), ly = lgef(y); |
int fl; |
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; |
if (lx < ly) return -1; |
if (lx < ly) return -1; |
for (i=lx-1; i>1; i--) |
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; |
int (*old)(GEN,GEN) = polcmp_coeff_cmp; |
polcmp_coeff_cmp = cmp; |
polcmp_coeff_cmp = cmp; |
(void)sort_factor_gen(y,polcmp); |
(void)sort_factor_gen(y,&polcmp); |
polcmp_coeff_cmp = old; return y; |
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 |
GEN |
sort_factor_gen(GEN y, int (*cmp)(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; |
GEN a,b,A,B,w; |
a = (GEN)y[1]; n = lg(a); A = new_chunk(n); |
a = (GEN)y[1]; n = lg(a); A = new_chunk(n); |
b = (GEN)y[2]; B = 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; |
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 */ |
/* for internal use */ |
GEN |
GEN |
centermod_i(GEN x, GEN p, GEN ps2) |
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); |
if (!ps2) ps2 = shifti(p,-1); |
switch(typ(x)) |
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++) |
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; |
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; |
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; |
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 |
static GEN |
polidivis(GEN x, GEN y, GEN bound) |
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; |
GEN z,p1,y_lead; |
|
|
dy=degpol(y); |
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; |
if (gcmp1(y_lead)) y_lead = NULL; |
|
|
p1 = (GEN)x[dx]; |
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--) |
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++) |
for (j=i-dy+1; j<=i && j<=dz; j++) |
p1 = subii(p1, mulii((GEN)z[j],(GEN)y[i-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; |
if (absi_cmp(p1, bound) > 0) return NULL; |
p1 = gerepileupto(av,p1); |
p1 = gerepileupto(av,p1); |
z[i-dy] = (long)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; |
if (!i) break; |
} |
} |
z -= 2; |
z -= 2; |
z[1]=evalsigne(1) | evallgef(dz+3) | evalvarn(vx); |
z[1] = evalsigne(1) | evallgef(dz+3) | evalvarn(vx); |
return z; |
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) **/ |
/** QUADRATIC HENSEL LIFT (adapted from V. Shoup's NTL) **/ |
/** **/ |
/** **/ |
/***********************************************************************/ |
/***********************************************************************/ |
/* Setup for divide/conquer quadratic Hensel lift |
/* 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 |
* V = set of products of factors built as follows |
* 1) V[1..k] = initial a |
* 1) V[1..k] = initial a |
* 2) iterate: |
* 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] |
* link[i] = -j if V[i] = a[j] |
* j if V[i] = V[j] * V[j+1] |
* j if V[i] = V[j] * V[j+1] |
* Arrays (link, V, W) pre-allocated for 2k - 2 elements */ |
* Arrays (link, V, W) pre-allocated for 2k - 2 elements */ |
#if 0 /* restricted to Fp */ |
|
static void |
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) |
BuildTree(GEN link, GEN V, GEN W, GEN a, GEN T, GEN p) |
{ |
{ |
long k = lg(a)-1; |
long k = lg(a)-1; |
Line 340 BuildTree(GEN link, GEN V, GEN W, GEN a, GEN T, GEN p) |
|
Line 340 BuildTree(GEN link, GEN V, GEN W, GEN a, GEN T, GEN p) |
|
static void |
static void |
HenselLift(GEN V, GEN W, long j, GEN f, GEN T, GEN pd, GEN p0, int noinv) |
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); |
long space = lgef(f) * (lgefint(pd) + lgefint(p0) - 2); |
GEN a2,b2,g,z,s,t; |
GEN a2,b2,g,z,s,t; |
GEN a = (GEN)V[j], b = (GEN)V[j+1]; |
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, |
|
Line 351 HenselLift(GEN V, GEN W, long j, GEN f, GEN T, GEN pd, |
|
(void)new_chunk(space); /* HACK */ |
(void)new_chunk(space); /* HACK */ |
g = gadd(f, gneg_i(gmul(a,b))); |
g = gadd(f, gneg_i(gmul(a,b))); |
if (T) g = FpXQX_red(g, T, mulii(p0,pd)); |
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); |
z = FpXQX_mul(v,g, T,pd); |
t = FpXQX_divres(z,a, T,pd, &s); |
t = FpXQX_divres(z,a, T,pd, &s); |
t = gadd(gmul(u,g), gmul(t,b)); t = FpXQX_red(t, T, pd); |
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, |
|
Line 370 HenselLift(GEN V, GEN W, long j, GEN f, GEN T, GEN pd, |
|
g = gadd(gneg_i(g), gun); |
g = gadd(gneg_i(g), gun); |
|
|
if (T) g = FpXQX_red(g, T, mulii(p0,pd)); |
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); |
z = FpXQX_mul(v,g, T,pd); |
t = FpXQX_divres(z,a, T,pd, &s); |
t = FpXQX_divres(z,a, T,pd, &s); |
t = gadd(gmul(u,g), gmul(t,b)); t = FpXQX_red(t, T, pd); |
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 |
|
Line 422 MultiLift(GEN f, GEN a, GEN T, GEN p, long e0, int fla |
|
l = 0; E[++l] = e; |
l = 0; E[++l] = e; |
while (e > 1) { e = (e+1)/2; 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) |
if (flag != 2) |
{ |
{ |
|
|
polhensellift(GEN pol, GEN fct, GEN p, long exp) |
polhensellift(GEN pol, GEN fct, GEN p, long exp) |
{ |
{ |
GEN p1, p2; |
GEN p1, p2; |
long av = avma, i, j, l; |
long i, j, l; |
|
gpmem_t av = avma; |
|
|
/* we check the arguments */ |
/* we check the arguments */ |
if (typ(pol) != t_POL) |
if (typ(pol) != t_POL) |
err(talker, "not a polynomial in polhensellift"); |
err(talker, "not a polynomial in polhensellift"); |
if ((typ(fct) != t_COL && typ(fct) != t_VEC) || (lg(fct) < 3)) |
if ((typ(fct) != t_COL && typ(fct) != t_VEC) || (lg(fct) < 3)) |
err(talker, "not a factorization in polhensellift"); |
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"); |
err(talker, "not a prime number in polhensellift"); |
if (exp < 1) |
if (exp < 1) |
err(talker, "not a positive exponent in polhensellift"); |
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 (i = 1; i <= l; i++) |
for (j = 1; j < i; j++) |
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", |
err(talker, "polhensellift: factors %Z and %Z are not coprime", |
p1[i], p1[j]); |
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 |
static GEN |
two_factor_bound(GEN x) |
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; |
GEN *invbin, c, r = cgetr(3), z; |
|
|
x += 2; invbin = (GEN*)new_chunk(n+1); |
x += 2; invbin = (GEN*)new_chunk(n+1); |
Line 604 two_factor_bound(GEN x) |
|
Line 606 two_factor_bound(GEN x) |
|
} |
} |
#endif |
#endif |
|
|
|
/* A | S ==> |a_i| <= binom(d-1, i-1) || S ||_2 + binom(d-1, i) lc(S) */ |
static GEN |
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 */ |
/* all factors have coeffs less than the bound */ |
static GEN |
static GEN |
all_factor_bound(GEN x) |
all_factor_bound(GEN x) |
{ |
{ |
long i, n = lgef(x) - 3; |
long i, n = lgef(x) - 3; |
GEN t, z = gzero; |
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]); |
t = absi((GEN)x[n+2]); |
z = addii(t, addsi(1, racine(z))); |
z = addii(t, addsi(1, racine(z))); |
z = mulii(z, binome(stoi(n-1), n>>1)); |
z = mulii(z, binome(stoi(n-1), n>>1)); |
return shifti(mulii(t,z),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 |
/* Naive recombination of modular factors: combine up to maxK modular |
* factors, degree <= klim and divisible by hint |
* factors, degree <= klim and divisible by hint |
* |
* |
* target = polynomial we want to factor |
* target = polynomial we want to factor |
* famod = array of modular factors. Product should be congruent to |
* famod = array of modular factors. Product should be congruent to |
* target/lc(target) modulo p^a |
* 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 |
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 maxK, long klim,long hint) |
{ |
{ |
long K = 1, cnt = 1, i,j,k, curdeg, lfamod = lg(famod)-1; |
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 ind = cgetg(lfamod+1, t_VECSMALL); |
GEN degpol = cgetg(lfamod+1, t_VECSMALL); |
GEN degpol = cgetg(lfamod+1, t_VECSMALL); |
GEN degsofar = cgetg(lfamod+1, t_VECSMALL); |
GEN degsofar = cgetg(lfamod+1, t_VECSMALL); |
GEN listmod = cgetg(lfamod+1, t_COL); |
GEN listmod = cgetg(lfamod+1, t_COL); |
GEN fa = cgetg(lfamod+1, t_COL); |
GEN fa = cgetg(lfamod+1, t_COL); |
GEN res = cgetg(3, t_VEC); |
GEN res = cgetg(3, t_VEC); |
GEN bound = all_factor_bound(target); |
|
|
|
if (maxK < 0) maxK = lfamod-1; |
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); |
pa_bs2 = shifti(pa_b,-1); |
pb= gpowgs(p, b); pbs2 = shifti(pb,-1); |
pb= gpowgs(p, b); |
for (i=1; i <= lfamod; i++) |
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 */ |
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 */ |
* 1 <= ind[i] <= lfamod */ |
nextK: |
nextK: |
if (K > maxK || 2*K > lfamod) goto END; |
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)); |
fprintferr("\n### K = %d, %Z combinations\n", K,binome(stoi(lfamod), K)); |
setlg(ind, K+1); ind[1] = 1; |
setlg(ind, K+1); ind[1] = 1; |
|
Sbound = ((K+1)>>1); |
i = 1; curdeg = degpol[ind[1]]; |
i = 1; curdeg = degpol[ind[1]]; |
for(;;) |
for(;;) |
{ /* try all combinations of K factors */ |
{ /* try all combinations of K factors */ |
|
|
} |
} |
if (curdeg <= klim && curdeg % hint == 0) /* trial divide */ |
if (curdeg <= klim && curdeg % hint == 0) /* trial divide */ |
{ |
{ |
GEN y, q, cont, list; |
GEN y, q, list; |
ulong av; |
gpmem_t av; |
long t; |
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 (t > spa_bs2) t -= spa_b; |
if (labs(t) > ((K+1)>>1)) |
if (labs(t) > Sbound) |
{ |
{ |
if (DEBUGLEVEL>6) fprintferr("."); |
if (DEBUGLEVEL>6) fprintferr("."); |
goto NEXT; |
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 */ |
/* check trailing coeff */ |
y = lc; |
y = lc; |
for (i=1; i<=K; i++) |
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; |
avma = av; goto NEXT; |
} |
} |
y = lc; /* full computation */ |
y = lc; /* full computation */ |
for (i=1; i<=K; i++) |
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 */ |
/* 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; |
avma = av; goto NEXT; |
} |
} |
/* found a factor */ |
/* 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); |
list = cgetg(K+1, t_VEC); |
listmod[cnt] = (long)list; |
listmod[cnt] = (long)list; |
for (i=1; i<=K; i++) list[i] = famod[ind[i]]; |
for (i=1; i<=K; i++) list[i] = famod[ind[i]]; |
|
|
|
y = primpart(y); |
fa[cnt++] = (long)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++) |
for (i=j=k=1; i <= lfamod; i++) |
{ /* remove used factors */ |
{ /* remove used factors */ |
if (j <= K && i == ind[j]) j++; |
if (j <= K && i == ind[j]) j++; |
else |
else |
{ |
{ |
famod[k] = famod[i]; |
famod[k] = famod[i]; |
trace[k] = trace[i]; |
trace1[k] = trace1[i]; |
|
trace2[k] = trace2[i]; |
degpol[k] = degpol[i]; k++; |
degpol[k] = degpol[i]; k++; |
} |
} |
} |
} |
lfamod -= K; |
lfamod -= K; |
if (lfamod < 2*K) goto END; |
if (lfamod < 2*K) goto END; |
i = 1; curdeg = degpol[ind[1]]; |
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) |
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; |
continue; |
} |
} |
|
|
|
|
} |
} |
} |
} |
END: |
END: |
if (degpol(target) > 0) |
if (degpol(pol) > 0) |
{ /* leftover factor */ |
{ /* 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); |
setlg(famod, lfamod+1); |
listmod[cnt] = (long)dummycopy(famod); |
listmod[cnt] = (long)dummycopy(famod); |
fa[cnt++] = (long)target; |
fa[cnt++] = (long)pol; |
} |
} |
if (DEBUGLEVEL>6) fprintferr("\n"); |
if (DEBUGLEVEL>6) fprintferr("\n"); |
setlg(listmod, cnt); setlg(fa, cnt); |
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 |
long |
logint(GEN B, GEN y, GEN *ptq) |
logint(GEN B, GEN y, GEN *ptq) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
long e,i,fl; |
long e,i,fl; |
GEN q,pow2, r = y; |
GEN q,pow2, r = y; |
|
|
|
|
if (ptq) *ptq = gerepileuptoint(av, icopy(r)); else avma = av; |
if (ptq) *ptq = gerepileuptoint(av, icopy(r)); else avma = av; |
return e; |
return e; |
} |
} |
|
|
/* recombination of modular factors: van Hoeij's algorithm */ |
/* 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 |
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]); |
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++) |
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); |
x = y; y = shifti(y,1); |
} |
} |
for(;;) |
for(k=0; ; k++) |
{ |
{ |
z = shifti(addii(x,y), -1); |
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; |
y = z; |
else |
else |
x = z; |
x = z; |
Line 937 root_bound(GEN P) |
|
Line 963 root_bound(GEN P) |
|
return y; |
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 |
GEN |
special_pivot(GEN x) |
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]); |
long i,j, l = lg(H), h = lg(H[1]); |
for (i=1; i<h; i++) |
for (i=1; i<h; i++) |
{ |
{ |
Line 992 special_pivot(GEN x) |
|
Line 990 special_pivot(GEN x) |
|
return H; |
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 |
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 |
/* Recombination phase of Berlekamp-Zassenhaus algorithm using a variant of |
* van Hoeij's knapsack |
* van Hoeij's knapsack |
* |
* |
* P = monic squarefree in Z[X]. |
* P = squarefree in Z[X]. |
* famod = array of (lifted) modular factors mod p^a |
* 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) |
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); |
av = avma; lim = stack_lim(av, 1); |
TT = cgetg(n0+1, t_VEC); |
TT = cgetg(n0+1, t_VEC); |
T = cgetg(n0+1, t_MAT); |
Tra = cgetg(n0+1, t_MAT); |
for (i=1; i<=n0; i++) |
for (i=1; i<=n0; i++) |
{ |
{ |
TT[i] = 0; |
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); |
pa = gpowgs(p,a); |
famod = hensel_lift_fact(P,famod,NULL,p,pa,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++) |
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; |
TT[i] = (long)p2; |
p2 += 1+tmax; /* ignore traces number 0...tmax */ |
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; |
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); |
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++) |
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); |
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); |
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 */ |
/* use van Hoeij's knapsack algorithm */ |
static GEN |
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 */ |
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]; |
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) |
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 (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, <); |
} |
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 |
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 (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; |
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); |
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; |
d = 0; e = da; nfacp = 0; |
w = PolX; prime[2] = p; |
w = PolX; prime[2] = p; |
while (d < (e>>1)) |
while (d < (e>>1)) |
{ |
{ |
long lgg; |
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); |
g = u_FpX_gcd(z, u_FpX_sub(w, PolX, p), p); |
lgg = degpol(g); |
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); |
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", |
fprintferr("...tried prime %3ld (%-3ld factor%s). Time = %ld\n", |
p, nfacp, nfacp==1?"": "s", timer2()); |
p, nfacp, nfacp==1?"": "s", timer2()); |
if (min_deg(lbit-1,tabbit) > klim) { avma = av; return _col(a); } |
|
if (nfacp < nmax) |
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++; |
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) */ |
/* 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); |
long i, lx, tx = typ(x); |
GEN z; |
GEN z; |
if (is_scalar_t(tx)) return gcopy(x); |
if (is_scalar_t(tx)) return gcopy(x); |
|
if (d <= 0) err(talker,"need positive degree in gdeflate"); |
if (tx == t_POL) |
if (tx == t_POL) |
{ |
{ |
long vx = varn(x); |
long vx = varn(x); |
ulong av; |
gpmem_t av; |
if (vx < v) |
if (vx < v) |
{ |
{ |
lx = lgef(x); |
lx = lgef(x); |
Line 1499 polinflate(GEN x0, long d) |
|
Line 1556 polinflate(GEN x0, long d) |
|
return y; |
return y; |
} |
} |
|
|
|
/* Distinct Degree Factorization (deflating first) |
|
* Assume x squarefree, degree(x) > 0, x(0) != 0 */ |
GEN |
GEN |
squff2(GEN x, long hint) |
DDF2(GEN x, long hint) |
{ |
{ |
GEN L; |
GEN L; |
long m; |
long m; |
x = poldeflate(x, &m); |
x = poldeflate(x, &m); |
L = squff(x, hint); |
L = DDF(x, hint); |
if (m > 1) |
if (m > 1) |
{ |
{ |
GEN e, v, fa = decomp(stoi(m)); |
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); |
GEN L2 = cgetg(1,t_VEC); |
for (i=1; i < lg(L); i++) |
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; |
L = L2; |
} |
} |
} |
} |
return L; |
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 |
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 **/ |
/** 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); |
assign_or_fail((GEN)p1[1],p); |
t[3]=1; break; |
t[3]=1; break; |
case t_COMPLEX: |
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++) |
for (j=1; j<=2; j++) |
{ |
{ |
p2 = (GEN)p1[j]; |
p2 = (GEN)p1[j]; |
Line 1756 factor0(GEN x,long flag) |
|
Line 1828 factor0(GEN x,long flag) |
|
return NULL; /* not reached */ |
return NULL; /* not reached */ |
} |
} |
|
|
/* assume f and g coprime integer factorizations */ |
|
GEN |
GEN |
merge_factor_i(GEN f, GEN g) |
concat_factor(GEN f, GEN g) |
{ |
{ |
GEN h; |
GEN h; |
if (lg(f) == 1) return g; |
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 = cgetg(3,t_MAT); |
h[1] = (long)concatsp((GEN)f[1], (GEN)g[1]); |
h[1] = (long)concatsp((GEN)f[1], (GEN)g[1]); |
h[2] = (long)concatsp((GEN)f[2], (GEN)g[2]); |
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 |
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) |
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; |
GEN y,p,p1,p2,p3,p4,p5,pol; |
|
|
if (is_matvec_t(tx)) |
if (is_matvec_t(tx)) |
|
|
case t_INT: return decomp(x); |
case t_INT: return decomp(x); |
|
|
case t_FRACN: |
case t_FRACN: |
x = gred(x); /* fall through */ |
return gerepileupto(av, factor(gred(x))); |
case t_FRAC: |
case t_FRAC: |
p1 = decomp((GEN)x[1]); |
p1 = decomp((GEN)x[1]); |
p2 = decomp((GEN)x[2]); p2[2] = (long)gneg_i((GEN)p2[2]); |
p2 = decomp((GEN)x[2]); p2[2] = (long)gneg_i((GEN)p2[2]); |
|
|
|
|
case t_COMPLEX: y=cgetg(3,t_MAT); lx=lgef(x)-2; v=varn(x); |
case t_COMPLEX: y=cgetg(3,t_MAT); lx=lgef(x)-2; v=varn(x); |
av = avma; p1 = roots(x,pa); tetpil = avma; |
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); |
y[1]=lpile(av,tetpil,p2); |
p3=cgetg(lx,t_COL); for (i=1; i<lx; i++) p3[i] = un; |
p3=cgetg(lx,t_COL); for (i=1; i<lx; i++) p3[i] = un; |
y[2]=(long)p3; return y; |
y[2]=(long)p3; return y; |
|
|
p2[2] = (long)deg1pol_i((GEN)p1[2], (GEN)p1[1], v); |
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()); |
if (killv) setvarn(pol, fetch_var()); |
switch (typ2(tx)) |
switch (typ2(tx)) |
{ |
{ |
|
|
switch (typ1(tx)) |
switch (typ1(tx)) |
{ |
{ |
case t_POLMOD: |
case t_POLMOD: |
if (killv) delete_var(); |
if (killv) (void)delete_var(); |
return gerepileupto(av,p1); |
return gerepileupto(av,p1); |
case t_COMPLEX: p5 = gi; break; |
case t_COMPLEX: p5 = gi; break; |
case t_QUAD: p5=cgetg(4,t_QUAD); |
case t_QUAD: p5=cgetg(4,t_QUAD); |
|
|
if(typ(p4)==t_POLMOD) p3[j]=lsubst((GEN)p4[2],v,p5); |
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); |
tetpil=avma; y=cgetg(3,t_MAT); |
y[1]=lcopy(p2);y[2]=lcopy((GEN)p1[2]); |
y[1]=lcopy(p2);y[2]=lcopy((GEN)p1[2]); |
return gerepile(av,tetpil,y); |
return gerepile(av,tetpil,y); |
|
|
} |
} |
|
|
case t_RFRACN: |
case t_RFRACN: |
x=gred_rfrac(x); /* fall through */ |
return gerepileupto(av, factor(gred_rfrac(x))); |
case t_RFRAC: |
case t_RFRAC: |
p1=factor((GEN)x[1]); |
p1=factor((GEN)x[1]); |
p2=factor((GEN)x[2]); p3=gneg_i((GEN)p2[2]); |
p2=factor((GEN)x[2]); p3=gneg_i((GEN)p2[2]); |
|
|
#undef typs |
#undef typs |
#undef tsh |
#undef tsh |
|
|
|
/* assume n != 0, t_INT. Compute x^|n| using left-right binary powering */ |
GEN |
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)) |
divide_conquer_prod(GEN x, GEN (*mul)(GEN,GEN)) |
{ |
{ |
long i,k,lx = lg(x); |
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 static_nf; |
|
|
static GEN |
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 |
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 |
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 |
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 |
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); |
x = cgetg(lx,t_VEC); |
for (l=1,k=1; k<lx; k++) |
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); |
setlg(x,l); |
return gerepileupto(av, divide_conquer_prod(x, _mul)); |
return gerepileupto(av, divide_conquer_prod(x, _mul)); |
} |
} |
|
|
GEN |
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) |
factorback(GEN fa, GEN nf) |
{ |
{ |
return factorback_i(fa,nf,0); |
return factorback_i(fa,nf,NULL,0); |
} |
} |
|
|
GEN |
GEN |
gisirreducible(GEN x) |
gisirreducible(GEN x) |
{ |
{ |
long av=avma, tx = typ(x),l,i; |
long tx = typ(x), l, i; |
|
gpmem_t av=avma; |
GEN y; |
GEN y; |
|
|
if (is_matvec_t(tx)) |
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 |
static GEN |
triv_cont_gcd(GEN x, GEN y) |
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]) |
GEN p1 = (typ(x)==t_COMPLEX)? ggcd((GEN)x[1],(GEN)x[2]) |
: ggcd((GEN)x[2],(GEN)x[3]); |
: ggcd((GEN)x[2],(GEN)x[3]); |
tetpil=avma; return gerepile(av,tetpil,ggcd(p1,y)); |
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 |
static GEN |
cont_gcd(GEN x, GEN y) |
cont_gcd(GEN x, GEN y) |
{ |
{ |
long av = avma, tetpil; |
gpmem_t av = avma, tetpil; |
GEN p1 = content(x); |
GEN p1 = content(x); |
|
|
tetpil=avma; return gerepile(av,tetpil,ggcd(p1,y)); |
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); |
long v = ggval(x,(GEN)y[2]), w = valp(y); |
if (w < v) v = w; |
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)\ |
#define fix_frac(z) if (signe(z[2])<0)\ |
{\ |
{\ |
setsigne(z[1],-signe(z[1]));\ |
setsigne(z[1],-signe(z[1]));\ |
setsigne(z[2],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 |
GEN |
ggcd(GEN x, GEN y) |
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; |
GEN p1,z; |
|
|
if (tx>ty) { swap(x,y); lswap(tx,ty); } |
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]); |
for (i=1; i<l; i++) z[i]=lgcd(x,(GEN)y[i]); |
return z; |
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 (is_const_t(tx)) |
{ |
{ |
if (ty == t_FRACN) |
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); |
case t_INTMOD: z=cgetg(3,t_INTMOD); |
if (egalii((GEN)x[1],(GEN)y[1])) |
if (egalii((GEN)x[1],(GEN)y[1])) |
{ copyifstack(x[1],z[1]); } |
copyifstack(x[1],z[1]); |
else |
else |
z[1] = lmppgcd((GEN)x[1],(GEN)y[1]); |
z[1] = lmppgcd((GEN)x[1],(GEN)y[1]); |
if (gcmp1((GEN)z[1])) z[2] = zero; |
if (gcmp1((GEN)z[1])) z[2] = zero; |
Line 2118 ggcd(GEN x, GEN y) |
|
Line 2351 ggcd(GEN x, GEN y) |
|
return z; |
return z; |
|
|
case t_COMPLEX: |
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); |
return triv_cont_gcd(y,x); |
|
|
case t_PADIC: |
case t_PADIC: |
if (!egalii((GEN)x[2],(GEN)y[2])) return gun; |
if (!egalii((GEN)x[2],(GEN)y[2])) return gun; |
vx = valp(x); |
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: |
case t_QUAD: |
av=avma; p1=gdiv(x,y); |
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[1] = lmppgcd(x,(GEN)y[1]); |
z[2] = licopy((GEN)y[2]); return z; |
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); |
return triv_cont_gcd(y,x); |
|
|
case t_PADIC: |
case t_PADIC: |
Line 2185 ggcd(GEN x, GEN y) |
|
Line 2405 ggcd(GEN x, GEN y) |
|
{ |
{ |
case t_FRAC: |
case t_FRAC: |
av = avma; p1=mppgcd((GEN)x[1],(GEN)y[2]); avma = av; |
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); |
return ggcd((GEN)y[1], x); |
|
|
case t_COMPLEX: case t_QUAD: |
case t_COMPLEX: case t_QUAD: |
Line 2198 ggcd(GEN x, GEN y) |
|
Line 2418 ggcd(GEN x, GEN y) |
|
case t_FRAC: |
case t_FRAC: |
switch(ty) |
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); |
return triv_cont_gcd(y,x); |
|
|
case t_PADIC: |
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); |
case t_POLMOD: z=cgetg(3,t_POLMOD); |
if (gegal((GEN)x[1],(GEN)y[1])) |
if (gegal((GEN)x[1],(GEN)y[1])) |
{ copyifstack(x[1],z[1]); } |
copyifstack(x[1],z[1]); |
else |
else |
z[1] = lgcd((GEN)x[1],(GEN)y[1]); |
z[1] = lgcd((GEN)x[1],(GEN)y[1]); |
if (lgef(z[1])<=3) z[2]=zero; |
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: |
case t_RFRAC: |
av = avma; p1=ggcd((GEN)x[1],(GEN)y[2]); avma = av; |
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); |
return ggcd((GEN)y[1],x); |
|
|
case t_RFRACN: |
case t_RFRACN: |
Line 2267 ggcd(GEN x, GEN y) |
|
Line 2489 ggcd(GEN x, GEN y) |
|
return srgcd(x,y); |
return srgcd(x,y); |
|
|
case t_SER: |
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); |
case t_RFRAC: case t_RFRACN: av=avma; z=cgetg(3,ty); |
z[1]=lgcd(x,(GEN)y[1]); |
z[1]=lgcd(x,(GEN)y[1]); |
Line 2279 ggcd(GEN x, GEN y) |
|
Line 2501 ggcd(GEN x, GEN y) |
|
switch(ty) |
switch(ty) |
{ |
{ |
case t_SER: |
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: |
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; |
break; |
|
|
case t_RFRAC: case t_RFRACN: z=cgetg(3,ty); |
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])); |
p1 = gdiv((GEN)y[2], ggcd((GEN)x[2], (GEN)y[2])); |
tetpil = avma; |
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; |
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 */ |
return NULL; /* not reached */ |
} |
} |
|
|
Line 2305 GEN glcm0(GEN x, GEN y) |
|
Line 2527 GEN glcm0(GEN x, GEN y) |
|
GEN |
GEN |
glcm(GEN x, GEN y) |
glcm(GEN x, GEN y) |
{ |
{ |
long av,tx,ty,i,l; |
long tx, ty, i, l; |
|
gpmem_t av; |
GEN p1,p2,z; |
GEN p1,p2,z; |
|
|
ty = typ(y); |
ty = typ(y); |
Line 2340 glcm(GEN x, GEN y) |
|
Line 2563 glcm(GEN x, GEN y) |
|
return gerepileupto(av,p2); |
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 |
static GEN |
polgcdnun(GEN x, GEN y) |
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(;;) |
for(;;) |
{ |
{ |
av1 = avma; p1 = gres(x,y); |
av1 = avma; r = gres(x,y); |
if (gcmp0(p1)) |
if (pol_approx0(r, x, exact)) |
{ |
{ |
avma = av1; |
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); |
return (y==yorig)? gcopy(y): gerepileupto(av,y); |
} |
} |
x = y; y = p1; |
x = y; y = r; |
if (low_stack(lim,stack_lim(av,1))) |
if (low_stack(lim,stack_lim(av,1))) |
{ |
{ |
GEN *gptr[2]; x = gcopy(x); gptr[0]=&x; gptr[1]=&y; |
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 */ |
/* return 1 if coeff explosion is not possible */ |
static int |
static int |
issimplefield(GEN x) |
issimplefield(GEN x) |
{ |
{ |
long lx,i; |
|
switch(typ(x)) |
switch(typ(x)) |
{ |
{ |
case t_REAL: case t_INTMOD: case t_PADIC: case t_SER: |
case t_REAL: case t_INTMOD: case t_PADIC: case t_SER: |
return 1; |
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]); |
return issimplefield((GEN)x[1]) || issimplefield((GEN)x[2]); |
|
case t_POLMOD: |
|
return issimplepol((GEN)x[1]) || issimplepol((GEN)x[2]); |
} |
} |
return 0; |
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 |
static int |
can_use_modular_gcd(GEN x, GEN *mod, long *v) |
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; |
break; |
case t_POL: |
case t_POL: |
if (!isrational(p1)) return 0; |
if (!isrational(p1)) return 0; |
if (*v >= 0) |
if (*v >= 0) |
{ |
{ |
if (*v != varn(p1)) return 0; |
if (*v != varn(p1)) return 0; |
} else *v = varn(p1); |
} else *v = varn(p1); |
Line 2456 isinexactfield(GEN x) |
|
Line 2686 isinexactfield(GEN x) |
|
static GEN |
static GEN |
gcdmonome(GEN x, GEN y) |
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; |
GEN p1,p2; |
|
|
if (dx < e) e = dx; |
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)); |
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 |
static GEN |
polinvinexact(GEN x, GEN y) |
polinvinexact(GEN x, GEN y) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
long i,dx=degpol(x),dy=degpol(y),lz=dx+dy; |
long i,dx=degpol(x),dy=degpol(y),lz=dx+dy; |
GEN v,z; |
GEN v,z; |
|
|
Line 2490 polinvinexact(GEN x, GEN y) |
|
Line 2721 polinvinexact(GEN x, GEN y) |
|
static GEN |
static GEN |
polinvmod(GEN x, GEN y) |
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; |
GEN u,v,d,p1; |
|
|
while (vx!=vy) |
while (vx!=vy) |
|
|
content(GEN x) |
content(GEN x) |
{ |
{ |
long lx,i,tx=typ(x); |
long lx,i,tx=typ(x); |
ulong av,tetpil; |
gpmem_t av, tetpil; |
GEN p1,p2; |
GEN p1,p2; |
|
|
if (is_scalar_t(tx)) |
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; |
av = avma; |
switch(tx) |
switch(tx) |
{ |
{ |
|
|
} |
} |
|
|
GEN |
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; |
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 */ |
/* SUBRESULTANT */ |
Line 2647 gdivexact(GEN x, GEN y) |
|
Line 3043 gdivexact(GEN x, GEN y) |
|
case t_INTMOD: |
case t_INTMOD: |
case t_POLMOD: return gmul(x,ginv(y)); |
case t_POLMOD: return gmul(x,ginv(y)); |
case t_POL: |
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); |
lx = lgef(x); z = cgetg(lx,t_POL); |
for (i=2; i<lx; i++) z[i]=(long)gdivexact((GEN)x[i],y); |
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); |
tx = typ(x); ty = typ(y); |
if (is_scalar_t(tx) || is_scalar_t(ty)) |
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; |
return gun; |
} |
} |
if (tx!=t_POL || ty!=t_POL) err(typeer,"subresall"); |
if (tx!=t_POL || ty!=t_POL) err(typeer,"subresall"); |
if (varn(x)==varn(y)) return NULL; |
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 */ |
/* return coefficients s.t x = x_0 X^n + ... + x_n */ |
static GEN |
GEN |
revpol(GEN x) |
revpol(GEN x) |
{ |
{ |
long i,n = degpol(x); |
long i,n = degpol(x); |
|
|
return y; |
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)"); |
if (!signe(y)) err(talker,"euclidean division by zero (pseudorem)"); |
(void)new_chunk(2); |
(void)new_chunk(2); |
dx=degpol(x); x = revpol(x); |
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--; |
x[0] = lneg((GEN)x[0]); p--; |
for (i=1; i<=dy; i++) |
for (i=1; i<=dy; i++) |
|
{ |
x[i] = ladd(gmul((GEN)y[0], (GEN)x[i]), gmul((GEN)x[0],(GEN)y[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++) |
for ( ; i<=dx; i++) |
|
{ |
x[i] = lmul((GEN)y[0], (GEN)x[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])); |
do { x++; dx--; } while (dx >= 0 && gcmp0((GEN)x[0])); |
if (dx < dy) break; |
if (dx < dy) break; |
if (low_stack(lim,stack_lim(av2,1))) |
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[0]=evaltyp(t_POL) | evallg(lx); |
x[1]=evalsigne(1) | evalvarn(vx) | evallgef(lx); |
x[1]=evalsigne(1) | evalvarn(vx) | evallgef(lx); |
x = revpol(x) - 2; |
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 |
/* assume dx >= dy, y non constant |
* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */ |
* Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */ |
GEN |
GEN |
pseudodiv(GEN x, GEN y, GEN *ptr) |
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; |
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); |
(void)new_chunk(2); |
dx=degpol(x); x = revpol(x); |
dx=degpol(x); x = revpol(x); |
dy=degpol(y); y = revpol(y); dz=dx-dy; p = dz+1; |
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])); |
x[i] = ladd(gmul((GEN)y[0], (GEN)x[i]), gmul((GEN)x[0],(GEN)y[i])); |
for ( ; i<=dx; i++) |
for ( ; i<=dx; i++) |
x[i] = lmul((GEN)y[0], (GEN)x[i]); |
x[i] = lmul((GEN)y[0], (GEN)x[i]); |
x++; dx--; |
x++; dx--; |
while (dx >= dy && gcmp0((GEN)x[0])) { x++; dx--; z[iz++] = zero; } |
while (dx >= dy && gcmp0((GEN)x[0])) { x++; dx--; z[iz++] = zero; } |
if (dx < dy) break; |
if (dx < dy) break; |
if (low_stack(lim,stack_lim(av2,1))) |
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--; } |
while (dx >= 0 && gcmp0((GEN)x[0])) { x++; dx--; } |
if (dx < 0) |
if (dx < 0) |
x = zeropol(vx); |
x = zeropol(vx); |
else |
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[1] = evalsigne(1) | evalvarn(vx) | evallgef(lz); |
z = revpol(z) - 2; |
z = revpol(z) - 2; |
r = gmul(x, (GEN)ypow[p]); |
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; |
*ptr = r; return z; |
} |
} |
|
|
/* Return resultant(u,v). If sol != NULL: set *sol to the last non-zero |
/* 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 |
GEN |
subresall(GEN u, GEN v, GEN *sol) |
subresall(GEN u, GEN v, GEN *sol) |
{ |
{ |
ulong av,av2,lim; |
gpmem_t av, av2, lim; |
long degq,dx,dy,du,dv,dr,signh; |
long degq,dx,dy,du,dv,dr,signh; |
GEN z,g,h,r,p1,p2,cu,cv; |
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 ((r = init_resultant(u,v))) return r; |
|
|
if (isinexactfield(u) || isinexactfield(v)) return resultant2(u,v); |
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); |
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; |
av = avma; |
u = primitive_part(u, &cu); |
u = primitive_part(u, &cu); |
v = primitive_part(v, &cv); |
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); |
r = pseudorem(u,v); dr = lgef(r); |
if (dr == 2) |
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; |
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); |
u = v; p1 = g; g = leading_term(u); |
switch(degq) |
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); |
p1 = gmul(gpowgs(h,degq),p1); |
h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1)); |
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); |
v = gdivexact(r,p1); |
if (dr==3) break; |
if (dr==3) break; |
if (low_stack(lim,stack_lim(av2,1))) |
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); |
if(DEBUGMEM>1) err(warnmem,"subresall, dr = %ld",dr); |
gerepilemany(av2,gptr,4); |
gerepileall(av2,4, &u, &v, &g, &h); |
} |
} |
} |
} |
z = (GEN)v[2]; |
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)) */ |
if (signh < 0) z = gneg(z); /* z = resultant(ppart(x), ppart(y)) */ |
p2 = gun; |
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); |
z = gmul(z, p2); |
|
|
if (sol) u = gclone(u); |
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 |
static GEN |
scalar_res(GEN x, GEN y, GEN *U, GEN *V) |
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) */ |
/* compute U, V s.t Ux + Vy = resultant(x,y) */ |
GEN |
GEN |
subresext(GEN x, GEN y, GEN *U, GEN *V) |
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; |
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; } |
if (gcmp0(x) || gcmp0(y)) { *U = *V = gzero; return gzero; } |
tx=typ(x); ty=typ(y); |
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)) |
if (varn(x) != varn(y)) |
return (varn(x)<varn(y))? scalar_res(x,y,U,V) |
return (varn(x)<varn(y))? scalar_res(x,y,U,V) |
: scalar_res(y,x,V,U); |
: 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) |
if (dx < dy) |
{ |
{ |
pswap(U,V); lswap(dx,dy); swap(x,y); |
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]); |
*U = gzero; return gmul(*V,(GEN)y[2]); |
} |
} |
av = avma; |
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); |
q = pseudodiv(u,v, &r); dr = lgef(r); |
if (dr == 2) { *U = *V = gzero; avma = av; return gzero; } |
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 */ |
/* 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; |
um1 = uze; uze = p1; |
u = v; p1 = g; g = leading_term(u); |
u = v; p1 = g; g = leading_term(u); |
switch(degq) |
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); |
p1 = gmul(gpowgs(h,degq),p1); |
h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1)); |
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); |
v = gdivexact(r,p1); |
uze= gdivexact(uze,p1); |
uze= gdivexact(uze,p1); |
if (dr==3) break; |
if (dr==3) break; |
if (low_stack(lim,stack_lim(av2,1))) |
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); |
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]; |
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); |
z = gmul(z,p1); uze = gmul(uze,p1); |
} |
} |
if (signh < 0) { z = gneg_i(z); uze = gneg_i(uze); } |
if (signh < 0) { z = gneg_i(z); uze = gneg_i(uze); } |
p1 = gadd(z, gneg(gmul(uze,xprim))); |
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)), */ |
/* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */ |
p2 = gun; |
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; |
cu = cu? gdiv(p2,cu): p2; |
cv = cv? gdiv(p2,cv): p2; |
cv = cv? gdiv(p2,cv): p2; |
|
|
tetpil = avma; |
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 |
static GEN |
Line 2967 zero_bezout(GEN y, GEN *U, GEN *V) |
|
Line 3386 zero_bezout(GEN y, GEN *U, GEN *V) |
|
GEN |
GEN |
bezoutpol(GEN x, GEN y, GEN *U, GEN *V) |
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; |
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]; |
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)) |
if (varn(x)!=varn(y)) |
return (varn(x)<varn(y))? scalar_bezout(x,y,U,V) |
return (varn(x)<varn(y))? scalar_bezout(x,y,U,V) |
: scalar_bezout(y,x,V,U); |
: scalar_bezout(y,x,V,U); |
dx = lgef(x); dy = lgef(y); |
dx = degpol(x); dy = degpol(y); |
if (dx < dy) |
if (dx < dy) |
{ |
{ |
pswap(U,V); lswap(dx,dy); swap(x,y); |
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; |
u = primitive_part(x, &cu); xprim = u; |
v = primitive_part(y, &cv); yprim = v; |
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); |
q = pseudodiv(u,v, &r); dr = lgef(r); |
if (dr == 2) break; |
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; |
um1 = uze; uze = p1; |
u = v; p1 = g; g = leading_term(u); |
u = v; p1 = g; g = leading_term(u); |
switch(degq) |
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 (dr==3) break; |
if (low_stack(lim,stack_lim(av2,1))) |
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); |
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 (cu) uze = gdiv(uze,cu); |
if (cv) vze = gdiv(vze,cv); |
if (cv) vze = gdiv(vze,cv); |
p1 = ginv(content(v)); tetpil = avma; |
p1 = ginv(content(v)); |
|
|
|
tetpil = avma; |
uze = gmul(uze,p1); |
uze = gmul(uze,p1); |
vze = gmul(vze,p1); |
vze = gmul(vze,p1); |
z = gmul(v,p1); |
z = gmul(v,p1); |
gptr[0]=&uze; gptr[1]=&vze; gptr[2]=&z; |
gptr[0]=&uze; gptr[1]=&vze; gptr[2]=&z; |
gerepilemanysp(av,tetpil,gptr,3); |
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); |
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 |
static GEN |
nextSousResultant(GEN P, GEN Q, GEN Z, GEN s) |
nextSousResultant(GEN P, GEN Q, GEN Z, GEN s) |
{ |
{ |
GEN p0, q0, z0, H, A; |
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); |
z0 = leading_term(Z); |
p = degpol(P); p0 = leading_term(P); P = reductum(P); |
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)); |
A = gadd(A,gmul((GEN)P[j+2],H)); |
if (low_stack(lim,stack_lim(av,1))) |
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); |
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); |
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 |
GEN |
resultantducos(GEN P, GEN Q) |
resultantducos(GEN P, GEN Q) |
{ |
{ |
ulong av = avma, av2,lim; |
gpmem_t av = avma, av2, lim; |
long dP,dQ,delta; |
long dP,dQ,delta; |
GEN cP,cQ,Z,s; |
GEN cP,cQ,Z,s; |
|
|
Line 3134 resultantducos(GEN P, GEN Q) |
|
Line 3551 resultantducos(GEN P, GEN Q) |
|
if (degpol(Q) > 0) |
if (degpol(Q) > 0) |
{ |
{ |
av2 = avma; lim = stack_lim(av2,1); |
av2 = avma; lim = stack_lim(av2,1); |
s = gpuigs(leading_term(Q),delta); |
s = gpowgs(leading_term(Q),delta); |
Z = Q; |
Z = Q; |
Q = pseudorem(P, gneg(Q)); |
Q = pseudorem(P, gneg(Q)); |
P = Z; |
P = Z; |
Line 3142 resultantducos(GEN P, GEN Q) |
|
Line 3559 resultantducos(GEN P, GEN Q) |
|
{ |
{ |
if (low_stack(lim,stack_lim(av,1))) |
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)); |
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); |
delta = degpol(P) - degpol(Q); |
Z = Lazard2(Q, leading_term(Q), s, delta); |
Z = Lazard2(Q, leading_term(Q), s, delta); |
Line 3215 sylvestermatrix(GEN x, GEN y) |
|
Line 3631 sylvestermatrix(GEN x, GEN y) |
|
GEN |
GEN |
resultant2(GEN x, GEN y) |
resultant2(GEN x, GEN y) |
{ |
{ |
long av; |
gpmem_t av; |
GEN r; |
GEN r; |
if ((r = init_resultant(x,y))) return r; |
if ((r = init_resultant(x,y))) return r; |
av = avma; return gerepileupto(av,det(sylvestermatrix_i(x,y))); |
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 |
GEN |
polresultant0(GEN x, GEN y, long v, long flag) |
polresultant0(GEN x, GEN y, long v, long flag) |
{ |
{ |
long av = avma, m = 0; |
long m = 0; |
|
gpmem_t av = avma; |
|
|
if (v >= 0) |
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 */ |
/* GCD USING SUBRESULTANT */ |
/* */ |
/* */ |
/*******************************************************************/ |
/*******************************************************************/ |
extern GEN nfgcd(GEN P, GEN Q, GEN nf, GEN den); |
|
extern GEN to_polmod(GEN x, GEN mod); |
|
|
|
GEN |
GEN |
srgcd(GEN x, GEN y) |
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; |
GEN d,g,h,p1,p2,u,v,mod,cx,cy; |
|
|
if (!signe(y)) return gcopy(x); |
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 (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 |
else |
{ |
{ |
dx=lgef(x); dy=lgef(y); |
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); |
h = g = leading_term(u); |
break; |
break; |
default: |
default: |
v = gdiv(r,gmul(gpuigs(h,degq),g)); |
v = gdiv(r,gmul(gpowgs(h,degq),g)); |
g = leading_term(u); |
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))) |
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"); |
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); |
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); |
return gerepileupto(av,x); |
} |
} |
|
|
extern GEN qf_disc(GEN x, GEN y, GEN z); |
|
|
|
GEN |
GEN |
poldisc0(GEN x, long v) |
poldisc0(GEN x, long v) |
{ |
{ |
long av,tx=typ(x),i; |
long tx=typ(x), i; |
|
gpmem_t av; |
GEN z,p1,p2; |
GEN z,p1,p2; |
|
|
switch(tx) |
switch(tx) |
|
|
GEN |
GEN |
reduceddiscsmith(GEN pol) |
reduceddiscsmith(GEN pol) |
{ |
{ |
long av=avma,tetpil,i,j,n; |
long i, j, n; |
|
gpmem_t av=avma, tetpil; |
GEN polp,alpha,p1,m; |
GEN polp,alpha,p1,m; |
|
|
if (typ(pol)!=t_POL) err(typeer,"reduceddiscsmith"); |
if (typ(pol)!=t_POL) err(typeer,"reduceddiscsmith"); |
Line 3446 reduceddiscsmith(GEN pol) |
|
Line 3860 reduceddiscsmith(GEN pol) |
|
long |
long |
sturmpart(GEN x, GEN a, GEN b) |
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; |
GEN g,h,u,v; |
|
|
if (typ(x) != t_POL) err(typeer,"sturm"); |
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: |
case 1: |
p1 = gmul(h,p1); h = g; break; |
p1 = gmul(h,p1); h = g; break; |
default: |
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); |
v = gdivexact(r,p1); |
if (low_stack(lim,stack_lim(av,1))) |
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); |
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 |
GEN |
quadpoly0(GEN x, long v) |
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; |
GEN y,p1; |
|
|
if (is_matvec_t(tx)) |
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; |
return y; |
} |
} |
if (tx!=t_INT) err(arither1); |
if (tx != t_INT) err(arither1); |
if (v < 0) v = 0; |
if (v < 0) v = 0; |
sx=signe(x); |
sx = signe(x); |
if (!sx) err(talker,"zero discriminant in quadpoly"); |
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; |
y[2] = (long) p1; |
if (!res) y[3] = zero; |
if (!res) y[3] = zero; |
else |
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); |
y[3] = lnegi(gun); |
} |
} |
return y; |
y[4] = un; return y; |
} |
} |
|
|
GEN |
GEN |
Line 3609 newtonpoly(GEN x, GEN p) |
|
Line 4025 newtonpoly(GEN x, GEN p) |
|
{ |
{ |
GEN y; |
GEN y; |
long n,ind,a,b,c,u1,u2,r1,r2; |
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"); |
if (typ(x)!=t_POL) err(typeer,"newtonpoly"); |
n=degpol(x); if (n<=0) { y=cgetg(1,t_VEC); return y; } |
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; |
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 */ |
/* Factor polynomial a on the number field defined by polynomial t */ |
GEN |
GEN |
polfnf(GEN a, GEN t) |
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; |
long lx,i,k,e; |
int sqfree, tmonic; |
int sqfree, tmonic; |
|
|
Line 3655 polfnf(GEN a, GEN t) |
|
Line 4066 polfnf(GEN a, GEN t) |
|
a = fix_relative_pol(t,a,0); |
a = fix_relative_pol(t,a,0); |
alift = lift(a); |
alift = lift(a); |
p1 = content(alift); if (!gcmp1(p1)) { a = gdiv(a, p1); 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)); |
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); |
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 */ |
/* n guaranteed to be squarefree */ |
fa = squff2(n,0); lx = lg(fa); |
fa = DDF2(n,0); lx = lg(fa); |
y = cgetg(3,t_MAT); |
y = cgetg(3,t_MAT); |
p1 = cgetg(lx,t_COL); y[1] = (long)p1; |
p1 = cgetg(lx,t_COL); y[1] = (long)p1; |
p2 = cgetg(lx,t_COL); y[2] = (long)p2; |
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))); |
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); |
F = nfgcd(u, F, t, dent); |
if (typ(F) != t_POL || degpol(F) == 0) |
if (typ(F) != t_POL || degpol(F) == 0) |
err(talker,"reducible modulus in factornf"); |
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); |
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)); |
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 |
static GEN |
to_frac(GEN a, GEN b) |
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 |
GEN |
matratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN denom) |
matratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN denom) |
{ |
{ |
ulong ltop = avma; |
gpmem_t ltop = avma; |
GEN N, a; |
GEN N, a; |
long l2, l3; |
long l2, l3; |
long i, j; |
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 |
GEN |
polratlift(GEN P, GEN mod, GEN amax, GEN bmax, GEN denom) |
polratlift(GEN P, GEN mod, GEN amax, GEN bmax, GEN denom) |
{ |
{ |
ulong ltop = avma; |
gpmem_t ltop = avma; |
GEN Q,a; |
GEN Q,a; |
long l2, j; |
long l2, j; |
if (typ(P)!=t_POL) err(typeer,"polratlift"); |
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, |
* If not NULL, den must a a multiple of the denominator of the gcd, |
* for example the discriminant of nf. |
* 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 |
GEN |
nfgcd(GEN P, GEN Q, GEN nf, GEN den) |
nfgcd(GEN P, GEN Q, GEN nf, GEN den) |
{ |
{ |
ulong ltop = avma; |
gpmem_t ltop = avma; |
GEN sol, mod = NULL; |
GEN sol, mod = NULL; |
long x = varn(P); |
long x = varn(P); |
long y = varn(nf); |
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))); |
den = mulii(den, mppgcd(ZX_resultant(lP, nf), ZX_resultant(lQ, nf))); |
/*Modular GCD*/ |
/*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 p; |
long dM=0, dR; |
long dM=0, dR; |
GEN M, dsol, dens; |
GEN M, dsol; |
GEN R, ax, bo; |
GEN R, ax, bo; |
byteptr primepointer; |
byteptr primepointer; |
for (p = 27449, primepointer = diffptr + 3000; ; p += *(primepointer++)) |
for (p = 27449, primepointer = diffptr + 3000; ; ) |
{ |
{ |
if (!*primepointer) err(primer1); |
if (!*primepointer) err(primer1); |
if (!smodis(den, p)) |
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) |
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); |
dR = degpol(R); |
if (dR == 0) return scalarpol(gun, x); |
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); |
R = polpol_to_mat(R, d); |
/* previous primes divided Res(P/gcd, Q/gcd)? Discard them. */ |
/* 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))) |
if (low_stack(st_lim, stack_lim(btop, 1))) |
{ |
{ |
GEN *bptr[2]; bptr[0]=&M; bptr[1]=&mod; |
|
if (DEBUGMEM>1) err(warnmem,"nfgcd"); |
if (DEBUGMEM>1) err(warnmem,"nfgcd"); |
gerepilemany(btop, bptr, 2); |
gerepileall(btop, 2, &M, &mod); |
} |
} |
|
|
ax = gmulgs(mpinvmod(stoi(p), mod), p); |
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*/ |
/* I suspect it must be better to take amax > bmax*/ |
bo = racine(shifti(mod, -1)); |
bo = racine(shifti(mod, -1)); |
if ((sol = matratlift(M, mod, bo, bo, den)) == NULL) |
if ((sol = matratlift(M, mod, bo, bo, den)) == NULL) |
continue; |
goto repeat; |
dens = denom(sol); |
|
sol = mat_to_polpol(sol,x,y); |
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); |
return gerepilecopy(ltop, sol); |