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

Return to polarit1.c CVS log

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

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

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

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

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

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

#include "pari.h"

extern GEN get_bas_den(GEN bas);

extern GEN get_mul_table(GEN x,GEN bas,GEN invbas,GEN *T);

extern GEN get_mul_table(GEN x,GEN bas,GEN invbas);

extern GEN pol_to_monic(GEN pol, GEN *lead);

#define lswap(x,y) { long _t=x; x=y; y=_t; }

#define swap(x,y) { GEN _t=x; x=y; y=_t; }

/* see splitgen() for how to use these two */

GEN

setloop(GEN a)

{

a=icopy(a); new_chunk(2); /* dummy to get two cells of extra space */

a=icopy(a); (void)new_chunk(2); /* dummy to get two cells of extra space */

return a;

}

Line 86 incloop(GEN a)

Line 89 incloop(GEN a)

int

gdivise(GEN x, GEN y)

{

long av=avma;

gpmem_t av=avma;

x=gmod(x,y); avma=av; return gcmp0(x);

}

int

poldivis(GEN x, GEN y, GEN *z)

{

long av = avma;

gpmem_t av = avma;

GEN p1 = poldivres(x,y,ONLY_DIVIDES);

if (p1) { *z = p1; return 1; }

avma=av; return 0;

Line 108 poldivis(GEN x, GEN y, GEN *z)

Line 111 poldivis(GEN x, GEN y, GEN *z)

* if z = ONLY_REM return remainder, otherwise return quotient

* if z != NULL set *z to remainder

* *z is the last object on stack (and thus can be disposed of with cgiv

* instead of gerepile)

* instead of gerepile) */

GEN

poldivres(GEN x, GEN y, GEN *pr)

{

ulong avy,av,av1;

gpmem_t avy, av, av1;

long ty=typ(y),tx,vx,vy,dx,dy,dz,i,j,sx,lrem;

int remainder;

GEN z,p1,rem,y_lead,mod;

GEN (*f)(GEN,GEN);

if (pr == ONLY_DIVIDES_EXACT)

f = gdiv;

{ f = gdivexact; pr = ONLY_DIVIDES; }

else

f = gdiv;

if (is_scalar_t(ty))

{

if (pr == ONLY_REM) return gzero;

Line 149 poldivres(GEN x, GEN y, GEN *pr)

Line 147 poldivres(GEN x, GEN y, GEN *pr)

}

return f(x,y);

}

if (!signe(y)) err(talker,"euclidean division by zero (poldivres)");

dy=degpol(y); y_lead = (GEN)y[dy+2];

if (!signe(y)) err(talker,"division by zero in poldivrem");

dy = degpol(y);

y_lead = (GEN)y[dy+2];

if (gcmp0(y_lead)) /* normalize denominator if leading term is 0 */

{

err(warner,"normalizing a polynomial with 0 leading term");

Line 170 poldivres(GEN x, GEN y, GEN *pr)

Line 169 poldivres(GEN x, GEN y, GEN *pr)

}

return f(x, constant_term(y));

}

dx=degpol(x);

dx = degpol(x);

if (vx>vy || dx<dy)

{

if (pr)

Line 181 poldivres(GEN x, GEN y, GEN *pr)

Line 180 poldivres(GEN x, GEN y, GEN *pr)

}

return zeropol(vy);

}

dz=dx-dy; av=avma; /* to avoid gsub's later on */

/* x,y in R[X], y non constant */

dz = dx-dy; av = avma;

p1 = new_chunk(dy+3);

for (i=2; i<dy+3; i++)

{

GEN p2 = (GEN)y[i];

/* gneg to avoid gsub's later on */

p1[i] = isexactzero(p2)? 0: (long)gneg_i(p2);

}

y = p1;

Line 198 poldivres(GEN x, GEN y, GEN *pr)

Line 200 poldivres(GEN x, GEN y, GEN *pr)

default: if (gcmp1(y_lead)) y_lead = NULL;

mod = NULL;

}

avy=avma; z=cgetg(dz+3,t_POL);

avy=avma;

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

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]; remainder = (pr == ONLY_REM);

p1 = (GEN)x[dx];

z[dz]=y_lead? (long)f(p1,y_lead): lcopy(p1);

z[dz] = y_lead? (long)f(p1,y_lead): lcopy(p1);

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

{

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

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

if (y[i-j]) p1 = gadd(p1, gmul((GEN)z[j],(GEN)y[i-j]));

if (y[i-j] && z[j] != zero) p1 = gadd(p1, gmul((GEN)z[j],(GEN)y[i-j]));

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

if (!remainder) p1 = avma==av1? gcopy(p1): gerepileupto(av1,p1);

if (isexactzero(p1)) { avma=av1; p1 = gzero; }

else

p1 = avma==av1? gcopy(p1): gerepileupto(av1,p1);

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

}

if (!pr) return gerepileupto(av,z-2);

rem = (GEN)avma; av1 = (long)new_chunk(dx+3);

rem = (GEN)avma; av1 = (gpmem_t)new_chunk(dx+3);

for (sx=0; ; i--)

{

p1 = (GEN)x[i];

/* we always enter this loop at least once */

for (j=0; j<=i && j<=dz; j++)

if (y[i-j]) p1 = gadd(p1, gmul((GEN)z[j],(GEN)y[i-j]));

if (y[i-j] && z[j] != zero) p1 = gadd(p1, gmul((GEN)z[j],(GEN)y[i-j]));

if (mod && avma==av1) p1 = gmul(p1,mod);

if (!gcmp0(p1)) { sx = 1; break; } /* remainder is non-zero */

if (!isinexactreal(p1) && !isexactzero(p1)) break;

Line 231 poldivres(GEN x, GEN y, GEN *pr)

Line 237 poldivres(GEN x, GEN y, GEN *pr)

if (pr == ONLY_DIVIDES)

{

if (sx) { avma=av; return NULL; }

avma = (long)rem;

avma = (gpmem_t)rem;

return gerepileupto(av,z-2);

}

lrem=i+3; rem -= lrem;

if (avma==av1) { avma = (long)rem; p1 = gcopy(p1); }

if (avma==av1) { avma = (gpmem_t)rem; p1 = gcopy(p1); }

else p1 = gerepileupto((long)rem,p1);

else p1 = gerepileupto((gpmem_t)rem,p1);

rem[0]=evaltyp(t_POL) | evallg(lrem);

rem[1]=evalsigne(1) | evalvarn(vx) | evallgef(lrem);

rem += 2;

Line 245 poldivres(GEN x, GEN y, GEN *pr)

Line 251 poldivres(GEN x, GEN y, GEN *pr)

{

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

for (j=0; j<=i && j<=dz; j++)

if (y[i-j]) p1 = gadd(p1, gmul((GEN)z[j],(GEN)y[i-j]));

if (y[i-j] && z[j] != zero) p1 = gadd(p1, gmul((GEN)z[j],(GEN)y[i-j]));

if (mod && avma==av1) p1 = gmul(p1,mod);

rem[i]=avma==av1? lcopy(p1):lpileupto(av1,p1);

}

rem -= 2;

if (!sx) normalizepol_i(rem, lrem);

if (!sx) (void)normalizepol_i(rem, lrem);

if (remainder) return gerepileupto(av,rem);

if (pr == ONLY_REM) return gerepileupto(av,rem);

z -= 2;

{

GEN *gptr[2]; gptr[0]=&z; gptr[1]=&rem;

Line 286 factmod_init(GEN *F, GEN pp, long *p)

Line 292 factmod_init(GEN *F, GEN pp, long *p)

f = gmul(f, mod(gun,pp));

if (!signe(f)) err(zeropoler,"factmod");

f = lift_intern(f); d = lgef(f);

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

if (typ(f[i])!=t_INT) err(impl,"factormod for general polynomials");

*F = f; return d-3;

}

Line 352 GEN

Line 358 GEN

rootmod2(GEN f, GEN pp)

{

GEN g,y,ss,q,r, x_minus_s;

long p,av = avma,av1,d,i,nbrac;

long p, d, i, nbrac;

gpmem_t av = avma, av1;

if (!(d = factmod_init(&f, pp, &p))) { avma=av; return cgetg(1,t_COL); }

if (!p) err(talker,"prime too big in rootmod2");

Line 389 rootmod2(GEN f, GEN pp)

Line 396 rootmod2(GEN f, GEN pp)

free(y); return g;

}

/* assume x reduced mod p, monic, squarefree. Return one root, or NULL if

* irreducible */

static GEN

quadsolvemod(GEN x, GEN p, int unknown)

{

GEN b = (GEN)x[3], c = (GEN)x[2];

GEN s,u,D;

if (egalii(p, gdeux))

{

if (signe(c)) return NULL;

return gun;

}

D = subii(sqri(b), shifti(c,2));

D = resii(D,p);

if (unknown && kronecker(D,p) == -1) return NULL;

s = mpsqrtmod(D,p);

if (!s) err(talker,"not a prime in quadsolvemod");

u = addis(shifti(p,-1), 1); /* = 1/2 */

return modii(mulii(u, subii(s,b)), p);

}

static GEN

otherroot(GEN x, GEN r, GEN p)

{

GEN s = addii((GEN)x[3], r);

s = subii(p, s); if (signe(s) < 0) s = addii(s,p);

return s;

}

/* by splitting */

GEN

rootmod(GEN f, GEN p)

{

long av = avma,tetpil,n,i,j,la,lb;

long n, i, j, la, lb;

gpmem_t av = avma, tetpil;

GEN y,pol,a,b,q,pol0;

if (!factmod_init(&f, p, &i)) { avma=av; return cgetg(1,t_COL); }

i = p[lgefint(p)-1];

i = modBIL(p);

if ((i & 1) == 0) { avma = av; return root_mod_even(f,i); }

i=2; while (!signe(f[i])) i++;

if (i == 2) j = 1;

Line 437 rootmod(GEN f, GEN p)

Line 476 rootmod(GEN f, GEN p)

if (la) y[j+lb] = (long)FpX_normalize(a,p);

pol = gadd(polx[varn(f)], gun); pol0 = constant_term(pol);

while (j<=n)

{

{ /* cf FpX_split_berlekamp */

a=(GEN)y[j]; la=degpol(a);

a = (GEN)y[j]; la = degpol(a);

if (la==1)

y[j++] = lsubii(p, constant_term(a));

else if (la==2)

{

GEN d = subii(sqri((GEN)a[3]), shifti((GEN)a[2],2));

GEN r = quadsolvemod(a, p, 0);

GEN e = mpsqrtmod(d,p), u = addis(q, 1); /* u = 1/2 */

y[j++] = (long)r;

y[j++] = lmodii(mulii(u,subii(e,(GEN)a[3])), p);

y[j++] = (long)otherroot(a,r, p);

y[j++] = lmodii(mulii(u,negi(addii(e,(GEN)a[3]))), p);

}

else for (pol0[2]=1; ; pol0[2]++)

{

b = ZX_s_add(FpXQ_pow(pol,q, a,p), -1); /* pol^(p-1)/2 - 1 */

b = FpX_gcd(a,b, p); lb = degpol(b);

if (lb && lb<la)

if (lb && lb < la)

{

b = FpX_normalize(b, p);

y[j+lb] = (long)FpX_div(a,b, p);

y[j] = (long)b; break;

}

if (pol0[2] == 100 && !BSW_psp(p))

err(talker, "not a prime in polrootsmod");

}

tetpil = avma; y = gerepile(av,tetpil,sort(y));

Line 483 rootmod0(GEN f, GEN p, long flag)

Line 523 rootmod0(GEN f, GEN p, long flag)

/* FACTORISATION MODULO p */

/* */

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

#define FqX_div(x,y,T,p) FpXQX_divres((x),(y),(T),(p),NULL)

#define FqX_rem(x,y,T,p) FpXQX_divres((x),(y),(T),(p),ONLY_REM)

#define FqX_red FpXQX_red

#define FqX_sqr FpXQX_sqr

static GEN spec_FpXQ_pow(GEN x, GEN p, GEN S);

extern GEN pol_to_vec(GEN x, long N);

extern GEN FqXQ_pow(GEN x, GEN n, GEN S, GEN T, GEN p);

extern GEN FpXQX_from_Kronecker(GEN z, GEN pol, GEN p);

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

extern GEN u_normalizepol(GEN x, long lx);

extern GEN Fq_pow(GEN x, GEN n, GEN pol, GEN p);

/* Functions giving information on the factorisation. */

/* u in Z[X], return kernel of (Frob - Id) over Fp[X] / u */

static GEN

Berlekamp_ker(GEN u, GEN p)

FpM_Berlekamp_ker(GEN u, GEN p)

{

long i,j,d,N = degpol(u);

long j,N = degpol(u);

GEN vker,v,w,Q,p1,p2;

GEN vker,v,w,Q,p1;

if (DEBUGLEVEL > 7) timer2();

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

Q = cgetg(N+1,t_MAT); Q[1] = (long)zerocol(N);

w = v = FpXQ_pow(polx[varn(u)],p,u,p);

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

{

Q[j] = lgetg(N+1,t_COL); p1 = (GEN)Q[j];

p1 = pol_to_vec(w, N);

d = lgef(w)-1; p2 = w+1;

for (i=1; i<d ; i++) p1[i] = p2[i];

for ( ; i<=N; i++) p1[i] = zero;

p1[j] = laddis((GEN)p1[j], -1);

Q[j] = (long)p1;

if (j < N)

{

ulong av = avma;

gpmem_t av = avma;

w = gerepileupto(av, FpX_res(gmul(w,v), u, p));

}

Line 514 Berlekamp_ker(GEN u, GEN p)

Line 562 Berlekamp_ker(GEN u, GEN p)

return vker;

}

static GEN

FqM_Berlekamp_ker(GEN u, GEN T, GEN q, GEN p)

{

long j,N = degpol(u);

GEN vker,v,w,Q,p1;

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

Q = cgetg(N+1,t_MAT); Q[1] = (long)zerocol(N);

w = v = FqXQ_pow(polx[varn(u)], q, u, T, p);

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

{

p1 = pol_to_vec(w, N);

p1[j] = laddgs((GEN)p1[j], -1);

Q[j] = (long)p1;

if (j < N)

{

gpmem_t av = avma;

w = gerepileupto(av, FpXQX_divres(FpXQX_mul(w,v, T,p), u,T,p,ONLY_REM));

}

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

vker = FqM_ker(Q,T,p);

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

return vker;

}

/* f in ZZ[X] and p a prime number. */

long

FpX_is_squarefree(GEN f, GEN p)

{

long av = avma;

gpmem_t av = avma;

GEN z;

z = FpX_gcd(f,derivpol(f),p);

avma = av;

Line 533 FpX_is_squarefree(GEN f, GEN p)

Line 606 FpX_is_squarefree(GEN f, GEN p)

long

FpX_nbroots(GEN f, GEN p)

{

long av = avma, n=lgef(f);

long n=lgef(f);

gpmem_t av = avma;

GEN z;

if (n <= 4) return n-3;

f = FpX_red(f, p);

Line 545 FpX_nbroots(GEN f, GEN p)

Line 619 FpX_nbroots(GEN f, GEN p)

long

FpX_is_totally_split(GEN f, GEN p)

{

long av = avma, n=lgef(f);

long n=degpol(f);

gpmem_t av = avma;

GEN z;

if (n <= 4) return 1;

if (n <= 1) return 1;

if (!is_bigint(p) && n-3 > p[2]) return 0;

if (!is_bigint(p) && n > p[2]) return 0;

f = FpX_red(f, p);

z = FpXQ_pow(polx[varn(f)], p, f, p);

avma = av; return lgef(z)==4 && gcmp1((GEN)z[3]) && !signe(z[2]);

Line 559 FpX_is_totally_split(GEN f, GEN p)

Line 634 FpX_is_totally_split(GEN f, GEN p)

long

FpX_nbfact(GEN u, GEN p)

{

ulong av = avma;

gpmem_t av = avma;

GEN vker = Berlekamp_ker(u,p);

GEN vker = FpM_Berlekamp_ker(u,p);

avma = av; return lg(vker)-1;

}

Line 576 static GEN gsmul(GEN a,GEN b){return FpX_mul(a,b,modul

Line 651 static GEN gsmul(GEN a,GEN b){return FpX_mul(a,b,modul

GEN

FpV_roots_to_pol(GEN V, GEN p, long v)

{

ulong ltop=avma;

gpmem_t ltop=avma;

long i;

GEN g=cgetg(lg(V),t_VEC);

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

g[i]=(long)deg1pol(gun,negi((GEN)V[i]),v);

g[i]=(long)deg1pol(gun,modii(negi((GEN)V[i]),p),v);

modulo=p;

g=divide_conquer_prod(g,&gsmul);

return gerepileupto(ltop,g);

Line 595 trivfact(void)

Line 670 trivfact(void)

y[2]=lgetg(1,t_COL); return y;

}

static void

fqunclone(GEN x, GEN a, GEN p)

{

long i,j,lx = lgef(x);

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

{

GEN p1 = (GEN)x[i];

if (typ(p1) == t_POLMOD) { p1[1] = (long)a; p1 = (GEN)p1[2]; }

if (typ(p1) == t_INTMOD) p1[1] = (long)p;

else /* t_POL */

for (j = lgef(p1)-1; j > 1; j--)

{

GEN p2 = (GEN)p1[j];

if (typ(p2) == t_INTMOD) p2[1] = (long)p;

}

static GEN

try_pow(GEN w0, GEN pol, GEN p, GEN q, long r)

{

Line 639 try_pow(GEN w0, GEN pol, GEN p, GEN q, long r)

Line 696 try_pow(GEN w0, GEN pol, GEN p, GEN q, long r)

static void

split(long m, GEN *t, long d, GEN p, GEN q, long r, GEN S)

{

long ps,l,v,dv,av0,av;

long ps, l, v, dv;

gpmem_t av0, av;

GEN w,w0;

dv=degpol(*t); if (dv==d) return;

Line 648 split(long m, GEN *t, long d, GEN p, GEN q, long r, GE

Line 706 split(long m, GEN *t, long d, GEN p, GEN q, long r, GE

{

if (ps==2)

{

w0=w=gpuigs(polx[v],m-1); m+=2;

w0 = w = FpXQ_pow(polx[v], stoi(m-1), *t, gdeux); m += 2;

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

w = gadd(w0, spec_FpXQ_pow(w, p, S));

}

Line 672 split(long m, GEN *t, long d, GEN p, GEN q, long r, GE

Line 730 split(long m, GEN *t, long d, GEN p, GEN q, long r, GE

static void

splitgen(GEN m, GEN *t, long d, GEN p, GEN q, long r)

{

long l,v,dv,av;

long l, v, dv;

gpmem_t av;

GEN w;

dv=degpol(*t); if (dv==d) return;

Line 705 init_pow_p_mod_pT(GEN p, GEN T)

Line 764 init_pow_p_mod_pT(GEN p, GEN T)

S[1] = (long)FpXQ_pow(polx[v], p, T, p);

/* use as many squarings as possible */

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

{

p1 = gsqr((GEN)S[i>>1]);

S[i] = (long)FpX_res(p1, T, p);

if (i == n-1) break;

p1 = gmul((GEN)S[i], (GEN)S[1]);

S[i+1] = (long)FpX_res(p1, T, p);

}

return S;

}

/* compute x^p, S is as above */

static GEN

spec_FpXQ_pow(GEN x, GEN p, GEN S)

{

long av = avma, lim = stack_lim(av,1), i,dx = degpol(x);

long i, dx = degpol(x);

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

GEN x0 = x+2, z;

z = (GEN)x0[0];

if (dx < 0) err(talker, "zero polynomial in FpXQ_pow. %Z not prime", p);

Line 746 spec_FpXQ_pow(GEN x, GEN p, GEN S)

Line 806 spec_FpXQ_pow(GEN x, GEN p, GEN S)

static GEN

factcantor0(GEN f, GEN pp, long flag)

{

long i,j,k,d,e,vf,p,nbfact,tetpil,av = avma;

long i, j, k, d, e, vf, p, nbfact;

gpmem_t tetpil, av = avma;

GEN ex,y,f2,g,g1,u,v,pd,q;

GEN *t;

Line 861 simplefactmod(GEN f, GEN p)

Line 922 simplefactmod(GEN f, GEN p)

return factcantor0(f,p,1);

}

/* vector of polynomials (in v) whose coeffs are given by the columns of x */

GEN

mat_to_vecpol(GEN x, long v)

col_to_ff(GEN x, long v)

{

long i,j, lx = lg(x), lcol = lg(x[1]);

long i, k = lg(x);

GEN y = cgetg(lx, t_VEC);

GEN p;

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

while (--k && gcmp0((GEN)x[k]));

if (k <= 1) return k? (GEN)x[1]: gzero;

i = k+2; p = cgetg(i,t_POL);

p[1] = evalsigne(1) | evallgef(i) | evalvarn(v);

x--; for (k=2; k<i; k++) p[k] = x[k];

return p;

}

GEN

vec_to_pol(GEN x, long v)

{

long i, k = lg(x);

GEN p;

while (--k && gcmp0((GEN)x[k]));

if (!k) return zeropol(v);

i = k+2; p = cgetg(i,t_POL);

p[1] = evalsigne(1) | evallgef(i) | evalvarn(v);

x--; for (k=2; k<i; k++) p[k] = x[k];

return p;

}

GEN

u_vec_to_pol(GEN x)

{

long i, k = lg(x);

GEN p;

while (--k && !x[k]);

if (!k) return u_zeropol();

i = k+2; p = cgetg(i,t_POL);

p[1] = evalsigne(1) | evallgef(i) | evalvarn(0);

x--; for (k=2; k<i; k++) p[k] = x[k];

return p;

}

#if 0

static GEN

u_FpM_FpV_mul(GEN x, GEN y, ulong p)

{

long i,k,l,lx=lg(x), ly=lg(y);

GEN z;

if (lx != ly) err(operi,"* [mod p]",x,y);

if (lx==1) return cgetg(1,t_VECSMALL);

l = lg(x[1]);

z = cgetg(l,t_VECSMALL);

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

{

GEN p1, col = (GEN)x[j];

ulong p1 = 0;

long k = lcol;

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

{

p1 += coeff(x,i,k) * y[k];

if (p1 & HIGHBIT) p1 %= p;

}

z[i] = p1 % p;

}

return z;

}

#endif

while (k-- && gcmp0((GEN)col[k]));

/* u_vec_to_pol(u_FpM_FpV_mul(x, u_pol_to_vec(y), p)) */

i=k+2; p1=cgetg(i,t_POL);

GEN

p1[1] = evalsigne(1) | evallgef(i) | evalvarn(v);

u_FpM_FpX_mul(GEN x, GEN y, ulong p)

col--; for (k=2; k<i; k++) p1[k] = col[k];

{

y[j] = (long)p1;

long i,k,l, ly=lgef(y)-1;

GEN z;

if (ly==1) return u_zeropol();

l = lg(x[1]);

y++;

z = vecsmall_const(l,0) + 1;

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

{

GEN c;

if (!y[k]) continue;

c = (GEN)x[k];

if (y[k] == 1)

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

{

z[i] += c[i];

if (z[i] & HIGHBIT) z[i] %= p;

}

else

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

{

z[i] += c[i] * y[k];

if (z[i] & HIGHBIT) z[i] %= p;

}

for (i=1; i<l; i++) z[i] %= p;

while (--l && !z[l]);

if (!l) return u_zeropol();

*z-- = evalsigne(1) | evallgef(l+2) | evalvarn(0);

return z;

}

/* return the (N-dimensional) vector of coeffs of p */

GEN

pol_to_vec(GEN x, long N)

{

long i, l;

GEN z = cgetg(N+1,t_COL);

if (typ(x) != t_POL)

{

z[1] = (long)x;

for (i=2; i<=N; i++) z[i]=zero;

return z;

}

l = lgef(x)-1; x++;

for (i=1; i<l ; i++) z[i]=x[i];

for ( ; i<=N; i++) z[i]=zero;

return z;

}

/* return the (N-dimensional) vector of coeffs of p */

GEN

u_pol_to_vec(GEN x, long N)

{

long i, l;

GEN z = cgetg(N+1,t_VECSMALL);

if (typ(x) != t_VECSMALL) err(typeer,"u_pol_to_vec");

l = lgef(x)-1; x++;

for (i=1; i<l ; i++) z[i]=x[i];

for ( ; i<=N; i++) z[i]=0;

return z;

}

/* vector of polynomials (in v) whose coeffs are given by the columns of x */

GEN

mat_to_vecpol(GEN x, long v)

{

long j, lx = lg(x);

GEN y = cgetg(lx, t_VEC);

for (j=1; j<lx; j++) y[j] = (long)vec_to_pol((GEN)x[j], v);

return y;

}

/* matrix whose entries are given by the coeffs of the polynomials in

* vector v (considered as degree n-1 polynomials) */

GEN

vecpol_to_mat(GEN v, long n)

{

long i,j,d,N = lg(v);

long j, N = lg(v);

GEN p1,w, y = cgetg(N, t_MAT);

GEN y = cgetg(N, t_MAT);

if (typ(v) != t_VEC) err(typeer,"vecpol_to_mat");

n++;

for (j=1; j<N; j++) y[j] = (long)pol_to_vec((GEN)v[j], n);

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

{

p1 = cgetg(n,t_COL); y[j] = (long)p1;

w = (GEN)v[j];

if (typ(w) != t_POL) { p1[1] = (long)w; i=2; }

else

{

d=lgef(w)-1; w++;

for (i=1; i<d; i++) p1[i] = w[i];

}

for ( ; i<n; i++) p1[i] = zero;

}

return y;

}

/* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */

GEN

mat_to_polpol(GEN x, long v,long w)

{

long i,j, lx = lg(x), lcol = lg(x[1]);

long j, lx = lg(x);

GEN y = cgetg(lx+1, t_POL);

y[1]=evalsigne(1) | evallgef(lx+1) | evalvarn(v);

y++;

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

for (j=1; j<lx; j++) y[j] = (long)vec_to_pol((GEN)x[j], w);

{

return normalizepol_i(--y, lx+1);

GEN p1, col = (GEN)x[j];

long k;

i=lcol+1; p1=cgetg(i,t_POL);

p1[1] = evalsigne(1) | evallgef(i) | evalvarn(w);

col--; for (k=2; k<i; k++) p1[k] = col[k];

y[j] = (long)normalizepol_i(p1,i);

}

return normalizepol_i(--y,lx+1);

}

/* matrix whose entries are given by the coeffs of the polynomial v in

Line 930 mat_to_polpol(GEN x, long v,long w)

Line 1093 mat_to_polpol(GEN x, long v,long w)

GEN

polpol_to_mat(GEN v, long n)

{

long i,j,d,N = lgef(v)-1;

long j, N = lgef(v)-1;

GEN p1,w, y = cgetg(N, t_MAT);

GEN y = cgetg(N, t_MAT);

if (typ(v) != t_POL) err(typeer,"polpol_to_mat");

n++;v++;

v++;

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

for (j=1; j<N; j++) y[j] = (long)pol_to_vec((GEN)v[j], n);

return y;

}

/* P(X,Y) --> P(Y,X), n-1 is the degree in Y */

GEN

swap_polpol(GEN x, long n, long w)

{

long j, lx = lgef(x), ly = n+3;

long v=varn(x);

GEN y = cgetg(ly, t_POL);

y[1]=evalsigne(1) | evallgef(ly) | evalvarn(v);

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

{

p1 = cgetg(n,t_COL); y[j] = (long)p1;

long k;

w = (GEN)v[j];

GEN p1=cgetg(lx,t_POL);

if (typ(w) != t_POL) { p1[1] = (long)w; i=2; }

p1[1] = evalsigne(1) | evallgef(lx) | evalvarn(w);

else

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

{

if( j<lgef(x[k]))

d=lgef(w)-1; w++;

p1[k] = mael(x,k,j);

for (i=1; i<d; i++) p1[i] = w[i];

else

}

p1[k] = zero;

for ( ; i<n; i++) p1[i] = zero;

y[j] = (long)normalizepol_i(p1,lx);

}

return y;

return normalizepol_i(y,ly);

}

/* set x <-- x + c*y mod p */

static void

split_berlekamp_addmul(GEN x, GEN y, long c, long p)

u_FpX_addmul(GEN x, GEN y, long c, long p)

{

long i,lx,ly,l;

if (!c) return;

lx = lgef(x); ly = lgef(y); l = min(lx,ly);

lx = lgef(x);

ly = lgef(y); l = min(lx,ly);

if (p & ~MAXHALFULONG)

{

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

for ( ; i<ly; i++) x[i] = mulssmod(c,y[i],p);

if (l == lx)

for ( ; i<ly; i++) x[i] = mulssmod(c,y[i],p);

}

else

{

for (i=2; i<l; i++) x[i] = ((ulong)x[i] + (ulong)(c*y[i])) % p;

for ( ; i<ly; i++) x[i] = (c*y[i]) % p;

if (l == lx)

for ( ; i<ly; i++) x[i] = (c*y[i]) % p;

}

do i--; while (i>1 && !x[i]);

(void)u_normalizepol(x,i);

if (i==1) setsigne(x,0); else { setsigne(x,1); setlgef(x,i+1); }

}

static long

small_rand(long p)

{

return (p==2)? ((mymyrand() & 0x1000) == 0): mymyrand() % p;

}

GEN

FpX_rand(long d1, long v, GEN p)

{

long i, d = d1+2;

GEN y;

y = cgetg(d,t_POL); y[1] = evalsigne(1) | evalvarn(v);

for (i=2; i<d; i++) y[i] = (long)genrand(p);

(void)normalizepol_i(y,d); return y;

}

/* return a random polynomial in F_q[v], degree < d1 */

GEN

FqX_rand(long d1, long v, GEN T, GEN p)

{

long i, d = d1+2, k = degpol(T), w = varn(T);

GEN y;

y = cgetg(d,t_POL); y[1] = evalsigne(1) | evalvarn(v);

for (i=2; i<d; i++) y[i] = (long)FpX_rand(k, w, p);

(void)normalizepol_i(y,d); return y;

}

#define set_irred(i) { if ((i)>ir) swap(t[i],t[ir]); ir++;}

long

split_berlekamp(GEN *t, GEN pp, GEN pps2)

FpX_split_berlekamp(GEN *t, GEN p)

{

GEN u = *t, p1, p2, vker,pol;

GEN u = *t, a,b,po2,vker,pol;

long av,N = degpol(u), d,i,kk,l1,l2,p, vu = varn(u);

long N = degpol(u), d, i, ir, L, la, lb, ps, vu = varn(u);

ulong av0 = avma;

gpmem_t av0 = avma;

vker = Berlekamp_ker(u,pp);

vker = FpM_Berlekamp_ker(u,p);

vker = mat_to_vecpol(vker,vu);

d = lg(vker)-1;

p = is_bigint(pp)? 0: pp[2];

ps = is_bigint(p)? 0: p[2];

if (p)

if (ps)

{

avma = av0; p1 = cgetg(d+1, t_VEC); /* hack: hidden gerepile */

avma = av0; a = cgetg(d+1, t_VEC); /* hack: hidden gerepile */

for (i=1; i<=d; i++) p1[i] = (long)pol_to_small((GEN)vker[i]);

for (i=1; i<=d; i++) a[i] = (long)pol_to_small((GEN)vker[i]);

vker = p1;

vker = a;

}

po2 = shifti(p, -1); /* (p-1) / 2 */

pol = cgetg(N+3,t_POL);

ir = 0;

for (kk=1; kk<d; )

/* t[i] irreducible for i <= ir, still to be treated for i < L */

for (L=1; L<d; )

{

GEN polt;

if (p)

if (ps)

{

if (p==2)

pol[2] = small_rand(ps); /* vker[1] = 1 */

{

pol[1] = evallgef(pol[2]? 3: 2);

pol[2] = ((mymyrand() & 0x1000) == 0);

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

pol[1] = evallgef(pol[2]? 3: 2);

u_FpX_addmul(pol,(GEN)vker[i],small_rand(ps), ps);

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

split_berlekamp_addmul(pol,(GEN)vker[i],(mymyrand()&0x1000)?0:1, p);

}

else

{

pol[2] = mymyrand()%p; /* vker[1] = 1 */

pol[1] = evallgef(pol[2]? 3: 2);

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

split_berlekamp_addmul(pol,(GEN)vker[i],mymyrand()%p, p);

}

polt = small_to_pol(pol,vu);

}

else

{

pol[2] = (long)genrand(pp);

pol[2] = (long)genrand(p);

pol[1] = evallgef(signe(pol[2])? 3: 2) | evalvarn(vu);

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

pol = gadd(pol, gmul((GEN)vker[i], genrand(pp)));

pol = gadd(pol, gmul((GEN)vker[i], genrand(p)));

polt = FpX_red(pol,pp);

polt = FpX_red(pol,p);

}

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

for (i=ir; i<L && L<d; i++)

{

p1=t[i-1]; l1=degpol(p1);

a = t[i]; la = degpol(a);

if (l1>1)

if (la == 1) { set_irred(i); }

else if (la == 2)

{

av = avma; p2 = FpX_res(polt, p1, pp);

GEN r = quadsolvemod(a,p,1);

if (lgef(p2) <= 3) { avma=av; continue; }

if (r)

p2 = FpXQ_pow(p2,pps2, p1,pp);

if (!signe(p2)) err(talker,"%Z not a prime in split_berlekamp",pp);

p2 = ZX_s_add(p2, -1);

p2 = FpX_gcd(p1,p2, pp); l2=degpol(p2);

if (l2>0 && l2<l1)

{

p2 = FpX_normalize(p2, pp);

t[i] = deg1pol_i(gun, subii(p,r), vu); r = otherroot(a,r,p);

t[i-1] = p2; kk++;

t[L] = deg1pol_i(gun, subii(p,r), vu); L++;

t[kk-1] = FpX_div(p1,p2,pp);

if (DEBUGLEVEL > 7) msgtimer("new factor");

}

set_irred(i);

}

else

{

gpmem_t av = avma;

b = FpX_res(polt, a, p);

if (degpol(b) == 0) { avma=av; continue; }

b = ZX_s_add(FpXQ_pow(b,po2, a,p), -1);

b = FpX_gcd(a,b, p); lb = degpol(b);

if (lb && lb < la)

{

b = FpX_normalize(b, p);

t[L] = FpX_div(a,b,p);

t[i]= b; L++;

}

else avma = av;

}

Line 1044 split_berlekamp(GEN *t, GEN pp, GEN pps2)

Line 1250 split_berlekamp(GEN *t, GEN pp, GEN pps2)

return d;

}

static GEN

FqX_deriv(GEN f, GEN T, GEN p)

{

return FqX_red(derivpol(f), T, p);

}

GEN

FqX_gcd(GEN P, GEN Q, GEN T, GEN p)

{

GEN g = T? FpXQX_safegcd(P,Q,T,p): FpX_gcd(P,Q,p);

if (!g) err(talker,"factmod9: %Z is reducible mod p!", T);

return g;

}

long

FqX_is_squarefree(GEN P, GEN T, GEN p)

{

gpmem_t av = avma;

GEN z = FqX_gcd(P, derivpol(P), T, p);

avma = av;

return degpol(z)==0;

}

long

FqX_split_berlekamp(GEN *t, GEN q, GEN T, GEN p)

{

GEN u = *t, a,b,qo2,vker,pol;

long N = degpol(u), vu = varn(u), vT = varn(T), dT = degpol(T);

long d, i, ir, L, la, lb;

vker = FqM_Berlekamp_ker(u,T,q,p);

vker = mat_to_vecpol(vker,vu);

d = lg(vker)-1;

qo2 = shifti(q, -1); /* (q-1) / 2 */

pol = cgetg(N+3,t_POL);

ir = 0;

/* t[i] irreducible for i < ir, still to be treated for i < L */

for (L=1; L<d; )

{

GEN polt;

pol[2] = (long)FpX_rand(dT,vT,p);

pol[1] = evallgef(signe(pol[2])? 3: 2) | evalvarn(vu);

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

pol = gadd(pol, gmul((GEN)vker[i], FpX_rand(dT,vT,p)));

polt = FqX_red(pol,T,p);

for (i=ir; i<L && L<d; i++)

{

a = t[i]; la = degpol(a);

if (la == 1) { set_irred(i); }

else

{

gpmem_t av = avma;

b = FqX_rem(polt, a, T,p);

if (!degpol(b)) { avma=av; continue; }

b = FqXQ_pow(b,qo2, a,T,p);

if (!degpol(b)) { avma=av; continue; }

b[2] = ladd((GEN)b[2], gun);

b = FqX_gcd(a,b, T,p); lb = degpol(b);

if (lb && lb < la)

{

b = FpXQX_normalize(b, T,p);

t[L] = FqX_div(a,b,T,p);

t[i]= b; L++;

}

else avma = av;

}

return d;

}

GEN

factmod0(GEN f, GEN pp)

{

long i,j,k,e,p,N,nbfact,av = avma,tetpil,d;

long i, j, k, e, p, N, nbfact, d;

gpmem_t av = avma, tetpil;

GEN pps2,ex,y,f2,p1,g1,u, *t;

if (!(d = factmod_init(&f, pp, &p))) { avma=av; return trivfact(); }

Line 1065 factmod0(GEN f, GEN pp)

Line 1342 factmod0(GEN f, GEN pp)

{

k++; if (p && !(k%p)) { k++; f2 = FpX_div(f2,g1,pp); }

p1 = FpX_gcd(f2,g1, pp); u = g1; g1 = p1;

if (lgef(p1)!=3)

if (degpol(p1))

{

u = FpX_div( u,p1,pp);

f2= FpX_div(f2,p1,pp);

}

/* u is square-free (product of irred. of multiplicity e * k) */

N = degpol(u);

if (N)

{

/* here u is square-free (product of irred. of multiplicity e * k) */

t[nbfact] = FpX_normalize(u,pp);

d = (N==1)? 1: split_berlekamp(t+nbfact, pp, pps2);

d = (N==1)? 1: FpX_split_berlekamp(t+nbfact, pp);

for (j=0; j<d; j++) ex[nbfact+j] = e*k;

nbfact += d;

}

Line 1097 factmod0(GEN f, GEN pp)

Line 1374 factmod0(GEN f, GEN pp)

GEN

factmod(GEN f, GEN pp)

{

long tetpil,av=avma;

gpmem_t tetpil, av=avma;

long nbfact;

long j;

GEN y,u,v;

Line 1135 static GEN

Line 1412 static GEN

padic_pol_to_int(GEN f)

{

long i, l = lgef(f);

f = gdiv(f,content(f));

GEN c = content(f);

if (gcmp0(c)) /* O(p^n) can occur */

{

if (typ(c) != t_PADIC) err(typeer,"padic_pol_to_int");

f = gdiv(f, gpowgs((GEN)c[2], valp(c)));

}

else

f = gdiv(f,c);

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

switch(typ(f[i]))

{

Line 1146 padic_pol_to_int(GEN f)

Line 1430 padic_pol_to_int(GEN f)

return f;

}

/* return invlead * (x + O(pr)), x in Z or Z_p, pr = p^r */

/* return invlead * (x + O(pr)), x in Z or Z_p, pr = p^r. May return gzero */

static GEN

int_to_padic(GEN x, GEN p, GEN pr, long r, GEN invlead)

{

GEN p1,y;

long v,sx, av = avma;

long v, sx;

gpmem_t av = avma;

if (typ(x) == t_PADIC)

{

Line 1171 int_to_padic(GEN x, GEN p, GEN pr, long r, GEN invlead

Line 1456 int_to_padic(GEN x, GEN p, GEN pr, long r, GEN invlead

{

y[4] = lmodii(p1,pr); r -= v;

}

else

{

y[4] = zero; v = r; r = 0;

}

y[3] = (long)pr;

y[2] = (long)p;

y[1] = evalprecp(r)|evalvalp(v);

return invlead? gerepileupto(av, gmul(invlead,y)): y;

}

/* as above. always return a t_PADIC */

static GEN

strict_int_to_padic(GEN x, GEN p, GEN pr, long r, GEN invlead)

{

GEN y = int_to_padic(x, p, pr, r, invlead);

if (y == gzero) y = ggrandocp(p,r);

return y;

}

/* return (x + O(p^r)) normalized (multiply by a unit such that leading coeff

* is a power of p), x in Z[X] (or Z_p[X]) */

static GEN

Line 1192 pol_to_padic(GEN x, GEN pr, GEN p, long r)

Line 1486 pol_to_padic(GEN x, GEN pr, GEN p, long r)

if (gcmp1(lead)) lead = NULL;

else

{

long av = avma;

gpmem_t av = avma;

v = ggval(lead,p);

if (v) lead = gdiv(lead, gpowgs(p,v));

lead = int_to_padic(lead,p,pr,r,NULL);

Line 1206 pol_to_padic(GEN x, GEN pr, GEN p, long r)

Line 1500 pol_to_padic(GEN x, GEN pr, GEN p, long r)

static GEN

padic_trivfact(GEN x, GEN p, long r)

{

GEN p1, y = cgetg(3,t_MAT);

GEN y = cgetg(3,t_MAT);

p1=cgetg(2,t_COL); y[1]=(long)p1; p = icopy(p);

p = icopy(p);

p1[1]=(long)pol_to_padic(x,gpowgs(p,r),p,r);

y[1]=(long)_col(pol_to_padic(x,gpowgs(p,r),p,r));

p1=cgetg(2,t_COL); y[2]=(long)p1;

y[2]=(long)_col(gun); return y;

p1[1]=un; return y;

}

/* a etant un p-adique, retourne le vecteur des racines p-adiques de f

/* Assume x reduced mod p. q = p^e (simplified version of int_to_padic).

* congrues a a modulo p dans le cas ou on suppose f(a) congru a 0 modulo p

* If p = 2, is defined (and reduced) mod 4 [from rootmod] */

* (ou a 4 si p=2).

GEN

apprgen(GEN f, GEN a)

Fp_to_Zp(GEN x, GEN p, GEN q, long e)

{

GEN fp,p1,p,P,pro,x,x2,u,ip;

GEN y = cgetg(5, t_PADIC);

long av=avma,tetpil,v,Ps,i,j,k,lu,n,fl2;

if (egalii(p, x)) /* implies p = x = 2 */

{

x = gun; q = shifti(q, -1);

y[1] = evalprecp(e-1) | evalvalp(1);

}

else if (!signe(x)) y[1] = evalprecp(0) | evalvalp(e);

else y[1] = evalprecp(e) | evalvalp(0);

y[2] = (long)p;

y[3] = (long)q;

y[4] = (long)x; return y;

}

if (typ(f)!=t_POL) err(notpoler,"apprgen");

/* f != 0 in Z[X], primitive, a t_PADIC s.t f(a) = 0 mod p */

if (gcmp0(f)) err(zeropoler,"apprgen");

static GEN

if (typ(a) != t_PADIC) err(rootper1);

apprgen_i(GEN f, GEN a)

{

GEN fp,u,p,q,P,res,a0,rac;

long prec,i,j,k;

prec = gcmp0(a)? valp(a): precp(a);

if (prec <= 1) return _vec(a);

fp = derivpol(f); u = ggcd(f,fp);

if (degpol(u) > 0) { f = gdeuc(f,u); fp = derivpol(f); }

p = (GEN)a[2];

P = egalii(p,gdeux)? stoi(4): p;

a0= gmod(a, P);

#if 0 /* assumption */

if (!gcmp0(FpX_eval(f,a0,P))) err(rootper2);

#endif

if (!gcmp0(FpX_eval(fp,a0,p))) /* simple zero */

{

res = rootpadiclift(f, a0, p, prec);

res = strict_int_to_padic(res,p,gpowgs(p,prec),prec,NULL);

return _vec(res);

}

f = poleval(f, gadd(a, gmul(P,polx[varn(f)])));

f = padic_pol_to_int(f);

fp=derivpol(f); p1=ggcd(f,fp);

f = gdiv(f, gpowgs(p,ggval(f, p)));

if (lgef(p1)>3) { f=gdeuc(f,p1); fp=derivpol(f); }

p=(GEN)a[2]; p1=poleval(f,a);

res = cgetg(degpol(f)+1,t_VEC);

v=ggval(p1,p); if (v <= 0) err(rootper2);

q = gpowgs(p,prec);

fl2=egalii(p,gdeux);

rac = lift_intern(rootmod(f, P));

if (fl2 && v==1) err(rootper2);

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

v=ggval(poleval(fp,a),p);

if (!v) /* simple zero */

{

while (!gcmp0(p1))

u = apprgen_i(f, Fp_to_Zp((GEN)rac[i], p,q,prec));

{

for (k=1; k<lg(u); k++) res[j++] = ladd(a, gmul(P,(GEN)u[k]));

a = gsub(a,gdiv(p1,poleval(fp,a)));

p1 = poleval(f,a);

}

tetpil=avma; pro=cgetg(2,t_VEC); pro[1]=lcopy(a);

return gerepile(av,tetpil,pro);

}

n=degpol(f); pro=cgetg(n+1,t_VEC);

setlg(res,j); return res;

}

if (is_bigint(p)) err(impl,"apprgen for p>=2^31");

/* a t_PADIC, return vector of p-adic roots of f equal to a (mod p)

x = ggrandocp(p, valp(a) | precp(a));

* We assume 1) f(a) = 0 mod p (mod 4 if p=2).

if (fl2)

* 2) leading coeff prime to p. */

GEN

apprgen(GEN f, GEN a)

{

gpmem_t av = avma;

if (typ(f) != t_POL) err(notpoler,"apprgen");

if (gcmp0(f)) err(zeropoler,"apprgen");

if (typ(a) != t_PADIC) err(rootper1);

return gerepilecopy(av, apprgen_i(padic_pol_to_int(f), a));

}

static GEN

rootpadic_i(GEN f, GEN p, long prec)

{

GEN fp,y,z,q,rac;

long lx,i,j,k;

fp = derivpol(f); z = ggcd(f,fp);

if (degpol(z) > 0) { f = gdeuc(f,z); fp = derivpol(f); }

rac = rootmod(f, (egalii(p,gdeux) && prec>=2)? stoi(4): p);

rac = lift_intern(rac); lx = lg(rac);

if (prec==1)

{

x2=ggrandocp(p,2); P = stoi(4);

y = cgetg(lx,t_COL);

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

y[i] = (long)Fp_to_Zp((GEN)rac[i],p,p,1);

return y;

}

else

y = cgetg(degpol(f)+1,t_COL);

q = gpowgs(p,prec);

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

{

x2=ggrandocp(p,1); P = p;

z = apprgen_i(f, Fp_to_Zp((GEN)rac[i], p,q,prec));

for (k=1; k<lg(z); k++,j++) y[j] = z[k];

}

f = poleval(f, gadd(a,gmul(P,polx[varn(f)])));

setlg(y,j); return y;

if (!gcmp0(f)) f = gdiv(f,gpuigs(p,ggval(f, p)));

}

Ps = itos(P);

for (j=0,i=0; i<Ps; i++)

/* reverse x in place */

static void

polreverse(GEN x)

{

long i, j;

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

for (i=2, j=lgef(x)-1; i<j; i++, j--) lswap(x[i], x[j]);

(void)normalizepol(x);

}

static GEN

pnormalize(GEN f, GEN p, long prec, long n, GEN *plead, long *pprec, int *prev)

{

*plead = leading_term(f);

*pprec = prec;

*prev = 0;

if (!gcmp1(*plead))

{

ip=stoi(i);

long val = ggval(*plead,p), val1 = ggval(constant_term(f),p);

if (gcmp0(poleval(f,gadd(ip,x2))))

if (val1 < val)

{

u=apprgen(f,gadd(x,ip)); lu=lg(u);

*prev = 1; polreverse(f);

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

/* take care of loss of precision from leading coeff of factor

{

* (whose valuation is <= val) */

j++; pro[j]=ladd(a,gmul(P,(GEN)u[k]));

*pprec += val;

}

val = val1;

}

*pprec += val * n;

}

setlg(pro,j+1); return gerepilecopy(av,pro);

return pol_to_monic(f, plead);

}

/* Retourne le vecteur des racines p-adiques de f en precision r */

/* return p-adic roots of f, precision prec */

GEN

rootpadic(GEN f, GEN p, long r)

rootpadic(GEN f, GEN p, long prec)

{

GEN fp,y,z,p1,pr,rac;

gpmem_t av = avma;

long lx,i,j,k,n,av=avma,tetpil,fl2;

GEN lead,y;

long PREC,i,k;

int reverse;

if (typ(f)!=t_POL) err(notpoler,"rootpadic");

if (gcmp0(f)) err(zeropoler,"rootpadic");

if (r<=0) err(rootper4);

if (prec <= 0) err(rootper4);

f = padic_pol_to_int(f);

fp=derivpol(f); p1=ggcd(f,fp);

f = pnormalize(f, p, prec, 1, &lead, &PREC, &reverse);

if (lgef(p1)>3) { f=gdeuc(f,p1); fp=derivpol(f); }

fl2=egalii(p,gdeux); rac=(fl2 && r>=2)? rootmod(f,stoi(4)): rootmod(f,p);

y = rootpadic_i(f,p,PREC);

lx=lg(rac); p=gclone(p);

k = lg(y);

if (r==1)

if (lead)

{

for (i=1; i<k; i++) y[i] = ldiv((GEN)y[i], lead);

tetpil=avma; y=cgetg(lx,t_COL);

if (reverse)

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

for (i=1; i<k; i++) y[i] = linv((GEN)y[i]);

{

return gerepilecopy(av, y);

z=cgetg(5,t_PADIC); y[i]=(long)z;

z[1] = evalprecp(1)|evalvalp(0);

z[2] = z[3] = (long)p;

z[4] = lcopy(gmael(rac,i,2));

}

return gerepile(av,tetpil,y);

}

n=degpol(f); y=cgetg(n+1,t_COL);

j=0; pr = NULL;

z = cgetg(5,t_PADIC);

z[2] = (long)p;

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

{

p1 = gmael(rac,i,2);

if (signe(p1))

{

if (!fl2 || mod2(p1))

{

z[1] = evalvalp(0)|evalprecp(r);

z[4] = (long)p1;

}

else

{

z[1] = evalvalp(1)|evalprecp(r);

z[4] = un;

}

if (!pr) pr=gpuigs(p,r);

z[3] = (long)pr;

}

else

{

z[1] = evalvalp(r);

z[3] = un;

z[4] = (long)p1;

}

p1 = apprgen(f,z);

for (k=1; k<lg(p1); k++) y[++j]=p1[k];

}

setlg(y,j+1); return gerepilecopy(av,y);

}

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

/* rootpadicfast */

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

/*lift accelerator. The author of the idea is unknown.*/

/* lift accelerator */

long hensel_lift_accel(long n, long *pmask)

{

long a,j;

Line 1365 long hensel_lift_accel(long n, long *pmask)

Line 1689 long hensel_lift_accel(long n, long *pmask)

/* STANDARD USE

There exists p a prime number and a>0 such that

q=p^a

f in ZZ[X], with leading term prime to p.

S must be a simple root mod p.

Line 1375 long hensel_lift_accel(long n, long *pmask)

Line 1699 long hensel_lift_accel(long n, long *pmask)

GEN

rootpadiclift(GEN T, GEN S, GEN p, long e)

{

ulong ltop=avma;

gpmem_t ltop=avma;

long x;

GEN qold, q, qm1;

GEN W, Tr, Sr, Wr = gzero;

Line 1387 rootpadiclift(GEN T, GEN S, GEN p, long e)

Line 1711 rootpadiclift(GEN T, GEN S, GEN p, long e)

W=FpX_eval(deriv(Tr, x),S,q);

W=mpinvmod(W,q);

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

{

qm1 = (mask&(1<<i))?sqri(qm1):mulii(qm1, q);

q = mulii(qm1, p);

Tr = FpX_red(T,q);

Line 1424 rootpadicliftroots(GEN f, GEN S, GEN q, long e)

Line 1748 rootpadicliftroots(GEN f, GEN S, GEN q, long e)

y[n-1]=(long) rootpadiclift(f, (GEN) S[n-1], q, e);

else/* distinct-->totally split-->use trace trick */

{

ulong av=avma;

gpmem_t av=avma;

GEN z;

z=(GEN)f[lgef(f)-2];/*-trace(roots)*/

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

Line 1441 rootpadicliftroots(GEN f, GEN S, GEN q, long e)

Line 1765 rootpadicliftroots(GEN f, GEN S, GEN q, long e)

f must have no multiple roots mod p.

return p-adics roots of f with prec pr, as integers (implicitly mod pr)

GEN

rootpadicfast(GEN f, GEN p, long e)

{

ulong ltop=avma;

gpmem_t ltop=avma;

GEN S,y;

S=lift(rootmod(f,p));/*no multiplicity*/

if (lg(S)==1)/*no roots*/

Line 1461 rootpadicfast(GEN f, GEN p, long e)

Line 1785 rootpadicfast(GEN f, GEN p, long e)

return y;

}

/* Same as rootpadiclift for the polynomial X^n-a,

* but here, n can be much larger.

* TODO: generalize to sparse polynomials.

GEN

padicsqrtnlift(GEN a, GEN n, GEN S, GEN p, long e)

{

ulong ltop=avma;

gpmem_t ltop=avma;

GEN qold, q, qm1;

GEN W, Sr, Wr = gzero;

long i, nb, mask;

Line 1476 padicsqrtnlift(GEN a, GEN n, GEN S, GEN p, long e)

Line 1800 padicsqrtnlift(GEN a, GEN n, GEN S, GEN p, long e)

W = modii(mulii(n,powmodulo(S,subii(n,gun),q)),q);

W = mpinvmod(W,q);

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

{

qm1 = (mask&(1<<i))?sqri(qm1):mulii(qm1, q);

q = mulii(qm1, p);

Sr = S;

Line 1512 getprec(GEN x, long prec, GEN *p)

Line 1836 getprec(GEN x, long prec, GEN *p)

return prec;

}

/* a appartenant a une extension finie de Q_p, retourne le vecteur des

/* a belongs to finite extension of Q_p, return all roots of f equal to a

* racines de f congrues a a modulo p dans le cas ou on suppose f(a) congru a

* mod p. We assume f(a) = 0 (mod p) [mod 4 if p=2] */

* 0 modulo p (ou a 4 si p=2).

GEN

apprgen9(GEN f, GEN a)

{

GEN fp,p1,p,pro,x,x2,u,ip,t,vecg;

GEN fp,p1,p,pro,x,x2,u,ip,T,vecg;

long av=avma,tetpil,v,ps_1,i,j,k,lu,n,prec,d,va,fl2;

long v, ps_1, i, j, k, n, prec, d, va, fl2;

gpmem_t av=avma, tetpil;

if (typ(f)!=t_POL) err(notpoler,"apprgen9");

if (gcmp0(f)) err(zeropoler,"apprgen9");

if (typ(a)==t_PADIC) return apprgen(f,a);

if (typ(a)!=t_POLMOD || typ(a[2])!=t_POL) err(rootper1);

fp=derivpol(f); p1=ggcd(f,fp);

if (lgef(p1)>3) { f=gdeuc(f,p1); fp=derivpol(f); }

fp = derivpol(f); u = ggcd(f,fp);

t=(GEN)a[1];

if (degpol(u) > 0) { f = gdeuc(f,u); fp = derivpol(f); }

T = (GEN)a[1];

prec = getprec((GEN)a[2], BIGINT, &p);

prec = getprec(t, prec, &p);

prec = getprec(T, prec, &p);

if (prec==BIGINT) err(rootper1);

p1=poleval(f,a); v=ggval(lift_intern(p1),p); if (v<=0) err(rootper2);

p1 = poleval(f,a); v = ggval(lift_intern(p1),p); if (v<=0) err(rootper2);

fl2=egalii(p,gdeux);

fl2 = egalii(p,gdeux);

if (fl2 && v==1) err(rootper2);

v=ggval(lift_intern(poleval(fp,a)), p);

if (!v)

Line 1547 apprgen9(GEN f, GEN a)

Line 1871 apprgen9(GEN f, GEN a)

tetpil=avma; pro=cgetg(2,t_COL); pro[1]=lcopy(a);

return gerepile(av,tetpil,pro);

}

n=degpol(f); pro=cgetg(n+1,t_COL); j=0;

n=degpol(f); pro=cgetg(n+1,t_COL);

if (is_bigint(p)) err(impl,"apprgen9 for p>=2^31");

x=gmodulcp(ggrandocp(p,prec), t);

x=gmodulcp(ggrandocp(p,prec), T);

if (fl2)

{

ps_1=3; x2=ggrandocp(p,2); p=stoi(4);

Line 1560 apprgen9(GEN f, GEN a)

Line 1884 apprgen9(GEN f, GEN a)

ps_1=itos(p)-1; x2=ggrandocp(p,1);

}

f = poleval(f,gadd(a,gmul(p,polx[varn(f)])));

if (!gcmp0(f)) f=gdiv(f,gpuigs(p,ggval(f,p)));

f = gdiv(f,gpowgs(p,ggval(f,p)));

d=degpol(t); vecg=cgetg(d+1,t_COL);

d=degpol(T); vecg=cgetg(d+1,t_COL);

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

vecg[i] = (long)setloop(gzero);

va=varn(t);

va=varn(T); j = 1;

for(;;) /* loop through F_q */

{

ip=gmodulcp(gtopoly(vecg,va),t);

ip=gmodulcp(gtopoly(vecg,va),T);

if (gcmp0(poleval(f,gadd(ip,x2))))

{

u=apprgen9(f,gadd(ip,x)); lu=lg(u);

u = apprgen9(f,gadd(ip,x));

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

for (k=1; k<lg(u); k++) pro[j++] = ladd(a, gmul(p,(GEN)u[k]));

{

j++; pro[j]=ladd(a,gmul(p,(GEN)u[k]));

}

for (i=d; i; i--)

{

Line 1584 apprgen9(GEN f, GEN a)

Line 1906 apprgen9(GEN f, GEN a)

}

if (!i) break;

}

setlg(pro,j+1); return gerepilecopy(av,pro);

setlg(pro,j); return gerepilecopy(av,pro);

}

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

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

/* Factorisation p-adique d'un polynome */

/* */

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

/* FACTORIZATION in Zp[X] */

/* */

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

extern GEN ZX_squff(GEN f, GEN *ex);

extern GEN fact_from_DDF(GEN fa, GEN e, long n);

int

cmp_padic(GEN x, GEN y)

{

Line 1603 cmp_padic(GEN x, GEN y)

Line 1930 cmp_padic(GEN x, GEN y)

return cmpii((GEN)x[4], (GEN)y[4]);

}

/* factorise le polynome T=nf[1] dans Zp avec la precision pr */

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

/* Using Buchmann--Lenstra */

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

/* factor T = nf[1] in Zp to precision k */

static GEN

padicff2(GEN nf,GEN p,long pr)

padicff2(GEN nf,GEN p,long k)

{

long N=degpol(nf[1]),i,j,d,l;

GEN z, mat, D, U,fa, pk, dec_p, Ui, mulx;

GEN mat,V,D,fa,p1,pk,dec_p,pke,a,theta;

long i, l;

pk=gpuigs(p,pr); dec_p=primedec(nf,p);

mulx = eltmul_get_table(nf, gmael(nf,8,2)); /* mul by (x mod T) */

l=lg(dec_p); fa=cgetg(l,t_COL);

dec_p = primedec(nf,p);

l = lg(dec_p); fa = cgetg(l,t_COL);

D = NULL; /* -Wall */

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

{

p1 = (GEN)dec_p[i];

GEN P = (GEN)dec_p[i];

pke = idealpows(nf,p1, pr * itos((GEN)p1[3]));

long e = itos((GEN)P[3]), ef = e * itos((GEN)P[4]);

p1=smith2(pke); V=(GEN)p1[3]; D=(GEN)p1[1];

D = smithall(idealpows(nf,P, k*e), &U, NULL);

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

Ui= ginv(U); setlg(Ui, ef+1); /* cf smithrel */

if (! egalii(gcoeff(V,d,d),pk)) break;

U = rowextract_i(U, 1, ef);

a=ginv(D); theta=gmael(nf,8,2); mat=cgetg(d,t_MAT);

mat = gmul(U, gmul(mulx, Ui));

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

fa[i] = (long)caradj(mat,0,NULL);

{

p1 = gmul(D, element_mul(nf,theta,(GEN)a[j]));

setlg(p1,d); mat[j]=(long)p1;

}

fa[i]=(long)caradj(mat,0,NULL);

}

a = cgetg(l,t_COL); pk = icopy(pk);

pk = gcoeff(D,1,1); /* = p^k */

z = cgetg(l,t_COL); pk = icopy(pk);

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

a[i] = (long)pol_to_padic((GEN)fa[i],pk,p,pr);

z[i] = (long)pol_to_padic((GEN)fa[i],pk,p,k);

return a;

return z;

}

static GEN

padicff(GEN x,GEN p,long pr)

{

GEN q,basden,bas,invbas,mul,dx,nf,mat;

GEN q,basden,bas,invbas,mul,dx,dK,nf,fa,g,e;

long n=degpol(x),av=avma;

long n=degpol(x);

gpmem_t av=avma;

nf=cgetg(10,t_VEC); nf[1]=(long)x; dx=discsr(x);

nf=cgetg(10,t_VEC); nf[1]=(long)x;

mat=cgetg(3,t_MAT); mat[1]=lgetg(3,t_COL); mat[2]=lgetg(3,t_COL);

fa = cgetg(3,t_MAT);

coeff(mat,1,1)=(long)p; coeff(mat,1,2)=lstoi(pvaluation(dx,p,&q));

g = cgetg(3,t_COL); fa[1] = (long)g;

coeff(mat,2,1)=(long)q; coeff(mat,2,2)=un;

e = cgetg(3,t_COL); fa[2] = (long)e;

bas=allbase4(x,(long)mat,(GEN*)(nf+3),NULL);

dx = ZX_disc(x);

if (!carrecomplet(divii(dx,(GEN)nf[3]),(GEN*)(nf+4)))

g[1] = (long)p; e[1] = lstoi(pvaluation(absi(dx),p,&q));

err(bugparier,"factorpadic2 (incorrect discriminant)");

g[2] = (long)q; e[2] = un;

bas = nfbasis(x, &dK, 0, fa);

nf[3] = (long)dK;

nf[4] = divise( diviiexact(dx, dK), p )? (long)p: un;

basden = get_bas_den(bas);

invbas = QM_inv(vecpol_to_mat(bas,n), gun);

mul = get_mul_table(x,basden,invbas,NULL);

mul = get_mul_table(x,basden,invbas);

nf[7]=(long)bas;

nf[8]=(long)invbas;

nf[9]=(long)mul; nf[2]=nf[5]=nf[6]=zero;

Line 1656 padicff(GEN x,GEN p,long pr)

Line 1990 padicff(GEN x,GEN p,long pr)

}

GEN

factorpadic2(GEN x, GEN p, long r)

factorpadic2(GEN f, GEN p, long prec)

{

long av=avma,av2,k,i,j,i1,f,nbfac;

gpmem_t av = avma;

GEN res,p1,p2,y,d,a,ap,t,v,w;

GEN fa,ex,y;

GEN *fa;

long n,i,l;

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

if (typ(f)!=t_POL) err(notpoler,"factorpadic");

if (gcmp0(x)) err(zeropoler,"factorpadic");

if (gcmp0(f)) err(zeropoler,"factorpadic");

if (r<=0) err(rootper4);

if (prec<=0) err(rootper4);

if (lgef(x)==3) return trivfact();

n = degpol(f);

if (!gcmp1(leading_term(x)))

if (n==0) return trivfact();

if (n==1) return padic_trivfact(f,p,prec);

if (!gcmp1(leading_term(f)))

err(impl,"factorpadic2 for non-monic polynomial");

if (lgef(x)==4) return padic_trivfact(x,p,r);

f = padic_pol_to_int(f);

y=cgetg(3,t_MAT);

fa = (GEN*)new_chunk(lgef(x)-2);

fa = ZX_squff(f, &ex);

d=content(x); a=gdiv(x,d);

l = lg(fa); n = 0;

ap=derivpol(a); t=ggcd(a,ap); v=gdeuc(a,t);

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

w=gdeuc(ap,t); j=0; f=1; nbfac=0;

while (f)

{

j++; w=gsub(w,derivpol(v)); f=signe(w);

fa[i] = (long)padicff((GEN)fa[i],p,prec);

if (f) { res=ggcd(v,w); v=gdeuc(v,res); w=gdeuc(w,res); }

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

else res=v;

fa[j]=(lgef(res)>3) ? padicff(res,p,r) : cgetg(1,t_COL);

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

}

av2=avma; y=cgetg(3,t_MAT);

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

y = fact_from_DDF(fa,ex,n);

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

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

for (i=1,k=0; i<=j; i++)

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

{

p1[++k]=lcopy((GEN)fa[i][i1]); p2[k]=lstoi(i);

}

y = gerepile(av,av2,y);

sort_factor(y, cmp_padic); return y;

}

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

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

/* */

/* Using ROUND 4 */

/* FACTORISATION P-adique avec ROUND 4 */

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

/* */

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

extern GEN Decomp(GEN p,GEN f,long mf,GEN theta,GEN chi,GEN nu,long r);

extern GEN nilord(GEN p, GEN fx, long mf, GEN gx, long flag);

extern GEN hensel_lift_fact(GEN pol, GEN Q, GEN T, GEN p, GEN pev, long e);

static GEN

static int

squarefree(GEN f, GEN *ex)

expo_is_squarefree(GEN e)

{

GEN T,V,W,A,B;

long i, l = lg(e);

long n,i,k;

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

if (e[i] != 1) return 0;

T=ggcd(derivpol(f),f); V=gdeuc(f,T);

return 1;

n=lgef(f)-2; A=cgetg(n,t_COL); B=cgetg(n,t_COL);

k=1; i=1;

{

W=ggcd(T,V); T=gdeuc(T,W);

if (lgef(V) != lgef(W))

{

A[i]=ldeuc(V,W); B[i]=k; i++;

}

k++; V=W;

}

while (lgef(V)>3);

setlg(A,i); *ex=B; return A;

}

#define swap(x,y) { long _t=x; x=y; y=_t; }

/* reverse x in place */

static void

polreverse(GEN x)

{

long i, j;

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

for (i=2, j=lgef(x)-1; i<j; i++, j--) swap(x[i], x[j]);

(void)normalizepol(x);

}

GEN

factorpadic4(GEN f,GEN p,long prec)

{

GEN w,g,poly,fx,y,p1,p2,ex,pols,exps,ppow,lead;

gpmem_t av = avma;

long v=varn(f),n=degpol(f),av,tetpil,mfx,i,k,j,r,pr;

GEN w,g,poly,y,p1,p2,ex,pols,exps,ppow,lead;

long v=varn(f),n=degpol(f),mfx,i,k,j,r,pr;

int reverse = 0;

if (typ(f)!=t_POL) err(notpoler,"factorpadic");

Line 1751 factorpadic4(GEN f,GEN p,long prec)

Line 2049 factorpadic4(GEN f,GEN p,long prec)

if (prec<=0) err(rootper4);

if (n==0) return trivfact();

av=avma; f = padic_pol_to_int(f);

if (n==1) return padic_trivfact(f,p,prec);

if (n==1) return gerepileupto(av, padic_trivfact(f,p,prec));

f = padic_pol_to_int(f);

lead = leading_term(f); pr = prec;

f = pnormalize(f, p, prec, n-1, &lead, &pr, &reverse);

if (!gcmp1(lead))

{

long val = ggval(lead,p), val1 = ggval(constant_term(f),p);

if (val1 < val)

{

reverse = 1; polreverse(f);

/* take care of loss of precision from leading coeff of factor

* (whose valuation is <= val) */

pr += val;

val = val1;

}

pr += val * (n-1);

}

f = pol_to_monic(f, &lead);

poly=squarefree(f,&ex);

poly = ZX_squff(f,&ex);

pols=cgetg(n+1,t_COL);

pols = cgetg(n+1,t_COL);

exps=cgetg(n+1,t_COL); n = lg(poly);

exps = cgetg(n+1,t_COL); n = lg(poly);

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

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

{

long av1 = avma;

gpmem_t av1 = avma;

fx=(GEN)poly[i]; mfx=ggval(discsr(fx),p);

GEN fx = (GEN)poly[i], fa = factmod0(fx,p);

w = (GEN)factmod(fx,p)[1];

w = (GEN)fa[1];

if (!mfx)

if (expo_is_squarefree((GEN)fa[2]))

{ /* no repeated factors: Hensel lift */

p1 = hensel_lift_fact(fx, lift_intern(w), NULL, p, gpowgs(p,pr), pr);

p1 = hensel_lift_fact(fx, w, NULL, p, gpowgs(p,pr), pr);

p2 = stoi(ex[i]);

for (k=1; k<lg(p1); k++,j++)

{

Line 1789 factorpadic4(GEN f,GEN p,long prec)

Line 2073 factorpadic4(GEN f,GEN p,long prec)

continue;

}

/* use Round 4 */

mfx = ggval(ZX_disc(fx),p);

r = lg(w)-1;

g = lift_intern((GEN)w[r]);

g = (GEN)w[r];

p2 = (r == 1)? nilord(p,fx,mfx,g,pr)

: Decomp(p,fx,mfx,polx[v],fx,g, (pr<=mfx)? mfx+1: pr);

if (p2)

{

p2 = gerepileupto(av1,p2);

p2 = gerepilecopy(av1,p2);

p1 = (GEN)p2[1];

p2 = (GEN)p2[2];

for (k=1; k<lg(p1); k++,j++)

{

pols[j]=p1[k];

pols[j] = p1[k];

exps[j]=lmulis((GEN)p2[k],ex[i]);

exps[j] = lmulis((GEN)p2[k],ex[i]);

}

else

{

avma=av1;

avma = av1;

pols[j]=(long)fx;

pols[j] = (long)fx;

exps[j]=lstoi(ex[i]); j++;

exps[j] = lstoi(ex[i]); j++;

}

if (lead)

{

p1 = gmul(polx[v],lead);

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

{

pols[i] = (long)primpart( unscale_pol((GEN)pols[i], lead) );

p2 = poleval((GEN)pols[i], p1);

y = cgetg(3,t_MAT);

pols[i] = ldiv(p2, content(p2));

p1 = cgetg(j,t_COL); p = icopy(p); ppow = gpowgs(p,prec);

}

tetpil=avma; y=cgetg(3,t_MAT);

p1 = cgetg(j,t_COL); ppow = gpowgs(p,prec); p = icopy(p);

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

{

if (reverse) polreverse((GEN)pols[i]);

p1[i] = (long)pol_to_padic((GEN)pols[i],ppow,p,prec);

}

y[1]=(long)p1; setlg(exps,j);

y[1] = (long)p1; setlg(exps,j);

y[2]=lcopy(exps); y = gerepile(av,tetpil,y);

y[2] = lcopy(exps);

sort_factor(y, cmp_padic); return y;

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

}

GEN

Line 1851 factorpadic0(GEN f,GEN p,long r,long flag)

Line 2129 factorpadic0(GEN f,GEN p,long r,long flag)

/* */

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

extern GEN to_Kronecker(GEN P, GEN Q);

extern GEN from_Kronecker(GEN z, GEN pol);

static GEN spec_Fq_pow_mod_pol(GEN x, GEN S, GEN T, GEN p);

static GEN spec_Fq_pow_mod_pol(GEN x, GEN p, GEN a, GEN S);

static GEN

to_fq(GEN x, GEN a, GEN p)

to_Fq(GEN x, GEN T, GEN p)

{

long i,lx = lgef(x);

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

GEN z = cgetg(3,t_POLMOD), pol = cgetg(lx,t_POL);

GEN z = cgetg(3,t_POLMOD), y;

pol[1] = x[1];

if (lx == 2) setsigne(pol, 0);

if (tx == t_INT)

y = mod(x,p);

else

for (i=2; i<lx; i++) pol[i] = (long)mod((GEN)x[i], p);

{

/* assume deg(pol) < deg(a) */

if (tx != t_POL) err(typeer,"to_Fq");

z[1] = (long)a;

lx = lgef(x);

z[2] = (long)pol; return z;

y = cgetg(lx,t_POL);

y[1] = x[1];

if (lx == 2) setsigne(y, 0);

else

for (i=2; i<lx; i++) y[i] = (long)mod((GEN)x[i], p);

}

/* assume deg(y) < deg(T) */

z[1] = (long)T;

z[2] = (long)y; return z;

}

/* x POLMOD over Fq, return lift(x^n) */

static GEN

Kronecker_powmod(GEN x, GEN mod, GEN n)

to_Fq_pol(GEN x, GEN T, GEN p)

{

long lim,av,av0 = avma, i,j,m,v = varn(x);

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

GEN y, p1, p = NULL, pol = NULL;

if (tx != t_POL) err(typeer,"to_Fq_pol");

lx = lgef(x);

for (i=lgef(mod)-1; i>1; i--)

for (i=2; i<lx; i++) x[i] = (long)to_Fq((GEN)x[i],T,p);

{

return x;

p1 = (GEN)mod[i];

if (typ(p1) == t_POLMOD) { pol = (GEN)p1[1] ; break; }

}

if (!pol) err(talker,"need POLMOD coeffs in Kronecker_powmod");

for (i=lgef(pol)-1; i>1; i--)

{

p1 = (GEN)pol[i];

if (typ(p1) == t_INTMOD) { p = (GEN)p1[1] ; break; }

}

if (!p) err(talker,"need Fq coeffs in Kronecker_powmod");

x = lift_intern(to_Kronecker(x,pol));

/* adapted from powgi */

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

p1 = n+2; m = *p1;

y=x; j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;

for (i=lgefint(n)-2;;)

{

for (; j; m<<=1,j--)

{

y = gsqr(y);

y = from_Kronecker(FpX(y,p), pol);

setvarn(y, v);

y = gres(y, mod);

y = lift_intern(to_Kronecker(y,pol));

if (m<0)

{

y = gmul(y,x);

y = from_Kronecker(FpX(y,p), pol);

setvarn(y, v);

y = gres(y, mod);

y = lift_intern(to_Kronecker(y,pol));

}

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

{

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

y = gerepilecopy(av, y);

}

if (--i == 0) break;

m = *++p1, j = BITS_IN_LONG;

}

y = from_Kronecker(FpX(y,p),pol);

setvarn(y, v); return gerepileupto(av0, y);

}

GEN

/* assume T a clone */

FpX_rand(long d1, long v, GEN p)

static GEN

copy_then_free(GEN T)

{

long i, d = d1+2;

GEN t = forcecopy(T); gunclone(T); return t;

GEN y;

y = cgetg(d,t_POL); y[1] = evalsigne(1) | evalvarn(v);

for (i=2; i<d; i++) y[i] = (long)genrand(p);

normalizepol_i(y,d); return y;

}

/* return a random polynomial in F_q[v], degree < d1 */

static GEN

FqX_rand(long d1, long v, GEN p, GEN a)

to_Fq_fact(GEN t, GEN ex, long nbf, int sort, GEN unfq, gpmem_t av)

{

long i,j, d = d1+2, k = lgef(a)-1;

GEN T = (GEN)unfq[1]/*clone*/, y = cgetg(3,t_MAT), u,v,p;

GEN y,t;

long j,k, l = lg(t);

y = cgetg(d,t_POL); y[1] = evalsigne(1) | evalvarn(v);

u = cgetg(nbf,t_COL); y[1] = (long)u;

t = cgetg(k,t_POL); t[1] = a[1];

v = cgetg(nbf,t_COL); y[2] = (long)v;

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

for (j=1,k=0; j<l; j++)

{

if (ex[j])

for (j=2; j<k; j++) t[j] = (long)genrand(p);

{

normalizepol_i(t,k); y[i]=(long)to_fq(t,a,p);

k++;

}

u[k] = (long)simplify((GEN)t[j]); /* may contain pols of degree 0 */

normalizepol_i(y,d); return y;

v[k] = lstoi(ex[j]);

}

y = gerepileupto(av, y);

if (sort) y = sort_factor(y,cmp_pol);

T = copy_then_free(T);

p = (GEN)leading_term(T)[1];

u = (GEN)y[1];

for (j=1; j<nbf; j++) u[j] = (long)to_Fq_pol((GEN)u[j], T,p);

return y;

}

/* split into r factors of degree d */

static void

split9(GEN *t, long d, GEN p, GEN q, GEN a, GEN S)

FqX_split(GEN *t, long d, GEN q, GEN S, GEN T, GEN p)

{

long l,v,av,is2,cnt, dt = degpol(*t), da = degpol(a);

long l, v, is2, cnt, dt = degpol(*t), dT = degpol(T);

gpmem_t av;

GEN w,w0;

if (dt == d) return;

v = varn(*t);

if (DEBUGLEVEL > 6) timer2();

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

av = avma; is2 = egalii(p, gdeux);

for(cnt = 1;;cnt++)

{ /* splits *t with probability ~ 1 - 2^(1-r) */

w = w0 = FqX_rand(dt,v, p,a);

w = w0 = FqX_rand(dt,v, T,p);

for (l=1; l<d; l++) /* sum_{0<i<d} w^(q^i), result in (F_q)^r */

w = gadd(w0, spec_Fq_pow_mod_pol(w, p, a, S));

w = gadd(w0, spec_Fq_pow_mod_pol(w, S, T, p));

w = FqX_red(w, T,p);

if (is2)

{

w0 = w;

for (l=1; l<da; l++) /* sum_{0<i<k} w^(2^i), result in (F_2)^r */

for (l=1; l<dT; l++) /* sum_{0<i<k} w^(2^i), result in (F_2)^r */

w = gadd(w0, gres(gsqr(w), *t));

{

w = FqX_rem(FqX_sqr(w,T,p), *t, T,p);

w = FqX_red(gadd(w0,w), NULL,p);

}

else

{

w = Kronecker_powmod(w, *t, shifti(q,-1));

w = FqXQ_pow(w, shifti(q,-1), *t, T, p);

/* w in {-1,0,1}^r */

if (lgef(w) == 3) continue;

if (!degpol(w)) continue;

w[2] = ladd((GEN)w[2], gun);

}

w = ggcd(*t,w); l = degpol(w);

w = FqX_gcd(*t,w, T,p); l = degpol(w);

if (l && l != dt) break;

avma = av;

}

w = gerepileupto(av,w);

if (DEBUGLEVEL > 6)

fprintferr("[split9] time for splitting: %ld (%ld trials)\n",timer2(),cnt);

fprintferr("[FqX_split] splitting time: %ld (%ld trials)\n",timer2(),cnt);

l /= d; t[l]=gdeuc(*t,w); *t=w;

l /= d; t[l] = FqX_div(*t,w, T,p); *t = w;

split9(t+l,d,p,q,a,S);

FqX_split(t+l,d,q,S,T,p);

split9(t ,d,p,q,a,S);

FqX_split(t ,d,q,S,T,p);

}

/* to "compare" (real) scalars and t_INTMODs */

Line 2019 cmp_pol(GEN x, GEN y)

Line 2267 cmp_pol(GEN x, GEN y)

return 0;

}

/* assume n > 1, X a POLMOD over Fq */

/* assume n > 1, X over FqX */

/* return S = [ X^q, X^2q, ... X^(n-1)q ] mod T (in Fq[X]) in Kronecker form */

/* return S = [ X^q, X^2q, ... X^(n-1)q ] mod u (in Fq[X]) in Kronecker form */

static GEN

init_pow_q_mod_pT(GEN Xmod, GEN q, GEN a, GEN T)

init_pow_q_mod_pT(GEN X, GEN q, GEN u, GEN T, GEN p)

{

long i, n = degpol(T);

long i, n = degpol(u);

GEN p1, S = cgetg(n, t_VEC);

S[1] = (long)Kronecker_powmod((GEN)Xmod[2], (GEN)Xmod[1], q);

S[1] = (long)FqXQ_pow(X, q, u, T, p);

#if 1 /* use as many squarings as possible */

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

{

p1 = gsqr((GEN)S[i>>1]);

S[i] = lres(p1, T);

S[i] = (long)FqX_rem(p1, u, T,p);

if (i == n-1) break;

p1 = gmul((GEN)S[i], (GEN)S[1]);

S[i+1] = lres(p1, T);

S[i+1] = (long)FqX_rem(p1, u, T,p);

}

#else

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

{

p1 = gmul((GEN)S[i-1], (GEN)S[1]);

S[i] = lres(p1, T);

S[i] = (long)FqX_rem(p1, u, T,p);

}

#endif

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

S[i] = (long)lift_intern(to_Kronecker((GEN)S[i], a));

S[i] = (long)to_Kronecker((GEN)S[i], T);

return S;

}

/* compute x^q, S is as above */

static GEN

spec_Fq_pow_mod_pol(GEN x, GEN p, GEN a, GEN S)

spec_Fq_pow_mod_pol(GEN x, GEN S, GEN T, GEN p)

{

long av = avma, lim = stack_lim(av,1), i,dx = degpol(x);

long i, dx = degpol(x);

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

GEN x0 = x+2, z,c;

c = (GEN)x0[0];

z = c = (GEN)x0[0];

z = lift_intern(lift(c));

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

{

GEN d;

c = (GEN)x0[i];

if (gcmp0(c)) continue;

d = (GEN)S[i];

if (!gcmp1(c)) d = gmul(lift_intern(lift(c)),d);

if (!gcmp1(c)) d = gmul(c,d);

z = gadd(z, d);

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

{

Line 2072 spec_Fq_pow_mod_pol(GEN x, GEN p, GEN a, GEN S)

Line 2320 spec_Fq_pow_mod_pol(GEN x, GEN p, GEN a, GEN S)

z = gerepileupto(av, z);

}

z = FpX(z, p);

z = FpXQX_from_Kronecker(z, T, p);

z = from_Kronecker(z, a);

setvarn(z, varn(x)); return gerepileupto(av, z);

}

static long isabsolutepol(GEN f, GEN pp, GEN a)

static long

isabsolutepol(GEN f)

{

int i,res=1;

int i, l = lgef(f);

GEN c;

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

for(i=2; i<lg(f); i++)

{

c = (GEN) f[i];

GEN c = (GEN)f[i];

switch(typ(c))

if (typ(c) == t_POL && degpol(c) > 0) return 0;

{

case t_INT: /* OK*/

break;

case t_INTMOD:

if (gcmp((GEN)c[1],pp))

err(typeer,"factmod9");

break;

case t_POLMOD:

if (gcmp((GEN)c[1],a))

err(typeer,"factmod9");

isabsolutepol((GEN)c[1],pp,gzero);

isabsolutepol((GEN)c[2],pp,gzero);

if (degpol(c[1])>0)

res = 0;

break;

case t_POL:

isabsolutepol(c,pp,gzero);

if (degpol(c)>0)

res = 0;

break;

default:

err(typeer,"factmod9");

}

return res;

return 1;

}

GEN

factmod9(GEN f, GEN pp, GEN a)

factmod9(GEN f, GEN p, GEN T)

{

long av = avma, tetpil,p,i,j,k,d,e,vf,va,nbfact,nbf,pk;

gpmem_t av = avma;

GEN S,ex,y,f2,f3,df1,df2,g,g1,xmod,u,v,qqd,qq,unfp,unfq, *t;

long pg, i, j, k, d, e, N, vf, va, nbfact, nbf, pk;

GEN S,ex,f2,f3,df1,df2,g1,u,v,q,unfp,unfq, *t;

GEN frobinv,X;

if (typ(a)!=t_POL || typ(f)!=t_POL || gcmp0(a)) err(typeer,"factmod9");

if (typ(T)!=t_POL || typ(f)!=t_POL || gcmp0(T)) err(typeer,"factmod9");

vf=varn(f); va=varn(a);

vf = varn(f);

if (va<=vf) err(talker,"polynomial variable must be of higher priority than finite field\nvariable in factorff");

va = varn(T);

if (isabsolutepol(f, pp, a))

if (va <= vf)

err(talker,"polynomial variable must have higher priority in factorff");

unfp = gmodulsg(1,p); T = gmul(unfp,T);

unfq = gmodulo(gmul(unfp,polun[va]), T);

f = gmul(unfq,f);

if (!signe(f)) err(zeropoler,"factmod9");

d = degpol(f); if (!d) { avma = av; return trivfact(); }

f = simplify(lift_intern(lift_intern(f)));

T = lift_intern(T);

if (isabsolutepol(f))

{

GEN z= Fp_factor_rel0(simplify(lift(lift(f))), pp, lift(a));

GEN z = Fp_factor_rel0(f, p, T);

GEN t=(GEN)z[1],ex=(GEN)z[2];

return to_Fq_fact((GEN)z[1],(GEN)z[2],lg(z[1]), 0, unfq,av);

unfp = gmodulsg(1,pp);

unfq = gmodulcp(gmul(unfp, polun[va]), gmul(unfp,a));

nbfact=lg(t);

tetpil=avma;

y=cgetg(3,t_MAT);

u=cgetg(nbfact,t_COL); y[1]=(long)u;

v=cgetg(nbfact,t_COL); y[2]=(long)v;

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

{

u[j] = lmul((GEN)t[j],unfq);

v[j] = lstoi(ex[j]);

}

return gerepile(av,tetpil,y);

}

p = is_bigint(pp)? 0: itos(pp);

unfp=gmodulsg(1,pp); a=gmul(unfp,a);

unfq=gmodulo(gmul(unfp,polun[va]), a); a = (GEN)unfq[1];

f = gmul(unfq,f); if (!signe(f)) err(zeropoler,"factmod9");

d = degpol(f); if (!d) { avma=av; gunclone(a); return trivfact(); }

pg = is_bigint(p)? 0: itos(p);

S = df2 = NULL; /* gcc -Wall */

pp = gmael(a,2,1); /* out of the stack */

t = (GEN*)cgetg(d+1,t_VEC); ex = new_chunk(d+1);

frobinv = gpowgs(pp, lgef(a)-4);

frobinv = gpowgs(p, degpol(T)-1);

xmod = cgetg(3,t_POLMOD);

X = gmul(polx[vf],unfq);

X = polx[vf];

xmod[2] = (long)X;

q = gpowgs(p, degpol(T));

qq=gpuigs(pp,degpol(a));

e = nbfact = 1;

pk=1; df1=derivpol(f); f3=NULL;

pk = 1;

f3 = NULL;

df1 = FqX_deriv(f, T, p);

for(;;)

{

long du,dg;

while (gcmp0(df1))

{ /* needs d >= pp: p = 0 can't happen */

{ /* needs d >= p: pg = 0 can't happen */

pk *= p; e=pk;

pk *= pg; e = pk;

j=(degpol(f))/p+3; setlg(f,j); setlgef(f,j);

j = (degpol(f))/pg+3; setlg(f,j); setlgef(f,j);

for (i=2; i<j; i++) f[i] = (long)powgi((GEN)f[p*(i-2)+2], frobinv);

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

df1=derivpol(f); f3=NULL;

f[i] = (long)Fq_pow((GEN)f[pg*(i-2)+2], frobinv, T,p);

df1 = FqX_deriv(f, T, p); f3 = NULL;

}

f2 = f3? f3: ggcd(f,df1);

f2 = f3? f3: FqX_gcd(f,df1, T,p);

if (lgef(f2)==3) u = f;

if (!degpol(f2)) u = f;

else

{

g1=gdeuc(f,f2); df2=derivpol(f2);

g1 = FqX_div(f,f2, T,p); df2 = FqX_deriv(f2, T,p);

if (gcmp0(df2)) { u=g1; f3=f2; }

if (gcmp0(df2)) { u = g1; f3 = f2; }

else

{

f3=ggcd(f2,df2);

f3 = FqX_gcd(f2,df2, T,p);

if (lgef(f3)==3) u=gdeuc(g1,f2);

if (!degpol(f3)) u = FqX_div(g1,f2, T,p);

else

u=gdeuc(g1,gdeuc(f2,f3));

u = FqX_div(g1, FqX_div(f2,f3, T,p), T,p);

}

/* u is square-free (product of irreducibles of multiplicity e) */

qqd=gun; xmod[1]=(long)u;

#if 0

du = degpol(u); v = X;

N = degpol(u);

if (du > 1) S = init_pow_q_mod_pT(xmod, qq, a, u);

if (N)

for (d=1; d <= du>>1; d++)

{

qqd=mulii(qqd,qq);

t[nbfact] = FpXQX_normalize(u, T,p);

v = spec_Fq_pow_mod_pol(v, pp, a, S);

d = (N==1)? 1: FqX_split_berlekamp(t+nbfact, q, T, p);

g = ggcd(gsub(v,X),u);

for (j=0; j<d; j++) ex[nbfact+j] = e;

nbfact += d;

}

#else

{

GEN qqd = gun, g;

long dg;

N = degpol(u); v = X;

if (N > 1) S = init_pow_q_mod_pT(X, q, u, T, p);

for (d=1; d <= N>>1; d++)

{

qqd = mulii(qqd,q);

v = spec_Fq_pow_mod_pol(v, S, T, p);

g = FqX_gcd(gsub(v,X),u, T,p);

dg = degpol(g);

if (dg <= 0) continue;

Line 2198 factmod9(GEN f, GEN pp, GEN a)

Line 2429 factmod9(GEN f, GEN pp, GEN a)

j = nbfact+dg/d;

t[nbfact] = g;

split9(t+nbfact,d,pp,qq,a,S);

FqX_split(t+nbfact,d,q,S,T,p);

for (; nbfact<j; nbfact++) ex[nbfact]=e;

for (; nbfact<j; nbfact++) ex[nbfact] = e;

du -= dg;

N -= dg;

u = gdeuc(u,g);

u = FqX_div(u,g, T,p);

v = gres(v,u); xmod[1] = (long)u;

v = FqX_rem(v,u, T,p);

}

if (du) { t[nbfact]=u; ex[nbfact++]=e; }

if (N) { t[nbfact] = u; ex[nbfact++] = e; }

if (lgef(f2) == 3) break;

}

#endif

if (!degpol(f2)) break;

f=f2; df1=df2; e += pk;

f = f2; df1 = df2; e += pk;

}

nbf=nbfact; tetpil=avma; y=cgetg(3,t_MAT);

nbf = nbfact; setlg(t, nbfact);

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

{

t[j]=gdiv((GEN)t[j],leading_term(t[j]));

t[j] = FpXQX_normalize((GEN)t[j], T,p);

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

if (ex[k] && gegal(t[j],t[k]))

{

Line 2221 factmod9(GEN f, GEN pp, GEN a)

Line 2454 factmod9(GEN f, GEN pp, GEN a)

nbf--; break;

}

u=cgetg(nbf,t_COL); y[1]=(long)u;

return to_Fq_fact((GEN)t,ex,nbf, 1, unfq,av);

v=cgetg(nbf,t_COL); y[2]=(long)v;

for (j=1,k=0; j<nbfact; j++)

if (ex[j])

{

k++;

u[k]=(long)t[j];

v[k]=lstoi(ex[j]);

}

y = gerepile(av,tetpil,y);

u=(GEN)y[1];

{ /* put a back on the stack */

GEN tokill = a;

a = forcecopy(a);

gunclone(tokill);

}

pp = (GEN)leading_term(a)[1];

for (j=1; j<nbf; j++) fqunclone((GEN)u[j], a, pp);

(void)sort_factor(y, cmp_pol); return y;

}

/* See also: Isomorphisms between finite field and relative

* factorization in polarit3.c */

Line 2258 GEN zrhqr(GEN a,long PREC);

Line 2473 GEN zrhqr(GEN a,long PREC);

GEN

rootsold(GEN x, long l)

{

long av1=avma,i,j,f,g,gg,fr,deg,l0,l1,l2,l3,l4,ln;

long i, j, f, g, gg, fr, deg, ln;

gpmem_t av=avma, av0, av1, av2, av3;

long exc,expmin,m,deg0,k,ti,h,ii,e,e1,emax,v;

GEN y,xc,xd0,xd,xdabs,p1,p2,p3,p4,p5,p6,p7,p8;

GEN p9,p10,p11,p12,p14,p15,pa,pax,pb,pp,pq,ps;

Line 2283 rootsold(GEN x, long l)

Line 2499 rootsold(GEN x, long l)

if (gsigne(p2)>0) g=0;

} else if (ti != t_INT && ti != t_INTMOD && !is_frac_t(ti)) g=0;

}

l1=avma; p2=cgetg(3,t_COMPLEX);

av1=avma; p2=cgetg(3,t_COMPLEX);

p2[1]=lmppi(DEFAULTPREC); p2[2]=ldivrs((GEN)p2[1],10);

p11=cgetg(4,t_POL); p11[1]=evalsigne(1)+evallgef(4);

setvarn(p11,v); p11[3]=un;

Line 2317 rootsold(GEN x, long l)

Line 2533 rootsold(GEN x, long l)

deg=degpol(ps);

if (deg)

{

l3=avma; e=gexpo((GEN)ps[deg+2]); emax=e;

av3=avma; e=gexpo((GEN)ps[deg+2]); emax=e;

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

{

p3=(GEN)(ps[i]);

e1=gexpo(p3); if (e1>emax) emax=e1;

}

e=emax-e; if (e<0) e=0; avma=l3; if (ps!=pax) xd0=deriv(ps,v);

e=emax-e; if (e<0) e=0; avma=av3; if (ps!=pax) xd0=deriv(ps,v);

xdabs=cgetg(deg+2,t_POL); xdabs[1]=xd0[1];

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

{

l3=avma; p3=(GEN)xd0[i];

av3=avma; p3=(GEN)xd0[i];

p4=gabs(greal(p3),l);

p5=gabs(gimag(p3),l); l4=avma;

p5=gabs(gimag(p3),l);

xdabs[i]=lpile(l3,l4,gadd(p4,p5));

xdabs[i]=lpileupto(av3, gadd(p4,p5));

}

l0=avma; xc=gcopy(ps); xd=gcopy(xd0); l2=avma;

av0=avma; xc=gcopy(ps); xd=gcopy(xd0); av2=avma;

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

{

if (i==deg)

Line 2346 rootsold(GEN x, long l)

Line 2562 rootsold(GEN x, long l)

while (exc >= -20)

{

p7 = gneg_i(gdiv(p4, poleval(xd,p3)));

l3 = avma;

av3 = avma;

if (gcmp0(p5)) exc = -32;

else exc = expo(gnorm(p7))-expo(gnorm(p3));

avma = l3;

avma = av3;

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

{

p8=gadd(p3,p7); p9=poleval(xc,p8); p10=gnorm(p9);

Line 2358 rootsold(GEN x, long l)

Line 2574 rootsold(GEN x, long l)

GEN *gptr[3];

p3=p8; p4=p9; p5=p10;

gptr[0]=&p3; gptr[1]=&p4; gptr[2]=&p5;

gerepilemanysp(l2,l3,gptr,3);

gerepilemanysp(av2,av3,gptr,3);

break;

}

gshiftz(p7,-2,p7); avma=l3;

gshiftz(p7,-2,p7); avma=av3;

}

if (j > 10)

{

avma=av1;

avma=av;

if (DEBUGLEVEL)

{

fprintferr("too many iterations in rootsold(): ");

Line 2375 rootsold(GEN x, long l)

Line 2591 rootsold(GEN x, long l)

}

p1=(GEN)y[k+m*i]; setlg(p1[1],3); setlg(p1[2],3); gaffect(p3,p1);

avma=l2; p14=(GEN)p1[1]; p15=(GEN)p1[2];

avma=av2; p14=(GEN)p1[1]; p15=(GEN)p1[2];

for (ln=4; ln<=l; ln=(ln<<1)-2)

{

setlg(p14,ln); setlg(p15,ln);

Line 2384 rootsold(GEN x, long l)

Line 2600 rootsold(GEN x, long l)

p4=poleval(xc,p1);

p5=poleval(xd,p1); p6=gneg_i(gdiv(p4,p5));

settyp(p14,t_REAL); settyp(p15,t_REAL);

gaffect(gadd(p1,p6),p1); avma=l2;

gaffect(gadd(p1,p6),p1); avma=av2;

}

setlg(p14,l); setlg(p15,l);

Line 2404 rootsold(GEN x, long l)

Line 2620 rootsold(GEN x, long l)

p6=gdiv(p4,poleval(xdabs,gabs(p7,l)));

if (gexpo(p6)>=expmin)

{

avma=av1;

avma=av;

if (DEBUGLEVEL)

{

fprintferr("internal error in rootsold(): using roots2()\n");

Line 2412 rootsold(GEN x, long l)

Line 2628 rootsold(GEN x, long l)

}

return roots2(x,l);

}

avma=l2;

avma=av2;

if (expo(p1[2])<expmin && g)

{

gaffect(gzero,(GEN)p1[2]);

for (j=1; j<m; j++) gaffect(p1,(GEN)y[k+(i-1)*m+j]);

p11[2]=lneg((GEN)p1[1]);

l4=avma; xc=gerepile(l0,l4,gdeuc(xc,p11));

xc=gerepileupto(av0, gdeuc(xc,p11));

}

else

{

Line 2428 rootsold(GEN x, long l)

Line 2644 rootsold(GEN x, long l)

p1=gconj(p1);

for (j=1; j<=m; j++) gaffect(p1,(GEN)y[k+i*m+j]);

i++;

p12[2]=lnorm(p1); p12[3]=lmulsg(-2,(GEN)p1[1]); l4=avma;

p12[2]=lnorm(p1); p12[3]=lmulsg(-2,(GEN)p1[1]);

xc=gerepile(l0,l4,gdeuc(xc,p12));

xc=gerepileupto(av0, gdeuc(xc,p12));

}

else

{

p11[2]=lneg(p1); l4=avma;

p11[2]=lneg(p1);

xc=gerepile(l0,l4,gdeuc(xc,p11));

xc=gerepileupto(av0, gdeuc(xc,p11));

}

xd=deriv(xc,v); l2=avma;

xd=deriv(xc,v); av2=avma;

}

k += deg*m;

}

Line 2445 rootsold(GEN x, long l)

Line 2661 rootsold(GEN x, long l)

}

while (k!=deg0);

}

avma=l1;

avma=av1;

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

{

p1 = (GEN)y[j];

Line 2466 rootsold(GEN x, long l)

Line 2682 rootsold(GEN x, long l)

GEN

roots2(GEN pol,long PREC)

{

ulong av = avma;

gpmem_t av = avma;

long N,flagexactpol,flagrealpol,flagrealrac,ti,i,j;

long nbpol,k,av1,multiqol,deg,nbroot,fr,f;

long nbpol, k, multiqol, deg, nbroot, fr, f;

gpmem_t av1;

GEN p1,p2,rr,EPS,qol,qolbis,x,b,c,*ad,v,tabqol;

if (typ(pol)!=t_POL) err(typeer,"roots2");

Line 2481 roots2(GEN pol,long PREC)

Line 2698 roots2(GEN pol,long PREC)

p2=gneg_i(gdiv((GEN)pol[2],p1));

return gerepilecopy(av,p2);

}

EPS=realun(3); setexpo(EPS, 12 - bit_accuracy(PREC));

EPS = real2n(12 - bit_accuracy(PREC), 3);

flagrealpol=1; flagexactpol=1;

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

{

Line 2611 roots2(GEN pol,long PREC)

Line 2828 roots2(GEN pol,long PREC)

static GEN

laguer(GEN pol,long N,GEN y0,GEN EPS,long PREC)

{

long av = avma, av1,MAXIT,iter,i,j;

long MAXIT, iter, j;

gpmem_t av = avma, av1;

GEN rac,erre,I,x,abx,abp,abm,dx,x1,b,d,f,g,h,sq,gp,gm,g2,*ffrac;

MAXIT=MR*MT; rac=cgetg(3,t_COMPLEX);

MAXIT = MR*MT; rac = cgetc(PREC);

rac[1]=lgetr(PREC); rac[2]=lgetr(PREC);

av1 = avma;

I=cgetg(3,t_COMPLEX); I[1]=un; I[2]=un;

ffrac=(GEN*)new_chunk(MR+1); for (i=0; i<=MR; i++) ffrac[i]=cgetr(PREC);

ffrac = (GEN*)new_chunk(MR+1);

affrr(dbltor(0.0), ffrac[0]); affrr(dbltor(0.5), ffrac[1]);

ffrac[0] = dbltor(0.0);

affrr(dbltor(0.25),ffrac[2]); affrr(dbltor(0.75),ffrac[3]);

ffrac[1] = dbltor(0.5);

affrr(dbltor(0.13),ffrac[4]); affrr(dbltor(0.38),ffrac[5]);

ffrac[2] = dbltor(0.25);

affrr(dbltor(0.62),ffrac[6]); affrr(dbltor(0.88),ffrac[7]);

ffrac[3] = dbltor(0.75);

affrr(dbltor(1.0),ffrac[8]);

ffrac[4] = dbltor(0.13);

ffrac[5] = dbltor(0.38);

ffrac[6] = dbltor(0.62);

ffrac[7] = dbltor(0.88);

ffrac[8] = dbltor(1.0);

x=y0;

for (iter=1; iter<=MAXIT; iter++)

{

Line 2631 laguer(GEN pol,long N,GEN y0,GEN EPS,long PREC)

Line 2852 laguer(GEN pol,long N,GEN y0,GEN EPS,long PREC)

d=gzero; f=gzero; abx=QuickNormL1(x,PREC);

for (j=N-1; j>=0; j--)

{

f=gadd(gmul(x,f),d); d=gadd(gmul(x,d),b);

f = gadd(gmul(x,f),d);

b=gadd(gmul(x,b),(GEN)pol[j+2]);

d = gadd(gmul(x,d),b);

erre=gadd(QuickNormL1(b,PREC),gmul(abx,erre));

b = gadd(gmul(x,b), (GEN)pol[j+2]);

erre = gadd(QuickNormL1(b,PREC), gmul(abx,erre));

}

erre=gmul(erre,EPS);

erre = gmul(erre,EPS);

if (gcmp(QuickNormL1(b,PREC),erre)<=0)

{

gaffect(x,rac); avma = av1; return rac;

Line 2674 laguer(GEN pol,long N,GEN y0,GEN EPS,long PREC)

Line 2896 laguer(GEN pol,long N,GEN y0,GEN EPS,long PREC)

static GEN

balanc(GEN x)

{

ulong av = avma;

gpmem_t av = avma;

long last,i,j, sqrdx = (RADIX<<1), n = lg(x);

GEN r,c,cofgen,a;

Line 2728 hqr(GEN mat) /* find all the eigenvalues of the matrix

Line 2950 hqr(GEN mat) /* find all the eigenvalues of the matrix

anorm=fabs(a[1][1]);

for (i=2; i<=n; i++) for (j=(i-1); j<=n; j++) anorm+=fabs(a[i][j]);

nn=n; t=0.0;

nn=n; t=0.;

if (DEBUGLEVEL>3) { fprintferr("* Finding eigenvalues\n"); flusherr(); }

while (nn>=1)

{

Line 2737 hqr(GEN mat) /* find all the eigenvalues of the matrix

Line 2959 hqr(GEN mat) /* find all the eigenvalues of the matrix

{

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

{

s=fabs(a[l-1][l-1])+fabs(a[l][l]); if (s==0.0) s=anorm;

s=fabs(a[l-1][l-1])+fabs(a[l][l]); if (s==0.) s=anorm;

if ((double)(fabs(a[l][l-1])+s)==s) break;

}

x=a[nn][nn];

if (l==nn){ wr[nn]=x+t; wi[nn--]=0.0; }

if (l==nn){ wr[nn]=x+t; wi[nn--]=0.; }

else

{

y=a[nn-1][nn-1];

Line 2749 hqr(GEN mat) /* find all the eigenvalues of the matrix

Line 2971 hqr(GEN mat) /* find all the eigenvalues of the matrix

if (l == nn-1)

{

p=0.5*(y-x); q=p*p+w; z=sqrt(fabs(q)); x+=t;

if (q>=0.0 || fabs(q)<=eps)

if (q>=0. || fabs(q)<=eps)

{

z=p+SIGN(z,p); wr[nn-1]=wr[nn]=x+z;

if (fabs(z)>eps) wr[nn]=x-w/z;

wi[nn-1]=wi[nn]=0.0;

wi[nn-1]=wi[nn]=0.;

}

else{ wr[nn-1]=wr[nn]=x+p; wi[nn-1]=-(wi[nn]=z); }

nn-=2;

}

else

{

p = q = r = 0.0; /* for lint */

p = q = r = 0.; /* for lint */

if (its==30) err(talker,"too many iterations in hqr");

if (its==10 || its==20)

{

Line 2780 hqr(GEN mat) /* find all the eigenvalues of the matrix

Line 3002 hqr(GEN mat) /* find all the eigenvalues of the matrix

v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1]));

if ((double)(u+v)==v) break;

}

for (i=m+2; i<=nn; i++){ a[i][i-2]=0.0; if (i!=(m+2)) a[i][i-3]=0.0; }

for (i=m+2; i<=nn; i++){ a[i][i-2]=0.; if (i!=(m+2)) a[i][i-3]=0.; }

for (k=m; k<=nn-1; k++)

{

if (k!=m)

{

p=a[k][k-1]; q=a[k+1][k-1];

r = (k != nn-1)? a[k+2][k-1]: 0.0;

r = (k != nn-1)? a[k+2][k-1]: 0.;

x = fabs(p)+fabs(q)+fabs(r);

if (x != 0.0) { p/=x; q/=x; r/=x; }

if (x != 0.) { p/=x; q/=x; r/=x; }

}

s = SIGN(sqrt(p*p+q*q+r*r),p);

if (s == 0.0) continue;

if (s == 0.) continue;

if (k==m)

{ if (l!=m) a[k][k-1] = -a[k][k-1]; }

Line 2819 hqr(GEN mat) /* find all the eigenvalues of the matrix

Line 3041 hqr(GEN mat) /* find all the eigenvalues of the matrix

}

for (j=2; j<=n; j++) /* ordering the roots */

{

x=wr[j]; y=wi[j]; if (y) flj=1; else flj=0;

x=wr[j]; y=wi[j]; if (y != 0.) flj=1; else flj=0;

for (k=j-1; k>=1; k--)

{

if (wi[k]) flk=1; else flk=0;

if (wi[k] != 0.) flk=1; else flk=0;

if (flk<flj) break;

if (!flk && !flj && wr[k]<=x) break;

if (flk&&flj&& wr[k]<x) break;

Line 2835 hqr(GEN mat) /* find all the eigenvalues of the matrix

Line 3057 hqr(GEN mat) /* find all the eigenvalues of the matrix

for (i=1; i<=n; i++) free(a[i]); free(a); eig=cgetg(n+1,t_COL);

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

{

if (wi[i])

if (wi[i] != 0.)

{

GEN p1 = cgetg(3,t_COMPLEX);

eig[i]=(long)p1;

Line 2854 hqr(GEN mat) /* find all the eigenvalues of the matrix

Line 3076 hqr(GEN mat) /* find all the eigenvalues of the matrix

GEN

zrhqr(GEN a,long prec)

{

ulong av = avma;

gpmem_t av = avma;

long i,j,prec2, n = degpol(a), ex = -bit_accuracy(prec);

GEN aa,b,p1,rt,rr,hess,x,dx,y,newval,oldval;

hess = cgetg(n+1,t_MAT);

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

{

p1 = cgetg(n+1,t_COL); hess[j] = (long)p1;

p1[1] = lneg(gdiv((GEN)a[n-j+2],(GEN)a[n+2]));

for (i=2; i<=n; i++) p1[i] = (i==(j+1))? un: zero;

}

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