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, |
|
Line 21 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
/***********************************************************************/ |
/***********************************************************************/ |
#include "pari.h" |
#include "pari.h" |
extern GEN get_bas_den(GEN bas); |
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); |
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 */ |
/* see splitgen() for how to use these two */ |
GEN |
GEN |
setloop(GEN a) |
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; |
return a; |
} |
} |
|
|
|
|
int |
int |
gdivise(GEN x, GEN y) |
gdivise(GEN x, GEN y) |
{ |
{ |
long av=avma; |
gpmem_t av=avma; |
x=gmod(x,y); avma=av; return gcmp0(x); |
x=gmod(x,y); avma=av; return gcmp0(x); |
} |
} |
|
|
int |
int |
poldivis(GEN x, GEN y, GEN *z) |
poldivis(GEN x, GEN y, GEN *z) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
GEN p1 = poldivres(x,y,ONLY_DIVIDES); |
GEN p1 = poldivres(x,y,ONLY_DIVIDES); |
if (p1) { *z = p1; return 1; } |
if (p1) { *z = p1; return 1; } |
avma=av; return 0; |
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 = ONLY_REM return remainder, otherwise return quotient |
* if z != NULL set *z to remainder |
* if z != NULL set *z to remainder |
* *z is the last object on stack (and thus can be disposed of with cgiv |
* *z is the last object on stack (and thus can be disposed of with cgiv |
* instead of gerepile) |
* instead of gerepile) */ |
*/ |
|
GEN |
GEN |
poldivres(GEN x, GEN y, GEN *pr) |
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; |
long ty=typ(y),tx,vx,vy,dx,dy,dz,i,j,sx,lrem; |
int remainder; |
|
GEN z,p1,rem,y_lead,mod; |
GEN z,p1,rem,y_lead,mod; |
GEN (*f)(GEN,GEN); |
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 (is_scalar_t(ty)) |
{ |
{ |
if (pr == ONLY_REM) return gzero; |
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); |
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 */ |
if (gcmp0(y_lead)) /* normalize denominator if leading term is 0 */ |
{ |
{ |
err(warner,"normalizing a polynomial with 0 leading term"); |
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)); |
return f(x, constant_term(y)); |
} |
} |
dx=degpol(x); |
dx = degpol(x); |
if (vx>vy || dx<dy) |
if (vx>vy || dx<dy) |
{ |
{ |
if (pr) |
if (pr) |
Line 181 poldivres(GEN x, GEN y, GEN *pr) |
|
Line 180 poldivres(GEN x, GEN y, GEN *pr) |
|
} |
} |
return zeropol(vy); |
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); |
p1 = new_chunk(dy+3); |
for (i=2; i<dy+3; i++) |
for (i=2; i<dy+3; i++) |
{ |
{ |
GEN p2 = (GEN)y[i]; |
GEN p2 = (GEN)y[i]; |
|
/* gneg to avoid gsub's later on */ |
p1[i] = isexactzero(p2)? 0: (long)gneg_i(p2); |
p1[i] = isexactzero(p2)? 0: (long)gneg_i(p2); |
} |
} |
y = p1; |
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; |
default: if (gcmp1(y_lead)) y_lead = NULL; |
mod = 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; |
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--) |
for (i=dx-1; i>=dy; i--) |
{ |
{ |
av1=avma; p1=(GEN)x[i]; |
av1=avma; p1=(GEN)x[i]; |
for (j=i-dy+1; j<=i && j<=dz; j++) |
for (j=i-dy+1; j<=i && j<=dz; j++) |
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 (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; |
z[i-dy] = (long)p1; |
} |
} |
if (!pr) return gerepileupto(av,z-2); |
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--) |
for (sx=0; ; i--) |
{ |
{ |
p1 = (GEN)x[i]; |
p1 = (GEN)x[i]; |
/* we always enter this loop at least once */ |
/* we always enter this loop at least once */ |
for (j=0; j<=i && j<=dz; j++) |
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 (mod && avma==av1) p1 = gmul(p1,mod); |
if (!gcmp0(p1)) { sx = 1; break; } /* remainder is non-zero */ |
if (!gcmp0(p1)) { sx = 1; break; } /* remainder is non-zero */ |
if (!isinexactreal(p1) && !isexactzero(p1)) break; |
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 (pr == ONLY_DIVIDES) |
{ |
{ |
if (sx) { avma=av; return NULL; } |
if (sx) { avma=av; return NULL; } |
avma = (long)rem; |
avma = (gpmem_t)rem; |
return gerepileupto(av,z-2); |
return gerepileupto(av,z-2); |
} |
} |
lrem=i+3; rem -= lrem; |
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[0]=evaltyp(t_POL) | evallg(lrem); |
rem[1]=evalsigne(1) | evalvarn(vx) | evallgef(lrem); |
rem[1]=evalsigne(1) | evalvarn(vx) | evallgef(lrem); |
rem += 2; |
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]; |
av1=avma; p1 = (GEN)x[i]; |
for (j=0; j<=i && j<=dz; j++) |
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 (mod && avma==av1) p1 = gmul(p1,mod); |
rem[i]=avma==av1? lcopy(p1):lpileupto(av1,p1); |
rem[i]=avma==av1? lcopy(p1):lpileupto(av1,p1); |
} |
} |
rem -= 2; |
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; |
z -= 2; |
{ |
{ |
GEN *gptr[2]; gptr[0]=&z; gptr[1]=&rem; |
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)); |
f = gmul(f, mod(gun,pp)); |
if (!signe(f)) err(zeropoler,"factmod"); |
if (!signe(f)) err(zeropoler,"factmod"); |
f = lift_intern(f); d = lgef(f); |
f = lift_intern(f); d = lgef(f); |
for (i=2; i <d; i++) |
for (i=2; i <d; i++) |
if (typ(f[i])!=t_INT) err(impl,"factormod for general polynomials"); |
if (typ(f[i])!=t_INT) err(impl,"factormod for general polynomials"); |
*F = f; return d-3; |
*F = f; return d-3; |
} |
} |
|
|
rootmod2(GEN f, GEN pp) |
rootmod2(GEN f, GEN pp) |
{ |
{ |
GEN g,y,ss,q,r, x_minus_s; |
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 (!(d = factmod_init(&f, pp, &p))) { avma=av; return cgetg(1,t_COL); } |
if (!p) err(talker,"prime too big in rootmod2"); |
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; |
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 */ |
/* by splitting */ |
GEN |
GEN |
rootmod(GEN f, GEN p) |
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; |
GEN y,pol,a,b,q,pol0; |
|
|
if (!factmod_init(&f, p, &i)) { avma=av; return cgetg(1,t_COL); } |
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); } |
if ((i & 1) == 0) { avma = av; return root_mod_even(f,i); } |
i=2; while (!signe(f[i])) i++; |
i=2; while (!signe(f[i])) i++; |
if (i == 2) j = 1; |
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); |
if (la) y[j+lb] = (long)FpX_normalize(a,p); |
pol = gadd(polx[varn(f)], gun); pol0 = constant_term(pol); |
pol = gadd(polx[varn(f)], gun); pol0 = constant_term(pol); |
while (j<=n) |
while (j<=n) |
{ |
{ /* cf FpX_split_berlekamp */ |
a=(GEN)y[j]; la=degpol(a); |
a = (GEN)y[j]; la = degpol(a); |
if (la==1) |
if (la==1) |
y[j++] = lsubii(p, constant_term(a)); |
y[j++] = lsubii(p, constant_term(a)); |
else if (la==2) |
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]++) |
else for (pol0[2]=1; ; pol0[2]++) |
{ |
{ |
b = ZX_s_add(FpXQ_pow(pol,q, a,p), -1); /* pol^(p-1)/2 - 1 */ |
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); |
b = FpX_gcd(a,b, p); lb = degpol(b); |
if (lb && lb<la) |
if (lb && lb < la) |
{ |
{ |
b = FpX_normalize(b, p); |
b = FpX_normalize(b, p); |
y[j+lb] = (long)FpX_div(a,b, p); |
y[j+lb] = (long)FpX_div(a,b, p); |
y[j] = (long)b; break; |
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)); |
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 */ |
/* 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); |
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. */ |
/* Functions giving information on the factorisation. */ |
|
|
/* u in Z[X], return kernel of (Frob - Id) over Fp[X] / u */ |
/* u in Z[X], return kernel of (Frob - Id) over Fp[X] / u */ |
static GEN |
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); |
Q = cgetg(N+1,t_MAT); Q[1] = (long)zerocol(N); |
w = v = FpXQ_pow(polx[varn(u)],p,u,p); |
w = v = FpXQ_pow(polx[varn(u)],p,u,p); |
for (j=2; j<=N; j++) |
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); |
p1[j] = laddis((GEN)p1[j], -1); |
|
Q[j] = (long)p1; |
if (j < N) |
if (j < N) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
w = gerepileupto(av, FpX_res(gmul(w,v), u, p)); |
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; |
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. */ |
/* f in ZZ[X] and p a prime number. */ |
long |
long |
FpX_is_squarefree(GEN f, GEN p) |
FpX_is_squarefree(GEN f, GEN p) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
GEN z; |
GEN z; |
z = FpX_gcd(f,derivpol(f),p); |
z = FpX_gcd(f,derivpol(f),p); |
avma = av; |
avma = av; |
Line 533 FpX_is_squarefree(GEN f, GEN p) |
|
Line 606 FpX_is_squarefree(GEN f, GEN p) |
|
long |
long |
FpX_nbroots(GEN f, GEN p) |
FpX_nbroots(GEN f, GEN p) |
{ |
{ |
long av = avma, n=lgef(f); |
long n=lgef(f); |
|
gpmem_t av = avma; |
GEN z; |
GEN z; |
if (n <= 4) return n-3; |
if (n <= 4) return n-3; |
f = FpX_red(f, p); |
f = FpX_red(f, p); |
Line 545 FpX_nbroots(GEN f, GEN p) |
|
Line 619 FpX_nbroots(GEN f, GEN p) |
|
long |
long |
FpX_is_totally_split(GEN f, GEN p) |
FpX_is_totally_split(GEN f, GEN p) |
{ |
{ |
long av = avma, n=lgef(f); |
long n=degpol(f); |
|
gpmem_t av = avma; |
GEN z; |
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); |
f = FpX_red(f, p); |
z = FpXQ_pow(polx[varn(f)], p, f, p); |
z = FpXQ_pow(polx[varn(f)], p, f, p); |
avma = av; return lgef(z)==4 && gcmp1((GEN)z[3]) && !signe(z[2]); |
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 |
long |
FpX_nbfact(GEN u, GEN p) |
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; |
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 |
GEN |
FpV_roots_to_pol(GEN V, GEN p, long v) |
FpV_roots_to_pol(GEN V, GEN p, long v) |
{ |
{ |
ulong ltop=avma; |
gpmem_t ltop=avma; |
long i; |
long i; |
GEN g=cgetg(lg(V),t_VEC); |
GEN g=cgetg(lg(V),t_VEC); |
for(i=1;i<lg(V);i++) |
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; |
modulo=p; |
g=divide_conquer_prod(g,&gsmul); |
g=divide_conquer_prod(g,&gsmul); |
return gerepileupto(ltop,g); |
return gerepileupto(ltop,g); |
|
|
y[2]=lgetg(1,t_COL); return y; |
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 |
static GEN |
try_pow(GEN w0, GEN pol, GEN p, GEN q, long r) |
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 |
static void |
split(long m, GEN *t, long d, GEN p, GEN q, long r, GEN S) |
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; |
GEN w,w0; |
|
|
dv=degpol(*t); if (dv==d) return; |
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) |
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++) |
for (l=1; l<d; l++) |
w = gadd(w0, spec_FpXQ_pow(w, p, S)); |
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 |
static void |
splitgen(GEN m, GEN *t, long d, GEN p, GEN q, long r) |
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; |
GEN w; |
|
|
dv=degpol(*t); if (dv==d) return; |
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); |
S[1] = (long)FpXQ_pow(polx[v], p, T, p); |
/* use as many squarings as possible */ |
/* use as many squarings as possible */ |
for (i=2; i < n; i+=2) |
for (i=2; i < n; i+=2) |
{ |
{ |
p1 = gsqr((GEN)S[i>>1]); |
p1 = gsqr((GEN)S[i>>1]); |
S[i] = (long)FpX_res(p1, T, p); |
S[i] = (long)FpX_res(p1, T, p); |
if (i == n-1) break; |
if (i == n-1) break; |
p1 = gmul((GEN)S[i], (GEN)S[1]); |
p1 = gmul((GEN)S[i], (GEN)S[1]); |
S[i+1] = (long)FpX_res(p1, T, p); |
S[i+1] = (long)FpX_res(p1, T, p); |
} |
} |
return S; |
return S; |
} |
} |
|
|
/* compute x^p, S is as above */ |
/* compute x^p, S is as above */ |
static GEN |
static GEN |
spec_FpXQ_pow(GEN x, GEN p, GEN S) |
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; |
GEN x0 = x+2, z; |
z = (GEN)x0[0]; |
z = (GEN)x0[0]; |
if (dx < 0) err(talker, "zero polynomial in FpXQ_pow. %Z not prime", p); |
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 |
static GEN |
factcantor0(GEN f, GEN pp, long flag) |
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 ex,y,f2,g,g1,u,v,pd,q; |
GEN *t; |
GEN *t; |
|
|
Line 861 simplefactmod(GEN f, GEN p) |
|
Line 922 simplefactmod(GEN f, GEN p) |
|
return factcantor0(f,p,1); |
return factcantor0(f,p,1); |
} |
} |
|
|
/* vector of polynomials (in v) whose coeffs are given by the columns of x */ |
|
GEN |
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; |
return y; |
} |
} |
|
|
/* matrix whose entries are given by the coeffs of the polynomials in |
/* matrix whose entries are given by the coeffs of the polynomials in |
* vector v (considered as degree n-1 polynomials) */ |
* vector v (considered as degree n-1 polynomials) */ |
GEN |
GEN |
vecpol_to_mat(GEN v, long n) |
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"); |
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; |
return y; |
} |
} |
|
|
/* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */ |
/* polynomial (in v) of polynomials (in w) whose coeffs are given by the columns of x */ |
GEN |
GEN |
mat_to_polpol(GEN x, long v,long w) |
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); |
GEN y = cgetg(lx+1, t_POL); |
y[1]=evalsigne(1) | evallgef(lx+1) | evalvarn(v); |
y[1]=evalsigne(1) | evallgef(lx+1) | evalvarn(v); |
y++; |
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 |
/* 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 |
GEN |
polpol_to_mat(GEN v, long n) |
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"); |
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 */ |
/* set x <-- x + c*y mod p */ |
static void |
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; |
long i,lx,ly,l; |
if (!c) return; |
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) |
if (p & ~MAXHALFULONG) |
{ |
{ |
for (i=2; i<l; i++) x[i] = ((ulong)x[i]+ (ulong)mulssmod(c,y[i],p)) % p; |
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 |
else |
{ |
{ |
for (i=2; i<l; i++) x[i] = ((ulong)x[i] + (ulong)(c*y[i])) % p; |
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 |
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); |
vker = mat_to_vecpol(vker,vu); |
d = lg(vker)-1; |
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); |
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; |
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); |
polt = small_to_pol(pol,vu); |
} |
} |
else |
else |
{ |
{ |
pol[2] = (long)genrand(pp); |
pol[2] = (long)genrand(p); |
pol[1] = evallgef(signe(pol[2])? 3: 2) | evalvarn(vu); |
pol[1] = evallgef(signe(pol[2])? 3: 2) | evalvarn(vu); |
for (i=2; i<=d; i++) |
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; |
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; |
return d; |
} |
} |
|
|
|
static GEN |
|
FqX_deriv(GEN f, GEN T, GEN p) |
|
{ |
|
return FqX_red(derivpol(f), T, p); |
|
} |
GEN |
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) |
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; |
GEN pps2,ex,y,f2,p1,g1,u, *t; |
|
|
if (!(d = factmod_init(&f, pp, &p))) { avma=av; return trivfact(); } |
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); } |
k++; if (p && !(k%p)) { k++; f2 = FpX_div(f2,g1,pp); } |
p1 = FpX_gcd(f2,g1, pp); u = g1; g1 = p1; |
p1 = FpX_gcd(f2,g1, pp); u = g1; g1 = p1; |
if (lgef(p1)!=3) |
if (degpol(p1)) |
{ |
{ |
u = FpX_div( u,p1,pp); |
u = FpX_div( u,p1,pp); |
f2= FpX_div(f2,p1,pp); |
f2= FpX_div(f2,p1,pp); |
} |
} |
|
/* u is square-free (product of irred. of multiplicity e * k) */ |
N = degpol(u); |
N = degpol(u); |
if (N) |
if (N) |
{ |
{ |
/* here u is square-free (product of irred. of multiplicity e * k) */ |
|
t[nbfact] = FpX_normalize(u,pp); |
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; |
for (j=0; j<d; j++) ex[nbfact+j] = e*k; |
nbfact += d; |
nbfact += d; |
} |
} |
Line 1097 factmod0(GEN f, GEN pp) |
|
Line 1374 factmod0(GEN f, GEN pp) |
|
GEN |
GEN |
factmod(GEN f, GEN pp) |
factmod(GEN f, GEN pp) |
{ |
{ |
long tetpil,av=avma; |
gpmem_t tetpil, av=avma; |
long nbfact; |
long nbfact; |
long j; |
long j; |
GEN y,u,v; |
GEN y,u,v; |
|
|
padic_pol_to_int(GEN f) |
padic_pol_to_int(GEN f) |
{ |
{ |
long i, l = lgef(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++) |
for (i=2; i<l; i++) |
switch(typ(f[i])) |
switch(typ(f[i])) |
{ |
{ |
Line 1146 padic_pol_to_int(GEN f) |
|
Line 1430 padic_pol_to_int(GEN f) |
|
return 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 |
static GEN |
int_to_padic(GEN x, GEN p, GEN pr, long r, GEN invlead) |
int_to_padic(GEN x, GEN p, GEN pr, long r, GEN invlead) |
{ |
{ |
GEN p1,y; |
GEN p1,y; |
long v,sx, av = avma; |
long v, sx; |
|
gpmem_t av = avma; |
|
|
if (typ(x) == t_PADIC) |
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; |
y[4] = lmodii(p1,pr); r -= v; |
} |
} |
else |
else |
{ |
{ |
y[4] = zero; v = r; r = 0; |
y[4] = zero; v = r; r = 0; |
} |
} |
y[3] = (long)pr; |
y[3] = (long)pr; |
y[2] = (long)p; |
y[2] = (long)p; |
y[1] = evalprecp(r)|evalvalp(v); |
y[1] = evalprecp(r)|evalvalp(v); |
return invlead? gerepileupto(av, gmul(invlead,y)): y; |
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 |
/* 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]) */ |
* is a power of p), x in Z[X] (or Z_p[X]) */ |
static GEN |
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; |
if (gcmp1(lead)) lead = NULL; |
else |
else |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
v = ggval(lead,p); |
v = ggval(lead,p); |
if (v) lead = gdiv(lead, gpowgs(p,v)); |
if (v) lead = gdiv(lead, gpowgs(p,v)); |
lead = int_to_padic(lead,p,pr,r,NULL); |
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 |
static GEN |
padic_trivfact(GEN x, GEN p, long r) |
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 |
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); |
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 |
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 (typ(f)!=t_POL) err(notpoler,"rootpadic"); |
if (gcmp0(f)) err(zeropoler,"rootpadic"); |
if (gcmp0(f)) err(zeropoler,"rootpadic"); |
if (r<=0) err(rootper4); |
if (prec <= 0) err(rootper4); |
f = padic_pol_to_int(f); |
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 */ |
/* rootpadicfast */ |
/*************************************************************************/ |
/*************************************************************************/ |
|
|
/*lift accelerator. The author of the idea is unknown.*/ |
/* lift accelerator */ |
long hensel_lift_accel(long n, long *pmask) |
long hensel_lift_accel(long n, long *pmask) |
{ |
{ |
long a,j; |
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 |
/* STANDARD USE |
There exists p a prime number and a>0 such that |
There exists p a prime number and a>0 such that |
q=p^a |
q=p^a |
f in ZZ[X], with leading term prime to p. |
f in ZZ[X], with leading term prime to p. |
S must be a simple root mod 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 |
GEN |
rootpadiclift(GEN T, GEN S, GEN p, long e) |
rootpadiclift(GEN T, GEN S, GEN p, long e) |
{ |
{ |
ulong ltop=avma; |
gpmem_t ltop=avma; |
long x; |
long x; |
GEN qold, q, qm1; |
GEN qold, q, qm1; |
GEN W, Tr, Sr, Wr = gzero; |
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=FpX_eval(deriv(Tr, x),S,q); |
W=mpinvmod(W,q); |
W=mpinvmod(W,q); |
for(i=0;i<nb;i++) |
for(i=0;i<nb;i++) |
{ |
{ |
qm1 = (mask&(1<<i))?sqri(qm1):mulii(qm1, q); |
qm1 = (mask&(1<<i))?sqri(qm1):mulii(qm1, q); |
q = mulii(qm1, p); |
q = mulii(qm1, p); |
Tr = FpX_red(T,q); |
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); |
y[n-1]=(long) rootpadiclift(f, (GEN) S[n-1], q, e); |
else/* distinct-->totally split-->use trace trick */ |
else/* distinct-->totally split-->use trace trick */ |
{ |
{ |
ulong av=avma; |
gpmem_t av=avma; |
GEN z; |
GEN z; |
z=(GEN)f[lgef(f)-2];/*-trace(roots)*/ |
z=(GEN)f[lgef(f)-2];/*-trace(roots)*/ |
for(i=1; i<n-1;i++) |
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. |
f must have no multiple roots mod p. |
|
|
return p-adics roots of f with prec pr, as integers (implicitly mod pr) |
return p-adics roots of f with prec pr, as integers (implicitly mod pr) |
|
|
*/ |
*/ |
GEN |
GEN |
rootpadicfast(GEN f, GEN p, long e) |
rootpadicfast(GEN f, GEN p, long e) |
{ |
{ |
ulong ltop=avma; |
gpmem_t ltop=avma; |
GEN S,y; |
GEN S,y; |
S=lift(rootmod(f,p));/*no multiplicity*/ |
S=lift(rootmod(f,p));/*no multiplicity*/ |
if (lg(S)==1)/*no roots*/ |
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; |
return y; |
} |
} |
/* Same as rootpadiclift for the polynomial X^n-a, |
/* Same as rootpadiclift for the polynomial X^n-a, |
* but here, n can be much larger. |
* but here, n can be much larger. |
* TODO: generalize to sparse polynomials. |
* TODO: generalize to sparse polynomials. |
*/ |
*/ |
GEN |
GEN |
padicsqrtnlift(GEN a, GEN n, GEN S, GEN p, long e) |
padicsqrtnlift(GEN a, GEN n, GEN S, GEN p, long e) |
{ |
{ |
ulong ltop=avma; |
gpmem_t ltop=avma; |
GEN qold, q, qm1; |
GEN qold, q, qm1; |
GEN W, Sr, Wr = gzero; |
GEN W, Sr, Wr = gzero; |
long i, nb, mask; |
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 = modii(mulii(n,powmodulo(S,subii(n,gun),q)),q); |
W = mpinvmod(W,q); |
W = mpinvmod(W,q); |
for(i=0;i<nb;i++) |
for(i=0;i<nb;i++) |
{ |
{ |
qm1 = (mask&(1<<i))?sqri(qm1):mulii(qm1, q); |
qm1 = (mask&(1<<i))?sqri(qm1):mulii(qm1, q); |
q = mulii(qm1, p); |
q = mulii(qm1, p); |
Sr = S; |
Sr = S; |
Line 1512 getprec(GEN x, long prec, GEN *p) |
|
Line 1836 getprec(GEN x, long prec, GEN *p) |
|
return prec; |
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 |
GEN |
apprgen9(GEN f, GEN a) |
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 (typ(f)!=t_POL) err(notpoler,"apprgen9"); |
if (gcmp0(f)) err(zeropoler,"apprgen9"); |
if (gcmp0(f)) err(zeropoler,"apprgen9"); |
if (typ(a)==t_PADIC) return apprgen(f,a); |
if (typ(a)==t_PADIC) return apprgen(f,a); |
if (typ(a)!=t_POLMOD || typ(a[2])!=t_POL) err(rootper1); |
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((GEN)a[2], BIGINT, &p); |
prec = getprec(t, prec, &p); |
prec = getprec(T, prec, &p); |
if (prec==BIGINT) err(rootper1); |
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); |
if (fl2 && v==1) err(rootper2); |
v=ggval(lift_intern(poleval(fp,a)), p); |
v=ggval(lift_intern(poleval(fp,a)), p); |
if (!v) |
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); |
tetpil=avma; pro=cgetg(2,t_COL); pro[1]=lcopy(a); |
return gerepile(av,tetpil,pro); |
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"); |
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) |
if (fl2) |
{ |
{ |
ps_1=3; x2=ggrandocp(p,2); p=stoi(4); |
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); |
ps_1=itos(p)-1; x2=ggrandocp(p,1); |
} |
} |
f = poleval(f,gadd(a,gmul(p,polx[varn(f)]))); |
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++) |
for (i=1; i<=d; i++) |
vecg[i] = (long)setloop(gzero); |
vecg[i] = (long)setloop(gzero); |
va=varn(t); |
va=varn(T); j = 1; |
for(;;) /* loop through F_q */ |
for(;;) /* loop through F_q */ |
{ |
{ |
ip=gmodulcp(gtopoly(vecg,va),t); |
ip=gmodulcp(gtopoly(vecg,va),T); |
if (gcmp0(poleval(f,gadd(ip,x2)))) |
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--) |
for (i=d; i; i--) |
{ |
{ |
Line 1584 apprgen9(GEN f, GEN a) |
|
Line 1906 apprgen9(GEN f, GEN a) |
|
} |
} |
if (!i) break; |
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 |
int |
cmp_padic(GEN x, GEN y) |
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]); |
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 |
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++) |
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++) |
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 |
static GEN |
padicff(GEN x,GEN p,long pr) |
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); |
basden = get_bas_den(bas); |
invbas = QM_inv(vecpol_to_mat(bas,n), gun); |
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[7]=(long)bas; |
nf[8]=(long)invbas; |
nf[8]=(long)invbas; |
nf[9]=(long)mul; nf[2]=nf[5]=nf[6]=zero; |
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 |
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"); |
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 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 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); |
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; |
|
do |
|
{ |
|
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 |
GEN |
factorpadic4(GEN f,GEN p,long prec) |
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; |
int reverse = 0; |
|
|
if (typ(f)!=t_POL) err(notpoler,"factorpadic"); |
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 (prec<=0) err(rootper4); |
|
|
if (n==0) return trivfact(); |
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 */ |
{ /* 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]); |
p2 = stoi(ex[i]); |
for (k=1; k<lg(p1); k++,j++) |
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; |
continue; |
} |
} |
/* use Round 4 */ |
/* use Round 4 */ |
|
mfx = ggval(ZX_disc(fx),p); |
r = lg(w)-1; |
r = lg(w)-1; |
g = lift_intern((GEN)w[r]); |
g = (GEN)w[r]; |
p2 = (r == 1)? nilord(p,fx,mfx,g,pr) |
p2 = (r == 1)? nilord(p,fx,mfx,g,pr) |
: Decomp(p,fx,mfx,polx[v],fx,g, (pr<=mfx)? mfx+1: pr); |
: Decomp(p,fx,mfx,polx[v],fx,g, (pr<=mfx)? mfx+1: pr); |
if (p2) |
if (p2) |
{ |
{ |
p2 = gerepileupto(av1,p2); |
p2 = gerepilecopy(av1,p2); |
p1 = (GEN)p2[1]; |
p1 = (GEN)p2[1]; |
p2 = (GEN)p2[2]; |
p2 = (GEN)p2[2]; |
for (k=1; k<lg(p1); k++,j++) |
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 |
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) |
if (lead) |
{ |
|
p1 = gmul(polx[v],lead); |
|
for (i=1; i<j; i++) |
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++) |
for (i=1; i<j; i++) |
{ |
{ |
if (reverse) polreverse((GEN)pols[i]); |
if (reverse) polreverse((GEN)pols[i]); |
p1[i] = (long)pol_to_padic((GEN)pols[i],ppow,p,prec); |
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 |
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 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 |
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 |
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 |
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 |
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 */ |
/* split into r factors of degree d */ |
static void |
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; |
GEN w,w0; |
|
|
if (dt == d) return; |
if (dt == d) return; |
v = varn(*t); |
v = varn(*t); |
if (DEBUGLEVEL > 6) timer2(); |
if (DEBUGLEVEL > 6) (void)timer2(); |
av = avma; is2 = egalii(p, gdeux); |
av = avma; is2 = egalii(p, gdeux); |
for(cnt = 1;;cnt++) |
for(cnt = 1;;cnt++) |
{ /* splits *t with probability ~ 1 - 2^(1-r) */ |
{ /* 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 */ |
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) |
if (is2) |
{ |
{ |
w0 = w; |
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 |
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 */ |
/* w in {-1,0,1}^r */ |
if (lgef(w) == 3) continue; |
if (!degpol(w)) continue; |
w[2] = ladd((GEN)w[2], gun); |
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; |
if (l && l != dt) break; |
avma = av; |
avma = av; |
} |
} |
w = gerepileupto(av,w); |
w = gerepileupto(av,w); |
if (DEBUGLEVEL > 6) |
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 */ |
/* 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; |
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 |
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); |
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 */ |
#if 1 /* use as many squarings as possible */ |
for (i=2; i < n; i+=2) |
for (i=2; i < n; i+=2) |
{ |
{ |
p1 = gsqr((GEN)S[i>>1]); |
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; |
if (i == n-1) break; |
p1 = gmul((GEN)S[i], (GEN)S[1]); |
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 |
#else |
for (i=2; i < n; i++) |
for (i=2; i < n; i++) |
{ |
{ |
p1 = gmul((GEN)S[i-1], (GEN)S[1]); |
p1 = gmul((GEN)S[i-1], (GEN)S[1]); |
S[i] = lres(p1, T); |
S[i] = (long)FqX_rem(p1, u, T,p); |
} |
} |
#endif |
#endif |
for (i=1; i < n; i++) |
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; |
return S; |
} |
} |
|
|
/* compute x^q, S is as above */ |
/* compute x^q, S is as above */ |
static GEN |
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; |
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++) |
for (i = 1; i <= dx; i++) |
{ |
{ |
GEN d; |
GEN d; |
c = (GEN)x0[i]; |
c = (GEN)x0[i]; |
if (gcmp0(c)) continue; |
if (gcmp0(c)) continue; |
d = (GEN)S[i]; |
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); |
z = gadd(z, d); |
if (low_stack(lim, stack_lim(av,1))) |
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 = 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); |
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 |
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; |
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 */ |
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); |
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; |
e = nbfact = 1; |
pk=1; df1=derivpol(f); f3=NULL; |
pk = 1; |
|
f3 = NULL; |
|
df1 = FqX_deriv(f, T, p); |
for(;;) |
for(;;) |
{ |
{ |
long du,dg; |
|
while (gcmp0(df1)) |
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 |
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 |
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 |
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) */ |
/* 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); |
dg = degpol(g); |
if (dg <= 0) continue; |
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; |
j = nbfact+dg/d; |
|
|
t[nbfact] = g; |
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++) |
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++) |
for (k=1; k<j; k++) |
if (ex[k] && gegal(t[j],t[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; |
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 |
/* See also: Isomorphisms between finite field and relative |
* factorization in polarit3.c */ |
* factorization in polarit3.c */ |
Line 2258 GEN zrhqr(GEN a,long PREC); |
|
Line 2473 GEN zrhqr(GEN a,long PREC); |
|
GEN |
GEN |
rootsold(GEN x, long l) |
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; |
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 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; |
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; |
if (gsigne(p2)>0) g=0; |
} else if (ti != t_INT && ti != t_INTMOD && !is_frac_t(ti)) 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); |
p2[1]=lmppi(DEFAULTPREC); p2[2]=ldivrs((GEN)p2[1],10); |
p11=cgetg(4,t_POL); p11[1]=evalsigne(1)+evallgef(4); |
p11=cgetg(4,t_POL); p11[1]=evalsigne(1)+evallgef(4); |
setvarn(p11,v); p11[3]=un; |
setvarn(p11,v); p11[3]=un; |
Line 2317 rootsold(GEN x, long l) |
|
Line 2533 rootsold(GEN x, long l) |
|
deg=degpol(ps); |
deg=degpol(ps); |
if (deg) |
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++) |
for (i=2; i<deg+2; i++) |
{ |
{ |
p3=(GEN)(ps[i]); |
p3=(GEN)(ps[i]); |
e1=gexpo(p3); if (e1>emax) emax=e1; |
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]; |
xdabs=cgetg(deg+2,t_POL); xdabs[1]=xd0[1]; |
for (i=2; i<deg+2; i++) |
for (i=2; i<deg+2; i++) |
{ |
{ |
l3=avma; p3=(GEN)xd0[i]; |
av3=avma; p3=(GEN)xd0[i]; |
p4=gabs(greal(p3),l); |
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++) |
for (i=1; i<=deg; i++) |
{ |
{ |
if (i==deg) |
if (i==deg) |
Line 2346 rootsold(GEN x, long l) |
|
Line 2562 rootsold(GEN x, long l) |
|
while (exc >= -20) |
while (exc >= -20) |
{ |
{ |
p7 = gneg_i(gdiv(p4, poleval(xd,p3))); |
p7 = gneg_i(gdiv(p4, poleval(xd,p3))); |
l3 = avma; |
av3 = avma; |
if (gcmp0(p5)) exc = -32; |
if (gcmp0(p5)) exc = -32; |
else exc = expo(gnorm(p7))-expo(gnorm(p3)); |
else exc = expo(gnorm(p7))-expo(gnorm(p3)); |
avma = l3; |
avma = av3; |
for (j=1; j<=10; j++) |
for (j=1; j<=10; j++) |
{ |
{ |
p8=gadd(p3,p7); p9=poleval(xc,p8); p10=gnorm(p9); |
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]; |
GEN *gptr[3]; |
p3=p8; p4=p9; p5=p10; |
p3=p8; p4=p9; p5=p10; |
gptr[0]=&p3; gptr[1]=&p4; gptr[2]=&p5; |
gptr[0]=&p3; gptr[1]=&p4; gptr[2]=&p5; |
gerepilemanysp(l2,l3,gptr,3); |
gerepilemanysp(av2,av3,gptr,3); |
break; |
break; |
} |
} |
gshiftz(p7,-2,p7); avma=l3; |
gshiftz(p7,-2,p7); avma=av3; |
} |
} |
if (j > 10) |
if (j > 10) |
{ |
{ |
avma=av1; |
avma=av; |
if (DEBUGLEVEL) |
if (DEBUGLEVEL) |
{ |
{ |
fprintferr("too many iterations in rootsold(): "); |
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); |
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) |
for (ln=4; ln<=l; ln=(ln<<1)-2) |
{ |
{ |
setlg(p14,ln); setlg(p15,ln); |
setlg(p14,ln); setlg(p15,ln); |
Line 2384 rootsold(GEN x, long l) |
|
Line 2600 rootsold(GEN x, long l) |
|
p4=poleval(xc,p1); |
p4=poleval(xc,p1); |
p5=poleval(xd,p1); p6=gneg_i(gdiv(p4,p5)); |
p5=poleval(xd,p1); p6=gneg_i(gdiv(p4,p5)); |
settyp(p14,t_REAL); settyp(p15,t_REAL); |
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); |
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))); |
p6=gdiv(p4,poleval(xdabs,gabs(p7,l))); |
if (gexpo(p6)>=expmin) |
if (gexpo(p6)>=expmin) |
{ |
{ |
avma=av1; |
avma=av; |
if (DEBUGLEVEL) |
if (DEBUGLEVEL) |
{ |
{ |
fprintferr("internal error in rootsold(): using roots2()\n"); |
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); |
return roots2(x,l); |
} |
} |
avma=l2; |
avma=av2; |
if (expo(p1[2])<expmin && g) |
if (expo(p1[2])<expmin && g) |
{ |
{ |
gaffect(gzero,(GEN)p1[2]); |
gaffect(gzero,(GEN)p1[2]); |
for (j=1; j<m; j++) gaffect(p1,(GEN)y[k+(i-1)*m+j]); |
for (j=1; j<m; j++) gaffect(p1,(GEN)y[k+(i-1)*m+j]); |
p11[2]=lneg((GEN)p1[1]); |
p11[2]=lneg((GEN)p1[1]); |
l4=avma; xc=gerepile(l0,l4,gdeuc(xc,p11)); |
xc=gerepileupto(av0, gdeuc(xc,p11)); |
} |
} |
else |
else |
{ |
{ |
Line 2428 rootsold(GEN x, long l) |
|
Line 2644 rootsold(GEN x, long l) |
|
p1=gconj(p1); |
p1=gconj(p1); |
for (j=1; j<=m; j++) gaffect(p1,(GEN)y[k+i*m+j]); |
for (j=1; j<=m; j++) gaffect(p1,(GEN)y[k+i*m+j]); |
i++; |
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 |
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; |
k += deg*m; |
} |
} |
Line 2445 rootsold(GEN x, long l) |
|
Line 2661 rootsold(GEN x, long l) |
|
} |
} |
while (k!=deg0); |
while (k!=deg0); |
} |
} |
avma=l1; |
avma=av1; |
for (j=2; j<=deg0; j++) |
for (j=2; j<=deg0; j++) |
{ |
{ |
p1 = (GEN)y[j]; |
p1 = (GEN)y[j]; |
Line 2466 rootsold(GEN x, long l) |
|
Line 2682 rootsold(GEN x, long l) |
|
GEN |
GEN |
roots2(GEN pol,long PREC) |
roots2(GEN pol,long PREC) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
long N,flagexactpol,flagrealpol,flagrealrac,ti,i,j; |
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; |
GEN p1,p2,rr,EPS,qol,qolbis,x,b,c,*ad,v,tabqol; |
|
|
if (typ(pol)!=t_POL) err(typeer,"roots2"); |
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)); |
p2=gneg_i(gdiv((GEN)pol[2],p1)); |
return gerepilecopy(av,p2); |
return gerepilecopy(av,p2); |
} |
} |
EPS=realun(3); setexpo(EPS, 12 - bit_accuracy(PREC)); |
EPS = real2n(12 - bit_accuracy(PREC), 3); |
flagrealpol=1; flagexactpol=1; |
flagrealpol=1; flagexactpol=1; |
for (i=2; i<=N+2; i++) |
for (i=2; i<=N+2; i++) |
{ |
{ |
Line 2611 roots2(GEN pol,long PREC) |
|
Line 2828 roots2(GEN pol,long PREC) |
|
static GEN |
static GEN |
laguer(GEN pol,long N,GEN y0,GEN EPS,long PREC) |
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; |
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; |
av1 = avma; |
I=cgetg(3,t_COMPLEX); I[1]=un; I[2]=un; |
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; |
x=y0; |
for (iter=1; iter<=MAXIT; iter++) |
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); |
d=gzero; f=gzero; abx=QuickNormL1(x,PREC); |
for (j=N-1; j>=0; j--) |
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) |
if (gcmp(QuickNormL1(b,PREC),erre)<=0) |
{ |
{ |
gaffect(x,rac); avma = av1; return rac; |
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 |
static GEN |
balanc(GEN x) |
balanc(GEN x) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
long last,i,j, sqrdx = (RADIX<<1), n = lg(x); |
long last,i,j, sqrdx = (RADIX<<1), n = lg(x); |
GEN r,c,cofgen,a; |
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]); |
anorm=fabs(a[1][1]); |
for (i=2; i<=n; i++) for (j=(i-1); j<=n; j++) anorm+=fabs(a[i][j]); |
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(); } |
if (DEBUGLEVEL>3) { fprintferr("* Finding eigenvalues\n"); flusherr(); } |
while (nn>=1) |
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--) |
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; |
if ((double)(fabs(a[l][l-1])+s)==s) break; |
} |
} |
x=a[nn][nn]; |
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 |
else |
{ |
{ |
y=a[nn-1][nn-1]; |
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) |
if (l == nn-1) |
{ |
{ |
p=0.5*(y-x); q=p*p+w; z=sqrt(fabs(q)); x+=t; |
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; |
z=p+SIGN(z,p); wr[nn-1]=wr[nn]=x+z; |
if (fabs(z)>eps) wr[nn]=x-w/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); } |
else{ wr[nn-1]=wr[nn]=x+p; wi[nn-1]=-(wi[nn]=z); } |
nn-=2; |
nn-=2; |
} |
} |
else |
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==30) err(talker,"too many iterations in hqr"); |
if (its==10 || its==20) |
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])); |
v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1])); |
if ((double)(u+v)==v) break; |
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++) |
for (k=m; k<=nn-1; k++) |
{ |
{ |
if (k!=m) |
if (k!=m) |
{ |
{ |
p=a[k][k-1]; q=a[k+1][k-1]; |
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); |
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); |
s = SIGN(sqrt(p*p+q*q+r*r),p); |
if (s == 0.0) continue; |
if (s == 0.) continue; |
|
|
if (k==m) |
if (k==m) |
{ if (l!=m) a[k][k-1] = -a[k][k-1]; } |
{ 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 */ |
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--) |
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) break; |
if (!flk && !flj && wr[k]<=x) break; |
if (!flk && !flj && wr[k]<=x) 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++) free(a[i]); free(a); eig=cgetg(n+1,t_COL); |
for (i=1; i<=n; i++) |
for (i=1; i<=n; i++) |
{ |
{ |
if (wi[i]) |
if (wi[i] != 0.) |
{ |
{ |
GEN p1 = cgetg(3,t_COMPLEX); |
GEN p1 = cgetg(3,t_COMPLEX); |
eig[i]=(long)p1; |
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 |
GEN |
zrhqr(GEN a,long prec) |
zrhqr(GEN a,long prec) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
long i,j,prec2, n = degpol(a), ex = -bit_accuracy(prec); |
long i,j,prec2, n = degpol(a), ex = -bit_accuracy(prec); |
GEN aa,b,p1,rt,rr,hess,x,dx,y,newval,oldval; |
GEN aa,b,p1,rt,rr,hess,x,dx,y,newval,oldval; |
|
|
hess = cgetg(n+1,t_MAT); |
hess = cgetg(n+1,t_MAT); |
for (j=1; j<=n; j++) |
for (j=1; j<=n; j++) |
{ |
{ |
p1 = cgetg(n+1,t_COL); hess[j] = (long)p1; |
p1 = cgetg(n+1,t_COL); hess[j] = (long)p1; |
p1[1] = lneg(gdiv((GEN)a[n-j+2],(GEN)a[n+2])); |
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; |
for (i=2; i<=n; i++) p1[i] = (i==(j+1))? un: zero; |
} |
} |