=================================================================== RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/polarit3.c,v retrieving revision 1.1 retrieving revision 1.2 diff -u -p -r1.1 -r1.2 --- OpenXM_contrib/pari-2.2/src/basemath/Attic/polarit3.c 2001/10/02 11:17:05 1.1 +++ OpenXM_contrib/pari-2.2/src/basemath/Attic/polarit3.c 2002/09/11 07:26:52 1.2 @@ -1,4 +1,4 @@ -/* $Id: polarit3.c,v 1.1 2001/10/02 11:17:05 noro Exp $ +/* $Id: polarit3.c,v 1.2 2002/09/11 07:26:52 noro Exp $ Copyright (C) 2000 The PARI group. @@ -33,6 +33,11 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston, #endif #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) */ @@ -80,7 +85,7 @@ specpol(GEN x, long nx) /* multiplication in Fp[X], p small */ -static GEN +GEN u_normalizepol(GEN x, long lx) { long i; @@ -120,14 +125,14 @@ u_FpX_gcd_i(GEN a, GEN b, ulong p) GEN 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)); } int 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); avma = av; return (degpol(d) == 0); } @@ -225,28 +230,53 @@ u_FpX_addshift(GEN x, GEN y, ulong p, long d) *--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_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) +GEN +u_zeropol() +{ + GEN z = cgetg(2, t_VECSMALL); + z[1] = evallgef(2) | evalvarn(0); return z; +} + 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[1] = evalsigne((d >= 0)) | evallgef(d+3) | evalvarn(0); 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 -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; } 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>=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 will be much slower with big unreduced coefficient @@ -869,7 +945,7 @@ FpX_Fp_mul(GEN y,GEN x,GEN p) GEN FpXQ_inv(GEN x,GEN T,GEN p) { - ulong av; + gpmem_t av; GEN U; if (!T) return mpinvmod(x,p); @@ -880,9 +956,16 @@ FpXQ_inv(GEN x,GEN T,GEN p) } 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; GEN z = FpX_Fp_mul((GEN)V[1],(GEN)W[1],NULL); for(i=2;ipol, 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 */ GEN u_FpXQ_pow(GEN x, GEN n, GEN pol, ulong p) { - ulong av = avma; - GEN p1 = n+2, y = x; - long m,i,j; - m = *p1; - j = 1+bfffo(m); m <<= j; j = BITS_IN_LONG-j; - for (i=lgefint(n)-2;;) - { - 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; - } + gpmem_t av = avma; + u_FpX_muldata D; + GEN y; + D.pol = pol; + D.p = p; + y = leftright_pow(x, n, (void*)&D, &_u_sqr, &_u_mul); return gerepileupto(av, y); } @@ -1093,9 +1199,10 @@ u_FpXQ_pow(GEN x, GEN n, GEN pol, ulong p) GEN FpXQ_pow(GEN x, GEN n, GEN pol, GEN p) { - ulong av, ltop = avma, lim=stack_lim(avma,1); - long m,i,j, vx = varn(x); - GEN p1 = n+2, y; + gpmem_t av = avma; + FpX_muldata D; + long vx = varn(x); + GEN y; if (!signe(n)) return polun[vx]; if (signe(n) < 0) { @@ -1107,41 +1214,28 @@ FpXQ_pow(GEN x, GEN n, GEN pol, GEN p) if (OK_ULONG(p)) { ulong pp = p[2]; - pol = u_Fp_FpX(pol,0, pp); - x = u_Fp_FpX(x,0, pp); + pol = u_Fp_FpX(pol, pp); + x = u_Fp_FpX(x, pp); y = u_FpXQ_pow(x, n, pol, pp); y = small_to_pol(y,vx); } else { av = avma; - m = *p1; y = x; - j = 1+bfffo(m); m <<= j; j = BITS_IN_LONG-j; - for (i=lgefint(n)-2;;) - { - 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; - } + D.pol = pol; + D.p = p; + y = leftright_pow(x, n, (void*)&D, &_sqr, &_mul); } - 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] */ @@ -1152,25 +1246,25 @@ FpXQ_pow(GEN x, GEN n, GEN pol, GEN p) GEN FpXX_red(GEN z, GEN p) { - GEN res; - int i; - res=cgetg(lgef(z),t_POL); - res[1] = evalsigne(1) | evalvarn(varn(z)) | evallgef(lgef(z)); - for(i=2;iS, 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] */ /* */ /*******************************************************************/ @@ -1328,7 +1534,7 @@ FpXQX_safegcd(GEN P, GEN Q, GEN T, GEN p) #define FqX_sqr FpXQX_sqr #define FqX_red FpXQX_red 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) { @@ -1342,7 +1548,7 @@ static GEN fgmul(GEN a,GEN b){return FqX_mul(a,b,Tmodu GEN FqV_roots_to_pol(GEN V, GEN Tp, GEN p, long v) { - ulong ltop=avma; + gpmem_t ltop=avma; long k; GEN W=cgetg(lg(V),t_VEC); for(k=1;k=6) - fprintferr("FF l-Gen:next %Z",z); + z = FpX_add(z, monomial(stoi(u%pp),v,x), NULL); + u /= pp; v++; } + if ( DEBUGLEVEL>=6 ) fprintferr("FF l-Gen:next %Z\n",z); m1 = m = FpXQ_pow(z,r,T,p); + if (gcmp1(m)) { avma = av1; continue; } + for (i=1; i=1 - * pgcd(r,l)=1 - * m=y^(q/l) - * y n'est pas une puissance l-ieme - * m!=1 - * ouf! - */ +/* Solve x^l = a mod (p,T) + * l must be prime + * q = p^degpol(T)-1 = (l^e)*r, with e>=1 and pgcd(r,l)=1 + * m = y^(q/l) + * y not an l-th power [ m != 1 ] */ GEN 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; GEN p1,p2,u1,u2,v,w,z; - bezout(r,l,&u1,&u2); - v=FpXQ_pow(a,u2,T,p); - w=FpXQ_pow(a,modii(mulii(negi(u1),r),q),T,p); + if (gcmp1(a)) return gcopy(a); + + (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); while (!gcmp1(w)) { - /* if p is not prime, next loop will not end */ - k=0; - p1=w; - do - { - z=p1; - p1=FpXQ_pow(p1,l,T,p); + k = 0; + p1 = w; + do { /* if p is not prime, loop will not end */ + z = p1; + p1 = FpXQ_pow(p1,l,T,p); k++; - }while(!gcmp1(p1)); + } while (!gcmp1(p1)); if (k==e) { avma=av; return NULL; } 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 - giant step instead*/ - p1= FpXQ_pow(y,modii(mulsi(i,gpowgs(l,e-k-1)),q),T,p); + for (i=1; !gcmp1(p2); i++) p2 = FpXQ_mul(p2,m,T,p);/*TODO: BS/GS instead*/ + p1= FpXQ_pow(y, modii(mulsi(i,gpowgs(l,e-k-1)), q), T,p); m = FpXQ_pow(m,stoi(i),T,p); e = k; v = FpXQ_mul(p1,v,T,p); @@ -1436,95 +1630,71 @@ ffsqrtlmod(GEN a, GEN l, GEN T ,GEN p , GEN q, long e, w = FpXQ_mul(y,w,T,p); if (low_stack(lim, stack_lim(av,1))) { - GEN *gptr[4]; if(DEBUGMEM>1) err(warnmem,"ffsqrtlmod"); - gptr[0]=&y; gptr[1]=&v; gptr[2]=&w; gptr[3]=&m; - gerepilemany(av,gptr,4); + gerepileall(av,4, &y,&v,&w,&m); } } - 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 - -return a solution of - -x^n=a mod p - -1)If there is no solution return NULL and if zetan is not NULL set zetan to gzero. - -2) If there is solution there are exactly m=gcd(p-1,n) of them. - -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 -*/ +/* Solve x^n = a mod p: n integer, a in Fp[X]/(T) [ p prime, T irred. mod p ] + * + * 1) if no solution, return NULL and (if zetan != NULL) set zetan to gzero. + * + * 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 + * of solutions is {x*zetan^k;k=0 to m-1} + * + * If a = 0, return 0 and (if zetan != NULL) set zetan = gun */ 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; GEN m,u1,u2; 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) err(typeer,"ffsqrtnmod"); - if (lgef(T)==3) - err(constpoler,"ffsqrtnmod"); - if(!signe(n)) - err(talker,"1/0 exponent in ffsqrtnmod"); - 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))); - m=bezout(n,q,&u1,&u2); - if (gcmp(m,n)) - { - 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 (!degpol(T)) err(constpoler,"ffsqrtnmod"); + if (!signe(n)) err(talker,"1/0 exponent in ffsqrtnmod"); + 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))); + m = bezout(n,q,&u1,&u2); + if (!egalii(m,n)) a = FpXQ_pow(a, modii(u1,q), T,p); + if (zetan) z = polun[varn(T)]; + lim = stack_lim(ltop,1); if (!gcmp1(m)) { - m=decomp(m); - av1=avma; + m = decomp(m); av1 = avma; for (i = lg(m[1])-1; i; i--) { - l=gcoeff(m,i,1); j=itos(gcoeff(m,i,2)); - e=pvaluation(q,l,&r); - y=fflgen(l,e,r,T,p,&zeta); - if (zetan) z=FpXQ_mul(z,FpXQ_pow(y,gpowgs(l,e-j),T,p),T,p); - do + l = gcoeff(m,i,1); + j = itos(gcoeff(m,i,2)); + e = pvaluation(q,l,&r); + if(DEBUGLEVEL>=6) (void)timer2(); + 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); - if (!a){avma=ltop;return NULL;} - j--; - }while (j); + a = ffsqrtlmod(a,l,T,p,q,e,r,y,zeta); + if (!a) {avma=ltop; return NULL;} + } 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 (zetan) - { - z=gcopy(z); - gptr[0]=&a;gptr[1]=&z; - gerepilemanysp(av1,lbot,gptr,2); - } - else - a=gerepileupto(av1,a); - lbot=av1; + gerepileall(av1,zetan? 2: 1, &a,&z); } } } if (zetan) { - *zetan=gcopy(z); - gptr[0]=&a;gptr[1]=zetan; - gerepilemanysp(ltop,lbot,gptr,2); + *zetan = z; + gerepileall(ltop,2,&a,zetan); } else - a=gerepileupto(ltop,a); + a = gerepileupto(ltop, a); return a; } /*******************************************************************/ @@ -1541,7 +1711,7 @@ matrixpow(long n, long m, GEN y, GEN P,GEN l) GEN Fp_inv_isom(GEN S,GEN T, GEN p) { - ulong ltop = avma, lbot; + gpmem_t lbot, ltop = avma; GEN M, V; int n, i; long x; @@ -1557,25 +1727,9 @@ Fp_inv_isom(GEN S,GEN T, GEN p) V = gtopolyrev(V, x); return gerepile(ltop, lbot, V); } -#if 0 -/*Old, trivial, implementation*/ -GEN -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 + +/* Let M the matrix of the x^p Frobenius automorphism. + * Compute x^(p^i) for i=0..r */ static GEN polfrobenius(GEN M, GEN p, long r, long v) @@ -1585,61 +1739,96 @@ polfrobenius(GEN M, GEN p, long r, long v) V = cgetg(r+2,t_VEC); V[1] = (long) polx[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]; for (i = 3; i <= r+1; ++i) { W = FpM_FpV_mul(M,W,p); - V[i] = (long) gtopolyrev(W,v); + V[i] = (long) vec_to_pol(W,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 matpolfrobenius(GEN V, GEN P, GEN T, GEN p) { + gpmem_t btop; long i; long l=degpol(T); long v=varn(T); - GEN M,W; + GEN M,W,Mi; GEN PV=gtovec(P); + GEN *gptr[2]; + long lV=lg(V); PV=cgetg(degpol(P)+2,t_VEC); for(i=1;i=4) timer2(); - V=polfrobenius(MA,l,r,varn(U)); + if (DEBUGLEVEL>=4) (void)timer2(); + V=polfrobenius(MA,l,r,vu); if (DEBUGLEVEL>=4) msgtimer("pol[frobenius]"); - M=matpolfrobenius(V,lU,P,l); + M=matpolfrobenius(V,U,P,l); if (DEBUGLEVEL>=4) msgtimer("matrix cyclo"); A=FpM_ker(M,l); if (DEBUGLEVEL>=4) msgtimer("kernel"); @@ -1652,18 +1841,17 @@ 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 * 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[1]=A[1]; R[r]=(long)FpM_FpV_mul(MA,gmul((GEN)A[1],ib0),l); for(i=r-1;i>1;i--) 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); for(i=1;i=4) timer2(); 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 (DEBUGLEVEL>=4) msgtimer("matrixpow"); A=Ap=zeropol(vp); 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 don't. (Well,now that we don't in the other case also, it is more dubious to treat cases apart...)*/ @@ -1710,21 +1896,23 @@ 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); z=negi(lift((GEN)z[1])); ipg=stoi(pg); - if (DEBUGLEVEL>=4) timer2(); + if (DEBUGLEVEL>=4) (void)timer2(); A=FpM_ker(gaddmat(z, MA),l); if (lg(A)!=2) err(talker,"ZZ_%Z[%Z]/(%Z) is not a field in Fp_intersect" ,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); if (lg(B)!=2) err(talker,"ZZ_%Z[%Z]/(%Z) is not a field in Fp_intersect" ,l,polx[vq],Q); - B=gtopolyrev((GEN)B[1],vq); + B=vec_to_pol((GEN)B[1],vq); if (DEBUGLEVEL>=4) msgtimer("FpM_ker"); An=(GEN) FpXQ_pow(A,ipg,P,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); if ( !L ) err(talker,"Polynomials not irreducible in Fp_intersect"); @@ -1733,26 +1921,24 @@ Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN } else { - GEN L,An,Bn,ipg,U,lU,z; - U=gmael(factmod(cyclo(pg,MAXVARN),l),1,1); - lU=lift(U); + GEN L,An,Bn,ipg,U,z; + U=lift(gmael(factmod(cyclo(pg,MAXVARN),l),1,1)); ipg=stoi(pg); - A=intersect_ker(P, MA, l, U, lU); - B=intersect_ker(Q, MB, l, U, lU); - /*Somewhat ugly, but it is a proof that POLYMOD are useful and - powerful.*/ - if (DEBUGLEVEL>=4) timer2(); - An=lift(lift((GEN)lift(gpowgs(gmodulcp(A,P),pg))[2])); - Bn=lift(lift((GEN)lift(gpowgs(gmodulcp(B,Q),pg))[2])); + A=intersect_ker(P, MA, U, l); + B=intersect_ker(Q, MB, U, l); + if (DEBUGLEVEL>=4) (void)timer2(); + An=(GEN) FpXQYQ_pow(A,stoi(pg),U,P,l)[2]; + Bn=(GEN) FpXQYQ_pow(B,stoi(pg),U,Q,l)[2]; if (DEBUGLEVEL>=4) msgtimer("pows [P,Q]"); - z=FpXQ_inv(Bn,lU,l); - z=FpXQ_mul(An,z,lU,l); - L=ffsqrtnmod(z,ipg,lU,l,NULL); + z=FpXQ_inv(Bn,U,l); + z=FpXQ_mul(An,z,U,l); + L=ffsqrtnmod(z,ipg,U,l,NULL); if (DEBUGLEVEL>=4) msgtimer("ffsqrtn"); if ( !L ) err(talker,"Polynomials not irreducible in Fp_intersect"); - B=gsubst(lift(lift(gmul(B,L))),MAXVARN,gzero); - A=gsubst(lift(lift(A)),MAXVARN,gzero); + B=FpXQX_FpXQ_mul(B,L,U,l); + B=gsubst(B,MAXVARN,gzero); + A=gsubst(A,MAXVARN,gzero); } } if (e!=0) @@ -1784,7 +1970,7 @@ Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN for(;i<=np;i++) VP[i]=zero; } Ap=FpM_invimage(MA,VP,l); - Ap=gtopolyrev(Ap,vp); + Ap=vec_to_pol(Ap,vp); if (j) { By=FpXQ_mul(By,FpXQ_pow(Bp,lmun,Q,l),Q,l); @@ -1792,7 +1978,7 @@ Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN for(;i<=nq;i++) VQ[i]=zero; } Bp=FpM_invimage(MB,VQ,l); - Bp=gtopolyrev(Bp,vq); + Bp=vec_to_pol(Bp,vq); if (DEBUGLEVEL>=4) msgtimer("FpM_invimage"); } }/*FpX_add is not clean, so we must do it *before* lbot=avma*/ @@ -1811,7 +1997,7 @@ Fp_intersect(long n, GEN P, GEN Q, GEN l,GEN *SP, GEN * isomorphism. */ GEN Fp_isom(GEN P,GEN Q,GEN l) { - ulong ltop=avma; + gpmem_t ltop=avma; GEN SP,SQ,R; Fp_intersect(degpol(P),P,Q,l,&SP,&SQ,NULL,NULL); R=Fp_inv_isom(SP,P,l); @@ -1819,33 +2005,39 @@ GEN Fp_isom(GEN P,GEN Q,GEN l) return gerepileupto(ltop,R); } 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; - GEN R,V,ld,Tl; + gpmem_t ltop=avma; + GEN R,V,Tl,z,M; long n,k,m; long v; v=varn(P); n=degpol(P); 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); V=cgetg(m+1,t_VEC); V[1]=lpolx[w]; + z=pol_to_vec((GEN)V[1],n); 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); + if (DEBUGLEVEL>=4) msgtimer("Fp_factorgalois: pol"); 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 * a Frobenius cycle] */ GEN Fp_factor_irred(GEN P,GEN l, GEN Q) { - ulong ltop=avma,av; - GEN SP,SQ,MP,MQ,M,MF,E,V,IR,res; + gpmem_t ltop=avma, av; + GEN SP,SQ,MP,MQ,M,FP,FQ,E,V,IR,res; long np=degpol(P),nq=degpol(Q); long i,d=cgcd(np,nq); long vp=varn(P),vq=varn(Q); @@ -1855,10 +2047,13 @@ Fp_factor_irred(GEN P,GEN l, GEN Q) res[1]=lcopy(P); return res; } - MF=matrixpow(nq,nq,FpXQ_pow(polx[vq],l,Q,l),Q,l); - Fp_intersect(d,P,Q,l,&SP,&SQ,NULL,MF); + if (DEBUGLEVEL>=4) (void)timer2(); + 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; - E=Fp_factorgalois(P,l,d,vq); + E=Fp_factorgalois(P,l,d,vq,FP); E=polpol_to_mat(E,np); MP = matrixpow(np,d,SP,P,l); IR = (GEN)FpM_sindexrank(MP,l)[1]; @@ -1872,7 +2067,7 @@ Fp_factor_irred(GEN P,GEN l, GEN Q) V = cgetg(d+1,t_VEC); V[1]=(long)M; 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); for(i=1;i<=d;i++) res[i]=(long)mat_to_polpol((GEN)V[i],vp,vq); @@ -1880,7 +2075,7 @@ Fp_factor_irred(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; long n,nbfact=0,nmax=lgef(P)-2; long i; @@ -1909,7 +2104,7 @@ GEN Fp_factor_rel0(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 j; GEN y,u,v; @@ -2011,15 +2206,9 @@ FpV_red(GEN z, GEN p) GEN FpM_red(GEN z, GEN p) { - long i,j,l = lg(z), m = lg((GEN)z[1]); - GEN x,y; - x = cgetg(l,t_MAT); - for (i=1; i=0 && !c[i]) i--; c = u_normalizepol(c-2, i+3); - if (pr == ONLY_REM) { free(q); return c; } *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 */ GEN 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; if (!p) return poldivres(x,y,pr); 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 (pr) @@ -2341,16 +2610,17 @@ FpX_divres(GEN x, GEN y, GEN p, GEN *pr) if (OK_ULONG(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); + GEN a = u_Fp_FpX(x, pp); + GEN b = u_Fp_FpX(y, pp); + 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) { - rem = small_to_pol(*pr,vx); - free(*pr); *pr = rem; + setlg(*pr, lgef(*pr)); *pr = dummycopy(*pr); + *pr = small_to_pol_ip(*pr,vx); } - z = small_to_pol(zz,vx); - free(zz); free(a); free(b); return z; + return small_to_pol_ip(z,vx); } lead = gcmp1(lead)? NULL: gclone(mpinvmod(lead,p)); avma = av0; @@ -2370,7 +2640,7 @@ FpX_divres(GEN x, GEN y, GEN p, GEN *pr) } 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--) { p1 = (GEN)x[i]; @@ -2384,12 +2654,12 @@ FpX_divres(GEN x, GEN y, GEN p, GEN *pr) { if (lead) gunclone(lead); if (sx) { avma=av0; return NULL; } - avma = (long)rem; return z-2; + 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((long)rem,tetpil,p1); + p1 = gerepile((gpmem_t)rem,tetpil,p1); rem += 2; rem[i]=(long)p1; for (i--; i>=0; i--) { @@ -2400,7 +2670,7 @@ FpX_divres(GEN x, GEN y, GEN p, GEN *pr) } rem -= 2; 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); *pr = rem; return z-2; } @@ -2409,7 +2679,8 @@ FpX_divres(GEN x, GEN y, GEN p, GEN *pr) GEN 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; if (!p) return poldivres(x,y,pr); @@ -2443,17 +2714,6 @@ FpXQX_divres(GEN x, GEN y, GEN T, GEN p, GEN *pr) #if 0 /* FIXME: to be done */ if (OK_ULONG(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 lead = gcmp1(lead)? NULL: gclone(FpXQ_inv(lead,T,p)); @@ -2474,7 +2734,7 @@ FpXQX_divres(GEN x, GEN y, GEN T, GEN p, GEN *pr) } 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--) { p1 = (GEN)x[i]; @@ -2488,12 +2748,12 @@ FpXQX_divres(GEN x, GEN y, GEN T, GEN p, GEN *pr) { if (lead) gunclone(lead); if (sx) { avma=av0; return NULL; } - avma = (long)rem; return z-2; + 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((long)rem,tetpil,p1); + p1 = gerepile((gpmem_t)rem,tetpil,p1); rem += 2; rem[i]=(long)p1; for (i--; i>=0; i--) { @@ -2504,21 +2764,132 @@ FpXQX_divres(GEN x, GEN y, GEN T, GEN p, GEN *pr) } rem -= 2; 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); *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=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 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; (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); } - b = u_Fp_FpX(y,0, pp); + b = u_Fp_FpX(y, pp); a = u_FpX_gcd_i(a,b, pp); avma = av; return small_to_pol(a, varn(x)); } @@ -2528,7 +2899,7 @@ GEN FpX_gcd(GEN x, GEN y, GEN p) { GEN a,b,c; - long av0,av; + gpmem_t av0, av; if (OK_ULONG(p)) return FpX_gcd_long(x,y,p); av0=avma; @@ -2541,7 +2912,28 @@ FpX_gcd(GEN x, GEN y, GEN p) avma = av; return gerepileupto(av0, a); } +/*Return 1 if gcd can be computed + * else return a factor of p*/ + 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) { long i,lz,lx = lgef(x), ly = lgef(y); @@ -2564,7 +2956,7 @@ u_FpX_sub(GEN x, GEN y, ulong p) /* list of u_FpX in gptr, return then as GEN */ 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)); long i; @@ -2586,11 +2978,11 @@ u_FpX_extgcd(GEN a, GEN b, ulong p, GEN *ptu, GEN *ptv { GEN q,z,u,v, x = a, y = b; - u = u_zeropol(0); - v= u_scalarpol(1, 0); /* v = 1 */ + u = u_zeropol(); + v= u_scalarpol(1); /* v = 1 */ 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) */ z = u_FpX_sub(u, u_FpX_mul(q,v, p), p); u = v; v = z; /* (u,v) = (v, u - q v) */ @@ -2605,11 +2997,11 @@ static GEN FpX_extgcd_long(GEN x, GEN y, GEN p, GEN *ptu, GEN *ptv) { ulong pp = (ulong)p[2]; - long av = avma; + gpmem_t av = avma; GEN a, b, d; - a = u_Fp_FpX(x,0, pp); - b = u_Fp_FpX(y,0, pp); + a = u_Fp_FpX(x, pp); + b = u_Fp_FpX(y, pp); d = u_FpX_extgcd(a,b, pp, ptu,ptv); { GEN *gptr[3]; gptr[0] = &d; gptr[1] = ptu; gptr[2] = ptv; @@ -2620,11 +3012,12 @@ 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 * ux + vy = gcd (mod p) */ +/*TODO: Document the typ() of *ptu and *ptv*/ GEN FpX_extgcd(GEN x, GEN y, GEN p, GEN *ptu, GEN *ptv) { 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); ltop=avma; @@ -2658,7 +3051,7 @@ GEN FpXQX_extgcd(GEN x, GEN y, GEN T, GEN p, GEN *ptu, GEN *ptv) { GEN a,b,q,r,u,v,d,d1,v1; - long ltop,lbot; + gpmem_t ltop, lbot; #if 0 /* FIXME To be done...*/ if (OK_ULONG(p)) return FpXQX_extgcd_long(x,y,T,p,ptu,ptv); @@ -2689,20 +3082,23 @@ FpXQX_extgcd(GEN x, GEN y, GEN T, GEN p, GEN *ptu, GEN *ptu = u; *ptv = v; return d; } -extern GEN caractducos(GEN p, GEN x, int v); - +/*x must be reduced*/ GEN FpXQ_charpoly(GEN x, GEN T, GEN p) { - ulong ltop=avma; - GEN R=lift(caractducos(FpX(T,p),FpX(x,p),varn(T))); + gpmem_t ltop=avma; + 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); } GEN FpXQ_minpoly(GEN x, GEN T, GEN p) { - ulong ltop=avma; + gpmem_t ltop=avma; GEN R=FpXQ_charpoly(x, T, p); GEN G=FpX_gcd(R,derivpol(R),p); G=FpX_normalize(G,p); @@ -2714,7 +3110,8 @@ FpXQ_minpoly(GEN x, GEN T, GEN p) static GEN 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; if (b == amod) return NULL; @@ -2725,10 +3122,10 @@ u_chrem_coprime(GEN a, ulong b, GEN q, ulong p, ulong } /* centerlift(u mod p) */ -GEN +long u_center(ulong u, ulong p, ulong ps2) { - return stoi((long) (u > ps2)? u - p: u); + return (long) (u > ps2)? u - p: u; } GEN @@ -2738,7 +3135,7 @@ ZX_init_CRT(GEN Hp, ulong p, long v) GEN H = cgetg(l, t_POL); H[1] = evalsigne(1)|evallgef(l)|evalvarn(v); for (i=2; i 0) h = subii(h,qp); + x[i] = (long)h; + } + setlgef(x,lp); *ptH = H = x; + stable = 0; lp = l; + } for (i=2; i>=1; if (!n) return y; - z = mulssmod(z,z,p); - } -} - /* separate from u_FpX_divrem for maximal speed. Implicitly malloc = 0 */ GEN u_FpX_rem(GEN x, GEN y, ulong p) @@ -2873,16 +3264,16 @@ u_FpX_rem(GEN x, GEN y, ulong p) long dx,dy,dz,i,j; ulong p1,inv; - dy = degpol(y); if (!dy) return u_zeropol(0); + dy = degpol(y); if (!dy) return u_zeropol(); 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; y += 2; - z = u_allocpol(dz, 1) + 2; + z = u_mallocpol(dz) + 2; inv = y[dy]; if (inv != 1UL) inv = u_invmod(inv,p); - c = u_allocpol(dy,0) + 2; + c = u_getpol(dy) + 2; if (u_OK_ULONG(p)) { z[dz] = (inv*x[dx]) % p; @@ -2892,7 +3283,7 @@ u_FpX_rem(GEN x, GEN y, ulong p) for (j=i-dy+1; j<=i && j<=dz; j++) { p1 += z[j]*y[i-j]; - if (p1 & HIGHBIT) p1 = p1 % p; + if (p1 & HIGHBIT) p1 %= p; } p1 %= p; z[i-dy] = p1? ((p - p1)*inv) % p: 0; @@ -2903,7 +3294,7 @@ u_FpX_rem(GEN x, GEN y, ulong p) for (j=1; j<=i && j<=dz; 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); } @@ -2934,7 +3325,8 @@ ulong u_FpX_resultant(GEN a, GEN b, ulong p) { long da,db,dc,cnt; - ulong lb,av, res = 1UL; + ulong lb, res = 1UL; + gpmem_t av; GEN c; if (!signe(a) || !signe(b)) return 0; @@ -2943,7 +3335,7 @@ u_FpX_resultant(GEN a, GEN b, ulong p) if (db > da) { 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 */ cnt = 0; av = avma; @@ -2954,7 +3346,7 @@ u_FpX_resultant(GEN a, GEN b, ulong p) a = b; b = c; dc = degpol(c); 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 (++cnt == 4) { cnt = 0; avma = av; } da = db; /* = degpol(a) */ @@ -2963,12 +3355,60 @@ u_FpX_resultant(GEN a, GEN b, ulong 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 */ ulong u_FpX_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GEN *ptV) { 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; if (!signe(x) || !signe(y)) return 0; @@ -2978,29 +3418,29 @@ u_FpX_extresultant(GEN a, GEN b, ulong p, GEN *ptU, GE { swap(x,y); lswap(dx,dy); pswap(ptU, ptV); a = x; b = y; - if (dx & dy & 1) res = p-res; + if (both_odd(dx,dy)) res = p-res; } - u = u_zeropol(0); - v = u_scalarpol(1, 0); /* v = 1 */ + u = u_zeropol(); + v = u_scalarpol(1); /* v = 1 */ while (dy) { /* b u = x (a), b v = y (a) */ 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) */ dz = degpol(z); if (dz < 0) { avma = av; return 0; } z = u_FpX_sub(u, u_FpX_mul(q,v, p), p); 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); dx = dy; /* = degpol(x) */ dy = dz; /* = degpol(y) */ } res = mulssmod(res, powuumod(y[2], dx, 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; - 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 */ *ptU = u; *ptV = v; return res; @@ -3012,7 +3452,8 @@ ulong u_FpX_resultant_all(GEN a, GEN b, long *C0, long *C1, GEN dglist, ulong p) { long da,db,dc,cnt,ind; - ulong lb,av,cx = 1, res = 1UL; + ulong lb, cx = 1, res = 1UL; + gpmem_t av; GEN c; if (C0) { *C0 = 1; *C1 = 0; } @@ -3022,7 +3463,7 @@ u_FpX_resultant_all(GEN a, GEN b, long *C0, long *C1, if (db > da) { 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 */ if (!da) return 1; @@ -3039,7 +3480,7 @@ u_FpX_resultant_all(GEN a, GEN b, long *C0, long *C1, { /* check that Euclidean remainder sequence doesn't degenerate */ if (dc != dglist[ind]) { avma = av; return 0; } /* update resultant */ - if (da & db & 1) res = p-res; + if (both_odd(da,db)) res = p-res; if (lb != 1) { ulong t = powuumod(lb, da - dc, p); @@ -3089,7 +3530,9 @@ u_FpV_roots_to_pol(GEN a, ulong p) { p2 = cgetg(5,t_VECSMALL); p1[k++] = (long)p2; 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); } if (i < lx) @@ -3123,7 +3566,7 @@ static GEN u_pol_comp(GEN P, ulong u, ulong v, ulong p) { long i, l = lgef(P); - GEN y = u_allocpol(l-3, 0); + GEN y = u_getpol(l-3); for (i=2; i1; 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 */ GEN u_FpV_polint(GEN xa, GEN ya, ulong p) { GEN T,dP, P = NULL, Q = u_FpV_roots_to_pol(xa, p); long i, n = lg(xa); - ulong av, inv; - av = avma; (void)new_chunk(n+3); /* storage space for P */ + ulong inv; + gpmem_t av = avma; for (i=1; i1) err(warnmem,"FpV_polint"); + if (!P) avma = av; else P = gerepileupto(av, P); + } + } + return P? P: zeropol(0); +} + static void u_FpV_polint_all(GEN xa, GEN ya, GEN C0, GEN C1, ulong p) { @@ -3274,53 +3763,38 @@ u_FpV_polint_all(GEN xa, GEN ya, GEN C0, GEN C1, ulong T = u_FpX_div_by_X_x(Q, xa[i], p); inv = u_invmod(u_FpX_eval(T,xa[i], p), p); -#if 0 - if (i < n-1 && (ulong)(xa[i] + xa[i+1]) == p) + if (ya[i]) { - if (ya[i]) - dP = u_pol_comp(T, mulssmod(ya[i],inv,p), mulssmod(ya[i+1],inv,p), p); - 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} */ + dP = u_FpX_Fp_mul(T, mulssmod(ya[i],inv,p), p); + P = P ? u_FpX_add(P , dP , p): dP; } - else -#endif + if (C0[i]) { - if (ya[i]) - { - 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; - } + dP0= u_FpX_Fp_mul(T, mulssmod(C0[i],inv,p), p); + P0= P0? u_FpX_add(P0, dP0, p): dP0; } + 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)); - C0[1] = (long) (P0? P0: u_zeropol(0)); - C1[1] = (long) (P1? P1: u_zeropol(0)); + ya[1] = (long) (P ? P : u_zeropol()); + C0[1] = (long) (P0? P0: u_zeropol()); + C1[1] = (long) (P1? P1: u_zeropol()); } /* b a vector of polynomials representing B in Fp[X][Y], evaluate at X = x, * Return 0 in case of degree drop. */ static GEN -vec_u_FpX_eval(GEN b, ulong x, ulong p) +u_vec_FpX_eval(GEN b, ulong x, ulong p) { GEN z; long i, lb = lgef(b); 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>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; i1) 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= pp) + { + l = lgef(a); + a[0] = evaltyp(t_POL) | evallg(l); + a[1] = evalsigne(1)|evalvarn(vY)|evallgef(l); + for (i=2; i>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, ... */ #define next_lambda(a) (a>0 ? -a : 1-a) @@ -3372,8 +4131,9 @@ ZY_ZXY_ResBound(GEN A, GEN B) GEN ZY_ZXY_resultant_all(GEN A, GEN B0, long *lambda, GEN *LERS) { - int checksqfree = lambda? 1: 0, delete = 0, first = 1, stable; - ulong av = avma, av2, lim, bound; + int checksqfree = lambda? 1: 0, delvar = 0, first = 1, stable; + ulong bound; + gpmem_t av = avma, av2, lim; long i,n, lb, dres = degpol(A)*degpol(B0), nmax = (dres+1)>>1; 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; @@ -3392,7 +4152,7 @@ ZY_ZXY_resultant_all(GEN A, GEN B0, long *lambda, GEN y = cgetg(dres+2, t_VECSMALL); if (vY == MAXVARN) { - vY = fetch_var(); delete = 1; + vY = fetch_var(); delvar = 1; B0 = gsubst(B0, MAXVARN, polx[vY]); A = dummycopy(A); setvarn(A, vY); } @@ -3433,18 +4193,22 @@ INIT: goto END; } + /* make sure p large enough */ + while (p < (dres<<1)) NEXT_PRIME_VIADIFF(p,d); + H = H0 = H1 = NULL; - lb = lgef(B); b = u_allocpol(degpol(B), 0); - bound = ZY_ZXY_ResBound(A,B); + lb = lgef(B); b = u_getpol(degpol(B)); + bound = ZY_ZXY_ResBound(A, B); if (DEBUGLEVEL>4) fprintferr("bound for resultant coeffs: 2^%ld\n",bound); + check_theta(bound); 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; i1) 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: - 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 (LERS) { @@ -3562,9 +4319,9 @@ END: } 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]. @@ -3573,10 +4330,10 @@ ZY_ZXY_resultant(GEN A, GEN B0, long *lambda) GEN ZX_caract_sqf(GEN A, GEN B, long *lambda, long v) { - ulong av = avma; + gpmem_t av = avma; GEN B0, R, a; long dB; - int delete; + int delvar; if (v < 0) v = 0; switch (typ(B)) @@ -3587,10 +4344,10 @@ ZX_caract_sqf(GEN A, GEN B, long *lambda, long v) if (lambda) { B = scalarpol(B,varn(A)); dB = 0; break;} return gerepileupto(av, gpowgs(gsub(polx[v], B), degpol(A))); } - delete = 0; + delvar = 0; if (varn(A) == 0) { - long v0 = fetch_var(); delete = 1; + long v0 = fetch_var(); delvar = 1; A = dummycopy(A); setvarn(A,v0); B = dummycopy(B); setvarn(B,v0); } @@ -3598,30 +4355,49 @@ ZX_caract_sqf(GEN A, GEN B, long *lambda, long v) B0[2] = (long)gneg_i(B); B0[3] = un; R = ZY_ZXY_resultant(A, B0, lambda); - if (delete) delete_var(); + if (delvar) (void)delete_var(); setvarn(R, v); a = leading_term(A); if (!gcmp1(a)) R = gdiv(R, gpowgs(a, dB)); return gerepileupto(av, R); } + GEN 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 trivial_case(GEN A, GEN B) { + long d; 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; } +/* Res(A, B/dB), assuming the A,B in Z[X] and result is integer */ 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; GEN q, a, b, H; byteptr d = diffptr + 3000; @@ -3630,19 +4406,36 @@ ZX_resultant_all(GEN A, GEN B, ulong bound) if ((H = trivial_case(A,B)) || (H = trivial_case(B,A))) return H; q = H = NULL; 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); + check_theta(bound); + dp = 0; /* gcc -Wall */ for(;;) { - p += *d++; if (!*d) err(primer1); - a = u_Fp_FpX(A, 0, p); - b = u_Fp_FpX(B, 0, p); + NEXT_PRIME_VIADIFF_CHECK(p,d); + if (dB) { dp = smodis(dB, p); if (!dp) continue; } + + a = u_Fp_FpX(A, p); + b = u_Fp_FpX(B, p); Hp= u_FpX_resultant(a, b, p); + if (dp) Hp = mulssmod(Hp, u_powmod(dp, -degA, p), p); if (!H) { 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] */ { @@ -3664,15 +4457,29 @@ ZX_resultant_all(GEN A, GEN B, ulong bound) } 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 */ GEN ZX_disc_all(GEN x, ulong bound) { - ulong av = avma; - GEN l, d = ZX_resultant_all(x, derivpol(x), bound); - l = leading_term(x); if (!gcmp1(l)) d = divii(d,l); + gpmem_t av = avma; + GEN l, d = ZX_resultant_all(x, derivpol(x), NULL, bound); + l = leading_term(x); if (!gcmp1(l)) d = diviiexact(d,l); if (degpol(x) & 2) d = negi(d); return gerepileuptoint(av,d); } @@ -3681,41 +4488,37 @@ GEN ZX_disc(GEN x) { return ZX_disc_all(x,0); } int ZX_is_squarefree(GEN x) { - ulong av = avma; + gpmem_t av = avma; int d = (lgef(modulargcd(x,derivpol(x))) == 3); avma = av; return d; } -/* h integer, P in Z[X]. Return h^degpol(P) P(x / h) */ -GEN -ZX_rescale_pol(GEN P, GEN h) +static GEN +_gcd(GEN a, GEN b) { - long i, l = lgef(P); - GEN Q = cgetg(l,t_POL), hi = gun; - Q[l-1] = P[l-1]; - for (i=l-2; i>=2; i--) - { - hi = gmul(hi,h); Q[i] = lmul((GEN)P[i], hi); - } - Q[1] = P[1]; return Q; + if (!a) a = gun; + if (!b) b = gun; + return ggcd(a,b); } /* A0 and B0 in Q[X] */ GEN 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; - 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; if ((typ(A0) | typ(B0)) !=t_POL) err(notpoler,"modulargcd"); if (!signe(A0)) return gcopy(B0); if (!signe(B0)) return gcopy(A0); - A = content(A0); - B = content(B0); D = ggcd(A,B); - A = gcmp1(A)? A0: gdiv(A0,A); - B = gcmp1(B)? B0: gdiv(B0,B); + A = primitive_part(A0, &a); check_pol_int(A, "modulargcd"); + B = primitive_part(B0, &b); check_pol_int(B, "modulargcd"); + D = _gcd(a,b); + if (varn(A) != varn(B)) err(talker,"different variables in modulargcd"); + /* A, B in Z[X] */ g = mppgcd(leading_term(A), leading_term(B)); /* multiple of lead(gcd) */ if (degpol(A) < degpol(B)) swap(A, B); @@ -3723,35 +4526,35 @@ modulargcd(GEN A0, GEN B0) av2 = avma; H = NULL; d += 3000; /* 27449 = prime(3000) */ - for(p = 27449; ; p+= *d++) + for(p = 27449; ;) { if (!*d) err(primer1); - if (!umodiu(g,p)) continue; + if (!umodiu(g,p)) goto repeat; - a = u_Fp_FpX(A, 0, p); - b = u_Fp_FpX(B, 0, p); Hp = u_FpX_gcd_i(a,b, p); + a = u_Fp_FpX(A, p); + b = u_Fp_FpX(B, p); Hp = u_FpX_gcd_i(a,b, p); m = degpol(Hp); 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 */ Hp = u_FpX_normalize(Hp, p); else { 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) { /* First time or degree drop [all previous p were as above; restart]. */ 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); - if (ZX_incremental_CRT(H, Hp, q, qp, p)) + if (ZX_incremental_CRT(&H, Hp, q, qp, p)) { /* H stable: check divisibility */ - if (!is_pm1(g)) { p1 = content(H); if (!is_pm1(p1)) H = gdiv(H,p1); } - if (!signe(gres(A,H)) && !signe(gres(B,H))) break; /* DONE */ + if (!is_pm1(g)) H = primpart(H); + if (gcmp0(pseudorem(A,H)) && gcmp0(pseudorem(B,H))) break; /* DONE */ if (DEBUGLEVEL) fprintferr("modulargcd: trial division failed"); } @@ -3762,65 +4565,238 @@ modulargcd(GEN A0, GEN B0) if (DEBUGMEM>1) err(warnmem,"modulargcd"); gerepilemany(av2,gptr,2); } + repeat: + NEXT_PRIME_VIADIFF_CHECK(p,d); } 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 -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; 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; - if (typ(B0) != t_POL) err(notpoler,"ZX_invmod"); + if (typ(B0) != t_POL) err(notpoler,"QX_invmod"); if (typ(A0) != t_POL) { if (is_scalar_t(typ(A0))) return ginv(A0); - err(notpoler,"ZX_invmod"); + err(notpoler,"QX_invmod"); } - A = content(A0); D = A; - B = content(B0); - A = gcmp1(A)? A0: gdiv(A0,A); - B = gcmp1(B)? B0: gdiv(B0,B); + if (degpol(A0) < 15) return ginvmod(A0,B0); + A = primitive_part(A0, &D); + B = primpart(B0); /* A, B in Z[X] */ av2 = avma; U = NULL; d += 3000; /* 27449 = prime(3000) */ - for(p = 27449; ; p+= *d++) + for(p = 27449; ; ) { if (!*d) err(primer1); - a = u_Fp_FpX(A, 0, p); - b = u_Fp_FpX(B, 0, p); + a = u_Fp_FpX(A, p); + b = u_Fp_FpX(B, p); /* 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) { /* First time */ U = ZX_init_CRT(Up,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); - stable = ZX_incremental_CRT(U, Up, q,qp, p); - stable &= ZX_incremental_CRT(V, Vp, q,qp, p); + stable = ZX_incremental_CRT(&U, Up, q,qp, p); + stable &= ZX_incremental_CRT(&V, Vp, q,qp, p); if (stable) { /* all stable: check divisibility */ res = gadd(gmul(A,U), gmul(B,V)); 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; if (low_stack(avlim, stack_lim(av,1))) { 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); } + repeat: + NEXT_PRIME_VIADIFF_CHECK(p,d); } - D = gmul(D,res); + D = D? gmul(D,res): res; 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 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