version 1.1, 2001/10/02 11:17:05 |
version 1.2, 2002/09/11 07:26:52 |
Line 33 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
Line 33 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
#endif |
#endif |
#define OK_ULONG(p) (lgefint(p) == 3 && u_OK_ULONG(p[2])) |
#define OK_ULONG(p) (lgefint(p) == 3 && u_OK_ULONG(p[2])) |
|
|
|
#define both_odd(x,y) ((x)&(y)&1) |
|
extern GEN caractducos(GEN p, GEN x, int v); |
|
extern double mylog2(GEN z); |
|
extern GEN revpol(GEN x); |
|
|
/*******************************************************************/ |
/*******************************************************************/ |
/* */ |
/* */ |
/* KARATSUBA (for polynomials) */ |
/* KARATSUBA (for polynomials) */ |
Line 80 specpol(GEN x, long nx) |
|
Line 85 specpol(GEN x, long nx) |
|
|
|
/* multiplication in Fp[X], p small */ |
/* multiplication in Fp[X], p small */ |
|
|
static GEN |
GEN |
u_normalizepol(GEN x, long lx) |
u_normalizepol(GEN x, long lx) |
{ |
{ |
long i; |
long i; |
Line 120 u_FpX_gcd_i(GEN a, GEN b, ulong p) |
|
Line 125 u_FpX_gcd_i(GEN a, GEN b, ulong p) |
|
GEN |
GEN |
u_FpX_gcd(GEN a, GEN b, ulong p) |
u_FpX_gcd(GEN a, GEN b, ulong p) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
return gerepileupto(av, u_FpX_gcd_i(a,b,p)); |
return gerepileupto(av, u_FpX_gcd_i(a,b,p)); |
} |
} |
|
|
int |
int |
u_FpX_is_squarefree(GEN z, ulong p) |
u_FpX_is_squarefree(GEN z, ulong p) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
GEN d = u_FpX_gcd_i(z, u_FpX_deriv(z,p) , p); |
GEN d = u_FpX_gcd_i(z, u_FpX_deriv(z,p) , p); |
avma = av; return (degpol(d) == 0); |
avma = av; return (degpol(d) == 0); |
} |
} |
Line 225 u_FpX_addshift(GEN x, GEN y, ulong p, long d) |
|
Line 230 u_FpX_addshift(GEN x, GEN y, ulong p, long d) |
|
*--zd = evaltyp(t_VECSMALL) | evallg(lz); return zd; |
*--zd = evaltyp(t_VECSMALL) | evallg(lz); return zd; |
} |
} |
|
|
#define u_zeropol(malloc) u_allocpol(-1,malloc) |
|
#define u_FpX_mul(x,y, p) u_FpX_Kmul(x+2,y+2,p, lgef(x)-2,lgef(y)-2) |
#define u_FpX_mul(x,y, p) u_FpX_Kmul(x+2,y+2,p, lgef(x)-2,lgef(y)-2) |
#define u_FpX_sqr(x, p) u_FpX_Ksqr(x+2,p, lgef(x)-2) |
#define u_FpX_sqr(x, p) u_FpX_Ksqr(x+2,p, lgef(x)-2) |
#define u_FpX_add(x,y, p) u_FpX_addspec(x+2,y+2,p, lgef(x)-2,lgef(y)-2) |
#define u_FpX_add(x,y, p) u_FpX_addspec(x+2,y+2,p, lgef(x)-2,lgef(y)-2) |
|
|
|
GEN |
|
u_zeropol() |
|
{ |
|
GEN z = cgetg(2, t_VECSMALL); |
|
z[1] = evallgef(2) | evalvarn(0); return z; |
|
} |
|
|
static GEN |
static GEN |
u_allocpol(long d, int malloc) |
u_mallocpol(long d) |
{ |
{ |
GEN z = malloc? (GEN)gpmalloc((d+3) * sizeof(long)): new_chunk(d+3); |
GEN z = (GEN)gpmalloc((d+3) * sizeof(long)); |
z[0] = evaltyp(t_VECSMALL) | evallg(d+3); |
z[0] = evaltyp(t_VECSMALL) | evallg(d+3); |
z[1] = evalsigne((d >= 0)) | evallgef(d+3) | evalvarn(0); |
z[1] = evalsigne((d >= 0)) | evallgef(d+3) | evalvarn(0); |
return z; |
return z; |
} |
} |
|
static GEN |
|
u_getpol(long d) |
|
{ |
|
GEN z = cgetg(d + 3, t_VECSMALL); |
|
z[1] = evalsigne((d >= 0)) | evallgef(d+3) | evalvarn(0); |
|
return z; |
|
} |
|
|
static GEN |
static GEN |
u_scalarpol(ulong x, int malloc) |
u_scalarpol(ulong x) |
{ |
{ |
GEN z = u_allocpol(0,malloc); |
GEN z; |
|
if (!x) return u_zeropol(); |
|
z = u_getpol(0); |
z[2] = x; return z; |
z[2] = x; return z; |
} |
} |
|
|
static GEN |
static GEN |
|
u_FpX_neg(GEN x, ulong p) |
|
{ |
|
long i, l = lgef(x); |
|
GEN z = cgetg(l, t_VECSMALL); |
|
z[1] = x[1]; |
|
for (i=2; i<l; i++) z[i] = x[i]? p - x[i]: 0; |
|
return z; |
|
} |
|
|
|
static GEN |
u_FpX_neg_i(GEN x, ulong p) |
u_FpX_neg_i(GEN x, ulong p) |
{ |
{ |
long i, l = lgef(x); |
long i, l = lgef(x); |
Line 257 u_FpX_neg_i(GEN x, ulong p) |
|
Line 287 u_FpX_neg_i(GEN x, ulong p) |
|
|
|
/* shift polynomial + gerepile */ |
/* shift polynomial + gerepile */ |
static GEN |
static GEN |
u_FpX_shiftip(long av, GEN x, long v) |
u_FpX_shiftip(gpmem_t av, GEN x, long v) |
{ |
{ |
long i, lx = lgef(x), ly; |
long i, lx = lgef(x), ly; |
GEN y; |
GEN y; |
|
|
u_FpX_Kmul(GEN a, GEN b, ulong p, long na, long nb) |
u_FpX_Kmul(GEN a, GEN b, ulong p, long na, long nb) |
{ |
{ |
GEN a0,c,c0; |
GEN a0,c,c0; |
long av,n0,n0a,i, v = 0; |
long n0, n0a, i, v = 0; |
|
gpmem_t av; |
|
|
while (na && !a[0]) { a++; na--; v++; } |
while (na && !a[0]) { a++; na--; v++; } |
while (nb && !b[0]) { b++; nb--; v++; } |
while (nb && !b[0]) { b++; nb--; v++; } |
if (na < nb) swapspec(a,b, na,nb); |
if (na < nb) swapspec(a,b, na,nb); |
if (!nb) return u_zeropol(0); |
if (!nb) return u_zeropol(); |
|
|
av = avma; |
av = avma; |
if (nb < u_MUL_LIMIT) |
if (nb < u_MUL_LIMIT) |
Line 328 u_FpX_sqrpol(GEN x, ulong p, long nx) |
|
Line 359 u_FpX_sqrpol(GEN x, ulong p, long nx) |
|
ulong p1; |
ulong p1; |
GEN z; |
GEN z; |
|
|
if (!nx) return u_zeropol(0); |
if (!nx) return u_zeropol(); |
lz = (nx << 1) + 1, nz = lz-2; |
lz = (nx << 1) + 1, nz = lz-2; |
z = cgetg(lz,t_VECSMALL) + 2; |
z = cgetg(lz,t_VECSMALL) + 2; |
for (i=0; i<nx; i++) |
for (i=0; i<nx; i++) |
|
|
u_FpX_Ksqr(GEN a, ulong p, long na) |
u_FpX_Ksqr(GEN a, ulong p, long na) |
{ |
{ |
GEN a0,c,c0,c1; |
GEN a0,c,c0,c1; |
long av,n0,n0a,i, v = 0; |
long n0, n0a, i, v = 0; |
|
gpmem_t av; |
|
|
while (na && !a[0]) { a++; na--; v += 2; } |
while (na && !a[0]) { a++; na--; v += 2; } |
av = avma; |
av = avma; |
|
|
mulpol_limb(GEN x, GEN y, char *ynonzero, long a, long b) |
mulpol_limb(GEN x, GEN y, char *ynonzero, long a, long b) |
{ |
{ |
GEN p1 = NULL; |
GEN p1 = NULL; |
long i,av = avma; |
long i; |
|
gpmem_t av = avma; |
for (i=a; i<b; i++) |
for (i=a; i<b; i++) |
if (ynonzero[i]) |
if (ynonzero[i]) |
{ |
{ |
|
|
quickmul(GEN a, GEN b, long na, long nb) |
quickmul(GEN a, GEN b, long na, long nb) |
{ |
{ |
GEN a0,c,c0; |
GEN a0,c,c0; |
long av,n0,n0a,i, v = 0; |
long n0, n0a, i, v = 0; |
|
gpmem_t av; |
|
|
while (na && isexactzero((GEN)a[0])) { a++; na--; v++; } |
while (na && isexactzero((GEN)a[0])) { a++; na--; v++; } |
while (nb && isexactzero((GEN)b[0])) { b++; nb--; v++; } |
while (nb && isexactzero((GEN)b[0])) { b++; nb--; v++; } |
Line 597 quickmul(GEN a, GEN b, long na, long nb) |
|
Line 631 quickmul(GEN a, GEN b, long na, long nb) |
|
GEN |
GEN |
sqrpol(GEN x, long nx) |
sqrpol(GEN x, long nx) |
{ |
{ |
long av,i,j,l,lz,nz; |
long i, j, l, lz, nz; |
|
gpmem_t av; |
GEN p1,z; |
GEN p1,z; |
char *p2; |
char *p2; |
|
|
|
|
quicksqr(GEN a, long na) |
quicksqr(GEN a, long na) |
{ |
{ |
GEN a0,c,c0,c1; |
GEN a0,c,c0,c1; |
long av,n0,n0a,i, v = 0; |
long n0, n0a, i, v = 0; |
|
gpmem_t av; |
|
|
while (na && isexactzero((GEN)a[0])) { a++; na--; v += 2; } |
while (na && isexactzero((GEN)a[0])) { a++; na--; v += 2; } |
if (v) (void)new_chunk(v); |
if (v) (void)new_chunk(v); |
Line 663 There are clean and memory efficient. |
|
Line 699 There are clean and memory efficient. |
|
GEN |
GEN |
FpX_center(GEN T,GEN mod) |
FpX_center(GEN T,GEN mod) |
{/*OK centermod exists, but is not so clean*/ |
{/*OK centermod exists, but is not so clean*/ |
ulong av; |
gpmem_t av; |
long i, l=lg(T); |
long i, l=lg(T); |
GEN P,mod2; |
GEN P,mod2; |
P=cgetg(l,t_POL); |
P=cgetg(l,t_POL); |
Line 712 FpX_add(GEN x,GEN y,GEN p) |
|
Line 748 FpX_add(GEN x,GEN y,GEN p) |
|
for (i=2; i<ly; i++) z[i]=laddii((GEN)x[i],(GEN)y[i]); |
for (i=2; i<ly; i++) z[i]=laddii((GEN)x[i],(GEN)y[i]); |
for ( ; i<lx; i++) z[i]=licopy((GEN)x[i]); |
for ( ; i<lx; i++) z[i]=licopy((GEN)x[i]); |
(void)normalizepol_i(z, lx); |
(void)normalizepol_i(z, lx); |
if (lgef(z) == 2) { avma = (long)(z + lx); z = zeropol(varn(x)); } |
if (lgef(z) == 2) { avma = (gpmem_t)(z + lx); z = zeropol(varn(x)); } |
if (p) z= FpX_red(z, p); |
if (p) z= FpX_red(z, p); |
return z; |
return z; |
} |
} |
Line 738 FpX_sub(GEN x,GEN y,GEN p) |
|
Line 774 FpX_sub(GEN x,GEN y,GEN p) |
|
for ( ; i<ly; i++) z[i]=lnegi((GEN)y[i]); |
for ( ; i<ly; i++) z[i]=lnegi((GEN)y[i]); |
/*polynomial is always normalized*/ |
/*polynomial is always normalized*/ |
} |
} |
if (lgef(z) == 2) { avma = (long)(z + lz); z = zeropol(varn(x)); } |
if (lgef(z) == 2) { avma = (gpmem_t)(z + lz); z = zeropol(varn(x)); } |
if (p) z= FpX_red(z, p); |
if (p) z= FpX_red(z, p); |
return z; |
return z; |
} |
} |
Line 758 FpX_sqr(GEN x,GEN p) |
|
Line 794 FpX_sqr(GEN x,GEN p) |
|
if (!p) return z; |
if (!p) return z; |
return FpX_red(z, p); |
return FpX_red(z, p); |
} |
} |
#define u_FpX_div(x,y,p) u_FpX_divrem((x),(y),(p),(0),NULL) |
#define u_FpX_div(x,y,p) u_FpX_divrem((x),(y),(p),NULL) |
|
|
/* Product of y and x in Z/pZ[X]/(pol), as t_VECSMALL. Assume OK_ULONG(p) */ |
/* Product of y and x in Z/pZ[X]/(pol), as t_VECSMALL. Assume OK_ULONG(p) */ |
static GEN |
static GEN |
Line 784 FpXQ_invsafe(GEN x, GEN T, GEN p) |
|
Line 820 FpXQ_invsafe(GEN x, GEN T, GEN p) |
|
|
|
if (typ(x) != t_POL) return mpinvmod(x, p); |
if (typ(x) != t_POL) return mpinvmod(x, p); |
z = FpX_extgcd(x, T, p, &U, &V); |
z = FpX_extgcd(x, T, p, &U, &V); |
if (lgef(z) != 3) return NULL; |
if (degpol(z)) return NULL; |
z = mpinvmod((GEN)z[2], p); |
z = mpinvmod((GEN)z[2], p); |
return FpX_Fp_mul(U, z, p); |
return FpX_Fp_mul(U, z, p); |
} |
} |
|
|
/* Product of y and x in Z/pZ[X]/(T) |
/* Product of y and x in Z/pZ[X]/(T) |
* return lift(lift(Mod(x*y*Mod(1,p),T*Mod(1,p)))) */ |
* return lift(lift(Mod(x*y*Mod(1,p),T*Mod(1,p)))) */ |
|
/* x and y must be poynomials in the same var as T. |
|
* t_INT are not allowed. Use Fq_mul instead. |
|
*/ |
GEN |
GEN |
FpXQ_mul(GEN y,GEN x,GEN T,GEN p) |
FpXQ_mul(GEN y,GEN x,GEN T,GEN p) |
{ |
{ |
GEN z; |
GEN z; |
if (typ(x) == t_INT) |
if (typ(x) == t_INT || typ(y) == t_INT) |
{ |
err(bugparier,"FpXQ_mul, t_INT are absolutely forbidden"); |
if (typ(y) == t_INT) return modii(mulii(x,y), p); |
|
return FpX_Fp_mul(y, x, p); |
|
} |
|
if (typ(y) == t_INT) return FpX_Fp_mul(x, y, p); |
|
z = quickmul(y+2, x+2, lgef(y)-2, lgef(x)-2); setvarn(z,varn(y)); |
z = quickmul(y+2, x+2, lgef(y)-2, lgef(x)-2); setvarn(z,varn(y)); |
z = FpX_red(z, p); return FpX_res(z,T, p); |
z = FpX_red(z, p); return FpX_res(z,T, p); |
} |
} |
Line 838 ZX_s_add(GEN y,long x) |
|
Line 873 ZX_s_add(GEN y,long x) |
|
return y; |
return y; |
} |
} |
/* y is a polynomial in ZZ[X] and x an integer. |
/* y is a polynomial in ZZ[X] and x an integer. |
* If p is NULL, no reduction is perfomed and return x*y |
* If p is NULL, no reduction is performed and return x*y |
* |
* |
* else the result is lift(y*x*Mod(1,p)) |
* else the result is lift(y*x*Mod(1,p)) |
*/ |
*/ |
Line 860 FpX_Fp_mul(GEN y,GEN x,GEN p) |
|
Line 895 FpX_Fp_mul(GEN y,GEN x,GEN p) |
|
* End of unclean functions. * |
* End of unclean functions. * |
* * |
* * |
*****************************************************************/ |
*****************************************************************/ |
|
|
|
/* modular power */ |
|
ulong |
|
powuumod(ulong x, ulong n0, ulong p) |
|
{ |
|
ulong y, z, n; |
|
if (n0 <= 2) |
|
{ /* frequent special cases */ |
|
if (n0 == 2) return mulssmod(x,x,p); |
|
if (n0 == 1) return x; |
|
if (n0 == 0) return 1; |
|
} |
|
y = 1; z = x; n = n0; |
|
for(;;) |
|
{ |
|
if (n&1) y = mulssmod(y,z,p); |
|
n>>=1; if (!n) return y; |
|
z = mulssmod(z,z,p); |
|
} |
|
} |
|
|
|
GEN |
|
powgumod(GEN x, ulong n0, GEN p) |
|
{ |
|
GEN y, z; |
|
ulong n; |
|
if (n0 <= 2) |
|
{ /* frequent special cases */ |
|
if (n0 == 2) return resii(sqri(x),p); |
|
if (n0 == 1) return x; |
|
if (n0 == 0) return gun; |
|
} |
|
y = gun; z = x; n = n0; |
|
for(;;) |
|
{ |
|
if (n&1) y = resii(mulii(y,z),p); |
|
n>>=1; if (!n) return y; |
|
z = resii(sqri(z),p); |
|
} |
|
} |
|
|
/***************************************************************** |
/***************************************************************** |
Clean and with no reduced hypothesis. Beware that some operations |
Clean and with no reduced hypothesis. Beware that some operations |
will be much slower with big unreduced coefficient |
will be much slower with big unreduced coefficient |
Line 869 FpX_Fp_mul(GEN y,GEN x,GEN p) |
|
Line 945 FpX_Fp_mul(GEN y,GEN x,GEN p) |
|
GEN |
GEN |
FpXQ_inv(GEN x,GEN T,GEN p) |
FpXQ_inv(GEN x,GEN T,GEN p) |
{ |
{ |
ulong av; |
gpmem_t av; |
GEN U; |
GEN U; |
|
|
if (!T) return mpinvmod(x,p); |
if (!T) return mpinvmod(x,p); |
Line 880 FpXQ_inv(GEN x,GEN T,GEN p) |
|
Line 956 FpXQ_inv(GEN x,GEN T,GEN p) |
|
} |
} |
|
|
GEN |
GEN |
FpXV_FpV_dotproduct(GEN V, GEN W, GEN p) |
FpXQ_div(GEN x,GEN y,GEN T,GEN p) |
{ |
{ |
long ltop=avma; |
gpmem_t av = avma; |
|
return gerepileupto(av, FpXQ_mul(x,FpXQ_inv(y,T,p),T,p)); |
|
} |
|
|
|
GEN |
|
FpXV_FpV_innerprod(GEN V, GEN W, GEN p) |
|
{ |
|
gpmem_t ltop=avma; |
long i; |
long i; |
GEN z = FpX_Fp_mul((GEN)V[1],(GEN)W[1],NULL); |
GEN z = FpX_Fp_mul((GEN)V[1],(GEN)W[1],NULL); |
for(i=2;i<lg(V);i++) |
for(i=2;i<lg(V);i++) |
Line 946 long brent_kung_optpow(long d, long n) |
|
Line 1029 long brent_kung_optpow(long d, long n) |
|
return (d+l-1)/l; |
return (d+l-1)/l; |
} |
} |
|
|
/*Close to FpXV_FpV_dotproduct*/ |
/*Close to FpXV_FpV_innerprod*/ |
|
|
static GEN |
static GEN |
spec_compo_powers(GEN P, GEN V, long a, long n) |
spec_compo_powers(GEN P, GEN V, long a, long n) |
Line 969 spec_compo_powers(GEN P, GEN V, long a, long n) |
|
Line 1052 spec_compo_powers(GEN P, GEN V, long a, long n) |
|
GEN |
GEN |
FpX_FpXQV_compo(GEN P, GEN V, GEN T, GEN p) |
FpX_FpXQV_compo(GEN P, GEN V, GEN T, GEN p) |
{ |
{ |
ulong ltop=avma; |
gpmem_t ltop=avma; |
long l=lg(V)-1; |
long l=lg(V)-1; |
GEN z,u; |
GEN z,u; |
long d=degpol(P),cnt=0; |
long d=degpol(P),cnt=0; |
Line 1003 FpX_FpXQV_compo(GEN P, GEN V, GEN T, GEN p) |
|
Line 1086 FpX_FpXQV_compo(GEN P, GEN V, GEN T, GEN p) |
|
GEN |
GEN |
FpX_FpXQ_compo(GEN T,GEN x,GEN pol,GEN p) |
FpX_FpXQ_compo(GEN T,GEN x,GEN pol,GEN p) |
{ |
{ |
ulong ltop=avma; |
gpmem_t ltop=avma; |
GEN z; |
GEN z; |
long d=degpol(T),rtd; |
long d=degpol(T),rtd; |
if (!signe(T)) return zeropol(varn(T)); |
if (!signe(T)) return zeropol(varn(T)); |
Line 1020 FpX_FpXQ_compo(GEN T,GEN x,GEN pol,GEN p) |
|
Line 1103 FpX_FpXQ_compo(GEN T,GEN x,GEN pol,GEN p) |
|
GEN |
GEN |
FpX_eval(GEN x,GEN y,GEN p) |
FpX_eval(GEN x,GEN y,GEN p) |
{ |
{ |
ulong av; |
gpmem_t av; |
GEN p1,r,res; |
GEN p1,r,res; |
long i,j; |
long i,j; |
i=lgef(x)-1; |
i=lgef(x)-1; |
Line 1035 FpX_eval(GEN x,GEN y,GEN p) |
|
Line 1118 FpX_eval(GEN x,GEN y,GEN p) |
|
for (j=i; !signe((GEN)x[j]); j--) |
for (j=i; !signe((GEN)x[j]); j--) |
if (j==2) |
if (j==2) |
{ |
{ |
if (i!=j) y = powmodulo(y,stoi(i-j+1),p); |
if (i!=j) y = powgumod(y,i-j+1,p); |
p1=mulii(p1,y); |
p1=mulii(p1,y); |
goto fppoleval;/*sorry break(2) no implemented*/ |
goto fppoleval;/*sorry break(2) no implemented*/ |
} |
} |
r = (i==j)? y: powmodulo(y,stoi(i-j+1),p); |
r = (i==j)? y: powgumod(y,i-j+1,p); |
p1 = modii(addii(mulii(p1,r), (GEN)x[j]),p); |
p1 = modii(addii(mulii(p1,r), (GEN)x[j]),p); |
} |
} |
fppoleval: |
fppoleval: |
Line 1056 FpX_eval(GEN x,GEN y,GEN p) |
|
Line 1139 FpX_eval(GEN x,GEN y,GEN p) |
|
GEN |
GEN |
FpX_chinese_coprime(GEN x,GEN y,GEN Tx,GEN Ty,GEN Tz,GEN p) |
FpX_chinese_coprime(GEN x,GEN y,GEN Tx,GEN Ty,GEN Tz,GEN p) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
GEN ax,p1; |
GEN ax,p1; |
ax = FpX_mul(FpXQ_inv(Tx,Ty,p), Tx,p); |
ax = FpX_mul(FpXQ_inv(Tx,Ty,p), Tx,p); |
p1=FpX_mul(ax, FpX_sub(y,x,p),p); |
p1=FpX_mul(ax, FpX_sub(y,x,p),p); |
Line 1066 FpX_chinese_coprime(GEN x,GEN y,GEN Tx,GEN Ty,GEN Tz,G |
|
Line 1149 FpX_chinese_coprime(GEN x,GEN y,GEN Tx,GEN Ty,GEN Tz,G |
|
return gerepileupto(av,p1); |
return gerepileupto(av,p1); |
} |
} |
|
|
|
typedef struct { |
|
GEN pol, p; |
|
} FpX_muldata; |
|
typedef struct { |
|
GEN pol; |
|
ulong p; |
|
} u_FpX_muldata; |
|
|
|
static GEN |
|
_u_sqr(void *data, GEN x) |
|
{ |
|
u_FpX_muldata *D = (u_FpX_muldata*)data; |
|
return u_FpXQ_sqr(x, D->pol, D->p); |
|
} |
|
static GEN |
|
_u_mul(void *data, GEN x, GEN y) |
|
{ |
|
u_FpX_muldata *D = (u_FpX_muldata*)data; |
|
return u_FpXQ_mul(x,y, D->pol, D->p); |
|
} |
|
static GEN |
|
_sqr(void *data, GEN x) |
|
{ |
|
FpX_muldata *D = (FpX_muldata*)data; |
|
return FpXQ_sqr(x, D->pol, D->p); |
|
} |
|
static GEN |
|
_mul(void *data, GEN x, GEN y) |
|
{ |
|
FpX_muldata *D = (FpX_muldata*)data; |
|
return FpXQ_mul(x,y, D->pol, D->p); |
|
} |
|
|
/* assume n > 0 */ |
/* assume n > 0 */ |
GEN |
GEN |
u_FpXQ_pow(GEN x, GEN n, GEN pol, ulong p) |
u_FpXQ_pow(GEN x, GEN n, GEN pol, ulong p) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
GEN p1 = n+2, y = x; |
u_FpX_muldata D; |
long m,i,j; |
GEN y; |
m = *p1; |
D.pol = pol; |
j = 1+bfffo(m); m <<= j; j = BITS_IN_LONG-j; |
D.p = p; |
for (i=lgefint(n)-2;;) |
y = leftright_pow(x, n, (void*)&D, &_u_sqr, &_u_mul); |
{ |
|
for (; j; m<<=1,j--) |
|
{ |
|
y = u_FpXQ_sqr(y,pol,p); |
|
if (m<0) |
|
y = u_FpXQ_mul(y,x,pol,p); |
|
} |
|
if (--i == 0) break; |
|
m = *++p1, j = BITS_IN_LONG; |
|
} |
|
return gerepileupto(av, y); |
return gerepileupto(av, y); |
} |
} |
|
|
Line 1093 u_FpXQ_pow(GEN x, GEN n, GEN pol, ulong p) |
|
Line 1199 u_FpXQ_pow(GEN x, GEN n, GEN pol, ulong p) |
|
GEN |
GEN |
FpXQ_pow(GEN x, GEN n, GEN pol, GEN p) |
FpXQ_pow(GEN x, GEN n, GEN pol, GEN p) |
{ |
{ |
ulong av, ltop = avma, lim=stack_lim(avma,1); |
gpmem_t av = avma; |
long m,i,j, vx = varn(x); |
FpX_muldata D; |
GEN p1 = n+2, y; |
long vx = varn(x); |
|
GEN y; |
if (!signe(n)) return polun[vx]; |
if (!signe(n)) return polun[vx]; |
if (signe(n) < 0) |
if (signe(n) < 0) |
{ |
{ |
Line 1107 FpXQ_pow(GEN x, GEN n, GEN pol, GEN p) |
|
Line 1214 FpXQ_pow(GEN x, GEN n, GEN pol, GEN p) |
|
if (OK_ULONG(p)) |
if (OK_ULONG(p)) |
{ |
{ |
ulong pp = p[2]; |
ulong pp = p[2]; |
pol = u_Fp_FpX(pol,0, pp); |
pol = u_Fp_FpX(pol, pp); |
x = u_Fp_FpX(x,0, pp); |
x = u_Fp_FpX(x, pp); |
y = u_FpXQ_pow(x, n, pol, pp); |
y = u_FpXQ_pow(x, n, pol, pp); |
y = small_to_pol(y,vx); |
y = small_to_pol(y,vx); |
} |
} |
else |
else |
{ |
{ |
av = avma; |
av = avma; |
m = *p1; y = x; |
D.pol = pol; |
j = 1+bfffo(m); m <<= j; j = BITS_IN_LONG-j; |
D.p = p; |
for (i=lgefint(n)-2;;) |
y = leftright_pow(x, n, (void*)&D, &_sqr, &_mul); |
{ |
|
for (; j; m<<=1,j--) |
|
{ |
|
y = FpXQ_sqr(y,pol,p); |
|
if (low_stack(lim, stack_lim(av,1))) |
|
{ |
|
if(DEBUGMEM>1) err(warnmem,"[1]: FpXQ_pow"); |
|
y = gerepileupto(av, y); |
|
} |
|
if (m<0) |
|
y = FpXQ_mul(y,x,pol,p); |
|
if (low_stack(lim, stack_lim(av,1))) |
|
{ |
|
if(DEBUGMEM>1) err(warnmem,"[2]: FpXQ_pow"); |
|
y = gerepileupto(av, y); |
|
} |
|
} |
|
if (--i == 0) break; |
|
m = *++p1, j = BITS_IN_LONG; |
|
} |
|
} |
} |
return gerepileupto(ltop,y); |
return gerepileupto(av, y); |
} |
} |
|
|
|
GEN |
|
Fq_pow(GEN x, GEN n, GEN pol, GEN p) |
|
{ |
|
if (typ(x) == t_INT) return powmodulo(x,n,p); |
|
return FpXQ_pow(x,n,pol,p); |
|
} |
|
|
/*******************************************************************/ |
/*******************************************************************/ |
/* */ |
/* */ |
/* Fp[X][Y] */ |
/* Fp[X][Y] */ |
Line 1152 FpXQ_pow(GEN x, GEN n, GEN pol, GEN p) |
|
Line 1246 FpXQ_pow(GEN x, GEN n, GEN pol, GEN p) |
|
GEN |
GEN |
FpXX_red(GEN z, GEN p) |
FpXX_red(GEN z, GEN p) |
{ |
{ |
GEN res; |
GEN res; |
int i; |
int i; |
res=cgetg(lgef(z),t_POL); |
res=cgetg(lgef(z),t_POL); |
res[1] = evalsigne(1) | evalvarn(varn(z)) | evallgef(lgef(z)); |
res[1] = evalsigne(1) | evalvarn(varn(z)) | evallgef(lgef(z)); |
for(i=2;i<lgef(res);i++) |
for(i=2;i<lgef(res);i++) |
if (typ(z[i])!=t_INT) |
if (typ(z[i])!=t_INT) |
|
{ |
|
gpmem_t av=avma; |
|
res[i]=(long)FpX_red((GEN)z[i],p); |
|
if (lgef(res[i])<=3) |
{ |
{ |
ulong av=avma; |
if (lgef(res[i])==2) {avma=av;res[i]=zero;} |
res[i]=(long)FpX_red((GEN)z[i],p); |
else res[i]=lpilecopy(av,gmael(res,i,2)); |
if (lgef(res[i])<=3) |
|
{ |
|
if (lgef(res[i])==2) {avma=av;res[i]=zero;} |
|
else res[i]=lpilecopy(av,gmael(res,i,2)); |
|
} |
|
} |
} |
else |
} |
res[i]=lmodii((GEN)z[i],p); |
else |
res=normalizepol_i(res,lgef(res)); |
res[i]=lmodii((GEN)z[i],p); |
return res; |
res=normalizepol_i(res,lgef(res)); |
|
return res; |
} |
} |
|
|
/*******************************************************************/ |
/*******************************************************************/ |
Line 1178 FpXX_red(GEN z, GEN p) |
|
Line 1272 FpXX_red(GEN z, GEN p) |
|
/* (Fp[X]/(Q))[Y] */ |
/* (Fp[X]/(Q))[Y] */ |
/* */ |
/* */ |
/*******************************************************************/ |
/*******************************************************************/ |
|
|
extern GEN to_Kronecker(GEN P, GEN Q); |
extern GEN to_Kronecker(GEN P, GEN Q); |
GEN /*Somewhat copy-pasted...*/ |
|
/*Not malloc nor warn-clean.*/ |
/*Not malloc nor warn-clean.*/ |
FpXQX_from_Kronecker(GEN z, GEN pol, GEN p) |
GEN |
|
FpXQX_from_Kronecker(GEN Z, GEN T, GEN p) |
{ |
{ |
long i,j,lx,l = lgef(z), N = (degpol(pol)<<1) + 1; |
long i,j,lx,l = lgef(Z), N = (degpol(T)<<1) + 1; |
GEN x, t = cgetg(N,t_POL); |
GEN x, t = cgetg(N,t_POL), z = FpX_red(Z, p); |
t[1] = pol[1] & VARNBITS; |
t[1] = T[1] & VARNBITS; |
lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL); |
lx = (l-2) / (N-2); |
if (isonstack(pol)) pol = gcopy(pol); |
x = cgetg(lx+3,t_POL); |
for (i=2; i<lx+2; i++) |
for (i=2; i<lx+2; i++) |
{ |
{ |
for (j=2; j<N; j++) t[j] = z[j]; |
for (j=2; j<N; j++) t[j] = z[j]; |
z += (N-2); |
z += (N-2); |
x[i] = (long)FpX_res(normalizepol_i(t,N), pol, p); |
x[i] = (long)FpX_res(normalizepol_i(t,N), T,p); |
} |
} |
N = (l-2) % (N-2) + 2; |
N = (l-2) % (N-2) + 2; |
for (j=2; j<N; j++) t[j] = z[j]; |
for (j=2; j<N; j++) t[j] = z[j]; |
x[i] = (long)FpX_res(normalizepol_i(t,N), pol, p); |
x[i] = (long)FpX_res(normalizepol_i(t,N), T,p); |
return normalizepol_i(x, i+1); |
return normalizepol_i(x, i+1); |
} |
} |
/*Unused/untested*/ |
|
GEN |
GEN |
|
FpVQX_red(GEN z, GEN T, GEN p) |
|
{ |
|
GEN res; |
|
int i, l = lg(z); |
|
res=cgetg(l,typ(z)); |
|
for(i=1;i<l;i++) |
|
if (typ(z[i]) != t_INT) |
|
{ |
|
if (T) |
|
res[i]=(long)FpX_res((GEN)z[i],T,p); |
|
else |
|
res[i]=(long)FpX_red((GEN)z[i],p); |
|
} |
|
else |
|
res[i] = lmodii((GEN)z[i],p); |
|
return res; |
|
} |
|
GEN |
FpXQX_red(GEN z, GEN T, GEN p) |
FpXQX_red(GEN z, GEN T, GEN p) |
{ |
{ |
GEN res; |
GEN res; |
int i; |
int i, l = lgef(z); |
res=cgetg(lgef(z),t_POL); |
res=cgetg(l,t_POL); |
res[1] = evalsigne(1) | evalvarn(varn(z)) | evallgef(lgef(z)); |
res[1] = evalsigne(1) | evalvarn(varn(z)) | evallgef(lgef(z)); |
for(i=2;i<lgef(res);i++) |
for(i=2;i<l;i++) |
if (typ(z[i])!=t_INT) |
if (typ(z[i]) != t_INT) |
{ |
{ |
if (T) |
if (T) |
res[i]=(long)FpX_res((GEN)z[i],T,p); |
res[i]=(long)FpX_res((GEN)z[i],T,p); |
Line 1225 FpXQX_red(GEN z, GEN T, GEN p) |
|
Line 1337 FpXQX_red(GEN z, GEN T, GEN p) |
|
GEN |
GEN |
FpXQX_mul(GEN x, GEN y, GEN T, GEN p) |
FpXQX_mul(GEN x, GEN y, GEN T, GEN p) |
{ |
{ |
ulong ltop; |
gpmem_t ltop; |
GEN z,kx,ky; |
GEN z,kx,ky; |
long vx; |
long vx; |
if (!T) return FpX_mul(x,y,p); |
if (!T) return FpX_mul(x,y,p); |
Line 1234 FpXQX_mul(GEN x, GEN y, GEN T, GEN p) |
|
Line 1346 FpXQX_mul(GEN x, GEN y, GEN T, GEN p) |
|
kx= to_Kronecker(x,T); |
kx= to_Kronecker(x,T); |
ky= to_Kronecker(y,T); |
ky= to_Kronecker(y,T); |
z = quickmul(ky+2, kx+2, lgef(ky)-2, lgef(kx)-2); |
z = quickmul(ky+2, kx+2, lgef(ky)-2, lgef(kx)-2); |
z = FpX_red(z,p); |
|
z = FpXQX_from_Kronecker(z,T,p); |
z = FpXQX_from_Kronecker(z,T,p); |
setvarn(z,vx);/*quickmul and Fq_from_Kronecker are not varn-clean*/ |
setvarn(z,vx);/*quickmul and Fq_from_Kronecker are not varn-clean*/ |
return gerepileupto(ltop,z); |
return gerepileupto(ltop,z); |
Line 1242 FpXQX_mul(GEN x, GEN y, GEN T, GEN p) |
|
Line 1353 FpXQX_mul(GEN x, GEN y, GEN T, GEN p) |
|
GEN/*Unused/untested*/ |
GEN/*Unused/untested*/ |
FpXQX_sqr(GEN x, GEN T, GEN p) |
FpXQX_sqr(GEN x, GEN T, GEN p) |
{ |
{ |
ulong ltop=avma; |
gpmem_t ltop=avma; |
GEN z,kx; |
GEN z,kx; |
long vx=varn(x); |
long vx=varn(x); |
kx= to_Kronecker(x,T); |
kx= to_Kronecker(x,T); |
z = quicksqr(kx+2, lgef(kx)-2); |
z = quicksqr(kx+2, lgef(kx)-2); |
z = FpX_red(z,p); |
|
z = FpXQX_from_Kronecker(z,T,p); |
z = FpXQX_from_Kronecker(z,T,p); |
setvarn(z,vx);/*quickmul and Fq_from_Kronecker are nor varn-clean*/ |
setvarn(z,vx);/*quickmul and Fq_from_Kronecker are nor varn-clean*/ |
return gerepileupto(ltop,z); |
return gerepileupto(ltop,z); |
} |
} |
|
|
GEN |
GEN |
FpXQX_FpXQ_mul(GEN P, GEN U, GEN T, GEN p) |
FpXQX_FpXQ_mul(GEN P, GEN U, GEN T, GEN p) |
{ |
{ |
Line 1259 FpXQX_FpXQ_mul(GEN P, GEN U, GEN T, GEN p) |
|
Line 1370 FpXQX_FpXQ_mul(GEN P, GEN U, GEN T, GEN p) |
|
GEN res = cgetg(lP,t_POL); |
GEN res = cgetg(lP,t_POL); |
res[1] = evalsigne(1) | evalvarn(varn(P)) | evallgef(lP); |
res[1] = evalsigne(1) | evalvarn(varn(P)) | evallgef(lP); |
for(i=2; i<lP; i++) |
for(i=2; i<lP; i++) |
res[i] = (long)FpXQ_mul(U,(GEN)P[i], T,p); |
res[i] = (long)Fq_mul(U,(GEN)P[i], T,p); |
return normalizepol_i(res,lgef(res)); |
return normalizepol_i(res,lgef(res)); |
} |
} |
|
|
Line 1280 monomial(GEN a, int degpol, int v) |
|
Line 1391 monomial(GEN a, int degpol, int v) |
|
GEN |
GEN |
FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p) |
FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p) |
{ |
{ |
ulong ltop = avma, btop, st_lim; |
gpmem_t btop, ltop = avma, st_lim; |
long dg, vx = varn(P); |
long dg, vx = varn(P); |
GEN U, q; |
GEN U, q; |
P = FpXX_red(P, p); |
P = FpXX_red(P, p); btop = avma; |
Q = FpXX_red(Q, p); |
Q = FpXX_red(Q, p); |
if (!signe(P)) return gerepileupto(ltop, Q); |
if (!signe(P)) return gerepileupto(ltop, Q); |
if (!signe(Q)) { avma = (ulong)P; return P; } |
if (!signe(Q)) { avma = btop; return P; } |
T = FpX_red(T, p); |
T = FpX_red(T, p); |
|
|
btop = avma; st_lim = stack_lim(btop, 1); |
btop = avma; st_lim = stack_lim(btop, 1); |
Line 1298 FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p) |
|
Line 1409 FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p) |
|
if (!U) { avma = ltop; return NULL; } |
if (!U) { avma = ltop; return NULL; } |
do /* set P := P % Q */ |
do /* set P := P % Q */ |
{ |
{ |
q = FpXQ_mul(U, gneg(leading_term(P)), T, p); |
q = Fq_mul(U, gneg(leading_term(P)), T, p); |
P = gadd(P, FpXQX_mul(monomial(q, dg, vx), Q, T, p)); |
P = gadd(P, FpXQX_mul(monomial(q, dg, vx), Q, T, p)); |
P = FpXQX_red(P, T, p); /* wasteful, but negligible */ |
P = FpXQX_red(P, T, p); /* wasteful, but negligible */ |
dg = lgef(P)-lgef(Q); |
dg = lgef(P)-lgef(Q); |
Line 1319 FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p) |
|
Line 1430 FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p) |
|
|
|
/*******************************************************************/ |
/*******************************************************************/ |
/* */ |
/* */ |
|
/* (Fp[X]/T(X))[Y] / S(Y) */ |
|
/* */ |
|
/*******************************************************************/ |
|
/* TODO: merge these 2 (I don't understand the first one) -- KB */ |
|
|
|
/*Preliminary implementation to speed up Fp_isom*/ |
|
typedef struct { |
|
GEN S, T, p; |
|
} FpXQYQ_muldata; |
|
|
|
static GEN |
|
FpXQYQ_redswap(GEN x, GEN S, GEN T, GEN p) |
|
{ |
|
gpmem_t ltop=avma; |
|
long n=degpol(S); |
|
long m=degpol(T); |
|
long v=varn(T),w=varn(S); |
|
GEN V = swap_polpol(x,n,w); |
|
setvarn(T,w); |
|
V = FpXQX_red(V,T,p); |
|
setvarn(T,v); |
|
V = swap_polpol(V,m,w); |
|
return gerepilecopy(ltop,V); |
|
} |
|
static GEN |
|
FpXQYQ_sqr(void *data, GEN x) |
|
{ |
|
FpXQYQ_muldata *D = (FpXQYQ_muldata*)data; |
|
return FpXQYQ_redswap(FpXQX_sqr(x, D->S, D->p),D->S,D->T,D->p); |
|
|
|
} |
|
static GEN |
|
FpXQYQ_mul(void *data, GEN x, GEN y) |
|
{ |
|
FpXQYQ_muldata *D = (FpXQYQ_muldata*)data; |
|
return FpXQYQ_redswap(FpXQX_mul(x,y, D->S, D->p),D->S,D->T,D->p); |
|
} |
|
|
|
/* x in Z[X,Y], S in Z[X] over Fq = Z[Y]/(p,T); compute lift(x^n mod (S,T,p)) */ |
|
GEN |
|
FpXQYQ_pow(GEN x, GEN n, GEN S, GEN T, GEN p) |
|
{ |
|
gpmem_t av = avma; |
|
FpXQYQ_muldata D; |
|
GEN y; |
|
D.S = S; |
|
D.T = T; |
|
D.p = p; |
|
y = leftright_pow(x, n, (void*)&D, &FpXQYQ_sqr, &FpXQYQ_mul); |
|
return gerepileupto(av, y); |
|
} |
|
|
|
typedef struct { |
|
GEN T, p, S; |
|
long v; |
|
} kronecker_muldata; |
|
|
|
static GEN |
|
_FqXQ_red(void *data, GEN x) |
|
{ |
|
kronecker_muldata *D = (kronecker_muldata*)data; |
|
GEN t = FpXQX_from_Kronecker(x, D->T,D->p); |
|
setvarn(t, D->v); |
|
t = FpXQX_divres(t, D->S,D->T,D->p, ONLY_REM); |
|
return to_Kronecker(t,D->T); |
|
} |
|
static GEN |
|
_FqXQ_mul(void *data, GEN x, GEN y) { |
|
return _FqXQ_red(data, FpX_mul(x,y,NULL)); |
|
} |
|
static GEN |
|
_FqXQ_sqr(void *data, GEN x) { |
|
return _FqXQ_red(data, FpX_sqr(x,NULL)); |
|
} |
|
|
|
/* x over Fq, return lift(x^n) mod S */ |
|
GEN |
|
FqXQ_pow(GEN x, GEN n, GEN S, GEN T, GEN p) |
|
{ |
|
gpmem_t av0 = avma; |
|
long v = varn(x); |
|
GEN y; |
|
kronecker_muldata D; |
|
|
|
x = to_Kronecker(x,T); |
|
D.S = S; |
|
D.T = T; |
|
D.p = p; |
|
D.v = v; |
|
y = leftright_pow(x, n, (void*)&D, &_FqXQ_sqr, &_FqXQ_mul); |
|
y = FpXQX_from_Kronecker(y, T,p); |
|
setvarn(y, v); return gerepileupto(av0, y); |
|
} |
|
/*******************************************************************/ |
|
/* */ |
/* Fq[X] */ |
/* Fq[X] */ |
/* */ |
/* */ |
/*******************************************************************/ |
/*******************************************************************/ |
Line 1328 FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p) |
|
Line 1534 FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p) |
|
#define FqX_sqr FpXQX_sqr |
#define FqX_sqr FpXQX_sqr |
#define FqX_red FpXQX_red |
#define FqX_red FpXQX_red |
static GEN |
static GEN |
Fq_neg(GEN x, GEN T, GEN p)/*T is not used but it is for consistency*/ |
Fq_neg(GEN x, GEN T/*unused*/, GEN p) |
{ |
{ |
switch(typ(x)==t_POL) |
switch(typ(x)==t_POL) |
{ |
{ |
Line 1342 static GEN fgmul(GEN a,GEN b){return FqX_mul(a,b,Tmodu |
|
Line 1548 static GEN fgmul(GEN a,GEN b){return FqX_mul(a,b,Tmodu |
|
GEN |
GEN |
FqV_roots_to_pol(GEN V, GEN Tp, GEN p, long v) |
FqV_roots_to_pol(GEN V, GEN Tp, GEN p, long v) |
{ |
{ |
ulong ltop=avma; |
gpmem_t ltop=avma; |
long k; |
long k; |
GEN W=cgetg(lg(V),t_VEC); |
GEN W=cgetg(lg(V),t_VEC); |
for(k=1;k<lg(V);k++) |
for(k=1;k<lg(V);k++) |
Line 1360 FqV_roots_to_pol(GEN V, GEN Tp, GEN p, long v) |
|
Line 1566 FqV_roots_to_pol(GEN V, GEN Tp, GEN p, long v) |
|
/*NO clean malloc*/ |
/*NO clean malloc*/ |
static GEN fflgen(GEN l, long e, GEN r, GEN T ,GEN p, GEN *zeta) |
static GEN fflgen(GEN l, long e, GEN r, GEN T ,GEN p, GEN *zeta) |
{ |
{ |
ulong av1; |
gpmem_t av1 = avma; |
GEN z,m,m1; |
GEN z, m, m1; |
long x=varn(T),k,u,v,pp,i; |
const long pp = is_bigint(p)? VERYBIGINT: itos(p); |
if (is_bigint(p)) |
long x=varn(T),k,u,v,i; |
pp=VERYBIGINT; |
|
else |
|
pp=itos(p); |
|
z=(lgef(T)==4)?polun[x]:polx[x]; |
|
|
|
av1 = avma; |
for (k=0; ; k++) |
for (k=1; ; k++) |
|
{ |
{ |
u=k;v=0; |
z = (degpol(T)==1)? polun[x]: polx[x]; |
while (u%pp==0){u/=pp;v++;} |
u = k/pp; v=2; /* FpX_Fp_add modify y */ |
if(!v) |
z = gaddgs(z, k%pp); |
z=gadd(z,gun); |
while(u) |
else |
|
{ |
{ |
z=gadd(z,gpowgs(polx[x],v)); |
z = FpX_add(z, monomial(stoi(u%pp),v,x), NULL); |
if (DEBUGLEVEL>=6) |
u /= pp; v++; |
fprintferr("FF l-Gen:next %Z",z); |
|
} |
} |
|
if ( DEBUGLEVEL>=6 ) fprintferr("FF l-Gen:next %Z\n",z); |
m1 = m = FpXQ_pow(z,r,T,p); |
m1 = m = FpXQ_pow(z,r,T,p); |
|
if (gcmp1(m)) { avma = av1; continue; } |
|
|
for (i=1; i<e; i++) |
for (i=1; i<e; i++) |
if (gcmp1(m=FpXQ_pow(m,l,T,p))) break; |
if (gcmp1(m = FpXQ_pow(m,l,T,p))) break; |
if (i==e) break; |
if (i==e) break; |
avma = av1; |
avma = av1; |
} |
} |
*zeta=m; |
*zeta = m; return m1; |
return m1; |
|
} |
} |
/* resoud x^l=a mod (p,T) |
/* Solve x^l = a mod (p,T) |
* l doit etre premier |
* l must be prime |
* q=p^degpol(T)-1 |
* q = p^degpol(T)-1 = (l^e)*r, with e>=1 and pgcd(r,l)=1 |
* q=(l^e)*r |
* m = y^(q/l) |
* e>=1 |
* y not an l-th power [ m != 1 ] */ |
* pgcd(r,l)=1 |
|
* m=y^(q/l) |
|
* y n'est pas une puissance l-ieme |
|
* m!=1 |
|
* ouf! |
|
*/ |
|
GEN |
GEN |
ffsqrtlmod(GEN a, GEN l, GEN T ,GEN p , GEN q, long e, GEN r, GEN y, GEN m) |
ffsqrtlmod(GEN a, GEN l, GEN T ,GEN p , GEN q, long e, GEN r, GEN y, GEN m) |
{ |
{ |
ulong av = avma,lim; |
gpmem_t av = avma, lim; |
long i,k; |
long i,k; |
GEN p1,p2,u1,u2,v,w,z; |
GEN p1,p2,u1,u2,v,w,z; |
|
|
bezout(r,l,&u1,&u2); |
if (gcmp1(a)) return gcopy(a); |
v=FpXQ_pow(a,u2,T,p); |
|
w=FpXQ_pow(a,modii(mulii(negi(u1),r),q),T,p); |
(void)bezout(r,l,&u1,&u2); /* result is 1 */ |
|
v = FpXQ_pow(a,u2,T,p); |
|
w = FpXQ_pow(a, modii(mulii(negi(u1),r),q), T,p); |
lim = stack_lim(av,1); |
lim = stack_lim(av,1); |
while (!gcmp1(w)) |
while (!gcmp1(w)) |
{ |
{ |
/* if p is not prime, next loop will not end */ |
k = 0; |
k=0; |
p1 = w; |
p1=w; |
do { /* if p is not prime, loop will not end */ |
do |
z = p1; |
{ |
p1 = FpXQ_pow(p1,l,T,p); |
z=p1; |
|
p1=FpXQ_pow(p1,l,T,p); |
|
k++; |
k++; |
}while(!gcmp1(p1)); |
} while (!gcmp1(p1)); |
if (k==e) { avma=av; return NULL; } |
if (k==e) { avma=av; return NULL; } |
p2 = FpXQ_mul(z,m,T,p); |
p2 = FpXQ_mul(z,m,T,p); |
for(i=1; !gcmp1(p2); i++) p2 = FpXQ_mul(p2,m,T,p);/*should be a baby step |
for (i=1; !gcmp1(p2); i++) p2 = FpXQ_mul(p2,m,T,p);/*TODO: BS/GS instead*/ |
giant step instead*/ |
p1= FpXQ_pow(y, modii(mulsi(i,gpowgs(l,e-k-1)), q), T,p); |
p1= FpXQ_pow(y,modii(mulsi(i,gpowgs(l,e-k-1)),q),T,p); |
|
m = FpXQ_pow(m,stoi(i),T,p); |
m = FpXQ_pow(m,stoi(i),T,p); |
e = k; |
e = k; |
v = FpXQ_mul(p1,v,T,p); |
v = FpXQ_mul(p1,v,T,p); |
Line 1436 ffsqrtlmod(GEN a, GEN l, GEN T ,GEN p , GEN q, long e, |
|
Line 1630 ffsqrtlmod(GEN a, GEN l, GEN T ,GEN p , GEN q, long e, |
|
w = FpXQ_mul(y,w,T,p); |
w = FpXQ_mul(y,w,T,p); |
if (low_stack(lim, stack_lim(av,1))) |
if (low_stack(lim, stack_lim(av,1))) |
{ |
{ |
GEN *gptr[4]; |
|
if(DEBUGMEM>1) err(warnmem,"ffsqrtlmod"); |
if(DEBUGMEM>1) err(warnmem,"ffsqrtlmod"); |
gptr[0]=&y; gptr[1]=&v; gptr[2]=&w; gptr[3]=&m; |
gerepileall(av,4, &y,&v,&w,&m); |
gerepilemany(av,gptr,4); |
|
} |
} |
} |
} |
return gerepilecopy(av,v); |
return gerepilecopy(av, v); |
} |
} |
/* n is an integer, a is in Fp[X]/(T), p is prime, T is irreducible mod p |
/* Solve x^n = a mod p: n integer, a in Fp[X]/(T) [ p prime, T irred. mod p ] |
|
* |
return a solution of |
* 1) if no solution, return NULL and (if zetan != NULL) set zetan to gzero. |
|
* |
x^n=a mod p |
* 2) If there is a solution, there are exactly m=gcd(p-1,n) of them. |
|
* If zetan != NULL, it is set to a primitive mth root of unity so that the set |
1)If there is no solution return NULL and if zetan is not NULL set zetan to gzero. |
* of solutions is {x*zetan^k;k=0 to m-1} |
|
* |
2) If there is solution there are exactly m=gcd(p-1,n) of them. |
* If a = 0, return 0 and (if zetan != NULL) set zetan = gun */ |
|
|
If zetan is not NULL, zetan is set to a primitive mth root of unity so that |
|
the set of solutions is {x*zetan^k;k=0 to m-1} |
|
|
|
If a=0 ,return 0 and if zetan is not NULL zetan is set to gun |
|
*/ |
|
GEN ffsqrtnmod(GEN a, GEN n, GEN T, GEN p, GEN *zetan) |
GEN ffsqrtnmod(GEN a, GEN n, GEN T, GEN p, GEN *zetan) |
{ |
{ |
ulong ltop=avma,lbot=0,av1,lim; |
gpmem_t ltop=avma, av1, lim; |
long i,j,e; |
long i,j,e; |
GEN m,u1,u2; |
GEN m,u1,u2; |
GEN q,r,zeta,y,l,z; |
GEN q,r,zeta,y,l,z; |
GEN *gptr[2]; |
|
if (typ(a) != t_POL || typ(n) != t_INT || typ(T) != t_POL || typ(p)!=t_INT) |
if (typ(a) != t_POL || typ(n) != t_INT || typ(T) != t_POL || typ(p)!=t_INT) |
err(typeer,"ffsqrtnmod"); |
err(typeer,"ffsqrtnmod"); |
if (lgef(T)==3) |
if (!degpol(T)) err(constpoler,"ffsqrtnmod"); |
err(constpoler,"ffsqrtnmod"); |
if (!signe(n)) err(talker,"1/0 exponent in ffsqrtnmod"); |
if(!signe(n)) |
if (gcmp1(n)) {if (zetan) *zetan=gun;return gcopy(a);} |
err(talker,"1/0 exponent in ffsqrtnmod"); |
if (gcmp0(a)) {if (zetan) *zetan=gun;return gzero;} |
if(gcmp1(n)) {if (zetan) *zetan=gun;return gcopy(a);} |
|
if(gcmp0(a)) {if (zetan) *zetan=gun;return gzero;} |
q = addsi(-1, gpowgs(p,degpol(T))); |
q=addsi(-1,gpowgs(p,degpol(T))); |
m = bezout(n,q,&u1,&u2); |
m=bezout(n,q,&u1,&u2); |
if (!egalii(m,n)) a = FpXQ_pow(a, modii(u1,q), T,p); |
if (gcmp(m,n)) |
if (zetan) z = polun[varn(T)]; |
{ |
lim = stack_lim(ltop,1); |
GEN b=modii(u1,q); |
|
lbot=avma; |
|
a=FpXQ_pow(a,b,T,p); |
|
} |
|
if (zetan) z=polun[varn(T)]; |
|
lim=stack_lim(ltop,1); |
|
if (!gcmp1(m)) |
if (!gcmp1(m)) |
{ |
{ |
m=decomp(m); |
m = decomp(m); av1 = avma; |
av1=avma; |
|
for (i = lg(m[1])-1; i; i--) |
for (i = lg(m[1])-1; i; i--) |
{ |
{ |
l=gcoeff(m,i,1); j=itos(gcoeff(m,i,2)); |
l = gcoeff(m,i,1); |
e=pvaluation(q,l,&r); |
j = itos(gcoeff(m,i,2)); |
y=fflgen(l,e,r,T,p,&zeta); |
e = pvaluation(q,l,&r); |
if (zetan) z=FpXQ_mul(z,FpXQ_pow(y,gpowgs(l,e-j),T,p),T,p); |
if(DEBUGLEVEL>=6) (void)timer2(); |
do |
y = fflgen(l,e,r,T,p,&zeta); |
|
if(DEBUGLEVEL>=6) msgtimer("fflgen"); |
|
if (zetan) z = FpXQ_mul(z, FpXQ_pow(y,gpowgs(l,e-j),T,p), T,p); |
|
for (; j; j--) |
{ |
{ |
lbot=avma; |
a = ffsqrtlmod(a,l,T,p,q,e,r,y,zeta); |
a=ffsqrtlmod(a,l,T,p,q,e,r,y,zeta); |
if (!a) {avma=ltop; return NULL;} |
if (!a){avma=ltop;return NULL;} |
} |
j--; |
|
}while (j); |
|
if (low_stack(lim, stack_lim(ltop,1))) |
if (low_stack(lim, stack_lim(ltop,1))) |
/* n can have lots of prime factors*/ |
{ /* n can have lots of prime factors */ |
{ |
|
if(DEBUGMEM>1) err(warnmem,"ffsqrtnmod"); |
if(DEBUGMEM>1) err(warnmem,"ffsqrtnmod"); |
if (zetan) |
gerepileall(av1,zetan? 2: 1, &a,&z); |
{ |
|
z=gcopy(z); |
|
gptr[0]=&a;gptr[1]=&z; |
|
gerepilemanysp(av1,lbot,gptr,2); |
|
} |
|
else |
|
a=gerepileupto(av1,a); |
|
lbot=av1; |
|
} |
} |
} |
} |
} |
} |
if (zetan) |
if (zetan) |
{ |
{ |
*zetan=gcopy(z); |
*zetan = z; |
gptr[0]=&a;gptr[1]=zetan; |
gerepileall(ltop,2,&a,zetan); |
gerepilemanysp(ltop,lbot,gptr,2); |
|
} |
} |
else |
else |
a=gerepileupto(ltop,a); |
a = gerepileupto(ltop, a); |
return a; |
return a; |
} |
} |
/*******************************************************************/ |
/*******************************************************************/ |
Line 1541 matrixpow(long n, long m, GEN y, GEN P,GEN l) |
|
Line 1711 matrixpow(long n, long m, GEN y, GEN P,GEN l) |
|
GEN |
GEN |
Fp_inv_isom(GEN S,GEN T, GEN p) |
Fp_inv_isom(GEN S,GEN T, GEN p) |
{ |
{ |
ulong ltop = avma, lbot; |
gpmem_t lbot, ltop = avma; |
GEN M, V; |
GEN M, V; |
int n, i; |
int n, i; |
long x; |
long x; |
Line 1557 Fp_inv_isom(GEN S,GEN T, GEN p) |
|
Line 1727 Fp_inv_isom(GEN S,GEN T, GEN p) |
|
V = gtopolyrev(V, x); |
V = gtopolyrev(V, x); |
return gerepile(ltop, lbot, V); |
return gerepile(ltop, lbot, V); |
} |
} |
#if 0 |
|
/*Old, trivial, implementation*/ |
/* Let M the matrix of the x^p Frobenius automorphism. |
GEN |
* Compute x^(p^i) for i=0..r |
intersect_ker(GEN P, GEN MA, GEN l, GEN U, GEN lU) |
|
{ |
|
ulong ltop=avma; |
|
long vp=varn(P), np=degpol(P); |
|
GEN A; |
|
A=FqM_ker(gaddmat(gneg(polx[MAXVARN]), *MA),lU,l); |
|
if (lg(A)!=2) |
|
err(talker,"ZZ_%Z[%Z]/(%Z) is not a field in Fp_intersect" |
|
,l,polx[vp],P); |
|
A=gmul((GEN)A[1],gmodulcp(gmodulcp(gun,l),U)); |
|
return gerepileupto(ltop,gtopolyrev(A,vp)); |
|
} |
|
#endif |
|
/* Let P a polynomial != 0 and M the matrix of the x^p Frobenius automorphism in |
|
* FFp[X]/T. Compute P(M) |
|
* not stack clean |
|
*/ |
*/ |
static GEN |
static GEN |
polfrobenius(GEN M, GEN p, long r, long v) |
polfrobenius(GEN M, GEN p, long r, long v) |
Line 1585 polfrobenius(GEN M, GEN p, long r, long v) |
|
Line 1739 polfrobenius(GEN M, GEN p, long r, long v) |
|
V = cgetg(r+2,t_VEC); |
V = cgetg(r+2,t_VEC); |
V[1] = (long) polx[v]; |
V[1] = (long) polx[v]; |
if (r == 0) return V; |
if (r == 0) return V; |
V[2] = (long) gtopolyrev((GEN)M[2],v); |
V[2] = (long) vec_to_pol((GEN)M[2],v); |
W = (GEN) M[2]; |
W = (GEN) M[2]; |
for (i = 3; i <= r+1; ++i) |
for (i = 3; i <= r+1; ++i) |
{ |
{ |
W = FpM_FpV_mul(M,W,p); |
W = FpM_FpV_mul(M,W,p); |
V[i] = (long) gtopolyrev(W,v); |
V[i] = (long) vec_to_pol(W,v); |
} |
} |
return V; |
return V; |
} |
} |
|
|
|
/* Let P a polynomial != 0 and M the matrix of the x^p Frobenius automorphism in |
|
* FFp[X]/T. Compute P(M) |
|
* V=polfrobenius(M, p, degpol(P), v) |
|
* not stack clean |
|
*/ |
|
|
static GEN |
static GEN |
matpolfrobenius(GEN V, GEN P, GEN T, GEN p) |
matpolfrobenius(GEN V, GEN P, GEN T, GEN p) |
{ |
{ |
|
gpmem_t btop; |
long i; |
long i; |
long l=degpol(T); |
long l=degpol(T); |
long v=varn(T); |
long v=varn(T); |
GEN M,W; |
GEN M,W,Mi; |
GEN PV=gtovec(P); |
GEN PV=gtovec(P); |
|
GEN *gptr[2]; |
|
long lV=lg(V); |
PV=cgetg(degpol(P)+2,t_VEC); |
PV=cgetg(degpol(P)+2,t_VEC); |
for(i=1;i<lg(PV);i++) |
for(i=1;i<lg(PV);i++) |
PV[i]=P[1+i]; |
PV[i]=P[1+i]; |
M=cgetg(l+1,t_VEC); |
M=cgetg(l+1,t_VEC); |
M[1]=(long)scalarpol(poleval(P,gun),v); |
M[1]=(long)scalarpol(poleval(P,gun),v); |
M[2]=(long)FpXV_FpV_dotproduct(V,PV,p); |
M[2]=(long)FpXV_FpV_innerprod(V,PV,p); |
W=cgetg(lg(V),t_VEC); |
btop=avma; |
for(i=1;i<lg(W);i++) |
gptr[0]=&Mi; |
|
gptr[1]=&W; |
|
W=cgetg(lV,t_VEC); |
|
for(i=1;i<lV;i++) |
W[i]=V[i]; |
W[i]=V[i]; |
for(i=3;i<=l;i++) |
for(i=3;i<=l;i++) |
{ |
{ |
long j; |
long j; |
for(j=1;j<lg(W);j++) |
gpmem_t bbot; |
W[j]=(long)FpXQ_mul((GEN)W[j],(GEN)V[j],T,p); |
GEN W2=cgetg(lV,t_VEC); |
M[i]=(long)FpXV_FpV_dotproduct(W,PV,p); |
for(j=1;j<lV;j++) |
|
W2[j]=(long)FpXQ_mul((GEN)W[j],(GEN)V[j],T,p); |
|
bbot=avma; |
|
Mi=FpXV_FpV_innerprod(W2,PV,p); |
|
W=gcopy(W2); |
|
gerepilemanysp(btop,bbot,gptr,2); |
|
btop=(gpmem_t)W; |
|
M[i]=(long)Mi; |
} |
} |
return vecpol_to_mat(M,l); |
return vecpol_to_mat(M,l); |
} |
} |
|
|
|
/* Let M the matrix of the Frobenius automorphism of Fp[X]/(T). |
|
* Compute M^d |
|
* TODO: use left-right binary (tricky!) |
|
*/ |
|
GEN |
|
FpM_frobenius_pow(GEN M, long d, GEN T, GEN p) |
|
{ |
|
gpmem_t ltop=avma; |
|
long i,l=degpol(T); |
|
GEN R, W = (GEN) M[2]; |
|
for (i = 2; i <= d; ++i) |
|
W = FpM_FpV_mul(M,W,p); |
|
R=matrixpow(l,l,vec_to_pol(W,varn(T)),T,p); |
|
return gerepilecopy(ltop,R); |
|
} |
|
|
/* Essentially we want to compute |
/* Essentially we want to compute |
* FqM_ker(MA-polx[MAXVARN],lU,l) |
* FqM_ker(MA-polx[MAXVARN],U,l) |
* To avoid use of matrix in Fq we procede as follows: |
* To avoid use of matrix in Fq we procede as follows: |
* We compute FpM_ker(U(MA)) and then we recover |
* We compute FpM_ker(U(MA),l) and then we recover |
* the eigen value by Galois action, see formula. |
* the eigen value by Galois action, see formula. |
*/ |
*/ |
static GEN |
static GEN |
intersect_ker(GEN P, GEN MA, GEN l, GEN U, GEN lU) |
intersect_ker(GEN P, GEN MA, GEN U, GEN l) |
{ |
{ |
ulong ltop=avma; |
gpmem_t ltop=avma; |
long vp=varn(P); |
long vp=varn(P); |
long vu=varn(lU), r=degpol(lU); |
long vu=varn(U), r=degpol(U); |
long i; |
long i; |
GEN A, R, M, ib0, V; |
GEN A, R, M, ib0, V; |
if (DEBUGLEVEL>=4) timer2(); |
if (DEBUGLEVEL>=4) (void)timer2(); |
V=polfrobenius(MA,l,r,varn(U)); |
V=polfrobenius(MA,l,r,vu); |
if (DEBUGLEVEL>=4) msgtimer("pol[frobenius]"); |
if (DEBUGLEVEL>=4) msgtimer("pol[frobenius]"); |
M=matpolfrobenius(V,lU,P,l); |
M=matpolfrobenius(V,U,P,l); |
if (DEBUGLEVEL>=4) msgtimer("matrix cyclo"); |
if (DEBUGLEVEL>=4) msgtimer("matrix cyclo"); |
A=FpM_ker(M,l); |
A=FpM_ker(M,l); |
if (DEBUGLEVEL>=4) msgtimer("kernel"); |
if (DEBUGLEVEL>=4) msgtimer("kernel"); |
Line 1652 intersect_ker(GEN P, GEN MA, GEN l, GEN U, GEN lU) |
|
Line 1841 intersect_ker(GEN P, GEN MA, GEN l, GEN U, GEN lU) |
|
* a_{i-1}=\phi(a_i)+b_ia_{r-1} i=r-1 to 1 |
* a_{i-1}=\phi(a_i)+b_ia_{r-1} i=r-1 to 1 |
* Where a_0=A[1] and b_i=U[i+2] |
* Where a_0=A[1] and b_i=U[i+2] |
*/ |
*/ |
ib0=negi(mpinvmod((GEN)lU[2],l)); |
ib0=negi(mpinvmod((GEN)U[2],l)); |
R=cgetg(r+1,t_MAT); |
R=cgetg(r+1,t_MAT); |
R[1]=A[1]; |
R[1]=A[1]; |
R[r]=(long)FpM_FpV_mul(MA,gmul((GEN)A[1],ib0),l); |
R[r]=(long)FpM_FpV_mul(MA,gmul((GEN)A[1],ib0),l); |
for(i=r-1;i>1;i--) |
for(i=r-1;i>1;i--) |
R[i]=(long)FpV_red(gadd(FpM_FpV_mul(MA,(GEN)R[i+1],l), |
R[i]=(long)FpV_red(gadd(FpM_FpV_mul(MA,(GEN)R[i+1],l), |
gmul((GEN)lU[i+2],(GEN)R[r])),l); |
gmul((GEN)U[i+2],(GEN)R[r])),l); |
R=gtrans_i(R); |
R=gtrans_i(R); |
for(i=1;i<lg(R);i++) |
for(i=1;i<lg(R);i++) |
R[i]=(long)gtopolyrev((GEN)R[i],vu); |
R[i]=(long)vec_to_pol((GEN)R[i],vu); |
A=gtopolyrev(R,vp); |
A=gtopolyrev(R,vp); |
A=gmul(A,gmodulcp(gmodulcp(gun,l),U)); |
|
return gerepileupto(ltop,A); |
return gerepileupto(ltop,A); |
} |
} |
|
|
Line 1680 intersect_ker(GEN P, GEN MA, GEN l, GEN U, GEN lU) |
|
Line 1868 intersect_ker(GEN P, GEN MA, GEN l, GEN U, GEN lU) |
|
void |
void |
Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN *SQ, GEN MA, GEN MB) |
Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN *SQ, GEN MA, GEN MB) |
{ |
{ |
ulong ltop=avma,lbot; |
gpmem_t lbot, ltop=avma; |
long vp,vq,np,nq,e,pg; |
long vp,vq,np,nq,e,pg; |
GEN q; |
GEN q; |
GEN A,B,Ap,Bp; |
GEN A,B,Ap,Bp; |
Line 1692 Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN |
|
Line 1880 Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN |
|
e=pvaluation(stoi(n),l,&q); |
e=pvaluation(stoi(n),l,&q); |
pg=itos(q); |
pg=itos(q); |
avma=ltop; |
avma=ltop; |
if (DEBUGLEVEL>=4) timer2(); |
|
if(!MA) MA=matrixpow(np,np,FpXQ_pow(polx[vp],l,P,l),P,l); |
if(!MA) MA=matrixpow(np,np,FpXQ_pow(polx[vp],l,P,l),P,l); |
if(!MB) MB=matrixpow(nq,nq,FpXQ_pow(polx[vq],l,Q,l),Q,l); |
if(!MB) MB=matrixpow(nq,nq,FpXQ_pow(polx[vq],l,Q,l),Q,l); |
if (DEBUGLEVEL>=4) msgtimer("matrixpow"); |
|
A=Ap=zeropol(vp); |
A=Ap=zeropol(vp); |
B=Bp=zeropol(vq); |
B=Bp=zeropol(vq); |
if (pg>1) |
if (pg > 1) |
{ |
{ |
if (gcmp0(modis(addis(l,-1),pg))) |
if (smodis(l,pg) == 1) |
/*We do not need to use relative extension in this setting, so |
/*We do not need to use relative extension in this setting, so |
we don't. (Well,now that we don't in the other case also, it is more |
we don't. (Well,now that we don't in the other case also, it is more |
dubious to treat cases apart...)*/ |
dubious to treat cases apart...)*/ |
Line 1710 Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN |
|
Line 1896 Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN |
|
if (lg(z)<2) err(talker,"%Z is not a prime in Fp_intersect",l); |
if (lg(z)<2) err(talker,"%Z is not a prime in Fp_intersect",l); |
z=negi(lift((GEN)z[1])); |
z=negi(lift((GEN)z[1])); |
ipg=stoi(pg); |
ipg=stoi(pg); |
if (DEBUGLEVEL>=4) timer2(); |
if (DEBUGLEVEL>=4) (void)timer2(); |
A=FpM_ker(gaddmat(z, MA),l); |
A=FpM_ker(gaddmat(z, MA),l); |
if (lg(A)!=2) |
if (lg(A)!=2) |
err(talker,"ZZ_%Z[%Z]/(%Z) is not a field in Fp_intersect" |
err(talker,"ZZ_%Z[%Z]/(%Z) is not a field in Fp_intersect" |
,l,polx[vp],P); |
,l,polx[vp],P); |
A=gtopolyrev((GEN)A[1],vp); |
A=vec_to_pol((GEN)A[1],vp); |
B=FpM_ker(gaddmat(z, MB),l); |
B=FpM_ker(gaddmat(z, MB),l); |
if (lg(B)!=2) |
if (lg(B)!=2) |
err(talker,"ZZ_%Z[%Z]/(%Z) is not a field in Fp_intersect" |
err(talker,"ZZ_%Z[%Z]/(%Z) is not a field in Fp_intersect" |
,l,polx[vq],Q); |
,l,polx[vq],Q); |
B=gtopolyrev((GEN)B[1],vq); |
B=vec_to_pol((GEN)B[1],vq); |
if (DEBUGLEVEL>=4) msgtimer("FpM_ker"); |
if (DEBUGLEVEL>=4) msgtimer("FpM_ker"); |
An=(GEN) FpXQ_pow(A,ipg,P,l)[2]; |
An=(GEN) FpXQ_pow(A,ipg,P,l)[2]; |
Bn=(GEN) FpXQ_pow(B,ipg,Q,l)[2]; |
Bn=(GEN) FpXQ_pow(B,ipg,Q,l)[2]; |
z=modii(mulii(An,mpinvmod(Bn,l)),l); |
if (!invmod(Bn,l,&z)) |
|
err(talker,"Polynomials not irreducible in Fp_intersect"); |
|
z=modii(mulii(An,z),l); |
L=mpsqrtnmod(z,ipg,l,NULL); |
L=mpsqrtnmod(z,ipg,l,NULL); |
if ( !L ) |
if ( !L ) |
err(talker,"Polynomials not irreducible in Fp_intersect"); |
err(talker,"Polynomials not irreducible in Fp_intersect"); |
Line 1733 Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN |
|
Line 1921 Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN |
|
} |
} |
else |
else |
{ |
{ |
GEN L,An,Bn,ipg,U,lU,z; |
GEN L,An,Bn,ipg,U,z; |
U=gmael(factmod(cyclo(pg,MAXVARN),l),1,1); |
U=lift(gmael(factmod(cyclo(pg,MAXVARN),l),1,1)); |
lU=lift(U); |
|
ipg=stoi(pg); |
ipg=stoi(pg); |
A=intersect_ker(P, MA, l, U, lU); |
A=intersect_ker(P, MA, U, l); |
B=intersect_ker(Q, MB, l, U, lU); |
B=intersect_ker(Q, MB, U, l); |
/*Somewhat ugly, but it is a proof that POLYMOD are useful and |
if (DEBUGLEVEL>=4) (void)timer2(); |
powerful.*/ |
An=(GEN) FpXQYQ_pow(A,stoi(pg),U,P,l)[2]; |
if (DEBUGLEVEL>=4) timer2(); |
Bn=(GEN) FpXQYQ_pow(B,stoi(pg),U,Q,l)[2]; |
An=lift(lift((GEN)lift(gpowgs(gmodulcp(A,P),pg))[2])); |
|
Bn=lift(lift((GEN)lift(gpowgs(gmodulcp(B,Q),pg))[2])); |
|
if (DEBUGLEVEL>=4) msgtimer("pows [P,Q]"); |
if (DEBUGLEVEL>=4) msgtimer("pows [P,Q]"); |
z=FpXQ_inv(Bn,lU,l); |
z=FpXQ_inv(Bn,U,l); |
z=FpXQ_mul(An,z,lU,l); |
z=FpXQ_mul(An,z,U,l); |
L=ffsqrtnmod(z,ipg,lU,l,NULL); |
L=ffsqrtnmod(z,ipg,U,l,NULL); |
if (DEBUGLEVEL>=4) msgtimer("ffsqrtn"); |
if (DEBUGLEVEL>=4) msgtimer("ffsqrtn"); |
if ( !L ) |
if ( !L ) |
err(talker,"Polynomials not irreducible in Fp_intersect"); |
err(talker,"Polynomials not irreducible in Fp_intersect"); |
B=gsubst(lift(lift(gmul(B,L))),MAXVARN,gzero); |
B=FpXQX_FpXQ_mul(B,L,U,l); |
A=gsubst(lift(lift(A)),MAXVARN,gzero); |
B=gsubst(B,MAXVARN,gzero); |
|
A=gsubst(A,MAXVARN,gzero); |
} |
} |
} |
} |
if (e!=0) |
if (e!=0) |
Line 1784 Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN |
|
Line 1970 Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN |
|
for(;i<=np;i++) VP[i]=zero; |
for(;i<=np;i++) VP[i]=zero; |
} |
} |
Ap=FpM_invimage(MA,VP,l); |
Ap=FpM_invimage(MA,VP,l); |
Ap=gtopolyrev(Ap,vp); |
Ap=vec_to_pol(Ap,vp); |
if (j) |
if (j) |
{ |
{ |
By=FpXQ_mul(By,FpXQ_pow(Bp,lmun,Q,l),Q,l); |
By=FpXQ_mul(By,FpXQ_pow(Bp,lmun,Q,l),Q,l); |
Line 1792 Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN |
|
Line 1978 Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN |
|
for(;i<=nq;i++) VQ[i]=zero; |
for(;i<=nq;i++) VQ[i]=zero; |
} |
} |
Bp=FpM_invimage(MB,VQ,l); |
Bp=FpM_invimage(MB,VQ,l); |
Bp=gtopolyrev(Bp,vq); |
Bp=vec_to_pol(Bp,vq); |
if (DEBUGLEVEL>=4) msgtimer("FpM_invimage"); |
if (DEBUGLEVEL>=4) msgtimer("FpM_invimage"); |
} |
} |
}/*FpX_add is not clean, so we must do it *before* lbot=avma*/ |
}/*FpX_add is not clean, so we must do it *before* lbot=avma*/ |
Line 1811 Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN |
|
Line 1997 Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN |
|
* isomorphism. */ |
* isomorphism. */ |
GEN Fp_isom(GEN P,GEN Q,GEN l) |
GEN Fp_isom(GEN P,GEN Q,GEN l) |
{ |
{ |
ulong ltop=avma; |
gpmem_t ltop=avma; |
GEN SP,SQ,R; |
GEN SP,SQ,R; |
Fp_intersect(degpol(P),P,Q,l,&SP,&SQ,NULL,NULL); |
Fp_intersect(degpol(P),P,Q,l,&SP,&SQ,NULL,NULL); |
R=Fp_inv_isom(SP,P,l); |
R=Fp_inv_isom(SP,P,l); |
Line 1819 GEN Fp_isom(GEN P,GEN Q,GEN l) |
|
Line 2005 GEN Fp_isom(GEN P,GEN Q,GEN l) |
|
return gerepileupto(ltop,R); |
return gerepileupto(ltop,R); |
} |
} |
GEN |
GEN |
Fp_factorgalois(GEN P,GEN l, long d, long w) |
Fp_factorgalois(GEN P,GEN l, long d, long w, GEN MP) |
{ |
{ |
ulong ltop=avma; |
gpmem_t ltop=avma; |
GEN R,V,ld,Tl; |
GEN R,V,Tl,z,M; |
long n,k,m; |
long n,k,m; |
long v; |
long v; |
v=varn(P); |
v=varn(P); |
n=degpol(P); |
n=degpol(P); |
m=n/d; |
m=n/d; |
ld=gpowgs(l,d); |
if (DEBUGLEVEL>=4) (void)timer2(); |
|
M=FpM_frobenius_pow(MP,d,P,l); |
|
if (DEBUGLEVEL>=4) msgtimer("Fp_factorgalois: frobenius power"); |
Tl=gcopy(P); setvarn(Tl,w); |
Tl=gcopy(P); setvarn(Tl,w); |
V=cgetg(m+1,t_VEC); |
V=cgetg(m+1,t_VEC); |
V[1]=lpolx[w]; |
V[1]=lpolx[w]; |
|
z=pol_to_vec((GEN)V[1],n); |
for(k=2;k<=m;k++) |
for(k=2;k<=m;k++) |
V[k]=(long)FpXQ_pow((GEN)V[k-1],ld,Tl,l); |
{ |
|
z=FpM_FpV_mul(M,z,l); |
|
V[k]=(long)vec_to_pol(z,w); |
|
} |
|
if (DEBUGLEVEL>=4) msgtimer("Fp_factorgalois: roots"); |
R=FqV_roots_to_pol(V,Tl,l,v); |
R=FqV_roots_to_pol(V,Tl,l,v); |
|
if (DEBUGLEVEL>=4) msgtimer("Fp_factorgalois: pol"); |
return gerepileupto(ltop,R); |
return gerepileupto(ltop,R); |
} |
} |
extern GEN mat_to_polpol(GEN x, long v,long w); |
|
extern GEN polpol_to_mat(GEN v, long n); |
|
/* P,Q irreducible over F_l. Factor P over FF_l[X] / Q [factors are ordered as |
/* P,Q irreducible over F_l. Factor P over FF_l[X] / Q [factors are ordered as |
* a Frobenius cycle] */ |
* a Frobenius cycle] */ |
GEN |
GEN |
Fp_factor_irred(GEN P,GEN l, GEN Q) |
Fp_factor_irred(GEN P,GEN l, GEN Q) |
{ |
{ |
ulong ltop=avma,av; |
gpmem_t ltop=avma, av; |
GEN SP,SQ,MP,MQ,M,MF,E,V,IR,res; |
GEN SP,SQ,MP,MQ,M,FP,FQ,E,V,IR,res; |
long np=degpol(P),nq=degpol(Q); |
long np=degpol(P),nq=degpol(Q); |
long i,d=cgcd(np,nq); |
long i,d=cgcd(np,nq); |
long vp=varn(P),vq=varn(Q); |
long vp=varn(P),vq=varn(Q); |
Line 1855 Fp_factor_irred(GEN P,GEN l, GEN Q) |
|
Line 2047 Fp_factor_irred(GEN P,GEN l, GEN Q) |
|
res[1]=lcopy(P); |
res[1]=lcopy(P); |
return res; |
return res; |
} |
} |
MF=matrixpow(nq,nq,FpXQ_pow(polx[vq],l,Q,l),Q,l); |
if (DEBUGLEVEL>=4) (void)timer2(); |
Fp_intersect(d,P,Q,l,&SP,&SQ,NULL,MF); |
FP=matrixpow(np,np,FpXQ_pow(polx[vp],l,P,l),P,l); |
|
FQ=matrixpow(nq,nq,FpXQ_pow(polx[vq],l,Q,l),Q,l); |
|
if (DEBUGLEVEL>=4) msgtimer("matrixpows"); |
|
Fp_intersect(d,P,Q,l,&SP,&SQ,FP,FQ); |
av=avma; |
av=avma; |
E=Fp_factorgalois(P,l,d,vq); |
E=Fp_factorgalois(P,l,d,vq,FP); |
E=polpol_to_mat(E,np); |
E=polpol_to_mat(E,np); |
MP = matrixpow(np,d,SP,P,l); |
MP = matrixpow(np,d,SP,P,l); |
IR = (GEN)FpM_sindexrank(MP,l)[1]; |
IR = (GEN)FpM_sindexrank(MP,l)[1]; |
Line 1872 Fp_factor_irred(GEN P,GEN l, GEN Q) |
|
Line 2067 Fp_factor_irred(GEN P,GEN l, GEN Q) |
|
V = cgetg(d+1,t_VEC); |
V = cgetg(d+1,t_VEC); |
V[1]=(long)M; |
V[1]=(long)M; |
for(i=2;i<=d;i++) |
for(i=2;i<=d;i++) |
V[i]=(long)FpM_mul(MF,(GEN)V[i-1],l); |
V[i]=(long)FpM_mul(FQ,(GEN)V[i-1],l); |
res=cgetg(d+1,t_COL); |
res=cgetg(d+1,t_COL); |
for(i=1;i<=d;i++) |
for(i=1;i<=d;i++) |
res[i]=(long)mat_to_polpol((GEN)V[i],vp,vq); |
res[i]=(long)mat_to_polpol((GEN)V[i],vp,vq); |
Line 1880 Fp_factor_irred(GEN P,GEN l, GEN Q) |
|
Line 2075 Fp_factor_irred(GEN P,GEN l, GEN Q) |
|
} |
} |
GEN Fp_factor_rel0(GEN P,GEN l, GEN Q) |
GEN Fp_factor_rel0(GEN P,GEN l, GEN Q) |
{ |
{ |
ulong ltop=avma,tetpil; |
gpmem_t ltop=avma, tetpil; |
GEN V,ex,F,y,R; |
GEN V,ex,F,y,R; |
long n,nbfact=0,nmax=lgef(P)-2; |
long n,nbfact=0,nmax=lgef(P)-2; |
long i; |
long i; |
Line 1909 GEN Fp_factor_rel0(GEN P,GEN l, GEN Q) |
|
Line 2104 GEN Fp_factor_rel0(GEN P,GEN l, GEN Q) |
|
} |
} |
GEN Fp_factor_rel(GEN P, GEN l, GEN Q) |
GEN Fp_factor_rel(GEN P, GEN l, GEN Q) |
{ |
{ |
long tetpil,av=avma; |
gpmem_t tetpil, av=avma; |
long nbfact; |
long nbfact; |
long j; |
long j; |
GEN y,u,v; |
GEN y,u,v; |
Line 2011 FpV_red(GEN z, GEN p) |
|
Line 2206 FpV_red(GEN z, GEN p) |
|
GEN |
GEN |
FpM_red(GEN z, GEN p) |
FpM_red(GEN z, GEN p) |
{ |
{ |
long i,j,l = lg(z), m = lg((GEN)z[1]); |
long i, l = lg(z); |
GEN x,y; |
GEN x = cgetg(l,t_MAT); |
x = cgetg(l,t_MAT); |
for (i=1; i<l; i++) x[i] = (long)FpV_red((GEN)z[i], p); |
for (i=1; i<l; i++) |
|
{ |
|
x[i]=lgetg(m,t_MAT);y=(GEN)x[i]; |
|
for(j=1; j<m ; j++) |
|
y[j] = lmodii(gmael(z,i,j),p); |
|
} |
|
return x; |
return x; |
} |
} |
|
|
Line 2042 FpXQX_normalize(GEN z, GEN T, GEN p) |
|
Line 2231 FpXQX_normalize(GEN z, GEN T, GEN p) |
|
} |
} |
|
|
GEN |
GEN |
|
small_to_vec(GEN z) |
|
{ |
|
long i, l = lg(z); |
|
GEN x = cgetg(l,t_VEC); |
|
for (i=1; i<l; i++) x[i] = lstoi(z[i]); |
|
return x; |
|
} |
|
|
|
GEN |
|
small_to_col(GEN z) |
|
{ |
|
long i, l = lg(z); |
|
GEN x = cgetg(l,t_COL); |
|
for (i=1; i<l; i++) x[i] = lstoi(z[i]); |
|
return x; |
|
} |
|
|
|
GEN |
small_to_mat(GEN z) |
small_to_mat(GEN z) |
{ |
{ |
long i, l = lg(z); |
long i, l = lg(z); |
GEN x = cgetg(l,t_MAT); |
GEN x = cgetg(l,t_MAT); |
for (i=1; i<l; i++) { x[i] = (long)gtovec((GEN)z[i]); settyp(x[i], t_COL); } |
for (i=1; i<l; i++) |
|
x[i] = (long)small_to_col((GEN)z[i]); |
return x; |
return x; |
} |
} |
|
|
Line 2066 small_to_pol(GEN z, long v) |
|
Line 2274 small_to_pol(GEN z, long v) |
|
GEN x = small_to_pol_i(z, lgef(z)); |
GEN x = small_to_pol_i(z, lgef(z)); |
setvarn(x,v); return x; |
setvarn(x,v); return x; |
} |
} |
|
/* same. In place */ |
|
GEN |
|
small_to_pol_ip(GEN z, long v) |
|
{ |
|
long i, l = lgef(z); |
|
for (i=2; i<l; i++) z[i] = lstoi(z[i]); |
|
settyp(z, t_POL); setvarn(z, v); return z; |
|
} |
|
|
GEN |
GEN |
pol_to_small(GEN x) |
pol_to_small(GEN x) |
{ |
{ |
long i, lx = lgef(x); |
long i, lx = lgef(x); |
GEN a = u_allocpol(lx-3, 0); |
GEN a = u_getpol(lx-3); |
for (i=2; i<lx; i++) a[i] = itos((GEN)x[i]); |
for (i=2; i<lx; i++) a[i] = itos((GEN)x[i]); |
return a; |
return a; |
} |
} |
|
|
/* z in Z[X,Y] representing an elt in F_p[X,Y] mod pol(Y) i.e F_q[X]) |
/* z in Z[X,Y] representing an elt in F_p[X,Y] mod T(Y) i.e F_q[X]) |
* in Kronecker form. Recover the "real" z, normalized |
* in Kronecker form. Recover the "real" z, normalized |
* If p = NULL, use generic functions and the coeff. ring implied by the |
* If p = NULL, use generic functions and the coeff. ring implied by the |
* coefficients of z */ |
* coefficients of z */ |
GEN |
GEN |
FqX_from_Kronecker(GEN z, GEN pol, GEN p) |
FqX_from_Kronecker(GEN z, GEN T, GEN p) |
{ |
{ |
long i,j,lx,l = lgef(z), N = (degpol(pol)<<1) + 1; |
long i,j,lx,l = lgef(z), N = (degpol(T)<<1) + 1; |
GEN a,x, t = cgetg(N,t_POL); |
GEN a,x, t = cgetg(N,t_POL); |
t[1] = pol[1] & VARNBITS; |
t[1] = T[1] & VARNBITS; |
lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL); |
lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL); |
if (isonstack(pol)) pol = gcopy(pol); |
if (isonstack(T)) T = gcopy(T); |
for (i=2; i<lx+2; i++) |
for (i=2; i<lx+2; i++) |
{ |
{ |
a = cgetg(3,t_POLMOD); x[i] = (long)a; |
a = cgetg(3,t_POLMOD); x[i] = (long)a; |
a[1] = (long)pol; |
a[1] = (long)T; |
for (j=2; j<N; j++) t[j] = z[j]; |
for (j=2; j<N; j++) t[j] = z[j]; |
z += (N-2); |
z += (N-2); |
a[2] = (long)FpX_res(normalizepol_i(t,N), pol,p); |
a[2] = (long)FpX_res(normalizepol_i(t,N), T,p); |
} |
} |
a = cgetg(3,t_POLMOD); x[i] = (long)a; |
a = cgetg(3,t_POLMOD); x[i] = (long)a; |
a[1] = (long)pol; |
a[1] = (long)T; |
N = (l-2) % (N-2) + 2; |
N = (l-2) % (N-2) + 2; |
for (j=2; j<N; j++) t[j] = z[j]; |
for (j=2; j<N; j++) t[j] = z[j]; |
a[2] = (long)FpX_res(normalizepol_i(t,N), pol,p); |
a[2] = (long)FpX_res(normalizepol_i(t,N), T,p); |
return normalizepol_i(x, i+1); |
return normalizepol_i(x, i+1); |
} |
} |
|
|
|
#if 0 |
|
/* z in Kronecker form representing a polynomial in FqX. Reduce mod (p,T) */ |
|
GEN |
|
FqX_Kronecker_red(GEN z, GEN T, GEN p) |
|
{ |
|
long i,j,lx,l = lgef(z), lT = lgef(T), N = ((dT-3)<<1) + 1; |
|
GEN a,x,y, t = cgetg(N,t_POL); |
|
t[1] = T[1] & VARNBITS; |
|
lx = (l-2) / (N-2); x = cgetg(lx+3,t_POL); |
|
x = y = z = FpX_red(z, p); |
|
for (i=2; i<lx+2; i++) |
|
{ |
|
for (j=2; j<N; j++) t[j] = z[j]; |
|
a = FpX_res(normalizepol_i(t,N), T,p); |
|
for (j=2; j<lT; j++) y[j] = a[j]; |
|
for ( ; j<N; j++) y[j] = zero; |
|
z += (N-2); |
|
y += (N-2); |
|
} |
|
N = (l-2) % (N-2) + 2; |
|
for (j=2; j<N; j++) t[j] = z[j]; |
|
a = FpX_res(normalizepol_i(t,N), T,p); |
|
for (j=2; j<lT; j++) y[j] = a[j]; |
|
for ( ; j<N; j++) y[j] = zero; |
|
return normalizepol_i(x, l); |
|
} |
|
#endif |
|
|
/* z in ?[X,Y] mod Q(Y) in Kronecker form ((subst(lift(z), x, y^(2deg(z)-1))) |
/* z in ?[X,Y] mod Q(Y) in Kronecker form ((subst(lift(z), x, y^(2deg(z)-1))) |
* Recover the "real" z, normalized */ |
* Recover the "real" z, normalized */ |
GEN |
GEN |
Line 2117 from_Kronecker(GEN z, GEN pol) |
|
Line 2361 from_Kronecker(GEN z, GEN pol) |
|
/* MODULAR GCD */ |
/* MODULAR GCD */ |
/* */ |
/* */ |
/*******************************************************************/ |
/*******************************************************************/ |
ulong xgcduu(ulong d, ulong d1, int f, ulong* v, ulong* v1, long *s); |
extern ulong xgcduu(ulong d, ulong d1, int f, ulong* v, ulong* v1, long *s); |
/* 1 / Mod(x,p) , or 0 if inverse doesn't exist */ |
/* 1 / Mod(x,p) , or 0 if inverse doesn't exist */ |
ulong |
ulong |
u_invmod(ulong x, ulong p) |
u_invmod(ulong x, ulong p) |
Line 2129 u_invmod(ulong x, ulong p) |
|
Line 2373 u_invmod(ulong x, ulong p) |
|
return xv; |
return xv; |
} |
} |
|
|
|
int |
|
u_val(ulong n, ulong p) |
|
{ |
|
ulong dummy; |
|
return svaluation(n,p,&dummy); |
|
} |
|
|
|
/* assume p^k is SMALL */ |
|
int |
|
u_pow(int p, int k) |
|
{ |
|
int i, pk; |
|
|
|
if (!k) return 1; |
|
if (p == 2) return 1<<k; |
|
pk = p; for (i=2; i<=k; i++) pk *= p; |
|
return pk; |
|
} |
|
|
#if 0 |
#if 0 |
static ulong |
static ulong |
umodratu(GEN a, ulong p) |
umodratu(GEN a, ulong p) |
Line 2145 umodratu(GEN a, ulong p) |
|
Line 2408 umodratu(GEN a, ulong p) |
|
|
|
/* return x[0 .. dx] mod p as t_VECSMALL. Assume x a t_POL/VECSMALL/INT */ |
/* return x[0 .. dx] mod p as t_VECSMALL. Assume x a t_POL/VECSMALL/INT */ |
GEN |
GEN |
u_Fp_FpX(GEN x, int malloc, ulong p) |
u_Fp_FpX(GEN x, ulong p) |
{ |
{ |
long i, lx; |
long i, lx; |
GEN a; |
GEN a; |
Line 2153 u_Fp_FpX(GEN x, int malloc, ulong p) |
|
Line 2416 u_Fp_FpX(GEN x, int malloc, ulong p) |
|
switch (typ(x)) |
switch (typ(x)) |
{ |
{ |
case t_VECSMALL: return x; |
case t_VECSMALL: return x; |
case t_INT: a = u_allocpol(0, malloc); |
case t_INT: a = u_getpol(0); |
a[2] = umodiu(x, p); return a; |
a[2] = umodiu(x, p); return a; |
} |
} |
lx = lgef(x); a = u_allocpol(lx-3, malloc); |
lx = lgef(x); a = u_getpol(lx-3); |
for (i=2; i<lx; i++) a[i] = (long)umodiu((GEN)x[i], p); |
for (i=2; i<lx; i++) a[i] = (long)umodiu((GEN)x[i], p); |
return u_normalizepol(a,lx); |
return u_normalizepol(a,lx); |
} |
} |
|
|
GEN |
GEN |
|
u_Fp_FpV(GEN x, ulong p) |
|
{ |
|
long i, n = lg(x); |
|
GEN y = cgetg(n,t_VECSMALL); |
|
for (i=1; i<n; i++) y[i] = (long)umodiu((GEN)x[i], p); |
|
return y; |
|
} |
|
|
|
GEN |
u_Fp_FpM(GEN x, ulong p) |
u_Fp_FpM(GEN x, ulong p) |
{ |
{ |
long i,j,m,n = lg(x); |
long j,n = lg(x); |
GEN y = cgetg(n,t_MAT); |
GEN y = cgetg(n,t_MAT); |
if (n == 1) return y; |
if (n == 1) return y; |
m = lg(x[1]); |
for (j=1; j<n; j++) y[j] = (long)u_Fp_FpV((GEN)x[j], p); |
for (j=1; j<n; j++) |
|
{ |
|
y[j] = (long)cgetg(m,t_VECSMALL); |
|
for (i=1; i<m; i++) coeff(y,i,j) = (long)umodiu(gcoeff(x,i,j), p); |
|
} |
|
return y; |
return y; |
} |
} |
|
|
static GEN |
static GEN |
u_FpX_Fp_mul(GEN y, ulong x,ulong p, int malloc) |
u_FpX_Fp_mul(GEN y, ulong x,ulong p) |
{ |
{ |
GEN z; |
GEN z; |
int i, l; |
int i, l; |
if (!x) return u_zeropol(malloc); |
if (!x) return u_zeropol(); |
l = lgef(y); z = u_allocpol(l-3, malloc); |
l = lgef(y); z = u_getpol(l-3); |
if (HIGHWORD(x | p)) |
if (HIGHWORD(x | p)) |
for(i=2; i<l; i++) z[i] = mulssmod(y[i], x, p); |
for(i=2; i<l; i++) z[i] = mulssmod(y[i], x, p); |
else |
else |
Line 2196 u_FpX_normalize(GEN z, ulong p) |
|
Line 2463 u_FpX_normalize(GEN z, ulong p) |
|
long l = lgef(z)-1; |
long l = lgef(z)-1; |
ulong p1 = z[l]; /* leading term */ |
ulong p1 = z[l]; /* leading term */ |
if (p1 == 1) return z; |
if (p1 == 1) return z; |
return u_FpX_Fp_mul(z, u_invmod(p1,p), p, 0); |
return u_FpX_Fp_mul(z, u_invmod(p1,p), p); |
} |
} |
|
|
static GEN |
static GEN |
u_copy(GEN x, int malloc) |
u_copy(GEN x) |
{ |
{ |
long i, l = lgef(x); |
long i, l = lgef(x); |
GEN z = u_allocpol(l-3, malloc); |
GEN z = u_getpol(l-3); |
for (i=2; i<l; i++) z[i] = x[i]; |
for (i=2; i<l; i++) z[i] = x[i]; |
return z; |
return z; |
} |
} |
|
|
/* as FpX_divres but working only on ulong types. ASSUME pr != ONLY_DIVIDES */ |
/* as FpX_divres but working only on ulong types. ASSUME pr != ONLY_DIVIDES |
|
* if relevant, *pr is the last object on stack */ |
GEN |
GEN |
u_FpX_divrem(GEN x, GEN y, ulong p, int malloc, GEN *pr) |
u_FpX_divrem(GEN x, GEN y, ulong p, GEN *pr) |
{ |
{ |
GEN z,q,c; |
GEN z,q,c; |
long dx,dy,dz,i,j; |
long dx,dy,dz,i,j; |
ulong p1,inv; |
ulong p1,inv; |
|
|
|
if (pr == ONLY_REM) return u_FpX_rem(x, y, p); |
dy = degpol(y); |
dy = degpol(y); |
if (!dy) |
if (!dy) |
{ |
{ |
if (pr) |
if (y[2] == 1UL) |
{ |
q = u_copy(x); |
if (pr == ONLY_REM) return u_zeropol(malloc); |
else |
*pr = u_zeropol(malloc); |
q = u_FpX_Fp_mul(x, u_invmod(y[2], p), p); |
} |
if (pr) *pr = u_zeropol(); |
if (y[2] == 1UL) return u_copy(x,malloc); |
return q; |
return u_FpX_Fp_mul(x, u_invmod(y[2], p), p, malloc); |
|
} |
} |
dx = degpol(x); |
dx = degpol(x); |
dz = dx-dy; |
dz = dx-dy; |
if (dz < 0) |
if (dz < 0) |
{ |
{ |
if (pr) |
q = u_zeropol(); |
{ |
if (pr) *pr = u_copy(x); |
c = u_copy(x, malloc); |
return q; |
if (pr == ONLY_REM) return c; |
|
*pr = c; |
|
} |
|
return u_zeropol(malloc); |
|
} |
} |
x += 2; |
x += 2; |
y += 2; |
y += 2; |
z = u_allocpol(dz, malloc || (pr == ONLY_REM)) + 2; |
z = u_getpol(dz) + 2; |
inv = y[dy]; |
inv = y[dy]; |
if (inv != 1UL) inv = u_invmod(inv,p); |
if (inv != 1UL) inv = u_invmod(inv,p); |
|
|
Line 2254 u_FpX_divrem(GEN x, GEN y, ulong p, int malloc, GEN *p |
|
Line 2518 u_FpX_divrem(GEN x, GEN y, ulong p, int malloc, GEN *p |
|
for (j=i-dy+1; j<=i && j<=dz; j++) |
for (j=i-dy+1; j<=i && j<=dz; j++) |
{ |
{ |
p1 += z[j]*y[i-j]; |
p1 += z[j]*y[i-j]; |
if (p1 & HIGHBIT) p1 = p1 % p; |
if (p1 & HIGHBIT) p1 %= p; |
} |
} |
p1 %= p; |
p1 %= p; |
z[i-dy] = p1? ((p - p1)*inv) % p: 0; |
z[i-dy] = p1? ((p - p1)*inv) % p: 0; |
Line 2274 u_FpX_divrem(GEN x, GEN y, ulong p, int malloc, GEN *p |
|
Line 2538 u_FpX_divrem(GEN x, GEN y, ulong p, int malloc, GEN *p |
|
q = u_normalizepol(z-2, dz+3); |
q = u_normalizepol(z-2, dz+3); |
if (!pr) return q; |
if (!pr) return q; |
|
|
c = u_allocpol(dy,malloc) + 2; |
c = u_getpol(dy) + 2; |
if (u_OK_ULONG(p)) |
if (u_OK_ULONG(p)) |
{ |
{ |
for (i=0; i<dy; i++) |
for (i=0; i<dy; i++) |
Line 2283 u_FpX_divrem(GEN x, GEN y, ulong p, int malloc, GEN *p |
|
Line 2547 u_FpX_divrem(GEN x, GEN y, ulong p, int malloc, GEN *p |
|
for (j=1; j<=i && j<=dz; j++) |
for (j=1; j<=i && j<=dz; j++) |
{ |
{ |
p1 += z[j]*y[i-j]; |
p1 += z[j]*y[i-j]; |
if (p1 & HIGHBIT) p1 = p1 % p; |
if (p1 & HIGHBIT) p1 %= p; |
} |
} |
c[i] = subssmod(x[i], p1%p, p); |
c[i] = subssmod(x[i], p1%p, p); |
} |
} |
Line 2300 u_FpX_divrem(GEN x, GEN y, ulong p, int malloc, GEN *p |
|
Line 2564 u_FpX_divrem(GEN x, GEN y, ulong p, int malloc, GEN *p |
|
} |
} |
i=dy-1; while (i>=0 && !c[i]) i--; |
i=dy-1; while (i>=0 && !c[i]) i--; |
c = u_normalizepol(c-2, i+3); |
c = u_normalizepol(c-2, i+3); |
if (pr == ONLY_REM) { free(q); return c; } |
|
*pr = c; return q; |
*pr = c; return q; |
} |
} |
|
|
|
/*FIXME: Unify the following 3 divrem routines. Treat the case x,y (lifted) in |
|
* R[X], y non constant. Given: (lifted) [inv(), mul()], (delayed) red() in R */ |
|
|
/* x and y in Z[X]. Possibly x in Z */ |
/* x and y in Z[X]. Possibly x in Z */ |
GEN |
GEN |
FpX_divres(GEN x, GEN y, GEN p, GEN *pr) |
FpX_divres(GEN x, GEN y, GEN p, GEN *pr) |
{ |
{ |
long vx,dx,dy,dz,i,j,av0,av,tetpil,sx,lrem; |
long vx, dx, dy, dz, i, j, sx, lrem; |
|
gpmem_t av0, av, tetpil; |
GEN z,p1,rem,lead; |
GEN z,p1,rem,lead; |
|
|
if (!p) return poldivres(x,y,pr); |
if (!p) return poldivres(x,y,pr); |
if (!signe(y)) err(talker,"division by zero in FpX_divres"); |
if (!signe(y)) err(talker,"division by zero in FpX_divres"); |
vx=varn(x); dy=degpol(y); dx=(typ(x)==t_INT)? 0: degpol(x); |
vx = varn(x); |
|
dy = degpol(y); |
|
dx = (typ(x)==t_INT)? 0: degpol(x); |
if (dx < dy) |
if (dx < dy) |
{ |
{ |
if (pr) |
if (pr) |
Line 2341 FpX_divres(GEN x, GEN y, GEN p, GEN *pr) |
|
Line 2610 FpX_divres(GEN x, GEN y, GEN p, GEN *pr) |
|
if (OK_ULONG(p)) |
if (OK_ULONG(p)) |
{ /* assume ab != 0 mod p */ |
{ /* assume ab != 0 mod p */ |
ulong pp = (ulong)p[2]; |
ulong pp = (ulong)p[2]; |
GEN a = u_Fp_FpX(x,1, pp); |
GEN a = u_Fp_FpX(x, pp); |
GEN b = u_Fp_FpX(y,1, pp); |
GEN b = u_Fp_FpX(y, pp); |
GEN zz= u_FpX_divrem(a,b,pp,1, pr); |
z = u_FpX_divrem(a,b,pp, pr); |
|
avma = av0; /* HACK: assume pr last on stack, then z */ |
|
setlg(z, lgef(z)); z = dummycopy(z); |
if (pr && pr != ONLY_DIVIDES && pr != ONLY_REM) |
if (pr && pr != ONLY_DIVIDES && pr != ONLY_REM) |
{ |
{ |
rem = small_to_pol(*pr,vx); |
setlg(*pr, lgef(*pr)); *pr = dummycopy(*pr); |
free(*pr); *pr = rem; |
*pr = small_to_pol_ip(*pr,vx); |
} |
} |
z = small_to_pol(zz,vx); |
return small_to_pol_ip(z,vx); |
free(zz); free(a); free(b); return z; |
|
} |
} |
lead = gcmp1(lead)? NULL: gclone(mpinvmod(lead,p)); |
lead = gcmp1(lead)? NULL: gclone(mpinvmod(lead,p)); |
avma = av0; |
avma = av0; |
Line 2370 FpX_divres(GEN x, GEN y, GEN p, GEN *pr) |
|
Line 2640 FpX_divres(GEN x, GEN y, GEN p, GEN *pr) |
|
} |
} |
if (!pr) { if (lead) gunclone(lead); return z-2; } |
if (!pr) { if (lead) gunclone(lead); return z-2; } |
|
|
rem = (GEN)avma; av = (long)new_chunk(dx+3); |
rem = (GEN)avma; av = (gpmem_t)new_chunk(dx+3); |
for (sx=0; ; i--) |
for (sx=0; ; i--) |
{ |
{ |
p1 = (GEN)x[i]; |
p1 = (GEN)x[i]; |
Line 2384 FpX_divres(GEN x, GEN y, GEN p, GEN *pr) |
|
Line 2654 FpX_divres(GEN x, GEN y, GEN p, GEN *pr) |
|
{ |
{ |
if (lead) gunclone(lead); |
if (lead) gunclone(lead); |
if (sx) { avma=av0; return NULL; } |
if (sx) { avma=av0; return NULL; } |
avma = (long)rem; return z-2; |
avma = (gpmem_t)rem; return z-2; |
} |
} |
lrem=i+3; rem -= lrem; |
lrem=i+3; rem -= lrem; |
rem[0]=evaltyp(t_POL) | evallg(lrem); |
rem[0]=evaltyp(t_POL) | evallg(lrem); |
rem[1]=evalsigne(1) | evalvarn(vx) | evallgef(lrem); |
rem[1]=evalsigne(1) | evalvarn(vx) | evallgef(lrem); |
p1 = gerepile((long)rem,tetpil,p1); |
p1 = gerepile((gpmem_t)rem,tetpil,p1); |
rem += 2; rem[i]=(long)p1; |
rem += 2; rem[i]=(long)p1; |
for (i--; i>=0; i--) |
for (i--; i>=0; i--) |
{ |
{ |
Line 2400 FpX_divres(GEN x, GEN y, GEN p, GEN *pr) |
|
Line 2670 FpX_divres(GEN x, GEN y, GEN p, GEN *pr) |
|
} |
} |
rem -= 2; |
rem -= 2; |
if (lead) gunclone(lead); |
if (lead) gunclone(lead); |
if (!sx) normalizepol_i(rem, lrem); |
if (!sx) (void)normalizepol_i(rem, lrem); |
if (pr == ONLY_REM) return gerepileupto(av0,rem); |
if (pr == ONLY_REM) return gerepileupto(av0,rem); |
*pr = rem; return z-2; |
*pr = rem; return z-2; |
} |
} |
Line 2409 FpX_divres(GEN x, GEN y, GEN p, GEN *pr) |
|
Line 2679 FpX_divres(GEN x, GEN y, GEN p, GEN *pr) |
|
GEN |
GEN |
FpXQX_divres(GEN x, GEN y, GEN T, GEN p, GEN *pr) |
FpXQX_divres(GEN x, GEN y, GEN T, GEN p, GEN *pr) |
{ |
{ |
long vx,dx,dy,dz,i,j,av0,av,tetpil,sx,lrem; |
long vx, dx, dy, dz, i, j, sx, lrem; |
|
gpmem_t av0, av, tetpil; |
GEN z,p1,rem,lead; |
GEN z,p1,rem,lead; |
|
|
if (!p) return poldivres(x,y,pr); |
if (!p) return poldivres(x,y,pr); |
Line 2443 FpXQX_divres(GEN x, GEN y, GEN T, GEN p, GEN *pr) |
|
Line 2714 FpXQX_divres(GEN x, GEN y, GEN T, GEN p, GEN *pr) |
|
#if 0 /* FIXME: to be done */ |
#if 0 /* FIXME: to be done */ |
if (OK_ULONG(p)) |
if (OK_ULONG(p)) |
{ /* assume ab != 0 mod p */ |
{ /* assume ab != 0 mod p */ |
ulong pp = (ulong)p[2]; |
|
GEN a = u_Fp_FpX(x,1, pp); |
|
GEN b = u_Fp_FpX(y,1, pp); |
|
GEN zz= u_FpX_divrem(a,b,pp,1, pr); |
|
if (pr && pr != ONLY_DIVIDES && pr != ONLY_REM) |
|
{ |
|
rem = small_to_pol(*pr,vx); |
|
free(*pr); *pr = rem; |
|
} |
|
z = small_to_pol(zz,vx); |
|
free(zz); free(a); free(b); return z; |
|
} |
} |
#endif |
#endif |
lead = gcmp1(lead)? NULL: gclone(FpXQ_inv(lead,T,p)); |
lead = gcmp1(lead)? NULL: gclone(FpXQ_inv(lead,T,p)); |
Line 2474 FpXQX_divres(GEN x, GEN y, GEN T, GEN p, GEN *pr) |
|
Line 2734 FpXQX_divres(GEN x, GEN y, GEN T, GEN p, GEN *pr) |
|
} |
} |
if (!pr) { if (lead) gunclone(lead); return z-2; } |
if (!pr) { if (lead) gunclone(lead); return z-2; } |
|
|
rem = (GEN)avma; av = (long)new_chunk(dx+3); |
rem = (GEN)avma; av = (gpmem_t)new_chunk(dx+3); |
for (sx=0; ; i--) |
for (sx=0; ; i--) |
{ |
{ |
p1 = (GEN)x[i]; |
p1 = (GEN)x[i]; |
Line 2488 FpXQX_divres(GEN x, GEN y, GEN T, GEN p, GEN *pr) |
|
Line 2748 FpXQX_divres(GEN x, GEN y, GEN T, GEN p, GEN *pr) |
|
{ |
{ |
if (lead) gunclone(lead); |
if (lead) gunclone(lead); |
if (sx) { avma=av0; return NULL; } |
if (sx) { avma=av0; return NULL; } |
avma = (long)rem; return z-2; |
avma = (gpmem_t)rem; return z-2; |
} |
} |
lrem=i+3; rem -= lrem; |
lrem=i+3; rem -= lrem; |
rem[0]=evaltyp(t_POL) | evallg(lrem); |
rem[0]=evaltyp(t_POL) | evallg(lrem); |
rem[1]=evalsigne(1) | evalvarn(vx) | evallgef(lrem); |
rem[1]=evalsigne(1) | evalvarn(vx) | evallgef(lrem); |
p1 = gerepile((long)rem,tetpil,p1); |
p1 = gerepile((gpmem_t)rem,tetpil,p1); |
rem += 2; rem[i]=(long)p1; |
rem += 2; rem[i]=(long)p1; |
for (i--; i>=0; i--) |
for (i--; i>=0; i--) |
{ |
{ |
Line 2504 FpXQX_divres(GEN x, GEN y, GEN T, GEN p, GEN *pr) |
|
Line 2764 FpXQX_divres(GEN x, GEN y, GEN T, GEN p, GEN *pr) |
|
} |
} |
rem -= 2; |
rem -= 2; |
if (lead) gunclone(lead); |
if (lead) gunclone(lead); |
if (!sx) normalizepol_i(rem, lrem); |
if (!sx) (void)normalizepol_i(rem, lrem); |
if (pr == ONLY_REM) return gerepileupto(av0,rem); |
if (pr == ONLY_REM) return gerepileupto(av0,rem); |
*pr = rem; return z-2; |
*pr = rem; return z-2; |
} |
} |
|
|
|
/* R = any base ring */ |
|
GEN |
|
RXQX_red(GEN P, GEN T) |
|
{ |
|
long i, l = lgef(P); |
|
GEN Q = cgetg(l, t_POL); |
|
Q[1] = P[1]; |
|
for (i=2; i<l; i++) Q[i] = lres((GEN)P[i], T); |
|
return Q; |
|
} |
|
|
|
/* R = any base ring */ |
|
GEN |
|
RXQV_red(GEN P, GEN T) |
|
{ |
|
long i, l = lg(P); |
|
GEN Q = cgetg(l, typ(P)); |
|
for (i=1; i<l; i++) Q[i] = lres((GEN)P[i], T); |
|
return Q; |
|
} |
|
|
|
/* x and y in (R[Y]/T)[X] (lifted), T in R[Y]. y preferably monic */ |
|
GEN |
|
RXQX_divrem(GEN x, GEN y, GEN T, GEN *pr) |
|
{ |
|
long vx, dx, dy, dz, i, j, sx, lrem; |
|
gpmem_t av0, av, tetpil; |
|
GEN z,p1,rem,lead; |
|
|
|
if (!signe(y)) err(talker,"division by zero in RXQX_divrem"); |
|
vx = varn(x); |
|
dx = degpol(x); |
|
dy = degpol(y); |
|
if (dx < dy) |
|
{ |
|
if (pr) |
|
{ |
|
av0 = avma; x = RXQX_red(x, T); |
|
if (pr == ONLY_DIVIDES) { avma=av0; return signe(x)? NULL: gzero; } |
|
if (pr == ONLY_REM) return x; |
|
*pr = x; |
|
} |
|
return zeropol(vx); |
|
} |
|
lead = leading_term(y); |
|
if (!dy) /* y is constant */ |
|
{ |
|
if (pr && pr != ONLY_DIVIDES) |
|
{ |
|
if (pr == ONLY_REM) return zeropol(vx); |
|
*pr = zeropol(vx); |
|
} |
|
if (gcmp1(lead)) return gcopy(x); |
|
av0 = avma; x = gmul(x, ginvmod(lead,T)); tetpil = avma; |
|
return gerepile(av0,tetpil,RXQX_red(x,T)); |
|
} |
|
av0 = avma; dz = dx-dy; |
|
lead = gcmp1(lead)? NULL: gclone(ginvmod(lead,T)); |
|
avma = av0; |
|
z = cgetg(dz+3,t_POL); |
|
z[1] = evalsigne(1) | evallgef(dz+3) | evalvarn(vx); |
|
x += 2; y += 2; z += 2; |
|
|
|
p1 = (GEN)x[dx]; av = avma; |
|
z[dz] = lead? lpileupto(av, gres(gmul(p1,lead), T)): lcopy(p1); |
|
for (i=dx-1; i>=dy; i--) |
|
{ |
|
av=avma; p1=(GEN)x[i]; |
|
for (j=i-dy+1; j<=i && j<=dz; j++) |
|
p1 = gsub(p1, gmul((GEN)z[j],(GEN)y[i-j])); |
|
if (lead) p1 = gmul(gres(p1, T), lead); |
|
tetpil=avma; z[i-dy] = lpile(av,tetpil, gres(p1, T)); |
|
} |
|
if (!pr) { if (lead) gunclone(lead); return z-2; } |
|
|
|
rem = (GEN)avma; av = (gpmem_t)new_chunk(dx+3); |
|
for (sx=0; ; i--) |
|
{ |
|
p1 = (GEN)x[i]; |
|
for (j=0; j<=i && j<=dz; j++) |
|
p1 = gsub(p1, gmul((GEN)z[j],(GEN)y[i-j])); |
|
tetpil=avma; p1 = gres(p1, T); if (signe(p1)) { sx = 1; break; } |
|
if (!i) break; |
|
avma=av; |
|
} |
|
if (pr == ONLY_DIVIDES) |
|
{ |
|
if (lead) gunclone(lead); |
|
if (sx) { avma=av0; return NULL; } |
|
avma = (gpmem_t)rem; return z-2; |
|
} |
|
lrem=i+3; rem -= lrem; |
|
rem[0]=evaltyp(t_POL) | evallg(lrem); |
|
rem[1]=evalsigne(1) | evalvarn(vx) | evallgef(lrem); |
|
p1 = gerepile((gpmem_t)rem,tetpil,p1); |
|
rem += 2; rem[i]=(long)p1; |
|
for (i--; i>=0; i--) |
|
{ |
|
av=avma; p1 = (GEN)x[i]; |
|
for (j=0; j<=i && j<=dz; j++) |
|
p1 = gsub(p1, gmul((GEN)z[j],(GEN)y[i-j])); |
|
tetpil=avma; rem[i]=lpile(av,tetpil, gres(p1, T)); |
|
} |
|
rem -= 2; |
|
if (lead) gunclone(lead); |
|
if (!sx) (void)normalizepol_i(rem, lrem); |
|
if (pr == ONLY_REM) return gerepileupto(av0,rem); |
|
*pr = rem; return z-2; |
|
} |
|
|
static GEN |
static GEN |
FpX_gcd_long(GEN x, GEN y, GEN p) |
FpX_gcd_long(GEN x, GEN y, GEN p) |
{ |
{ |
ulong pp = (ulong)p[2], av = avma; |
ulong pp = (ulong)p[2]; |
|
gpmem_t av = avma; |
GEN a,b; |
GEN a,b; |
|
|
(void)new_chunk((lgef(x) + lgef(y)) << 2); /* scratch space */ |
(void)new_chunk((lgef(x) + lgef(y)) << 2); /* scratch space */ |
a = u_Fp_FpX(x,0, pp); |
a = u_Fp_FpX(x, pp); |
if (!signe(a)) { avma = av; return FpX_red(y,p); } |
if (!signe(a)) { avma = av; return FpX_red(y,p); } |
b = u_Fp_FpX(y,0, pp); |
b = u_Fp_FpX(y, pp); |
a = u_FpX_gcd_i(a,b, pp); |
a = u_FpX_gcd_i(a,b, pp); |
avma = av; return small_to_pol(a, varn(x)); |
avma = av; return small_to_pol(a, varn(x)); |
} |
} |
|
|
FpX_gcd(GEN x, GEN y, GEN p) |
FpX_gcd(GEN x, GEN y, GEN p) |
{ |
{ |
GEN a,b,c; |
GEN a,b,c; |
long av0,av; |
gpmem_t av0, av; |
|
|
if (OK_ULONG(p)) return FpX_gcd_long(x,y,p); |
if (OK_ULONG(p)) return FpX_gcd_long(x,y,p); |
av0=avma; |
av0=avma; |
Line 2541 FpX_gcd(GEN x, GEN y, GEN p) |
|
Line 2912 FpX_gcd(GEN x, GEN y, GEN p) |
|
avma = av; return gerepileupto(av0, a); |
avma = av; return gerepileupto(av0, a); |
} |
} |
|
|
|
/*Return 1 if gcd can be computed |
|
* else return a factor of p*/ |
|
|
GEN |
GEN |
|
FpX_gcd_check(GEN x, GEN y, GEN p) |
|
{ |
|
GEN a,b,c; |
|
gpmem_t av=avma; |
|
|
|
a = FpX_red(x, p); |
|
b = FpX_red(y, p); |
|
while (signe(b)) |
|
{ |
|
GEN lead = leading_term(b); |
|
GEN g = mppgcd(lead,p); |
|
if (!is_pm1(g)) return gerepileupto(av,g); |
|
c = FpX_res(a,b,p); a=b; b=c; |
|
} |
|
avma = av; return gun; |
|
} |
|
|
|
GEN |
u_FpX_sub(GEN x, GEN y, ulong p) |
u_FpX_sub(GEN x, GEN y, ulong p) |
{ |
{ |
long i,lz,lx = lgef(x), ly = lgef(y); |
long i,lz,lx = lgef(x), ly = lgef(y); |
Line 2564 u_FpX_sub(GEN x, GEN y, ulong p) |
|
Line 2956 u_FpX_sub(GEN x, GEN y, ulong p) |
|
|
|
/* list of u_FpX in gptr, return then as GEN */ |
/* list of u_FpX in gptr, return then as GEN */ |
static void |
static void |
u_gerepilemany(long av, GEN* gptr[], long n, long v) |
u_gerepilemany(gpmem_t av, GEN* gptr[], long n, long v) |
{ |
{ |
GEN *l = (GEN*)gpmalloc(n*sizeof(GEN)); |
GEN *l = (GEN*)gpmalloc(n*sizeof(GEN)); |
long i; |
long i; |
Line 2586 u_FpX_extgcd(GEN a, GEN b, ulong p, GEN *ptu, GEN *ptv |
|
Line 2978 u_FpX_extgcd(GEN a, GEN b, ulong p, GEN *ptu, GEN *ptv |
|
{ |
{ |
GEN q,z,u,v, x = a, y = b; |
GEN q,z,u,v, x = a, y = b; |
|
|
u = u_zeropol(0); |
u = u_zeropol(); |
v= u_scalarpol(1, 0); /* v = 1 */ |
v= u_scalarpol(1); /* v = 1 */ |
while (signe(y)) |
while (signe(y)) |
{ |
{ |
q = u_FpX_divrem(x,y,p, 0,&z); |
q = u_FpX_divrem(x,y,p,&z); |
x = y; y = z; /* (x,y) = (y, x - q y) */ |
x = y; y = z; /* (x,y) = (y, x - q y) */ |
z = u_FpX_sub(u, u_FpX_mul(q,v, p), p); |
z = u_FpX_sub(u, u_FpX_mul(q,v, p), p); |
u = v; v = z; /* (u,v) = (v, u - q v) */ |
u = v; v = z; /* (u,v) = (v, u - q v) */ |
|
|
FpX_extgcd_long(GEN x, GEN y, GEN p, GEN *ptu, GEN *ptv) |
FpX_extgcd_long(GEN x, GEN y, GEN p, GEN *ptu, GEN *ptv) |
{ |
{ |
ulong pp = (ulong)p[2]; |
ulong pp = (ulong)p[2]; |
long av = avma; |
gpmem_t av = avma; |
GEN a, b, d; |
GEN a, b, d; |
|
|
a = u_Fp_FpX(x,0, pp); |
a = u_Fp_FpX(x, pp); |
b = u_Fp_FpX(y,0, pp); |
b = u_Fp_FpX(y, pp); |
d = u_FpX_extgcd(a,b, pp, ptu,ptv); |
d = u_FpX_extgcd(a,b, pp, ptu,ptv); |
{ |
{ |
GEN *gptr[3]; gptr[0] = &d; gptr[1] = ptu; gptr[2] = ptv; |
GEN *gptr[3]; gptr[0] = &d; gptr[1] = ptu; gptr[2] = ptv; |
Line 2620 FpX_extgcd_long(GEN x, GEN y, GEN p, GEN *ptu, GEN *pt |
|
Line 3012 FpX_extgcd_long(GEN x, GEN y, GEN p, GEN *ptu, GEN *pt |
|
|
|
/* x and y in Z[X], return lift(gcd(x mod p, y mod p)). Set u and v st |
/* x and y in Z[X], return lift(gcd(x mod p, y mod p)). Set u and v st |
* ux + vy = gcd (mod p) */ |
* ux + vy = gcd (mod p) */ |
|
/*TODO: Document the typ() of *ptu and *ptv*/ |
GEN |
GEN |
FpX_extgcd(GEN x, GEN y, GEN p, GEN *ptu, GEN *ptv) |
FpX_extgcd(GEN x, GEN y, GEN p, GEN *ptu, GEN *ptv) |
{ |
{ |
GEN a,b,q,r,u,v,d,d1,v1; |
GEN a,b,q,r,u,v,d,d1,v1; |
long ltop,lbot; |
gpmem_t ltop, lbot; |
|
|
if (OK_ULONG(p)) return FpX_extgcd_long(x,y,p,ptu,ptv); |
if (OK_ULONG(p)) return FpX_extgcd_long(x,y,p,ptu,ptv); |
ltop=avma; |
ltop=avma; |
|
|
FpXQX_extgcd(GEN x, GEN y, GEN T, GEN p, GEN *ptu, GEN *ptv) |
FpXQX_extgcd(GEN x, GEN y, GEN T, GEN p, GEN *ptu, GEN *ptv) |
{ |
{ |
GEN a,b,q,r,u,v,d,d1,v1; |
GEN a,b,q,r,u,v,d,d1,v1; |
long ltop,lbot; |
gpmem_t ltop, lbot; |
|
|
#if 0 /* FIXME To be done...*/ |
#if 0 /* FIXME To be done...*/ |
if (OK_ULONG(p)) return FpXQX_extgcd_long(x,y,T,p,ptu,ptv); |
if (OK_ULONG(p)) return FpXQX_extgcd_long(x,y,T,p,ptu,ptv); |
Line 2689 FpXQX_extgcd(GEN x, GEN y, GEN T, GEN p, GEN *ptu, GEN |
|
Line 3082 FpXQX_extgcd(GEN x, GEN y, GEN T, GEN p, GEN *ptu, GEN |
|
*ptu = u; *ptv = v; return d; |
*ptu = u; *ptv = v; return d; |
} |
} |
|
|
extern GEN caractducos(GEN p, GEN x, int v); |
/*x must be reduced*/ |
|
|
GEN |
GEN |
FpXQ_charpoly(GEN x, GEN T, GEN p) |
FpXQ_charpoly(GEN x, GEN T, GEN p) |
{ |
{ |
ulong ltop=avma; |
gpmem_t ltop=avma; |
GEN R=lift(caractducos(FpX(T,p),FpX(x,p),varn(T))); |
long v=varn(T); |
|
GEN R; |
|
T=gcopy(T); setvarn(T,MAXVARN); |
|
x=gcopy(x); setvarn(x,MAXVARN); |
|
R=FpY_FpXY_resultant(T,deg1pol(gun,FpX_neg(x,p),v),p); |
return gerepileupto(ltop,R); |
return gerepileupto(ltop,R); |
} |
} |
|
|
GEN |
GEN |
FpXQ_minpoly(GEN x, GEN T, GEN p) |
FpXQ_minpoly(GEN x, GEN T, GEN p) |
{ |
{ |
ulong ltop=avma; |
gpmem_t ltop=avma; |
GEN R=FpXQ_charpoly(x, T, p); |
GEN R=FpXQ_charpoly(x, T, p); |
GEN G=FpX_gcd(R,derivpol(R),p); |
GEN G=FpX_gcd(R,derivpol(R),p); |
G=FpX_normalize(G,p); |
G=FpX_normalize(G,p); |
Line 2714 FpXQ_minpoly(GEN x, GEN T, GEN p) |
|
Line 3110 FpXQ_minpoly(GEN x, GEN T, GEN p) |
|
static GEN |
static GEN |
u_chrem_coprime(GEN a, ulong b, GEN q, ulong p, ulong qinv, GEN pq) |
u_chrem_coprime(GEN a, ulong b, GEN q, ulong p, ulong qinv, GEN pq) |
{ |
{ |
ulong av = avma, d, amod = umodiu(a,p); |
ulong d, amod = umodiu(a, p); |
|
gpmem_t av = avma; |
GEN ax; |
GEN ax; |
|
|
if (b == amod) return NULL; |
if (b == amod) return NULL; |
Line 2725 u_chrem_coprime(GEN a, ulong b, GEN q, ulong p, ulong |
|
Line 3122 u_chrem_coprime(GEN a, ulong b, GEN q, ulong p, ulong |
|
} |
} |
|
|
/* centerlift(u mod p) */ |
/* centerlift(u mod p) */ |
GEN |
long |
u_center(ulong u, ulong p, ulong ps2) |
u_center(ulong u, ulong p, ulong ps2) |
{ |
{ |
return stoi((long) (u > ps2)? u - p: u); |
return (long) (u > ps2)? u - p: u; |
} |
} |
|
|
GEN |
GEN |
Line 2738 ZX_init_CRT(GEN Hp, ulong p, long v) |
|
Line 3135 ZX_init_CRT(GEN Hp, ulong p, long v) |
|
GEN H = cgetg(l, t_POL); |
GEN H = cgetg(l, t_POL); |
H[1] = evalsigne(1)|evallgef(l)|evalvarn(v); |
H[1] = evalsigne(1)|evallgef(l)|evalvarn(v); |
for (i=2; i<l; i++) |
for (i=2; i<l; i++) |
H[i] = (long)u_center(Hp[i], p, lim); |
H[i] = lstoi(u_center(Hp[i], p, lim)); |
return H; |
return H; |
} |
} |
|
|
Line 2753 ZM_init_CRT(GEN Hp, ulong p) |
|
Line 3150 ZM_init_CRT(GEN Hp, ulong p) |
|
cp = (GEN)Hp[j]; |
cp = (GEN)Hp[j]; |
c = cgetg(m, t_COL); |
c = cgetg(m, t_COL); |
H[j] = (long)c; |
H[j] = (long)c; |
for (i=1; i<l; i++) c[i] = (long)u_center(cp[i],p, lim); |
for (i=1; i<l; i++) c[i] = lstoi(u_center(cp[i],p, lim)); |
} |
} |
return H; |
return H; |
} |
} |
Line 2774 Z_incremental_CRT(GEN *H, ulong Hp, GEN q, GEN qp, ulo |
|
Line 3171 Z_incremental_CRT(GEN *H, ulong Hp, GEN q, GEN qp, ulo |
|
} |
} |
|
|
int |
int |
ZX_incremental_CRT(GEN H, GEN Hp, GEN q, GEN qp, ulong p) |
ZX_incremental_CRT(GEN *ptH, GEN Hp, GEN q, GEN qp, ulong p) |
{ |
{ |
GEN h, lim = shifti(qp,-1); |
GEN H = *ptH, h, lim = shifti(qp,-1); |
ulong qinv = u_invmod(umodiu(q,p), p); |
ulong qinv = u_invmod(umodiu(q,p), p); |
long i, l = lgef(H), lp = lgef(Hp); |
long i, l = lgef(H), lp = lgef(Hp); |
int stable = 1; |
int stable = 1; |
|
|
|
if (l < lp) |
|
{ /* degree increases */ |
|
GEN x = cgetg(lp, t_POL); |
|
for (i=1; i<l; i++) x[i] = H[i]; |
|
for ( ; i<lp; i++) |
|
{ |
|
h = stoi(Hp[i]); |
|
if (cmpii(h,lim) > 0) h = subii(h,qp); |
|
x[i] = (long)h; |
|
} |
|
setlgef(x,lp); *ptH = H = x; |
|
stable = 0; lp = l; |
|
} |
for (i=2; i<lp; i++) |
for (i=2; i<lp; i++) |
{ |
{ |
h = u_chrem_coprime((GEN)H[i],Hp[i],q,p,qinv,qp); |
h = u_chrem_coprime((GEN)H[i],Hp[i],q,p,qinv,qp); |
Line 2845 stopoly_gen(GEN m, GEN p, long v) |
|
Line 3256 stopoly_gen(GEN m, GEN p, long v) |
|
return y; |
return y; |
} |
} |
|
|
/* modular power */ |
|
ulong |
|
powuumod(ulong x, ulong n0, ulong p) |
|
{ |
|
ulong y, z, n; |
|
if (n0 <= 2) |
|
{ /* frequent special cases */ |
|
if (n0 == 2) return mulssmod(x,x,p); |
|
if (n0 == 1) return x; |
|
if (n0 == 0) return 1; |
|
} |
|
y = 1; z = x; n = n0; |
|
for(;;) |
|
{ |
|
if (n&1) y = mulssmod(y,z,p); |
|
n>>=1; if (!n) return y; |
|
z = mulssmod(z,z,p); |
|
} |
|
} |
|
|
|
/* separate from u_FpX_divrem for maximal speed. Implicitly malloc = 0 */ |
/* separate from u_FpX_divrem for maximal speed. Implicitly malloc = 0 */ |
GEN |
GEN |
u_FpX_rem(GEN x, GEN y, ulong p) |
u_FpX_rem(GEN x, GEN y, ulong p) |
Line 2873 u_FpX_rem(GEN x, GEN y, ulong p) |
|
Line 3264 u_FpX_rem(GEN x, GEN y, ulong p) |
|
long dx,dy,dz,i,j; |
long dx,dy,dz,i,j; |
ulong p1,inv; |
ulong p1,inv; |
|
|
dy = degpol(y); if (!dy) return u_zeropol(0); |
dy = degpol(y); if (!dy) return u_zeropol(); |
dx = degpol(x); |
dx = degpol(x); |
dz = dx-dy; if (dz < 0) return u_copy(x, 0); |
dz = dx-dy; if (dz < 0) return u_copy(x); |
x += 2; |
x += 2; |
y += 2; |
y += 2; |
z = u_allocpol(dz, 1) + 2; |
z = u_mallocpol(dz) + 2; |
inv = y[dy]; |
inv = y[dy]; |
if (inv != 1UL) inv = u_invmod(inv,p); |
if (inv != 1UL) inv = u_invmod(inv,p); |
|
|
c = u_allocpol(dy,0) + 2; |
c = u_getpol(dy) + 2; |
if (u_OK_ULONG(p)) |
if (u_OK_ULONG(p)) |
{ |
{ |
z[dz] = (inv*x[dx]) % p; |
z[dz] = (inv*x[dx]) % p; |
Line 2892 u_FpX_rem(GEN x, GEN y, ulong p) |
|
Line 3283 u_FpX_rem(GEN x, GEN y, ulong p) |
|
for (j=i-dy+1; j<=i && j<=dz; j++) |
for (j=i-dy+1; j<=i && j<=dz; j++) |
{ |
{ |
p1 += z[j]*y[i-j]; |
p1 += z[j]*y[i-j]; |
if (p1 & HIGHBIT) p1 = p1 % p; |
if (p1 & HIGHBIT) p1 %= p; |
} |
} |
p1 %= p; |
p1 %= p; |
z[i-dy] = p1? ((p - p1)*inv) % p: 0; |
z[i-dy] = p1? ((p - p1)*inv) % p: 0; |
Line 2903 u_FpX_rem(GEN x, GEN y, ulong p) |
|
Line 3294 u_FpX_rem(GEN x, GEN y, ulong p) |
|
for (j=1; j<=i && j<=dz; j++) |
for (j=1; j<=i && j<=dz; j++) |
{ |
{ |
p1 += z[j]*y[i-j]; |
p1 += z[j]*y[i-j]; |
if (p1 & HIGHBIT) p1 = p1 % p; |
if (p1 & HIGHBIT) p1 %= p; |
} |
} |
c[i] = subssmod(x[i], p1%p, p); |
c[i] = subssmod(x[i], p1%p, p); |
} |
} |
|
|
u_FpX_resultant(GEN a, GEN b, ulong p) |
u_FpX_resultant(GEN a, GEN b, ulong p) |
{ |
{ |
long da,db,dc,cnt; |
long da,db,dc,cnt; |
ulong lb,av, res = 1UL; |
ulong lb, res = 1UL; |
|
gpmem_t av; |
GEN c; |
GEN c; |
|
|
if (!signe(a) || !signe(b)) return 0; |
if (!signe(a) || !signe(b)) return 0; |
Line 2943 u_FpX_resultant(GEN a, GEN b, ulong p) |
|
Line 3335 u_FpX_resultant(GEN a, GEN b, ulong p) |
|
if (db > da) |
if (db > da) |
{ |
{ |
swapspec(a,b, da,db); |
swapspec(a,b, da,db); |
if (da & db & 1) res = p-res; |
if (both_odd(da,db)) res = p-res; |
} |
} |
if (!da) return 1; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */ |
if (!da) return 1; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */ |
cnt = 0; av = avma; |
cnt = 0; av = avma; |
Line 2954 u_FpX_resultant(GEN a, GEN b, ulong p) |
|
Line 3346 u_FpX_resultant(GEN a, GEN b, ulong p) |
|
a = b; b = c; dc = degpol(c); |
a = b; b = c; dc = degpol(c); |
if (dc < 0) { avma = av; return 0; } |
if (dc < 0) { avma = av; return 0; } |
|
|
if (da & db & 1) res = p - res; |
if (both_odd(da,db)) res = p - res; |
if (lb != 1) res = mulssmod(res, powuumod(lb, da - dc, p), p); |
if (lb != 1) res = mulssmod(res, powuumod(lb, da - dc, p), p); |
if (++cnt == 4) { cnt = 0; avma = av; } |
if (++cnt == 4) { cnt = 0; avma = av; } |
da = db; /* = degpol(a) */ |
da = db; /* = degpol(a) */ |
Line 2963 u_FpX_resultant(GEN a, GEN b, ulong p) |
|
Line 3355 u_FpX_resultant(GEN a, GEN b, ulong p) |
|
avma = av; return mulssmod(res, powuumod(b[2], da, p), p); |
avma = av; return mulssmod(res, powuumod(b[2], da, p), p); |
} |
} |
|
|
|
static GEN |
|
muliimod(GEN x, GEN y, GEN p) |
|
{ |
|
return modii(mulii(x,y), p); |
|
} |
|
|
|
#define FpX_rem(x,y,p) FpX_divres((x),(y),(p),ONLY_REM) |
|
|
|
/* Res(A,B) = Res(B,R) * lc(B)^(a-r) * (-1)^(ab+a+r), with R=A%B, a=deg(A) ...*/ |
|
GEN |
|
FpX_resultant(GEN a, GEN b, GEN p) |
|
{ |
|
long da,db,dc,cnt; |
|
gpmem_t av, lim; |
|
GEN c,lb, res = gun; |
|
|
|
if (!signe(a) || !signe(b)) return gzero; |
|
da = degpol(a); |
|
db = degpol(b); |
|
if (db > da) |
|
{ |
|
swapspec(a,b, da,db); |
|
if (both_odd(da,db)) res = subii(p, res); |
|
} |
|
if (!da) return gun; /* = res * a[2] ^ db, since 0 <= db <= da = 0 */ |
|
cnt = 0; av = avma; lim = stack_lim(av,2); |
|
while (db) |
|
{ |
|
lb = (GEN)b[db+2]; |
|
c = FpX_rem(a,b, p); |
|
a = b; b = c; dc = degpol(c); |
|
if (dc < 0) { avma = av; return 0; } |
|
|
|
if (both_odd(da,db)) res = subii(p, res); |
|
if (!gcmp1(lb)) res = muliimod(res, powgumod(lb, da - dc, p), p); |
|
if (low_stack(lim,stack_lim(av,2))) |
|
{ |
|
if (DEBUGMEM>1) err(warnmem,"FpX_resultant (da = %ld)",da); |
|
gerepileall(av,3, &a,&b,&res); |
|
} |
|
da = db; /* = degpol(a) */ |
|
db = dc; /* = degpol(b) */ |
|
} |
|
res = muliimod(res, powgumod((GEN)b[2], da, p), p); |
|
return gerepileuptoint(av, res); |
|
} |
|
|
/* If resultant is 0, *ptU and *ptU are not set */ |
/* If resultant is 0, *ptU and *ptU are not set */ |
ulong |
ulong |
u_FpX_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GEN *ptV) |
u_FpX_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GEN *ptV) |
{ |
{ |
GEN z,q,u,v, x = a, y = b; |
GEN z,q,u,v, x = a, y = b; |
ulong lb, av = avma, res = 1UL; |
ulong lb, res = 1UL; |
|
gpmem_t av = avma; |
long dx,dy,dz; |
long dx,dy,dz; |
|
|
if (!signe(x) || !signe(y)) return 0; |
if (!signe(x) || !signe(y)) return 0; |
Line 2978 u_FpX_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GE |
|
Line 3418 u_FpX_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GE |
|
{ |
{ |
swap(x,y); lswap(dx,dy); pswap(ptU, ptV); |
swap(x,y); lswap(dx,dy); pswap(ptU, ptV); |
a = x; b = y; |
a = x; b = y; |
if (dx & dy & 1) res = p-res; |
if (both_odd(dx,dy)) res = p-res; |
} |
} |
u = u_zeropol(0); |
u = u_zeropol(); |
v = u_scalarpol(1, 0); /* v = 1 */ |
v = u_scalarpol(1); /* v = 1 */ |
while (dy) |
while (dy) |
{ /* b u = x (a), b v = y (a) */ |
{ /* b u = x (a), b v = y (a) */ |
lb = y[dy+2]; |
lb = y[dy+2]; |
q = u_FpX_divrem(x,y, p, 0, &z); |
q = u_FpX_divrem(x,y, p, &z); |
x = y; y = z; /* (x,y) = (y, x - q y) */ |
x = y; y = z; /* (x,y) = (y, x - q y) */ |
dz = degpol(z); if (dz < 0) { avma = av; return 0; } |
dz = degpol(z); if (dz < 0) { avma = av; return 0; } |
z = u_FpX_sub(u, u_FpX_mul(q,v, p), p); |
z = u_FpX_sub(u, u_FpX_mul(q,v, p), p); |
u = v; v = z; /* (u,v) = (v, u - q v) */ |
u = v; v = z; /* (u,v) = (v, u - q v) */ |
|
|
if (dx & dy & 1) res = p - res; |
if (both_odd(dx,dy)) res = p - res; |
if (lb != 1) res = mulssmod(res, powuumod(lb, dx-dz, p), p); |
if (lb != 1) res = mulssmod(res, powuumod(lb, dx-dz, p), p); |
dx = dy; /* = degpol(x) */ |
dx = dy; /* = degpol(x) */ |
dy = dz; /* = degpol(y) */ |
dy = dz; /* = degpol(y) */ |
} |
} |
res = mulssmod(res, powuumod(y[2], dx, p), p); |
res = mulssmod(res, powuumod(y[2], dx, p), p); |
lb = mulssmod(res, u_invmod(y[2],p), p); |
lb = mulssmod(res, u_invmod(y[2],p), p); |
v = gerepileupto(av, u_FpX_Fp_mul(v, lb, p, 0)); |
v = gerepileupto(av, u_FpX_Fp_mul(v, lb, p)); |
av = avma; |
av = avma; |
u = u_FpX_sub(u_scalarpol(res,0), u_FpX_mul(b,v,p), p); |
u = u_FpX_sub(u_scalarpol(res), u_FpX_mul(b,v,p), p); |
u = gerepileupto(av, u_FpX_div(u,a,p)); /* = (res - b v) / a */ |
u = gerepileupto(av, u_FpX_div(u,a,p)); /* = (res - b v) / a */ |
*ptU = u; |
*ptU = u; |
*ptV = v; return res; |
*ptV = v; return res; |
|
|
u_FpX_resultant_all(GEN a, GEN b, long *C0, long *C1, GEN dglist, ulong p) |
u_FpX_resultant_all(GEN a, GEN b, long *C0, long *C1, GEN dglist, ulong p) |
{ |
{ |
long da,db,dc,cnt,ind; |
long da,db,dc,cnt,ind; |
ulong lb,av,cx = 1, res = 1UL; |
ulong lb, cx = 1, res = 1UL; |
|
gpmem_t av; |
GEN c; |
GEN c; |
|
|
if (C0) { *C0 = 1; *C1 = 0; } |
if (C0) { *C0 = 1; *C1 = 0; } |
Line 3022 u_FpX_resultant_all(GEN a, GEN b, long *C0, long *C1, |
|
Line 3463 u_FpX_resultant_all(GEN a, GEN b, long *C0, long *C1, |
|
if (db > da) |
if (db > da) |
{ |
{ |
swapspec(a,b, da,db); |
swapspec(a,b, da,db); |
if (da & db & 1) res = p-res; |
if (both_odd(da,db)) res = p-res; |
} |
} |
/* = res * a[2] ^ db, since 0 <= db <= da = 0 */ |
/* = res * a[2] ^ db, since 0 <= db <= da = 0 */ |
if (!da) return 1; |
if (!da) return 1; |
Line 3039 u_FpX_resultant_all(GEN a, GEN b, long *C0, long *C1, |
|
Line 3480 u_FpX_resultant_all(GEN a, GEN b, long *C0, long *C1, |
|
{ /* check that Euclidean remainder sequence doesn't degenerate */ |
{ /* check that Euclidean remainder sequence doesn't degenerate */ |
if (dc != dglist[ind]) { avma = av; return 0; } |
if (dc != dglist[ind]) { avma = av; return 0; } |
/* update resultant */ |
/* update resultant */ |
if (da & db & 1) res = p-res; |
if (both_odd(da,db)) res = p-res; |
if (lb != 1) |
if (lb != 1) |
{ |
{ |
ulong t = powuumod(lb, da - dc, p); |
ulong t = powuumod(lb, da - dc, p); |
Line 3089 u_FpV_roots_to_pol(GEN a, ulong p) |
|
Line 3530 u_FpV_roots_to_pol(GEN a, ulong p) |
|
{ |
{ |
p2 = cgetg(5,t_VECSMALL); p1[k++] = (long)p2; |
p2 = cgetg(5,t_VECSMALL); p1[k++] = (long)p2; |
p2[2] = mulssmod(a[i], a[i+1], p); |
p2[2] = mulssmod(a[i], a[i+1], p); |
p2[3] = (p<<1) - (a[i] + a[i+1]); |
p2[3] = a[i] + a[i+1]; |
|
if (p2[3] >= p) p2[3] -= p; |
|
if (p2[3]) p2[3] = p - p2[3]; /* - (a[i] + a[i+1]) mod p */ |
p2[4] = 1; p2[1] = evallgef(5); |
p2[4] = 1; p2[1] = evallgef(5); |
} |
} |
if (i < lx) |
if (i < lx) |
|
|
u_pol_comp(GEN P, ulong u, ulong v, ulong p) |
u_pol_comp(GEN P, ulong u, ulong v, ulong p) |
{ |
{ |
long i, l = lgef(P); |
long i, l = lgef(P); |
GEN y = u_allocpol(l-3, 0); |
GEN y = u_getpol(l-3); |
for (i=2; i<l; i++) |
for (i=2; i<l; i++) |
{ |
{ |
ulong t = P[i]; |
ulong t = P[i]; |
Line 3134 u_pol_comp(GEN P, ulong u, ulong v, ulong p) |
|
Line 3577 u_pol_comp(GEN P, ulong u, ulong v, ulong p) |
|
return u_normalizepol(y,l); |
return u_normalizepol(y,l); |
} |
} |
|
|
GEN roots_to_pol(GEN a, long v); |
extern GEN roots_to_pol(GEN a, long v); |
|
|
GEN |
GEN |
polint_triv(GEN xa, GEN ya) |
polint_triv(GEN xa, GEN ya) |
{ |
{ |
GEN P = NULL, Q = roots_to_pol(xa,0); |
GEN P = NULL, Q = roots_to_pol(xa,0); |
long i, n = lg(xa), av = avma, lim = stack_lim(av,2); |
long i, n = lg(xa); |
|
gpmem_t av = avma, lim = stack_lim(av, 2); |
for (i=1; i<n; i++) |
for (i=1; i<n; i++) |
{ |
{ |
GEN T,dP; |
GEN T,dP; |
|
|
u_FpX_div_by_X_x(GEN a, ulong x, ulong p) |
u_FpX_div_by_X_x(GEN a, ulong x, ulong p) |
{ |
{ |
long l = lgef(a), i; |
long l = lgef(a), i; |
GEN z = u_allocpol(l-4, 0), a0, z0; |
GEN z = u_getpol(l-4), a0, z0; |
a0 = a + l-1; |
a0 = a + l-1; |
z0 = z + l-2; *z0 = *a0--; |
z0 = z + l-2; *z0 = *a0--; |
if (u_OK_ULONG(p)) |
if (u_OK_ULONG(p)) |
Line 3234 u_FpX_div_by_X_x(GEN a, ulong x, ulong p) |
|
Line 3678 u_FpX_div_by_X_x(GEN a, ulong x, ulong p) |
|
return z; |
return z; |
} |
} |
|
|
|
static GEN |
|
FpX_div_by_X_x(GEN a, GEN x, GEN p) |
|
{ |
|
long l = lgef(a), i; |
|
GEN z = cgetg(l-1, t_POL), a0, z0; |
|
z[1] = evalsigne(1)|evalvarn(0)|evallgef(l-1); |
|
a0 = a + l-1; |
|
z0 = z + l-2; *z0 = *a0--; |
|
for (i=l-3; i>1; i--) /* z[i] = (a[i+1] + x*z[i+1]) % p */ |
|
{ |
|
GEN t = addii((GEN)*a0--, muliimod(x, (GEN)*z0--, p)); |
|
*z0 = (long)t; |
|
} |
|
return z; |
|
} |
|
|
/* xa, ya = t_VECSMALL */ |
/* xa, ya = t_VECSMALL */ |
GEN |
GEN |
u_FpV_polint(GEN xa, GEN ya, ulong p) |
u_FpV_polint(GEN xa, GEN ya, ulong p) |
{ |
{ |
GEN T,dP, P = NULL, Q = u_FpV_roots_to_pol(xa, p); |
GEN T,dP, P = NULL, Q = u_FpV_roots_to_pol(xa, p); |
long i, n = lg(xa); |
long i, n = lg(xa); |
ulong av, inv; |
ulong inv; |
av = avma; (void)new_chunk(n+3); /* storage space for P */ |
gpmem_t av = avma; |
for (i=1; i<n; i++) |
for (i=1; i<n; i++) |
{ |
{ |
if (!ya[i]) continue; |
if (!ya[i]) continue; |
Line 3253 u_FpV_polint(GEN xa, GEN ya, ulong p) |
|
Line 3713 u_FpV_polint(GEN xa, GEN ya, ulong p) |
|
i++; /* x_i = -x_{i+1} */ |
i++; /* x_i = -x_{i+1} */ |
} |
} |
else |
else |
dP = u_FpX_Fp_mul(T, mulssmod(ya[i],inv,p), p, 0); |
dP = u_FpX_Fp_mul(T, mulssmod(ya[i],inv,p), p); |
if (P) avma = av; |
|
P = P? u_FpX_add(P, dP, p): dP; |
P = P? u_FpX_add(P, dP, p): dP; |
} |
} |
return P? P: u_zeropol(0); |
return P? gerepileupto(av, P): u_zeropol(); |
} |
} |
|
|
|
GEN |
|
FpV_polint(GEN xa, GEN ya, GEN p) |
|
{ |
|
GEN inv,T,dP, P = NULL, Q = FpV_roots_to_pol(xa, p, 0); |
|
long i, n = lg(xa); |
|
gpmem_t av, lim; |
|
av = avma; lim = stack_lim(av,2); |
|
for (i=1; i<n; i++) |
|
{ |
|
if (!signe(ya[i])) continue; |
|
T = FpX_div_by_X_x(Q, (GEN)xa[i], p); |
|
inv = mpinvmod(FpX_eval(T,(GEN)xa[i], p), p); |
|
if (i < n-1 && egalii(addii((GEN)xa[i], (GEN)xa[i+1]), p)) |
|
{ |
|
dP = pol_comp(T, muliimod((GEN)ya[i], inv,p), |
|
muliimod((GEN)ya[i+1],inv,p)); |
|
i++; /* x_i = -x_{i+1} */ |
|
} |
|
else |
|
dP = FpX_Fp_mul(T, muliimod((GEN)ya[i],inv,p), p); |
|
P = P? FpX_add(P, dP, p): dP; |
|
if (low_stack(lim, stack_lim(av,2))) |
|
{ |
|
if (DEBUGMEM>1) err(warnmem,"FpV_polint"); |
|
if (!P) avma = av; else P = gerepileupto(av, P); |
|
} |
|
} |
|
return P? P: zeropol(0); |
|
} |
|
|
static void |
static void |
u_FpV_polint_all(GEN xa, GEN ya, GEN C0, GEN C1, ulong p) |
u_FpV_polint_all(GEN xa, GEN ya, GEN C0, GEN C1, ulong p) |
{ |
{ |
Line 3274 u_FpV_polint_all(GEN xa, GEN ya, GEN C0, GEN C1, ulong |
|
Line 3763 u_FpV_polint_all(GEN xa, GEN ya, GEN C0, GEN C1, ulong |
|
T = u_FpX_div_by_X_x(Q, xa[i], p); |
T = u_FpX_div_by_X_x(Q, xa[i], p); |
inv = u_invmod(u_FpX_eval(T,xa[i], p), p); |
inv = u_invmod(u_FpX_eval(T,xa[i], p), p); |
|
|
#if 0 |
if (ya[i]) |
if (i < n-1 && (ulong)(xa[i] + xa[i+1]) == p) |
|
{ |
{ |
if (ya[i]) |
dP = u_FpX_Fp_mul(T, mulssmod(ya[i],inv,p), p); |
dP = u_pol_comp(T, mulssmod(ya[i],inv,p), mulssmod(ya[i+1],inv,p), p); |
P = P ? u_FpX_add(P , dP , p): dP; |
if (C0[i]) |
|
dP0= u_pol_comp(T, mulssmod(C0[i],inv,p), mulssmod(C0[i+1],inv,p), p); |
|
if (C1[i]) |
|
dP1= u_pol_comp(T, mulssmod(C1[i],inv,p), mulssmod(C1[i+1],inv,p), p); |
|
i++; /* x_i = -x_{i+1} */ |
|
} |
} |
else |
if (C0[i]) |
#endif |
|
{ |
{ |
if (ya[i]) |
dP0= u_FpX_Fp_mul(T, mulssmod(C0[i],inv,p), p); |
{ |
P0= P0? u_FpX_add(P0, dP0, p): dP0; |
dP = u_FpX_Fp_mul(T, mulssmod(ya[i],inv,p), p, 0); |
|
P = P ? u_FpX_add(P , dP , p): dP; |
|
} |
|
if (C0[i]) |
|
{ |
|
dP0= u_FpX_Fp_mul(T, mulssmod(C0[i],inv,p), p, 0); |
|
P0= P0? u_FpX_add(P0, dP0, p): dP0; |
|
} |
|
if (C1[i]) |
|
{ |
|
dP1= u_FpX_Fp_mul(T, mulssmod(C1[i],inv,p), p, 0); |
|
P1= P1? u_FpX_add(P1, dP1, p): dP1; |
|
} |
|
} |
} |
|
if (C1[i]) |
|
{ |
|
dP1= u_FpX_Fp_mul(T, mulssmod(C1[i],inv,p), p); |
|
P1= P1? u_FpX_add(P1, dP1, p): dP1; |
|
} |
} |
} |
ya[1] = (long) (P ? P : u_zeropol(0)); |
ya[1] = (long) (P ? P : u_zeropol()); |
C0[1] = (long) (P0? P0: u_zeropol(0)); |
C0[1] = (long) (P0? P0: u_zeropol()); |
C1[1] = (long) (P1? P1: u_zeropol(0)); |
C1[1] = (long) (P1? P1: u_zeropol()); |
} |
} |
|
|
/* b a vector of polynomials representing B in Fp[X][Y], evaluate at X = x, |
/* b a vector of polynomials representing B in Fp[X][Y], evaluate at X = x, |
* Return 0 in case of degree drop. */ |
* Return 0 in case of degree drop. */ |
static GEN |
static GEN |
vec_u_FpX_eval(GEN b, ulong x, ulong p) |
u_vec_FpX_eval(GEN b, ulong x, ulong p) |
{ |
{ |
GEN z; |
GEN z; |
long i, lb = lgef(b); |
long i, lb = lgef(b); |
ulong leadz = u_FpX_eval((GEN)b[lb-1], x, p); |
ulong leadz = u_FpX_eval((GEN)b[lb-1], x, p); |
if (!leadz) return u_zeropol(0); |
if (!leadz) return u_zeropol(); |
|
|
z = u_allocpol(lb-3, 0); |
z = u_getpol(lb-3); |
for (i=2; i<lb-1; i++) |
for (i=2; i<lb-1; i++) |
z[i] = u_FpX_eval((GEN)b[i], x, p); |
z[i] = u_FpX_eval((GEN)b[i], x, p); |
z[i] = leadz; return z; |
z[i] = leadz; return z; |
Line 3328 vec_u_FpX_eval(GEN b, ulong x, ulong p) |
|
Line 3802 vec_u_FpX_eval(GEN b, ulong x, ulong p) |
|
|
|
/* as above, but don't care about degree drop */ |
/* as above, but don't care about degree drop */ |
static GEN |
static GEN |
vec_u_FpX_eval_gen(GEN b, ulong x, ulong p, int *drop) |
u_vec_FpX_eval_gen(GEN b, ulong x, ulong p, int *drop) |
{ |
{ |
GEN z; |
GEN z; |
long i, lb = lgef(b); |
long i, lb = lgef(b); |
z = u_allocpol(lb-3, 0); |
z = u_getpol(lb-3); |
for (i=2; i<lb; i++) |
for (i=2; i<lb; i++) |
z[i] = u_FpX_eval((GEN)b[i], x, p); |
z[i] = u_FpX_eval((GEN)b[i], x, p); |
z = u_normalizepol(z, lb); |
z = u_normalizepol(z, lb); |
Line 3340 vec_u_FpX_eval_gen(GEN b, ulong x, ulong p, int *drop) |
|
Line 3814 vec_u_FpX_eval_gen(GEN b, ulong x, ulong p, int *drop) |
|
return z; |
return z; |
} |
} |
|
|
|
static GEN |
|
vec_FpX_eval_gen(GEN b, GEN x, GEN p, int *drop) |
|
{ |
|
GEN z; |
|
long i, lb = lgef(b); |
|
z = cgetg(lb, t_POL); |
|
z[1] = b[1]; |
|
for (i=2; i<lb; i++) |
|
z[i] = (long)FpX_eval((GEN)b[i], x, p); |
|
z = normalizepol_i(z, lb); |
|
*drop = lb - lgef(z); |
|
return z; |
|
} |
|
|
/* Interpolate at roots of 1 and use Hadamard bound for univariate resultant: |
/* Interpolate at roots of 1 and use Hadamard bound for univariate resultant: |
* bound = N_2(A)^degpol B N_2(B)^degpol(A), where N_2(A) = sqrt(sum (N_1(Ai))^2) |
* bound = N_2(A)^degpol B N_2(B)^degpol(A), where |
* Return B such that bound < 2^B */ |
* N_2(A) = sqrt(sum (N_1(Ai))^2) |
|
* Return e such that Res(A, B) < 2^e */ |
ulong |
ulong |
ZY_ZXY_ResBound(GEN A, GEN B) |
ZY_ZXY_ResBound(GEN A, GEN B) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
GEN a = gzero, b = gzero, run = realun(DEFAULTPREC); |
GEN a = gzero, b = gzero, run = realun(DEFAULTPREC); |
long i , lA = lgef(A), lB = lgef(B); |
long i , lA = lgef(A), lB = lgef(B); |
for (i=2; i<lA; i++) a = addii(a, sqri((GEN)A[i])); |
for (i=2; i<lA; i++) a = addii(a, sqri((GEN)A[i])); |
Line 3356 ZY_ZXY_ResBound(GEN A, GEN B) |
|
Line 3845 ZY_ZXY_ResBound(GEN A, GEN B) |
|
if (typ(t) == t_POL) t = gnorml1(t, 0); |
if (typ(t) == t_POL) t = gnorml1(t, 0); |
b = addii(b, sqri(t)); |
b = addii(b, sqri(t)); |
} |
} |
b = gmul(gpowgs(mulir(a,run), degpol(B)), gpowgs(mulir(b,run), degpol(A))); |
a = mulir(a,run); |
|
b = mulir(b,run); |
|
b = gmul(gpowgs(a, degpol(B)), gpowgs(b, degpol(A))); |
avma = av; return 1 + (gexpo(b)>>1); |
avma = av; return 1 + (gexpo(b)>>1); |
} |
} |
|
|
|
/* return Res(a(Y), b(n,Y)) over Fp. la = leading_term(a) [for efficiency] */ |
|
static ulong |
|
u_FpX_resultant_after_eval(GEN a, GEN b, ulong n, ulong p, ulong la) |
|
{ |
|
int drop; |
|
GEN ev = u_vec_FpX_eval_gen(b, n, p, &drop); |
|
ulong r = u_FpX_resultant(a, ev, p); |
|
if (drop && la != 1) r = mulssmod(r, powuumod(la, drop,p),p); |
|
return r; |
|
} |
|
static GEN |
|
FpX_resultant_after_eval(GEN a, GEN b, GEN n, GEN p, GEN la) |
|
{ |
|
int drop; |
|
GEN ev = vec_FpX_eval_gen(b, n, p, &drop); |
|
GEN r = FpX_resultant(a, ev, p); |
|
if (drop && !gcmp1(la)) r = muliimod(r, powgumod(la, drop,p),p); |
|
return r; |
|
} |
|
|
|
/* assume dres := deg(Res_Y(a,b), X) <= deg(a,Y) * deg(b,X) < p */ |
|
static GEN |
|
u_FpY_FpXY_resultant(GEN a, GEN b, ulong p, long dres) |
|
{ |
|
ulong la; |
|
long i,n,nmax; |
|
GEN x,y; |
|
|
|
nmax = (dres+1)>>1; |
|
la = (ulong)leading_term(a); |
|
x = cgetg(dres+2, t_VECSMALL); |
|
y = cgetg(dres+2, t_VECSMALL); |
|
/* Evaluate at 0 (if dres even) and +/- n, so that P_n(X) = P_{-n}(-X), |
|
* where P_i is Lagrange polynomial: P_i(j) = 1 if i=j, 0 otherwise */ |
|
for (i=0,n = 1; n <= nmax; n++) |
|
{ |
|
i++; x[i] = n; y[i] = u_FpX_resultant_after_eval(a,b, x[i], p,la); |
|
i++; x[i] = p-n; y[i] = u_FpX_resultant_after_eval(a,b, x[i], p,la); |
|
} |
|
if (i == dres) |
|
{ |
|
i++; x[i] = 0; y[i] = u_FpX_resultant_after_eval(a,b, x[i], p,la); |
|
} |
|
return u_FpV_polint(x,y, p); |
|
} |
|
|
|
/* x^n mod p */ |
|
ulong |
|
u_powmod(ulong x, long n, ulong p) |
|
{ |
|
ulong y = 1, z; |
|
long m; |
|
|
|
if (n < 0) { n = -n; x = u_invmod(x, p); } |
|
m = n; |
|
z = x; |
|
for (;;) |
|
{ |
|
if (m&1) y = mulssmod(y,z, p); |
|
m >>= 1; if (!m) return y; |
|
z = mulssmod(z,z, p); |
|
} |
|
} |
|
|
|
/* x^n mod p, assume n > 0 */ |
|
static GEN |
|
u_FpX_pow(GEN x, long n, ulong p) |
|
{ |
|
GEN y = u_scalarpol(1), z; |
|
long m; |
|
m = n; |
|
z = x; |
|
for (;;) |
|
{ |
|
if (m&1) y = u_FpX_mul(y,z, p); |
|
m >>= 1; if (!m) return y; |
|
z = u_FpX_mul(z,z, p); |
|
} |
|
} |
|
|
|
static GEN |
|
u_FpXX_pseudorem(GEN x, GEN y, ulong p) |
|
{ |
|
long vx = varn(x), dx, dy, dz, i, lx, dp; |
|
gpmem_t av = avma, av2, lim; |
|
|
|
if (!signe(y)) err(talker,"euclidean division by zero (pseudorem)"); |
|
(void)new_chunk(2); |
|
dx=degpol(x); x = revpol(x); |
|
dy=degpol(y); y = revpol(y); dz=dx-dy; dp = dz+1; |
|
av2 = avma; lim = stack_lim(av2,1); |
|
for (;;) |
|
{ |
|
x[0] = (long)u_FpX_neg((GEN)x[0], p); dp--; |
|
for (i=1; i<=dy; i++) |
|
x[i] = (long)u_FpX_add( u_FpX_mul((GEN)y[0], (GEN)x[i], p), |
|
u_FpX_mul((GEN)x[0], (GEN)y[i], p), p ); |
|
for ( ; i<=dx; i++) |
|
x[i] = (long)u_FpX_mul((GEN)y[0], (GEN)x[i], p); |
|
do { x++; dx--; } while (dx >= 0 && !signe((GEN)x[0])); |
|
if (dx < dy) break; |
|
if (low_stack(lim,stack_lim(av2,1))) |
|
{ |
|
if(DEBUGMEM>1) err(warnmem,"pseudorem dx = %ld >= %ld",dx,dy); |
|
gerepilemanycoeffs(av2,x,dx+1); |
|
} |
|
} |
|
if (dx < 0) return u_zeropol(); |
|
lx = dx+3; x -= 2; |
|
x[0]=evaltyp(t_POL) | evallg(lx); |
|
x[1]=evalsigne(1) | evalvarn(vx) | evallgef(lx); |
|
x = revpol(x) - 2; |
|
if (dp) |
|
{ /* multiply by y[0]^dp [beware dummy vars from FpY_FpXY_resultant] */ |
|
GEN t = u_FpX_pow((GEN)y[0], dp, p); |
|
for (i=2; i<lx; i++) |
|
x[i] = (long)u_FpX_mul((GEN)x[i], t, p); |
|
} |
|
return gerepilecopy(av, x); |
|
} |
|
|
|
static GEN |
|
u_FpX_divexact(GEN x, GEN y, ulong p) |
|
{ |
|
long i, l; |
|
GEN z; |
|
if (degpol(y) == 0) |
|
{ |
|
ulong t = (ulong)y[2]; |
|
if (t == 1) return x; |
|
t = u_invmod(t, p); |
|
l = lgef(x); z = cgetg(l, t_POL); z[1] = x[1]; |
|
for (i=2; i<l; i++) z[i] = (long)u_FpX_Fp_mul((GEN)x[i],t,p); |
|
} |
|
else |
|
{ |
|
l = lgef(x); z = cgetg(l, t_POL); z[1] = x[1]; |
|
for (i=2; i<l; i++) z[i] = (long)u_FpX_div((GEN)x[i],y,p); |
|
} |
|
return z; |
|
} |
|
|
|
static GEN |
|
u_FpYX_subres(GEN u, GEN v, ulong p) |
|
{ |
|
gpmem_t av = avma, av2, lim; |
|
long degq,dx,dy,du,dv,dr,signh; |
|
GEN z,g,h,r,p1; |
|
|
|
dx=degpol(u); dy=degpol(v); signh=1; |
|
if (dx < dy) |
|
{ |
|
swap(u,v); lswap(dx,dy); |
|
if (both_odd(dx, dy)) signh = -signh; |
|
} |
|
if (dy < 0) return gzero; |
|
if (dy==0) return gerepileupto(av, u_FpX_pow((GEN)v[2],dx,p)); |
|
|
|
g = h = u_scalarpol(1); av2 = avma; lim = stack_lim(av2,1); |
|
for(;;) |
|
{ |
|
r = u_FpXX_pseudorem(u,v,p); dr = lgef(r); |
|
if (dr == 2) { avma = av; return gzero; } |
|
du = degpol(u); dv = degpol(v); degq = du-dv; |
|
u = v; p1 = g; g = leading_term(u); |
|
switch(degq) |
|
{ |
|
case 0: break; |
|
case 1: |
|
p1 = u_FpX_mul(h,p1, p); h = g; break; |
|
default: |
|
p1 = u_FpX_mul(u_FpX_pow(h,degq,p), p1, p); |
|
h = u_FpX_div(u_FpX_pow(g,degq,p), u_FpX_pow(h,degq-1,p), p); |
|
} |
|
if (both_odd(du,dv)) signh = -signh; |
|
v = u_FpX_divexact(r, p1, p); |
|
if (dr==3) break; |
|
if (low_stack(lim,stack_lim(av2,1))) |
|
{ |
|
if(DEBUGMEM>1) err(warnmem,"subresall, dr = %ld",dr); |
|
gerepileall(av2,4, &u, &v, &g, &h); |
|
} |
|
} |
|
z = (GEN)v[2]; |
|
if (dv > 1) z = u_FpX_div(u_FpX_pow(z,dv,p), u_FpX_pow(h,dv-1,p), p); |
|
if (signh < 0) z = u_FpX_neg(z,p); |
|
return gerepileupto(av, z); |
|
} |
|
|
|
/* return a t_POL (in dummy variable 0) whose coeffs are the coeffs of b, |
|
* in variable v. This is an incorrect PARI object if initially varn(b) < v. |
|
* We could return a vector of coeffs, but it is convenient to have degpol() |
|
* and friends available. Even in that case, it will behave nicely with all |
|
* functions treating a polynomial as a vector of coeffs (eg poleval). |
|
* FOR INTERNAL USE! */ |
|
GEN |
|
swap_vars(GEN b0, long v) |
|
{ |
|
long i, n = poldegree(b0, v); |
|
GEN b = cgetg(n+3, t_POL), x = b + 2; |
|
b[1] = evalsigne(1) | evallgef(n+3) | evalvarn(v); |
|
for (i=0; i<=n; i++) x[i] = (long)polcoeff_i(b0, i, v); |
|
return b; |
|
} |
|
|
|
/* assume varn(b) < varn(a) */ |
|
GEN |
|
FpY_FpXY_resultant(GEN a, GEN b0, GEN p) |
|
{ |
|
long i,n,dres,nmax, vX = varn(b0), vY = varn(a); |
|
GEN la,x,y,b = swap_vars(b0, vY); |
|
|
|
dres = degpol(a)*degpol(b0); |
|
if (OK_ULONG(p)) |
|
{ |
|
ulong pp = (ulong)p[2]; |
|
long l = lgef(b); |
|
|
|
a = u_Fp_FpX(a, pp); |
|
for (i=2; i<l; i++) |
|
b[i] = (long)u_Fp_FpX((GEN)b[i], pp); |
|
if (dres >= pp) |
|
{ |
|
l = lgef(a); |
|
a[0] = evaltyp(t_POL) | evallg(l); |
|
a[1] = evalsigne(1)|evalvarn(vY)|evallgef(l); |
|
for (i=2; i<l; i++) |
|
a[i] = (long)u_scalarpol(a[i]); |
|
x = u_FpYX_subres(a, b, pp); |
|
} |
|
else |
|
x = u_FpY_FpXY_resultant(a, b, pp, dres); |
|
return small_to_pol(x, vX); |
|
} |
|
|
|
nmax = (dres+1)>>1; |
|
la = leading_term(a); |
|
x = cgetg(dres+2, t_VEC); |
|
y = cgetg(dres+2, t_VEC); |
|
/* Evaluate at 0 (if dres even) and +/- n, so that P_n(X) = P_{-n}(-X), |
|
* where P_i is Lagrange polynomial: P_i(j) = 1 if i=j, 0 otherwise */ |
|
for (i=0,n = 1; n <= nmax; n++) |
|
{ |
|
i++; x[i] = lstoi(n); |
|
y[i] = (long)FpX_resultant_after_eval(a,b, (GEN)x[i], p,la); |
|
i++; x[i] = lsubis(p,n); |
|
y[i] = (long)FpX_resultant_after_eval(a,b, (GEN)x[i], p,la); |
|
} |
|
if (i == dres) |
|
{ |
|
i++; x[i] = zero; |
|
y[i] = (long)FpX_resultant_after_eval(a,b, (GEN)x[i], p,la); |
|
} |
|
x = FpV_polint(x,y, p); |
|
setvarn(x, vX); return x; |
|
} |
|
|
|
/* check that theta(maxprime) - theta(27448) >= 2^bound */ |
|
/* NB: theta(27449) ~ 27225.387, theta(x) > 0.98 x for x>7481 |
|
* (Schoenfeld, 1976 for x > 1155901 + direct calculations) */ |
|
static void |
|
check_theta(ulong bound) |
|
{ |
|
ulong c = (ulong)ceil((bound * LOG2 + 27225.388) / 0.98); |
|
if (maxprime() < c) |
|
err(talker,"not enough precalculated primes: need primelimit ~ %lu", c); |
|
} |
|
|
/* 0, 1, -1, 2, -2, ... */ |
/* 0, 1, -1, 2, -2, ... */ |
#define next_lambda(a) (a>0 ? -a : 1-a) |
#define next_lambda(a) (a>0 ? -a : 1-a) |
|
|
Line 3372 ZY_ZXY_ResBound(GEN A, GEN B) |
|
Line 4131 ZY_ZXY_ResBound(GEN A, GEN B) |
|
GEN |
GEN |
ZY_ZXY_resultant_all(GEN A, GEN B0, long *lambda, GEN *LERS) |
ZY_ZXY_resultant_all(GEN A, GEN B0, long *lambda, GEN *LERS) |
{ |
{ |
int checksqfree = lambda? 1: 0, delete = 0, first = 1, stable; |
int checksqfree = lambda? 1: 0, delvar = 0, first = 1, stable; |
ulong av = avma, av2, lim, bound; |
ulong bound; |
|
gpmem_t av = avma, av2, lim; |
long i,n, lb, dres = degpol(A)*degpol(B0), nmax = (dres+1)>>1; |
long i,n, lb, dres = degpol(A)*degpol(B0), nmax = (dres+1)>>1; |
long vX = varn(B0), vY = varn(A); /* assume vX < vY */ |
long vX = varn(B0), vY = varn(A); /* assume vX < vY */ |
GEN x,y,dglist,cB,B,q,a,b,ev,H,H0,H1,Hp,H0p,H1p,C0,C1; |
GEN x,y,dglist,cB,B,q,a,b,ev,H,H0,H1,Hp,H0p,H1p,C0,C1; |
Line 3392 ZY_ZXY_resultant_all(GEN A, GEN B0, long *lambda, GEN |
|
Line 4152 ZY_ZXY_resultant_all(GEN A, GEN B0, long *lambda, GEN |
|
y = cgetg(dres+2, t_VECSMALL); |
y = cgetg(dres+2, t_VECSMALL); |
if (vY == MAXVARN) |
if (vY == MAXVARN) |
{ |
{ |
vY = fetch_var(); delete = 1; |
vY = fetch_var(); delvar = 1; |
B0 = gsubst(B0, MAXVARN, polx[vY]); |
B0 = gsubst(B0, MAXVARN, polx[vY]); |
A = dummycopy(A); setvarn(A, vY); |
A = dummycopy(A); setvarn(A, vY); |
} |
} |
|
|
goto END; |
goto END; |
} |
} |
|
|
|
/* make sure p large enough */ |
|
while (p < (dres<<1)) NEXT_PRIME_VIADIFF(p,d); |
|
|
H = H0 = H1 = NULL; |
H = H0 = H1 = NULL; |
lb = lgef(B); b = u_allocpol(degpol(B), 0); |
lb = lgef(B); b = u_getpol(degpol(B)); |
bound = ZY_ZXY_ResBound(A,B); |
bound = ZY_ZXY_ResBound(A, B); |
if (DEBUGLEVEL>4) fprintferr("bound for resultant coeffs: 2^%ld\n",bound); |
if (DEBUGLEVEL>4) fprintferr("bound for resultant coeffs: 2^%ld\n",bound); |
|
check_theta(bound); |
|
|
for(;;) |
for(;;) |
{ |
{ |
p += *d++; if (!*d) err(primer1); |
NEXT_PRIME_VIADIFF_CHECK(p,d); |
|
|
a = u_Fp_FpX(A, 0, p); |
a = u_Fp_FpX(A, p); |
for (i=2; i<lb; i++) |
for (i=2; i<lb; i++) |
b[i] = (long)u_Fp_FpX((GEN)B[i], 0, p); |
b[i] = (long)u_Fp_FpX((GEN)B[i], p); |
if (LERS) |
if (LERS) |
{ |
{ |
if (!b[lb-1] || degpol(a) < degpol(A)) continue; /* p | lc(A)lc(B) */ |
if (!b[lb-1] || degpol(a) < degpol(A)) continue; /* p | lc(A)lc(B) */ |
|
|
setlg(dglist, 1); |
setlg(dglist, 1); |
for (n=0; n <= dres; n++) |
for (n=0; n <= dres; n++) |
{ |
{ |
ev = vec_u_FpX_eval(b, n, p); |
ev = u_vec_FpX_eval(b, n, p); |
(void)u_FpX_resultant_all(a, ev, NULL, NULL, dglist, p); |
(void)u_FpX_resultant_all(a, ev, NULL, NULL, dglist, p); |
if (lg(dglist)-1 == goal) break; |
if (lg(dglist)-1 == goal) break; |
} |
} |
|
|
|
|
for (i=0,n = 0; i <= dres; n++) |
for (i=0,n = 0; i <= dres; n++) |
{ |
{ |
i++; ev = vec_u_FpX_eval(b, n, p); |
i++; ev = u_vec_FpX_eval(b, n, p); |
x[i] = n; |
x[i] = n; |
y[i] = u_FpX_resultant_all(a, ev, C0+i, C1+i, dglist, p); |
y[i] = u_FpX_resultant_all(a, ev, C0+i, C1+i, dglist, p); |
if (!C1[i]) i--; /* C1(i) = 0. No way to recover C0(i) */ |
if (!C1[i]) i--; /* C1(i) = 0. No way to recover C0(i) */ |
|
|
H1p= (GEN)C1[1]; |
H1p= (GEN)C1[1]; |
} |
} |
else |
else |
{ |
{ /* cf u_FpXY_resultant */ |
ulong la = (ulong)leading_term(a); |
ulong la = (ulong)leading_term(a); |
int drop; |
|
/* Evaluate at 0 (if dres even) and +/- n, so that P_n(X) = P_{-n}(-X), |
/* Evaluate at 0 (if dres even) and +/- n, so that P_n(X) = P_{-n}(-X), |
* where P_i is Lagrange polynomial: P_i(j) = 1 if i=j, 0 otherwise */ |
* where P_i is Lagrange polynomial: P_i(j) = 1 if i=j, 0 otherwise */ |
for (i=0,n = 1; n <= nmax; n++) |
for (i=0,n = 1; n <= nmax; n++) |
{ |
{ |
ev = vec_u_FpX_eval_gen(b, n, p, &drop); |
i++; x[i] = n; y[i] = u_FpX_resultant_after_eval(a,b, x[i], p,la); |
i++; x[i] = n; y[i] = u_FpX_resultant(a, ev, p); |
i++; x[i] = p-n; y[i] = u_FpX_resultant_after_eval(a,b, x[i], p,la); |
if (drop && la != 1) y[i] = mulssmod(y[i], powuumod(la, drop,p),p); |
|
ev = vec_u_FpX_eval_gen(b, p-n, p, &drop); |
|
i++; x[i] = p-n; y[i] = u_FpX_resultant(a, ev, p); |
|
if (drop && la != 1) y[i] = mulssmod(y[i], powuumod(la, drop,p),p); |
|
} |
} |
if (i == dres) |
if (i == dres) |
{ |
{ |
ev = vec_u_FpX_eval_gen(b, 0, p, &drop); |
i++; x[i] = 0; y[i] = u_FpX_resultant_after_eval(a,b, x[i], p,la); |
i++; x[i] = 0; y[i] = u_FpX_resultant(a, ev, p); |
|
if (drop && la != 1) y[i] = mulssmod(y[i], powuumod(la, drop,p),p); |
|
} |
} |
Hp = u_FpV_polint(x,y, p); |
Hp = u_FpV_polint(x,y, p); |
} |
} |
|
|
else |
else |
{ |
{ |
GEN qp = muliu(q,p); |
GEN qp = muliu(q,p); |
stable = ZX_incremental_CRT(H, Hp, q,qp, p); |
stable = ZX_incremental_CRT(&H, Hp, q,qp, p); |
if (LERS) { |
if (LERS) { |
stable &= ZX_incremental_CRT(H0,H0p, q,qp, p); |
stable &= ZX_incremental_CRT(&H0,H0p, q,qp, p); |
stable &= ZX_incremental_CRT(H1,H1p, q,qp, p); |
stable &= ZX_incremental_CRT(&H1,H1p, q,qp, p); |
} |
} |
q = qp; |
q = qp; |
} |
} |
|
|
{ |
{ |
GEN *gptr[4]; gptr[0] = &H; gptr[1] = &q; gptr[2] = &H0; gptr[3] = &H1; |
GEN *gptr[4]; gptr[0] = &H; gptr[1] = &q; gptr[2] = &H0; gptr[3] = &H1; |
if (DEBUGMEM>1) err(warnmem,"ZY_ZXY_resultant"); |
if (DEBUGMEM>1) err(warnmem,"ZY_ZXY_resultant"); |
gerepilemany(av2,gptr,LERS? 4: 2); b = u_allocpol(degpol(B), 0); |
gerepilemany(av2,gptr,LERS? 4: 2); b = u_getpol(degpol(B)); |
} |
} |
} |
} |
END: |
END: |
setvarn(H, vX); if (delete) delete_var(); |
setvarn(H, vX); if (delvar) (void)delete_var(); |
if (cB) H = gmul(H, gpowgs(cB, degpol(A))); |
if (cB) H = gmul(H, gpowgs(cB, degpol(A))); |
if (LERS) |
if (LERS) |
{ |
{ |
|
|
} |
} |
|
|
GEN |
GEN |
ZY_ZXY_resultant(GEN A, GEN B0, long *lambda) |
ZY_ZXY_resultant(GEN A, GEN B, long *lambda) |
{ |
{ |
return ZY_ZXY_resultant_all(A, B0, lambda, NULL); |
return ZY_ZXY_resultant_all(A, B, lambda, NULL); |
} |
} |
|
|
/* If lambda = NULL, return caract(Mod(B, A)), A,B in Z[X]. |
/* If lambda = NULL, return caract(Mod(B, A)), A,B in Z[X]. |
Line 3573 ZY_ZXY_resultant(GEN A, GEN B0, long *lambda) |
|
Line 4330 ZY_ZXY_resultant(GEN A, GEN B0, long *lambda) |
|
GEN |
GEN |
ZX_caract_sqf(GEN A, GEN B, long *lambda, long v) |
ZX_caract_sqf(GEN A, GEN B, long *lambda, long v) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
GEN B0, R, a; |
GEN B0, R, a; |
long dB; |
long dB; |
int delete; |
int delvar; |
|
|
if (v < 0) v = 0; |
if (v < 0) v = 0; |
switch (typ(B)) |
switch (typ(B)) |
Line 3587 ZX_caract_sqf(GEN A, GEN B, long *lambda, long v) |
|
Line 4344 ZX_caract_sqf(GEN A, GEN B, long *lambda, long v) |
|
if (lambda) { B = scalarpol(B,varn(A)); dB = 0; break;} |
if (lambda) { B = scalarpol(B,varn(A)); dB = 0; break;} |
return gerepileupto(av, gpowgs(gsub(polx[v], B), degpol(A))); |
return gerepileupto(av, gpowgs(gsub(polx[v], B), degpol(A))); |
} |
} |
delete = 0; |
delvar = 0; |
if (varn(A) == 0) |
if (varn(A) == 0) |
{ |
{ |
long v0 = fetch_var(); delete = 1; |
long v0 = fetch_var(); delvar = 1; |
A = dummycopy(A); setvarn(A,v0); |
A = dummycopy(A); setvarn(A,v0); |
B = dummycopy(B); setvarn(B,v0); |
B = dummycopy(B); setvarn(B,v0); |
} |
} |
Line 3598 ZX_caract_sqf(GEN A, GEN B, long *lambda, long v) |
|
Line 4355 ZX_caract_sqf(GEN A, GEN B, long *lambda, long v) |
|
B0[2] = (long)gneg_i(B); |
B0[2] = (long)gneg_i(B); |
B0[3] = un; |
B0[3] = un; |
R = ZY_ZXY_resultant(A, B0, lambda); |
R = ZY_ZXY_resultant(A, B0, lambda); |
if (delete) delete_var(); |
if (delvar) (void)delete_var(); |
setvarn(R, v); a = leading_term(A); |
setvarn(R, v); a = leading_term(A); |
if (!gcmp1(a)) R = gdiv(R, gpowgs(a, dB)); |
if (!gcmp1(a)) R = gdiv(R, gpowgs(a, dB)); |
return gerepileupto(av, R); |
return gerepileupto(av, R); |
} |
} |
|
|
|
|
GEN |
GEN |
ZX_caract(GEN A, GEN B, long v) |
ZX_caract(GEN A, GEN B, long v) |
{ |
{ |
return ZX_caract_sqf(A, B, NULL, v); |
return (degpol(A) < 16) ? caractducos(A,B,v): ZX_caract_sqf(A,B, NULL, v); |
} |
} |
|
|
|
/* assume A integral, B in Q[v] */ |
|
GEN |
|
QX_caract(GEN A, GEN B, long v) |
|
{ |
|
gpmem_t av = avma; |
|
GEN cB, B0 = Q_primitive_part(B, &cB); |
|
GEN ch = ZX_caract(A, B0, v); |
|
if (cB) |
|
ch = gerepilecopy(av, rescale_pol(ch, cB)); |
|
return ch; |
|
} |
|
|
static GEN |
static GEN |
trivial_case(GEN A, GEN B) |
trivial_case(GEN A, GEN B) |
{ |
{ |
|
long d; |
if (typ(A) == t_INT) return gpowgs(A, degpol(B)); |
if (typ(A) == t_INT) return gpowgs(A, degpol(B)); |
if (degpol(A) == 0) return trivial_case((GEN)A[2],B); |
d = degpol(A); |
|
if (d == 0) return trivial_case((GEN)A[2],B); |
|
if (d < 0) return gzero; |
return NULL; |
return NULL; |
} |
} |
|
|
|
/* Res(A, B/dB), assuming the A,B in Z[X] and result is integer */ |
GEN |
GEN |
ZX_resultant_all(GEN A, GEN B, ulong bound) |
ZX_resultant_all(GEN A, GEN B, GEN dB, ulong bound) |
{ |
{ |
ulong av = avma, av2, lim, Hp; |
ulong Hp, dp; |
|
gpmem_t av = avma, av2, lim; |
|
long degA; |
int stable; |
int stable; |
GEN q, a, b, H; |
GEN q, a, b, H; |
byteptr d = diffptr + 3000; |
byteptr d = diffptr + 3000; |
Line 3630 ZX_resultant_all(GEN A, GEN B, ulong bound) |
|
Line 4406 ZX_resultant_all(GEN A, GEN B, ulong bound) |
|
if ((H = trivial_case(A,B)) || (H = trivial_case(B,A))) return H; |
if ((H = trivial_case(A,B)) || (H = trivial_case(B,A))) return H; |
q = H = NULL; |
q = H = NULL; |
av2 = avma; lim = stack_lim(av,2); |
av2 = avma; lim = stack_lim(av,2); |
if (!bound) bound = ZY_ZXY_ResBound(A,B); |
degA = degpol(A); |
|
if (!bound) |
|
{ |
|
bound = ZY_ZXY_ResBound(A, B); |
|
if (bound > 50000) |
|
{ |
|
GEN run = realun(MEDDEFAULTPREC); |
|
GEN Ar = gmul(A, run), Br = gmul(B, run); |
|
GEN R = subres(Ar,Br); |
|
bound = gexpo(R) + 1; |
|
} |
|
if (dB) bound -= (long)(mylog2(dB)*degA); |
|
} |
if (DEBUGLEVEL>4) fprintferr("bound for resultant: 2^%ld\n",bound); |
if (DEBUGLEVEL>4) fprintferr("bound for resultant: 2^%ld\n",bound); |
|
check_theta(bound); |
|
|
|
dp = 0; /* gcc -Wall */ |
for(;;) |
for(;;) |
{ |
{ |
p += *d++; if (!*d) err(primer1); |
NEXT_PRIME_VIADIFF_CHECK(p,d); |
a = u_Fp_FpX(A, 0, p); |
if (dB) { dp = smodis(dB, p); if (!dp) continue; } |
b = u_Fp_FpX(B, 0, p); |
|
|
a = u_Fp_FpX(A, p); |
|
b = u_Fp_FpX(B, p); |
Hp= u_FpX_resultant(a, b, p); |
Hp= u_FpX_resultant(a, b, p); |
|
if (dp) Hp = mulssmod(Hp, u_powmod(dp, -degA, p), p); |
if (!H) |
if (!H) |
{ |
{ |
stable = 0; q = utoi(p); |
stable = 0; q = utoi(p); |
H = u_center(Hp, p, p>>1); |
H = stoi(u_center(Hp, p, p>>1)); |
} |
} |
else /* could make it probabilistic ??? [e.g if stable twice, etc] */ |
else /* could make it probabilistic ??? [e.g if stable twice, etc] */ |
{ |
{ |
Line 3664 ZX_resultant_all(GEN A, GEN B, ulong bound) |
|
Line 4457 ZX_resultant_all(GEN A, GEN B, ulong bound) |
|
} |
} |
|
|
GEN |
GEN |
ZX_resultant(GEN A, GEN B) { return ZX_resultant_all(A,B,0); } |
ZX_resultant(GEN A, GEN B) { return ZX_resultant_all(A,B,NULL,0); } |
|
|
|
GEN |
|
ZX_QX_resultant(GEN A, GEN B) |
|
{ |
|
GEN c, d, n, R; |
|
gpmem_t av = avma; |
|
B = Q_primitive_part(B, &c); |
|
if (!c) return ZX_resultant(A,B); |
|
n = numer(c); |
|
d = denom(c); if (is_pm1(d)) d = NULL; |
|
R = ZX_resultant_all(A, B, d, 0); |
|
if (!is_pm1(n)) R = mulii(R, gpowgs(n, degpol(A))); |
|
return gerepileuptoint(av, R); |
|
} |
|
|
/* assume x has integral coefficients */ |
/* assume x has integral coefficients */ |
GEN |
GEN |
ZX_disc_all(GEN x, ulong bound) |
ZX_disc_all(GEN x, ulong bound) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
GEN l, d = ZX_resultant_all(x, derivpol(x), bound); |
GEN l, d = ZX_resultant_all(x, derivpol(x), NULL, bound); |
l = leading_term(x); if (!gcmp1(l)) d = divii(d,l); |
l = leading_term(x); if (!gcmp1(l)) d = diviiexact(d,l); |
if (degpol(x) & 2) d = negi(d); |
if (degpol(x) & 2) d = negi(d); |
return gerepileuptoint(av,d); |
return gerepileuptoint(av,d); |
} |
} |
Line 3681 GEN ZX_disc(GEN x) { return ZX_disc_all(x,0); } |
|
Line 4488 GEN ZX_disc(GEN x) { return ZX_disc_all(x,0); } |
|
int |
int |
ZX_is_squarefree(GEN x) |
ZX_is_squarefree(GEN x) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
int d = (lgef(modulargcd(x,derivpol(x))) == 3); |
int d = (lgef(modulargcd(x,derivpol(x))) == 3); |
avma = av; return d; |
avma = av; return d; |
} |
} |
|
|
/* h integer, P in Z[X]. Return h^degpol(P) P(x / h) */ |
static GEN |
GEN |
_gcd(GEN a, GEN b) |
ZX_rescale_pol(GEN P, GEN h) |
|
{ |
{ |
long i, l = lgef(P); |
if (!a) a = gun; |
GEN Q = cgetg(l,t_POL), hi = gun; |
if (!b) b = gun; |
Q[l-1] = P[l-1]; |
return ggcd(a,b); |
for (i=l-2; i>=2; i--) |
|
{ |
|
hi = gmul(hi,h); Q[i] = lmul((GEN)P[i], hi); |
|
} |
|
Q[1] = P[1]; return Q; |
|
} |
} |
|
|
/* A0 and B0 in Q[X] */ |
/* A0 and B0 in Q[X] */ |
GEN |
GEN |
modulargcd(GEN A0, GEN B0) |
modulargcd(GEN A0, GEN B0) |
{ |
{ |
GEN a,b,Hp,D,A,B,q,qp,H,g,p1; |
GEN a,b,Hp,D,A,B,q,qp,H,g; |
long m,n; |
long m,n; |
ulong p, av2, av = avma, avlim = stack_lim(av,1); |
ulong p; |
|
gpmem_t av2, av = avma, avlim = stack_lim(av, 1); |
byteptr d = diffptr; |
byteptr d = diffptr; |
|
|
if ((typ(A0) | typ(B0)) !=t_POL) err(notpoler,"modulargcd"); |
if ((typ(A0) | typ(B0)) !=t_POL) err(notpoler,"modulargcd"); |
if (!signe(A0)) return gcopy(B0); |
if (!signe(A0)) return gcopy(B0); |
if (!signe(B0)) return gcopy(A0); |
if (!signe(B0)) return gcopy(A0); |
A = content(A0); |
A = primitive_part(A0, &a); check_pol_int(A, "modulargcd"); |
B = content(B0); D = ggcd(A,B); |
B = primitive_part(B0, &b); check_pol_int(B, "modulargcd"); |
A = gcmp1(A)? A0: gdiv(A0,A); |
D = _gcd(a,b); |
B = gcmp1(B)? B0: gdiv(B0,B); |
if (varn(A) != varn(B)) err(talker,"different variables in modulargcd"); |
|
|
/* A, B in Z[X] */ |
/* A, B in Z[X] */ |
g = mppgcd(leading_term(A), leading_term(B)); /* multiple of lead(gcd) */ |
g = mppgcd(leading_term(A), leading_term(B)); /* multiple of lead(gcd) */ |
if (degpol(A) < degpol(B)) swap(A, B); |
if (degpol(A) < degpol(B)) swap(A, B); |
Line 3723 modulargcd(GEN A0, GEN B0) |
|
Line 4526 modulargcd(GEN A0, GEN B0) |
|
|
|
av2 = avma; H = NULL; |
av2 = avma; H = NULL; |
d += 3000; /* 27449 = prime(3000) */ |
d += 3000; /* 27449 = prime(3000) */ |
for(p = 27449; ; p+= *d++) |
for(p = 27449; ;) |
{ |
{ |
if (!*d) err(primer1); |
if (!*d) err(primer1); |
if (!umodiu(g,p)) continue; |
if (!umodiu(g,p)) goto repeat; |
|
|
a = u_Fp_FpX(A, 0, p); |
a = u_Fp_FpX(A, p); |
b = u_Fp_FpX(B, 0, p); Hp = u_FpX_gcd_i(a,b, p); |
b = u_Fp_FpX(B, p); Hp = u_FpX_gcd_i(a,b, p); |
m = degpol(Hp); |
m = degpol(Hp); |
if (m == 0) { H = polun[varn(A0)]; break; } /* coprime. DONE */ |
if (m == 0) { H = polun[varn(A0)]; break; } /* coprime. DONE */ |
if (m > n) continue; /* p | Res(A/G, B/G). Discard */ |
if (m > n) goto repeat; /* p | Res(A/G, B/G). Discard */ |
|
|
if (is_pm1(g)) /* make sure lead(H) = g mod p */ |
if (is_pm1(g)) /* make sure lead(H) = g mod p */ |
Hp = u_FpX_normalize(Hp, p); |
Hp = u_FpX_normalize(Hp, p); |
else |
else |
{ |
{ |
ulong t = mulssmod(umodiu(g, p), u_invmod(Hp[m+2],p), p); |
ulong t = mulssmod(umodiu(g, p), u_invmod(Hp[m+2],p), p); |
Hp = u_FpX_Fp_mul(Hp, t, p, 0); |
Hp = u_FpX_Fp_mul(Hp, t, p); |
} |
} |
if (m < n) |
if (m < n) |
{ /* First time or degree drop [all previous p were as above; restart]. */ |
{ /* First time or degree drop [all previous p were as above; restart]. */ |
H = ZX_init_CRT(Hp,p,varn(A0)); |
H = ZX_init_CRT(Hp,p,varn(A0)); |
q = utoi(p); n = m; continue; |
q = utoi(p); n = m; goto repeat; |
} |
} |
|
|
qp = muliu(q,p); |
qp = muliu(q,p); |
if (ZX_incremental_CRT(H, Hp, q, qp, p)) |
if (ZX_incremental_CRT(&H, Hp, q, qp, p)) |
{ /* H stable: check divisibility */ |
{ /* H stable: check divisibility */ |
if (!is_pm1(g)) { p1 = content(H); if (!is_pm1(p1)) H = gdiv(H,p1); } |
if (!is_pm1(g)) H = primpart(H); |
if (!signe(gres(A,H)) && !signe(gres(B,H))) break; /* DONE */ |
if (gcmp0(pseudorem(A,H)) && gcmp0(pseudorem(B,H))) break; /* DONE */ |
|
|
if (DEBUGLEVEL) fprintferr("modulargcd: trial division failed"); |
if (DEBUGLEVEL) fprintferr("modulargcd: trial division failed"); |
} |
} |
Line 3762 modulargcd(GEN A0, GEN B0) |
|
Line 4565 modulargcd(GEN A0, GEN B0) |
|
if (DEBUGMEM>1) err(warnmem,"modulargcd"); |
if (DEBUGMEM>1) err(warnmem,"modulargcd"); |
gerepilemany(av2,gptr,2); |
gerepilemany(av2,gptr,2); |
} |
} |
|
repeat: |
|
NEXT_PRIME_VIADIFF_CHECK(p,d); |
} |
} |
return gerepileupto(av, gmul(D,H)); |
return gerepileupto(av, gmul(D,H)); |
} |
} |
|
|
/* lift(1 / Mod(A,B)) */ |
/* lift(1 / Mod(A,B)). B0 a t_POL, A0 a scalar or a t_POL. Rational coeffs */ |
GEN |
GEN |
ZX_invmod(GEN A0, GEN B0) |
QX_invmod(GEN A0, GEN B0) |
{ |
{ |
GEN a,b,D,A,B,q,qp,Up,Vp,U,V,res; |
GEN a,b,D,A,B,q,qp,Up,Vp,U,V,res; |
long stable; |
long stable; |
ulong p, av2, av = avma, avlim = stack_lim(av,1); |
ulong p; |
|
gpmem_t av2, av = avma, avlim = stack_lim(av, 1); |
byteptr d = diffptr; |
byteptr d = diffptr; |
|
|
if (typ(B0) != t_POL) err(notpoler,"ZX_invmod"); |
if (typ(B0) != t_POL) err(notpoler,"QX_invmod"); |
if (typ(A0) != t_POL) |
if (typ(A0) != t_POL) |
{ |
{ |
if (is_scalar_t(typ(A0))) return ginv(A0); |
if (is_scalar_t(typ(A0))) return ginv(A0); |
err(notpoler,"ZX_invmod"); |
err(notpoler,"QX_invmod"); |
} |
} |
A = content(A0); D = A; |
if (degpol(A0) < 15) return ginvmod(A0,B0); |
B = content(B0); |
A = primitive_part(A0, &D); |
A = gcmp1(A)? A0: gdiv(A0,A); |
B = primpart(B0); |
B = gcmp1(B)? B0: gdiv(B0,B); |
|
/* A, B in Z[X] */ |
/* A, B in Z[X] */ |
av2 = avma; U = NULL; |
av2 = avma; U = NULL; |
d += 3000; /* 27449 = prime(3000) */ |
d += 3000; /* 27449 = prime(3000) */ |
for(p = 27449; ; p+= *d++) |
for(p = 27449; ; ) |
{ |
{ |
if (!*d) err(primer1); |
if (!*d) err(primer1); |
a = u_Fp_FpX(A, 0, p); |
a = u_Fp_FpX(A, p); |
b = u_Fp_FpX(B, 0, p); |
b = u_Fp_FpX(B, p); |
/* if p | Res(A/G, B/G), discard */ |
/* if p | Res(A/G, B/G), discard */ |
if (!u_FpX_extresultant(b,a,p, &Vp,&Up)) continue; |
if (!u_FpX_extresultant(b,a,p, &Vp,&Up)) goto repeat; |
|
|
if (!U) |
if (!U) |
{ /* First time */ |
{ /* First time */ |
U = ZX_init_CRT(Up,p,varn(A0)); |
U = ZX_init_CRT(Up,p,varn(A0)); |
V = ZX_init_CRT(Vp,p,varn(A0)); |
V = ZX_init_CRT(Vp,p,varn(A0)); |
q = utoi(p); continue; |
q = utoi(p); goto repeat; |
} |
} |
if (DEBUGLEVEL>5) msgtimer("ZX_invmod: mod %ld (bound 2^%ld)", p,expi(q)); |
if (DEBUGLEVEL>5) msgtimer("QX_invmod: mod %ld (bound 2^%ld)", p,expi(q)); |
qp = muliu(q,p); |
qp = muliu(q,p); |
stable = ZX_incremental_CRT(U, Up, q,qp, p); |
stable = ZX_incremental_CRT(&U, Up, q,qp, p); |
stable &= ZX_incremental_CRT(V, Vp, q,qp, p); |
stable &= ZX_incremental_CRT(&V, Vp, q,qp, p); |
if (stable) |
if (stable) |
{ /* all stable: check divisibility */ |
{ /* all stable: check divisibility */ |
res = gadd(gmul(A,U), gmul(B,V)); |
res = gadd(gmul(A,U), gmul(B,V)); |
if (degpol(res) == 0) break; /* DONE */ |
if (degpol(res) == 0) break; /* DONE */ |
if (DEBUGLEVEL) fprintferr("ZX_invmod: char 0 check failed"); |
if (DEBUGLEVEL) fprintferr("QX_invmod: char 0 check failed"); |
} |
} |
q = qp; |
q = qp; |
if (low_stack(avlim, stack_lim(av,1))) |
if (low_stack(avlim, stack_lim(av,1))) |
{ |
{ |
GEN *gptr[3]; gptr[0]=&q; gptr[1]=&U; gptr[2]=&V; |
GEN *gptr[3]; gptr[0]=&q; gptr[1]=&U; gptr[2]=&V; |
if (DEBUGMEM>1) err(warnmem,"ZX_invmod"); |
if (DEBUGMEM>1) err(warnmem,"QX_invmod"); |
gerepilemany(av2,gptr,3); |
gerepilemany(av2,gptr,3); |
} |
} |
|
repeat: |
|
NEXT_PRIME_VIADIFF_CHECK(p,d); |
} |
} |
D = gmul(D,res); |
D = D? gmul(D,res): res; |
return gerepileupto(av, gdiv(U,D)); |
return gerepileupto(av, gdiv(U,D)); |
} |
} |
|
|
|
/* irreducible (unitary) polynomial of degree n over Fp */ |
|
GEN |
|
ffinit_rand(GEN p,long n) |
|
{ |
|
gpmem_t av = avma; |
|
GEN pol; |
|
|
|
for(;; avma = av) |
|
{ |
|
pol = gadd(gpowgs(polx[0],n), FpX_rand(n-1,0, p)); |
|
if (FpX_is_irred(pol, p)) break; |
|
} |
|
return pol; |
|
} |
|
|
|
GEN |
|
FpX_direct_compositum(GEN A, GEN B, GEN p) |
|
{ |
|
GEN C,a,b,x; |
|
a = dummycopy(A); setvarn(a, MAXVARN); |
|
b = dummycopy(B); setvarn(b, MAXVARN); |
|
x = gadd(polx[0], polx[MAXVARN]); |
|
C = FpY_FpXY_resultant(a, poleval(b,x),p); |
|
return C; |
|
} |
|
|
|
GEN |
|
FpX_compositum(GEN A, GEN B, GEN p) |
|
{ |
|
GEN C, a,b; |
|
long k; |
|
|
|
a = dummycopy(A); setvarn(a, MAXVARN); |
|
b = dummycopy(B); setvarn(b, MAXVARN); |
|
for (k = 1;; k = next_lambda(k)) |
|
{ |
|
GEN x = gadd(polx[0], gmulsg(k, polx[MAXVARN])); |
|
C = FpY_FpXY_resultant(a, poleval(b,x),p); |
|
if (FpX_is_squarefree(C, p)) break; |
|
} |
|
return C; |
|
} |
|
|
|
/* return an extension of degree 2^l of F_2, assume l > 0 |
|
* using Adleman-Lenstra Algorithm. |
|
* Not stack clean. */ |
|
static GEN |
|
f2init(long l) |
|
{ |
|
long i; |
|
GEN a, q, T, S; |
|
|
|
if (l == 1) return cyclo(3, MAXVARN); |
|
|
|
a = gun; |
|
S = coefs_to_pol(4, gun,gun,gzero,gzero); /* #(#^2 + #) */ |
|
setvarn(S, MAXVARN); |
|
q = coefs_to_pol(3, gun,gun, S); /* X^2 + X + #(#^2+#) */ |
|
|
|
/* x^4+x+1, irred over F_2, minimal polynomial of a root of q */ |
|
T = coefs_to_pol(5, gun,gzero,gzero,gun,gun); |
|
|
|
for (i=2; i<l; i++) |
|
{ /* q = X^2 + X + a(#) irred. over K = F2[#] / (T(#)) |
|
* ==> X^2 + X + a(#) b irred. over K for any root b of q |
|
* ==> X^2 + X + (b^2+b)b */ |
|
setvarn(T, MAXVARN); |
|
T = FpY_FpXY_resultant(T, q, gdeux); |
|
/* T = minimal polynomial of b over F2 */ |
|
} |
|
return T; |
|
} |
|
|
|
/*Check if subcyclo(n,l,0) is irreducible modulo p*/ |
|
static long |
|
fpinit_check(GEN p, long n, long l) |
|
{ |
|
gpmem_t ltop=avma; |
|
long q,o; |
|
if (!isprime(stoi(n))) {avma=ltop; return 0;} |
|
q = smodis(p,n); |
|
if (!q) {avma=ltop; return 0;} |
|
o = itos(order(gmodulss(q,n))); |
|
avma = ltop; |
|
return ( cgcd((n-1)/o,l) == 1 ); |
|
} |
|
|
|
/* let k=2 if p%4==1, and k=4 else and assume k*p does not divide l. |
|
* Return an irreducible polynomial of degree l over F_p. |
|
* This a variant of an algorithm of Adleman and Lenstra |
|
* "Finding irreducible polynomials over finite fields", |
|
* ACM, 1986 (5) 350--355 |
|
* Not stack clean. |
|
*/ |
|
static GEN |
|
fpinit(GEN p, long l) |
|
{ |
|
ulong n = 1+l, k = 1; |
|
while (!fpinit_check(p,n,l)) { n += l; k++; } |
|
if (DEBUGLEVEL>=4) |
|
fprintferr("FFInit: using subcyclo(%ld, %ld)\n",n,l); |
|
return FpX_red(subcyclo(n,l,0),p); |
|
} |
|
|
|
GEN |
|
ffinit_fact(GEN p, long n) |
|
{ |
|
gpmem_t ltop=avma; |
|
GEN F; /* vecsmall */ |
|
GEN P; /* pol */ |
|
long i; |
|
F = decomp_primary_small(n); |
|
/* If n is even, then F[1] is 2^bfffo(n)*/ |
|
if (!(n&1) && egalii(p, gdeux)) |
|
P = f2init(vals(n)); |
|
else |
|
P = fpinit(p, F[1]); |
|
for (i = 2; i < lg(F); ++i) |
|
P = FpX_direct_compositum(fpinit(p, F[i]), P, p); |
|
return gerepileupto(ltop,FpX(P,p)); |
|
} |
|
|
|
GEN |
|
ffinit_nofact(GEN p, long n) |
|
{ |
|
gpmem_t av = avma; |
|
GEN P,Q=NULL; |
|
if (lgefint(p)==3) |
|
{ |
|
ulong lp=p[2], q; |
|
long v=svaluation(n,lp,&q); |
|
if (v>0) |
|
{ |
|
if (lp==2) |
|
Q=f2init(v); |
|
else |
|
Q=fpinit(p,n/q); |
|
n=q; |
|
} |
|
} |
|
if (n==1) |
|
P=Q; |
|
else |
|
{ |
|
P = fpinit(p, n); |
|
if (Q) P = FpX_direct_compositum(P, Q, p); |
|
} |
|
return gerepileupto(av, FpX(P,p)); |
|
} |
|
|
|
GEN |
|
ffinit(GEN p, long n, long v) |
|
{ |
|
gpmem_t ltop=avma; |
|
GEN P; |
|
if (n <= 0) err(talker,"non positive degree in ffinit"); |
|
if (typ(p) != t_INT) err(typeer, "ffinit"); |
|
if (v < 0) v = 0; |
|
if (n == 1) return FpX(polx[v],p); |
|
/*If we are in a easy case just use cyclo*/ |
|
if (fpinit_check(p, n + 1, n)) |
|
return gerepileupto(ltop,FpX(cyclo(n + 1, v),p)); |
|
if (lgefint(p)-2<BITS_IN_LONG-bfffo(n)) |
|
P=ffinit_fact(p,n); |
|
else |
|
P=ffinit_nofact(p,n); |
|
setvarn(P, v); |
|
return P; |
|
} |