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

Return to polarit3.c CVS log

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

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

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

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

Line 33 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) */

Line 80 specpol(GEN x, long nx)

Line 85 specpol(GEN x, long nx)

/* multiplication in Fp[X], p small */

static GEN

GEN

u_normalizepol(GEN x, long lx)

{

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

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);

}

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;

}

#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<l; i++) z[i] = x[i]? p - x[i]: 0;

return z;

}

static GEN

u_FpX_neg_i(GEN x, ulong p)

{

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 */

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;

GEN y;

Line 278 GEN

Line 308 GEN

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

{

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 (nb && !b[0]) { b++; nb--; v++; }

if (na < nb) swapspec(a,b, na,nb);

if (!nb) return u_zeropol(0);

if (!nb) return u_zeropol();

av = avma;

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;

GEN z;

if (!nx) return u_zeropol(0);

if (!nx) return u_zeropol();

lz = (nx << 1) + 1, nz = lz-2;

z = cgetg(lz,t_VECSMALL) + 2;

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

Line 378 GEN

Line 409 GEN

u_FpX_Ksqr(GEN a, ulong p, long na)

{

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; }

av = avma;

Line 434 GEN

Line 466 GEN

mulpol_limb(GEN x, GEN y, char *ynonzero, long a, long b)

{

GEN p1 = NULL;

long i,av = avma;

long i;

gpmem_t av = avma;

for (i=a; i<b; i++)

if (ynonzero[i])

{

Line 554 GEN

Line 587 GEN

quickmul(GEN a, GEN b, long na, long nb)

{

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 (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

sqrpol(GEN x, long nx)

{

long av,i,j,l,lz,nz;

long i, j, l, lz, nz;

gpmem_t av;

GEN p1,z;

char *p2;

Line 635 GEN

Line 670 GEN

quicksqr(GEN a, long na)

{

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; }

if (v) (void)new_chunk(v);

Line 663 There are clean and memory efficient.

Line 699 There are clean and memory efficient.

GEN

FpX_center(GEN T,GEN mod)

{/*OK centermod exists, but is not so clean*/

ulong av;

gpmem_t av;

long i, l=lg(T);

GEN P,mod2;

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<lx; i++) z[i]=licopy((GEN)x[i]);

(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);

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]);

/*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);

return z;

}

Line 758 FpX_sqr(GEN x,GEN p)

Line 794 FpX_sqr(GEN x,GEN p)

if (!p) return z;

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) */

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);

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);

return FpX_Fp_mul(U, z, p);

}

/* Product of y and x in Z/pZ[X]/(T)

* 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

FpXQ_mul(GEN y,GEN x,GEN T,GEN p)

{

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 = 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;

}

/* 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))

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. *

* *

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

/* 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

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

FpXQ_inv(GEN x,GEN T,GEN p)

{

ulong av;

gpmem_t av;

GEN U;

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

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;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;

}

/*Close to FpXV_FpV_dotproduct*/

/*Close to FpXV_FpV_innerprod*/

static GEN

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

FpX_FpXQV_compo(GEN P, GEN V, GEN T, GEN p)

{

ulong ltop=avma;

gpmem_t ltop=avma;

long l=lg(V)-1;

GEN z,u;

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

FpX_FpXQ_compo(GEN T,GEN x,GEN pol,GEN p)

{

ulong ltop=avma;

gpmem_t ltop=avma;

GEN z;

long d=degpol(T),rtd;

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

FpX_eval(GEN x,GEN y,GEN p)

{

ulong av;

gpmem_t av;

GEN p1,r,res;

long i,j;

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--)

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);

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);

}

fppoleval:

Line 1056 FpX_eval(GEN x,GEN y,GEN p)

Line 1139 FpX_eval(GEN x,GEN y,GEN p)

GEN

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;

ax = FpX_mul(FpXQ_inv(Tx,Ty,p), Tx,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);

}

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 */

GEN

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);

}

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

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) < 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))

{

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 = small_to_pol(y,vx);

}

else

{

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] */

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

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;i<lgef(res);i++)

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] */

/* */

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

extern GEN to_Kronecker(GEN P, GEN Q);

GEN /*Somewhat copy-pasted...*/

/*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 (j=2; j<N; j++) t[j] = z[j];

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;

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);

}

/*Unused/untested*/

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)

{

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));

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)

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

FpXQX_mul(GEN x, GEN y, GEN T, GEN p)

{

ulong ltop;

gpmem_t ltop;

GEN z,kx,ky;

long vx;

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);

ky= to_Kronecker(y,T);

z = quickmul(ky+2, kx+2, lgef(ky)-2, lgef(kx)-2);

z = FpX_red(z,p);

z = FpXQX_from_Kronecker(z,T,p);

setvarn(z,vx);/*quickmul and Fq_from_Kronecker are not varn-clean*/

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*/

FpXQX_sqr(GEN x, GEN T, GEN p)

{

ulong ltop=avma;

gpmem_t ltop=avma;

GEN z,kx;

long vx=varn(x);

kx= to_Kronecker(x,T);

z = quicksqr(kx+2, lgef(kx)-2);

z = FpX_red(z,p);

z = FpXQX_from_Kronecker(z,T,p);

setvarn(z,vx);/*quickmul and Fq_from_Kronecker are nor varn-clean*/

return gerepileupto(ltop,z);

}

GEN

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);

res[1] = evalsigne(1) | evalvarn(varn(P)) | evallgef(lP);

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));

}

Line 1280 monomial(GEN a, int degpol, int v)

Line 1391 monomial(GEN a, int degpol, int v)

GEN

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);

GEN U, q;

P = FpXX_red(P, p);

P = FpXX_red(P, p); btop = avma;

Q = FpXX_red(Q, p);

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);

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; }

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 = FpXQX_red(P, T, p); /* wasteful, but negligible */

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] */

/* */

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

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_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)

{

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

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<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*/

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);

if (gcmp1(m)) { avma = av1; continue; }

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;

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

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);

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);

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 */

z = p1;

{

p1 = FpXQ_pow(p1,l,T,p);

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

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);

e = k;

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);

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;

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)

{

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)

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))

{

m=decomp(m);

m = decomp(m); av1 = avma;

av1=avma;

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();

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)))

/* n can have lots of prime factors*/

{ /* n can have lots of prime factors */

{

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)

{

*zetan=gcopy(z);

*zetan = z;

gptr[0]=&a;gptr[1]=zetan;

gerepileall(ltop,2,&a,zetan);

gerepilemanysp(ltop,lbot,gptr,2);

}

else

a=gerepileupto(ltop,a);

a = gerepileupto(ltop, 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

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;

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);

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

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[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<lg(PV);i++)

PV[i]=P[1+i];

M=cgetg(l+1,t_VEC);

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];

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

{

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);

}

/* 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

* 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:

* 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.

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 vu=varn(lU), r=degpol(lU);

long vu=varn(U), r=degpol(U);

long i;

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]");

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");

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

* 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<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=gmul(A,gmodulcp(gmodulcp(gun,l),U));

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

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;

GEN q;

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);

pg=itos(q);

avma=ltop;

if (DEBUGLEVEL>=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...)*/

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);

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");

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

{

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);

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]");

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 ( !L )

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)

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;

}

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);

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;

}

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*/

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. */

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);

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);

}

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 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;

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 i,d=cgcd(np,nq);

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);

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;

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];

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[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);

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)

{

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;

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)

{

long tetpil,av=avma;

gpmem_t tetpil, av=avma;

long nbfact;

long j;

GEN y,u,v;

Line 2011 FpV_red(GEN z, GEN p)

Line 2206 FpV_red(GEN z, GEN p)

GEN

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;

}

Line 2042 FpXQX_normalize(GEN z, GEN T, GEN p)

Line 2231 FpXQX_normalize(GEN z, GEN T, GEN p)

}

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)

{

long i, l = lg(z);

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;

}

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));

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

pol_to_small(GEN 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]);

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

* If p = NULL, use generic functions and the coeff. ring implied by the

* coefficients of z */

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);

t[1] = pol[1] & VARNBITS;

t[1] = T[1] & VARNBITS;

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++)

{

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];

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[1] = (long)pol;

a[1] = (long)T;

N = (l-2) % (N-2) + 2;

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);

}

#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)))

* Recover the "real" z, normalized */

GEN

Line 2117 from_Kronecker(GEN z, GEN pol)

Line 2361 from_Kronecker(GEN z, GEN pol)

/* 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 */

ulong

u_invmod(ulong x, ulong p)

Line 2129 u_invmod(ulong x, ulong p)

Line 2373 u_invmod(ulong x, ulong p)

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

static ulong

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 */

GEN

u_Fp_FpX(GEN x, int malloc, ulong p)

u_Fp_FpX(GEN x, ulong p)

{

long i, lx;

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))

{

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;

}

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);

return u_normalizepol(a,lx);

}

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)

{

long i,j,m,n = lg(x);

long j,n = lg(x);

GEN y = cgetg(n,t_MAT);

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;

}

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;

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))

for(i=2; i<l; i++) z[i] = mulssmod(y[i], x, p);

else

Line 2196 u_FpX_normalize(GEN z, ulong p)

Line 2463 u_FpX_normalize(GEN z, ulong p)

long l = lgef(z)-1;

ulong p1 = z[l]; /* leading term */

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

u_copy(GEN x, int malloc)

u_copy(GEN 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];

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

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;

long dx,dy,dz,i,j;

ulong p1,inv;

if (pr == ONLY_REM) return u_FpX_rem(x, y, p);

dy = degpol(y);

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);

dz = dx-dy;

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;

y += 2;

z = u_allocpol(dz, malloc || (pr == ONLY_REM)) + 2;

z = u_getpol(dz) + 2;

inv = y[dy];

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++)

{

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;

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);

if (!pr) return q;

c = u_allocpol(dy,malloc) + 2;

c = u_getpol(dy) + 2;

if (u_OK_ULONG(p))

{

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++)

{

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);

}

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--;

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)

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))

{ /* assume ab != 0 mod p */

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)

{

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));

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; }

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];

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 (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--)

{

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;

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;

}

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

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);

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 (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));

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; }

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];

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 (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--)

{

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;

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<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

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));

}

Line 2528 GEN

Line 2899 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;

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);

}

/*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);

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 */

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;

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;

u = u_zeropol(0);

u = u_zeropol();

v= u_scalarpol(1, 0); /* v = 1 */

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) */

Line 2605 static GEN

Line 2997 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);

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);

{

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

* 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;

Line 2658 GEN

Line 3051 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);

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;

}

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;

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);

}

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);

Line 2714 FpXQ_minpoly(GEN x, GEN T, GEN p)

Line 3110 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;

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) */

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

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);

H[1] = evalsigne(1)|evallgef(l)|evalvarn(v);

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;

}

Line 2753 ZM_init_CRT(GEN Hp, ulong p)

Line 3150 ZM_init_CRT(GEN Hp, ulong p)

cp = (GEN)Hp[j];

c = cgetg(m, t_COL);

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;

}

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

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);

long i, l = lgef(H), lp = lgef(Hp);

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++)

{

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;

}

/* 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 */

GEN

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;

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;

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++)

{

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;

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++)

{

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);

}

Line 2934 ulong

Line 3325 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;

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)

{

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;

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);

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) */

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);

}

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;

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);

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)

{ /* 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;

Line 3012 ulong

Line 3452 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; }

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)

{

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;

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 */

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);

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[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)

Line 3123 static GEN

Line 3566 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; i<l; 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);

}

GEN roots_to_pol(GEN a, long v);

extern GEN roots_to_pol(GEN a, long v);

GEN

polint_triv(GEN xa, GEN ya)

{

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++)

{

GEN T,dP;

Line 3212 static GEN

Line 3656 static GEN

u_FpX_div_by_X_x(GEN a, ulong x, ulong p)

{

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;

z0 = z + l-2; *z0 = *a0--;

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;

}

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 */

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;

ulong inv;

av = avma; (void)new_chunk(n+3); /* storage space for P */

gpmem_t av = avma;

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

{

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} */

}

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;

}

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

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);

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,

* 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<lb-1; i++)

z[i] = u_FpX_eval((GEN)b[i], x, p);

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 */

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;

long i, lb = lgef(b);

z = u_allocpol(lb-3, 0);

z = u_getpol(lb-3);

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

z[i] = u_FpX_eval((GEN)b[i], x, p);

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;

}

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:

* 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

ZY_ZXY_ResBound(GEN A, GEN B)

{

ulong av = avma;

gpmem_t av = avma;

GEN a = gzero, b = gzero, run = realun(DEFAULTPREC);

long i , lA = lgef(A), lB = lgef(B);

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);

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);

}

/* 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, ... */

#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

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 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;

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);

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);

}

Line 3433 INIT:

Line 4193 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);

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);

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; 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 (!b[lb-1] || degpol(a) < degpol(A)) continue; /* p | lc(A)lc(B) */

Line 3455 INIT:

Line 4219 INIT:

setlg(dglist, 1);

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);

if (lg(dglist)-1 == goal) break;

}

Line 3473 INIT:

Line 4237 INIT:

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;

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) */

Line 3484 INIT:

Line 4248 INIT:

H1p= (GEN)C1[1];

}

else

{

{ /* cf u_FpXY_resultant */

ulong la = (ulong)leading_term(a);

int drop;

/* 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++)

{

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)

{

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);

}

Line 3525 INIT:

Line 4282 INIT:

else

{

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) {

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;

}

Line 3541 INIT:

Line 4298 INIT:

{

GEN *gptr[4]; gptr[0] = &H; gptr[1] = &q; gptr[2] = &H0; gptr[3] = &H1;

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:

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)

{

Line 3562 END:

Line 4319 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].

Line 3573 ZY_ZXY_resultant(GEN A, GEN B0, long *lambda)

Line 4330 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))

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;}

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);

}

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[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;

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;

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);

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);

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] */

{

Line 3664 ZX_resultant_all(GEN A, GEN B, ulong bound)

Line 4457 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;

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);

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

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) */

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] */

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);

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] */

g = mppgcd(leading_term(A), leading_term(B)); /* multiple of lead(gcd) */

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;

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);

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);

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 (!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");

}

Line 3762 modulargcd(GEN A0, GEN B0)

Line 4565 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;

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] */

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);

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 (!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(&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)

{ /* 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<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;

}

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>