Return to base1.c CVS log | Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath |
version 1.1, 2001/10/02 11:17:01 | version 1.2, 2002/09/11 07:26:48 | ||
---|---|---|---|
|
|
||
/**************************************************************/ | /**************************************************************/ | ||
#include "pari.h" | #include "pari.h" | ||
#include "parinf.h" | #include "parinf.h" | ||
GEN idealhermite_aux(GEN nf, GEN x); | |||
int new_galois_format = 0; | |||
extern GEN R_from_QR(GEN x, long prec); | |||
extern GEN to_polmod(GEN x, GEN mod); | |||
extern GEN roots_to_pol_r1r2(GEN a, long r1, long v); | |||
extern GEN idealhermite_aux(GEN nf, GEN x); | |||
extern GEN cauchy_bound(GEN p); | |||
extern GEN galoisbig(GEN x, long prec); | |||
extern GEN lllfp_marked(int M, GEN x, long D, long flag, long prec, int gram); | |||
extern GEN lllint_marked(long M, GEN x, long D, int g, GEN *h, GEN *fl, GEN *B); | |||
extern GEN mulmat_pol(GEN A, GEN x); | |||
extern ulong smulss(ulong x, ulong y, ulong *rem); | |||
void | void | ||
checkrnf(GEN rnf) | checkrnf(GEN rnf) | ||
{ | { | ||
if (typ(rnf)!=t_VEC) err(idealer1); | if (typ(rnf)!=t_VEC || lg(rnf)!=12) err(idealer1); | ||
if (lg(rnf)!=12) err(idealer1); | |||
} | } | ||
GEN | GEN | ||
checkbnf(GEN bnf) | _checkbnf(GEN bnf) | ||
{ | { | ||
if (typ(bnf)!=t_VEC) err(idealer1); | if (typ(bnf) == t_VEC) | ||
switch (lg(bnf)) | switch (lg(bnf)) | ||
{ | |||
case 11: return bnf; | |||
case 7: return checkbnf((GEN)bnf[1]); | |||
case 3: | |||
if (typ(bnf[2])==t_POLMOD) | |||
return checkbnf((GEN)bnf[1]); | |||
} | |||
return NULL; | |||
} | |||
GEN | |||
_checknf(GEN nf) | |||
{ | |||
if (typ(nf)==t_VEC) | |||
switch(lg(nf)) | |||
{ | |||
case 10: return nf; | |||
case 11: return checknf((GEN)nf[7]); | |||
case 7: return checknf((GEN)nf[1]); | |||
case 3: if (typ(nf[2]) == t_POLMOD) return checknf((GEN)nf[1]); | |||
} | |||
return NULL; | |||
} | |||
GEN | |||
checkbnf(GEN x) | |||
{ | |||
GEN bnf = _checkbnf(x); | |||
if (!bnf) | |||
{ | { | ||
case 11: return bnf; | if (_checknf(x)) err(talker,"please apply bnfinit first"); | ||
case 10: | err(idealer1); | ||
if (typ(bnf[1])==t_POL) | |||
err(talker,"please apply bnfinit first"); | |||
break; | |||
case 7: | |||
return checkbnf((GEN)bnf[1]); | |||
case 3: | |||
if (typ(bnf[2])==t_POLMOD) | |||
return checkbnf((GEN)bnf[1]); | |||
} | } | ||
err(idealer1); | return bnf; | ||
return NULL; /* not reached */ | |||
} | } | ||
GEN | GEN | ||
|
|
||
} | } | ||
GEN | GEN | ||
checknf(GEN nf) | checknf(GEN x) | ||
{ | { | ||
if (typ(nf)==t_POL) err(talker,"please apply nfinit first"); | GEN nf = _checknf(x); | ||
if (typ(nf)!=t_VEC) err(idealer1); | if (!nf) | ||
switch(lg(nf)) | |||
{ | { | ||
case 10: return nf; | if (typ(x)==t_POL) err(talker,"please apply nfinit first"); | ||
case 11: return checknf((GEN)nf[7]); | err(idealer1); | ||
case 7: return checknf((GEN)nf[1]); | |||
case 3: if (typ(nf[2]) == t_POLMOD) return checknf((GEN)nf[1]); | |||
} | } | ||
err(idealer1); | return nf; | ||
return NULL; /* not reached */ | |||
} | } | ||
void | void | ||
|
|
||
} | } | ||
void | void | ||
checkprhall(GEN prhall) | |||
{ | |||
if (typ(prhall) != t_VEC || lg(prhall) != 3 || typ(prhall[1]) != t_MAT) | |||
err(talker,"incorrect prhall format"); | |||
} | |||
void | |||
check_pol_int(GEN x, char *s) | check_pol_int(GEN x, char *s) | ||
{ | { | ||
long k = lgef(x)-1; | long k = lgef(x)-1; | ||
|
|
||
GEN | GEN | ||
tschirnhaus(GEN x) | tschirnhaus(GEN x) | ||
{ | { | ||
long a, v = varn(x), av = avma; | gpmem_t av = avma, av2; | ||
long a, v = varn(x); | |||
GEN u, p1 = cgetg(5,t_POL); | GEN u, p1 = cgetg(5,t_POL); | ||
if (typ(x)!=t_POL) err(notpoler,"tschirnhaus"); | if (typ(x)!=t_POL) err(notpoler,"tschirnhaus"); | ||
|
|
||
a = mymyrand() >> randshift; if (a==0) a =1; p1[4]=lstoi(a); | a = mymyrand() >> randshift; if (a==0) a =1; p1[4]=lstoi(a); | ||
a = mymyrand() >> (randshift-1); if (a>=4) a-=8; p1[3]=lstoi(a); | a = mymyrand() >> (randshift-1); if (a>=4) a-=8; p1[3]=lstoi(a); | ||
a = mymyrand() >> (randshift-1); if (a>=4) a-=8; p1[2]=lstoi(a); | a = mymyrand() >> (randshift-1); if (a>=4) a-=8; p1[2]=lstoi(a); | ||
u = caract2(x,p1,v); a=avma; | u = caract2(x,p1,v); av2=avma; | ||
} | } | ||
while (lgef(srgcd(u,derivpol(u))) > 3); /* while u not separable */ | while (lgef(srgcd(u,derivpol(u))) > 3); /* while u not separable */ | ||
if (DEBUGLEVEL>1) | if (DEBUGLEVEL>1) | ||
fprintferr("Tschirnhaus transform. New pol: %Z",u); | fprintferr("Tschirnhaus transform. New pol: %Z",u); | ||
avma=a; return gerepileupto(av,u); | avma=av2; return gerepileupto(av,u); | ||
} | } | ||
#undef randshift | #undef randshift | ||
|
|
||
return 0; | return 0; | ||
} | } | ||
/* assume pol in Z[X] */ | /* assume pol in Z[X]. Find C, L in Z such that POL = C pol(x / L) monic in Z[X]. | ||
* Return POL and set *ptlead = L */ | |||
GEN | GEN | ||
primitive_pol_to_monic(GEN pol, GEN *ptlead) | primitive_pol_to_monic(GEN pol, GEN *ptlead) | ||
{ | { | ||
|
|
||
return p1; | return p1; | ||
} | } | ||
GEN galoisbig(GEN x, long prec); | static GEN | ||
GEN roots_to_pol(GEN a, long v); | roots_to_ZX(GEN z, long *e) | ||
{ | |||
GEN a = roots_to_pol(z,0); | |||
GEN b = grndtoi(greal(a),e); | |||
long e1 = gexpo(gimag(a)); | |||
if (e1 > *e) *e = e1; | |||
return b; | |||
} | |||
static GEN | |||
_res(long n, long s, long k) | |||
{ | |||
GEN y = cgetg(4, t_VEC); | |||
y[1] = lstoi(n); | |||
y[2] = lstoi(s); | |||
if (!new_galois_format) k = (n == 6 && (k == 6 || k == 2))? 2: 1; | |||
y[3] = lstoi(k); return y; | |||
} | |||
GEN | GEN | ||
galois(GEN x, long prec) | galois(GEN x, long prec) | ||
{ | { | ||
long av=avma,av1,i,j,k,n,f,l,l2,e,e1,pr,ind; | gpmem_t av = avma, av1; | ||
GEN x1,p1,p2,p3,p4,p5,p6,w,y,z,ee; | long i,j,k,n,f,l,l2,e,e1,pr,ind; | ||
GEN x1,p1,p2,p3,p4,p5,w,z,ee; | |||
static int ind5[20]={2,5,3,4, 1,3,4,5, 1,5,2,4, 1,2,3,5, 1,4,2,3}; | static int ind5[20]={2,5,3,4, 1,3,4,5, 1,5,2,4, 1,2,3,5, 1,4,2,3}; | ||
static int ind6[60]={3,5,4,6, 2,6,4,5, 2,3,5,6, 2,4,3,6, 2,5,3,4, | static int ind6[60]={3,5,4,6, 2,6,4,5, 2,3,5,6, 2,4,3,6, 2,5,3,4, | ||
1,4,5,6, 1,5,3,6, 1,6,3,4, 1,3,4,5, 1,6,2,5, | 1,4,5,6, 1,5,3,6, 1,6,3,4, 1,3,4,5, 1,6,2,5, | ||
|
|
||
if (typ(x)!=t_POL) err(notpoler,"galois"); | if (typ(x)!=t_POL) err(notpoler,"galois"); | ||
n=degpol(x); if (n<=0) err(constpoler,"galois"); | n=degpol(x); if (n<=0) err(constpoler,"galois"); | ||
if (n>11) err(impl,"galois of degree higher than 11"); | if (n>11) err(impl,"galois of degree higher than 11"); | ||
x = gdiv(x,content(x)); | x = primpart(x); | ||
check_pol_int(x, "galois"); | check_pol_int(x, "galois"); | ||
if (gisirreducible(x) != gun) | if (gisirreducible(x) != gun) | ||
err(impl,"galois of reducible polynomial"); | err(impl,"galois of reducible polynomial"); | ||
if (n<4) | if (n<4) | ||
{ | { | ||
if (n<3) | if (n == 1) { avma = av; return _res(1, 1,1); } | ||
{ | if (n == 2) { avma = av; return _res(2,-1,1); } | ||
avma=av; y=cgetg(4,t_VEC); | /* n = 3 */ | ||
y[1] = (n==1)? un: deux; | f = carreparfait(ZX_disc(x)); | ||
y[2]=lnegi(gun); | avma = av; | ||
} | return f? _res(3,1,1): _res(6,-1,2); | ||
else /* n=3 */ | |||
{ | |||
f=carreparfait(discsr(x)); | |||
avma=av; y=cgetg(4,t_VEC); | |||
if (f) { y[1]=lstoi(3); y[2]=un; } | |||
else { y[1]=lstoi(6); y[2]=lnegi(gun); } | |||
} | |||
y[3]=un; return y; | |||
} | } | ||
x1 = x = primitive_pol_to_monic(x,NULL); av1=avma; | x1 = x = primitive_pol_to_monic(x,NULL); av1=avma; | ||
if (n>7) return galoisbig(x,prec); | if (n > 7) return galoisbig(x,prec); | ||
for(;;) | for(;;) | ||
{ | { | ||
GEN cb = cauchy_bound(x); | |||
switch(n) | switch(n) | ||
{ | { | ||
case 4: z = cgetg(7,t_VEC); | case 4: z = cgetg(7,t_VEC); | ||
prec = DEFAULTPREC + (long)((gexpo(cb)*18. / BITS_IN_LONG)); | |||
for(;;) | for(;;) | ||
{ | { | ||
p1=roots(x,prec); | p1=roots(x,prec); | ||
|
|
||
z[4] = (long)get_F4(transroot(p1,1,4)); | z[4] = (long)get_F4(transroot(p1,1,4)); | ||
z[5] = (long)get_F4(transroot(p1,2,3)); | z[5] = (long)get_F4(transroot(p1,2,3)); | ||
z[6] = (long)get_F4(transroot(p1,3,4)); | z[6] = (long)get_F4(transroot(p1,3,4)); | ||
p4 = roots_to_pol(z,0); | p5 = roots_to_ZX(z, &e); | ||
p5=grndtoi(greal(p4),&e); | |||
e1=gexpo(gimag(p4)); if (e1>e) e=e1; | |||
if (e <= -10) break; | if (e <= -10) break; | ||
prec = (prec<<1)-2; | prec = (prec<<1)-2; | ||
} | } | ||
p6 = modulargcd(p5, derivpol(p5)); | if (!ZX_is_squarefree(p5)) goto tchi; | ||
if (typ(p6)==t_POL && lgef(p6)>3) goto tchi; | |||
p2 = (GEN)factor(p5)[1]; | p2 = (GEN)factor(p5)[1]; | ||
switch(lg(p2)-1) | switch(lg(p2)-1) | ||
{ | { | ||
case 1: f=carreparfait(discsr(x)); avma=av; y=cgetg(4,t_VEC); | case 1: f = carreparfait(ZX_disc(x)); avma = av; | ||
y[3]=un; | return f? _res(12,1,4): _res(24,-1,5); | ||
if (f) { y[2]=un; y[1]=lstoi(12); return y; } | |||
y[2]=lnegi(gun); y[1]=lstoi(24); return y; | |||
case 2: avma=av; y=cgetg(4,t_VEC); | case 2: avma = av; return _res(8,-1,3); | ||
y[3]=un; y[2]=lnegi(gun); y[1]=lstoi(8); return y; | |||
case 3: avma=av; y=cgetg(4,t_VEC); | case 3: avma = av; | ||
y[1]=lstoi(4); y[3]=un; | return (degpol(p2[1])==2)? _res(4,1,2): _res(4,-1,1); | ||
y[2] = (lgef(p2[1])==5)? un: lnegi(gun); | |||
return y; | |||
default: err(bugparier,"galois (bug1)"); | default: err(bugparier,"galois (bug1)"); | ||
} | } | ||
case 5: z = cgetg(7,t_VEC); | case 5: z = cgetg(7,t_VEC); | ||
ee = new_chunk(7); w = new_chunk(7); | ee= cgetg(7,t_VECSMALL); | ||
w = cgetg(7,t_VECSMALL); | |||
prec = DEFAULTPREC + (long)((gexpo(cb)*21. / BITS_IN_LONG)); | |||
for(;;) | for(;;) | ||
{ | { | ||
for(;;) | for(;;) | ||
|
|
||
e1 = gexpo(gimag(p3)); if (e1>e) e=e1; | e1 = gexpo(gimag(p3)); if (e1>e) e=e1; | ||
ee[l]=e; z[l] = (long)p3; | ee[l]=e; z[l] = (long)p3; | ||
} | } | ||
p4 = roots_to_pol(z,0); | p5 = roots_to_ZX(z, &e); | ||
p5=grndtoi(greal(p4),&e); | |||
e1 = gexpo(gimag(p4)); if (e1>e) e=e1; | |||
if (e <= -10) break; | if (e <= -10) break; | ||
prec = (prec<<1)-2; | prec = (prec<<1)-2; | ||
} | } | ||
p6 = modulargcd(p5,derivpol(p5)); | if (!ZX_is_squarefree(p5)) goto tchi; | ||
if (typ(p6)==t_POL && lgef(p6)>3) goto tchi; | |||
p3=(GEN)factor(p5)[1]; | p3=(GEN)factor(p5)[1]; | ||
f=carreparfait(discsr(x)); | f=carreparfait(ZX_disc(x)); | ||
if (lg(p3)-1==1) | if (lg(p3)-1==1) | ||
{ | { | ||
avma=av; y=cgetg(4,t_VEC); y[3]=un; | avma = av; | ||
if (f) { y[2]=un; y[1]=lstoi(60); return y; } | return f? _res(60,1,4): _res(120,-1,5); | ||
else { y[2]=lneg(gun); y[1]=lstoi(120); return y; } | |||
} | } | ||
if (!f) | if (!f) { avma = av; return _res(20,-1,3); } | ||
{ | |||
avma=av; y=cgetg(4,t_VEC); | |||
y[3]=un; y[2]=lneg(gun); y[1]=lstoi(20); return y; | |||
} | |||
pr = - (bit_accuracy(prec) >> 1); | pr = - (bit_accuracy(prec) >> 1); | ||
for (l=1; l<=6; l++) | for (l=1; l<=6; l++) | ||
if (ee[l] <= pr && gcmp0(poleval(p5,(GEN)w[l]))) break; | if (ee[l] <= pr && gcmp0(poleval(p5,(GEN)w[l]))) break; | ||
|
|
||
p3=gzero; | p3=gzero; | ||
for (i=1; i<=5; i++) | for (i=1; i<=5; i++) | ||
{ | { | ||
j=(i%5)+1; | j = i+1; if (j>5) j -= 5; | ||
p3=gadd(p3,gmul(gmul((GEN)p2[i],(GEN)p2[j]), | p3 = gadd(p3,gmul(gmul((GEN)p2[i],(GEN)p2[j]), | ||
gsub((GEN)p2[j],(GEN)p2[i]))); | gsub((GEN)p2[j],(GEN)p2[i]))); | ||
} | } | ||
p5=gsqr(p3); p4=grndtoi(greal(p5),&e); | p5=gsqr(p3); p4=grndtoi(greal(p5),&e); | ||
e1 = gexpo(gimag(p5)); if (e1>e) e=e1; | e1 = gexpo(gimag(p5)); if (e1>e) e=e1; | ||
if (e <= -10) | if (e <= -10) | ||
{ | { | ||
if (gcmp0(p4)) goto tchi; | if (gcmp0(p4)) goto tchi; | ||
f=carreparfait(p4); avma=av; y=cgetg(4,t_VEC); | f = carreparfait(p4); avma = av; | ||
y[3]=y[2]=un; y[1]=lstoi(f?5:10); | return f? _res(5,1,1): _res(10,1,2); | ||
return y; | |||
} | } | ||
prec=(prec<<1)-2; | prec=(prec<<1)-2; | ||
} | } | ||
case 6: z = cgetg(7, t_VEC); | case 6: z = cgetg(7, t_VEC); | ||
prec = DEFAULTPREC + (long)((gexpo(cb)*42. / BITS_IN_LONG)); | |||
for(;;) | for(;;) | ||
{ | { | ||
for(;;) | for(;;) | ||
|
|
||
{ | { | ||
p5=gadd(gmul((GEN)p2[ind6[k]],(GEN)p2[ind6[k+1]]), | p5=gadd(gmul((GEN)p2[ind6[k]],(GEN)p2[ind6[k+1]]), | ||
gmul((GEN)p2[ind6[k+2]],(GEN)p2[ind6[k+3]])); | gmul((GEN)p2[ind6[k+2]],(GEN)p2[ind6[k+3]])); | ||
p3=gadd(p3,gmul(gsqr(gmul((GEN)p2[i],(GEN)p2[j])),p5)); k+=4; | p3=gadd(p3,gmul(gsqr(gmul((GEN)p2[i],(GEN)p2[j])),p5)); | ||
k += 4; | |||
} | } | ||
z[l] = (long)p3; | z[l] = (long)p3; | ||
} | } | ||
p4=roots_to_pol(z,0); | p5 = roots_to_ZX(z, &e); | ||
p5=grndtoi(greal(p4),&e); | |||
e1 = gexpo(gimag(p4)); if (e1>e) e=e1; | |||
if (e <= -10) break; | if (e <= -10) break; | ||
prec=(prec<<1)-2; | prec=(prec<<1)-2; | ||
} | } | ||
p6 = modulargcd(p5,derivpol(p5)); | if (!ZX_is_squarefree(p5)) goto tchi; | ||
if (typ(p6)==t_POL && lgef(p6)>3) goto tchi; | |||
p2=(GEN)factor(p5)[1]; | p2=(GEN)factor(p5)[1]; | ||
switch(lg(p2)-1) | switch(lg(p2)-1) | ||
{ | { | ||
|
|
||
gmul(gmul((GEN)p2[4],(GEN)p2[5]),(GEN)p2[6])); | gmul(gmul((GEN)p2[4],(GEN)p2[5]),(GEN)p2[6])); | ||
z[++ind] = (long)p3; | z[++ind] = (long)p3; | ||
} | } | ||
p4 = roots_to_pol(z,0); | p5 = roots_to_ZX(z, &e); | ||
p5=grndtoi(greal(p4),&e); | |||
e1 = gexpo(gimag(p4)); if (e1>e) e=e1; | |||
if (e <= -10) | if (e <= -10) | ||
{ | { | ||
p6 = modulargcd(p5,derivpol(p5)); | if (!ZX_is_squarefree(p5)) goto tchi; | ||
if (typ(p6)==t_POL && lgef(p6)>3) goto tchi; | p2 = (GEN)factor(p5)[1]; | ||
p2=(GEN)factor(p5)[1]; | f = carreparfait(ZX_disc(x)); | ||
f=carreparfait(discsr(x)); | avma = av; | ||
avma=av; y=cgetg(4,t_VEC); y[3]=un; | |||
if (lg(p2)-1==1) | if (lg(p2)-1==1) | ||
{ | return f? _res(360,1,15): _res(720,-1,16); | ||
if (f) { y[2]=un; y[1]=lstoi(360); } | |||
else { y[2]=lnegi(gun); y[1]=lstoi(720); } | |||
} | |||
else | else | ||
{ | return f? _res(36,1,10): _res(72,-1,13); | ||
if (f) { y[2]=un; y[1]=lstoi(36); } | |||
else { y[2]=lnegi(gun); y[1]=lstoi(72); } | |||
} | |||
return y; | |||
} | } | ||
prec=(prec<<1)-2; break; | prec=(prec<<1)-2; break; | ||
case 2: l2=degpol(p2[1]); if (l2>3) l2=6-l2; | case 2: l2=degpol(p2[1]); if (l2>3) l2=6-l2; | ||
switch(l2) | switch(l2) | ||
{ | { | ||
case 1: f=carreparfait(discsr(x)); | case 1: f = carreparfait(ZX_disc(x)); | ||
avma=av; y=cgetg(4,t_VEC); y[3]=un; | avma = av; | ||
if (f) { y[2]=un; y[1]=lstoi(60); } | return f? _res(60,1,12): _res(120,-1,14); | ||
else { y[2]=lneg(gun); y[1]=lstoi(120); } | case 2: f = carreparfait(ZX_disc(x)); | ||
return y; | if (f) { avma = av; return _res(24,1,7); } | ||
case 2: f=carreparfait(discsr(x)); | p3 = (degpol(p2[1])==2)? (GEN)p2[2]: (GEN)p2[1]; | ||
if (f) | f = carreparfait(ZX_disc(p3)); | ||
{ | avma = av; | ||
avma=av; y=cgetg(4,t_VEC); | return f? _res(24,-1,6): _res(48,-1,11); | ||
y[3]=y[2]=un; y[1]=lstoi(24); | case 3: f = carreparfait(ZX_disc((GEN)p2[1])) | ||
} | || carreparfait(ZX_disc((GEN)p2[2])); | ||
else | avma = av; | ||
{ | return f? _res(18,-1,5): _res(36,-1,9); | ||
p3=(lgef(p2[1])==5) ? (GEN)p2[2]:(GEN)p2[1]; | |||
f=carreparfait(discsr(p3)); | |||
avma=av; y=cgetg(4,t_VEC); y[2]=lneg(gun); | |||
if (f) { y[1]=lstoi(24); y[3]=deux; } | |||
else { y[1]=lstoi(48); y[3]=un; } | |||
} | |||
return y; | |||
case 3: f=carreparfait(discsr((GEN)p2[1])) | |||
|| carreparfait(discsr((GEN)p2[2])); | |||
avma=av; y=cgetg(4,t_VEC); | |||
y[3]=un; y[2]=lneg(gun); y[1]=lstoi(f? 18: 36); | |||
return y; | |||
} | } | ||
case 3: | case 3: | ||
for (l2=1; l2<=3; l2++) | for (l2=1; l2<=3; l2++) | ||
if (lgef(p2[l2])>=6) p3=(GEN)p2[l2]; | if (degpol(p2[l2]) >= 3) p3 = (GEN)p2[l2]; | ||
if (lgef(p3)==6) | if (degpol(p3) == 3) | ||
{ | { | ||
f=carreparfait(discsr(p3)); avma=av; y=cgetg(4,t_VEC); | f = carreparfait(ZX_disc(p3)); avma = av; | ||
y[2]=lneg(gun); y[1]=lstoi(f? 6: 12); | return f? _res(6,-1,1): _res(12,-1,3); | ||
} | } | ||
else | else | ||
{ | { | ||
f=carreparfait(discsr(x)); avma=av; y=cgetg(4,t_VEC); | f = carreparfait(ZX_disc(x)); avma = av; | ||
if (f) { y[2]=un; y[1]=lstoi(12); } | return f? _res(12,1,4): _res(24,-1,8); | ||
else { y[2]=lneg(gun); y[1]=lstoi(24); } | |||
} | } | ||
y[3]=un; return y; | case 4: avma = av; return _res(6,-1,2); | ||
case 4: avma=av; y=cgetg(4,t_VEC); | |||
y[1]=lstoi(6); y[2]=lneg(gun); y[3]=deux; return y; | |||
default: err(bugparier,"galois (bug3)"); | default: err(bugparier,"galois (bug3)"); | ||
} | } | ||
} | } | ||
case 7: z = cgetg(36,t_VEC); | case 7: z = cgetg(36,t_VEC); | ||
prec = DEFAULTPREC + (long)((gexpo(cb)*7. / BITS_IN_LONG)); | |||
for(;;) | for(;;) | ||
{ | { | ||
ind = 0; p1=roots(x,prec); p4=gun; | ind = 0; p1=roots(x,prec); | ||
for (i=1; i<=5; i++) | for (i=1; i<=5; i++) | ||
for (j=i+1; j<=6; j++) | for (j=i+1; j<=6; j++) | ||
{ | { | ||
p6 = gadd((GEN)p1[i],(GEN)p1[j]); | GEN t = gadd((GEN)p1[i],(GEN)p1[j]); | ||
for (k=j+1; k<=7; k++) | for (k=j+1; k<=7; k++) z[++ind] = ladd(t, (GEN)p1[k]); | ||
z[++ind] = (long) gadd(p6,(GEN)p1[k]); | |||
} | } | ||
p4 = roots_to_pol(z,0); | p5 = roots_to_ZX(z, &e); | ||
p5 = grndtoi(greal(p4),&e); | |||
e1 = gexpo(gimag(p4)); if (e1>e) e=e1; | |||
if (e <= -10) break; | if (e <= -10) break; | ||
prec = (prec<<1)-2; | prec = (prec<<1)-2; | ||
} | } | ||
p6 = modulargcd(p5,derivpol(p5)); | if (!ZX_is_squarefree(p5)) goto tchi; | ||
if (typ(p6)==t_POL && lgef(p6)>3) goto tchi; | |||
p2=(GEN)factor(p5)[1]; | p2=(GEN)factor(p5)[1]; | ||
switch(lg(p2)-1) | switch(lg(p2)-1) | ||
{ | { | ||
case 1: f=carreparfait(discsr(x)); avma=av; y=cgetg(4,t_VEC); y[3]=un; | case 1: f = carreparfait(ZX_disc(x)); avma = av; | ||
if (f) { y[2]=un; y[1]=lstoi(2520); } | return f? _res(2520,1,6): _res(5040,-1,7); | ||
else { y[2]=lneg(gun); y[1]=lstoi(5040); } | case 2: f = degpol(p2[1]); avma = av; | ||
return y; | return (f==7 || f==28)? _res(168,1,5): _res(42,-1,4); | ||
case 2: f=degpol(p2[1]); avma=av; y=cgetg(4,t_VEC); y[3]=un; | case 3: avma = av; return _res(21,1,3); | ||
if (f==7 || f==28) { y[2]=un; y[1]=lstoi(168); } | case 4: avma = av; return _res(14,-1,2); | ||
else { y[2]=lneg(gun); y[1]=lstoi(42); } | case 5: avma = av; return _res(7,1,1); | ||
return y; | |||
case 3: avma=av; y=cgetg(4,t_VEC); | |||
y[3]=y[2]=un; y[1]=lstoi(21); return y; | |||
case 4: avma=av; y=cgetg(4,t_VEC); | |||
y[3]=un; y[2]=lneg(gun); y[1]=lstoi(14); return y; | |||
case 5: avma=av; y=cgetg(4,t_VEC); | |||
y[3]=y[2]=un; y[1]=lstoi(7); return y; | |||
default: err(talker,"galois (bug2)"); | default: err(talker,"galois (bug2)"); | ||
} | } | ||
} | } | ||
tchi: avma=av1; x=tschirnhaus(x1); | tchi: avma = av1; x = tschirnhaus(x1); | ||
} | } | ||
} | } | ||
GEN | GEN | ||
galoisapply(GEN nf, GEN aut, GEN x) | galoisapply(GEN nf, GEN aut, GEN x) | ||
{ | { | ||
long av=avma,tetpil,lx,j,N; | gpmem_t av=avma,tetpil; | ||
long lx,j,N; | |||
GEN p,p1,y,pol; | GEN p,p1,y,pol; | ||
nf=checknf(nf); pol=(GEN)nf[1]; | nf=checknf(nf); pol=(GEN)nf[1]; | ||
|
|
||
int cmp_pol(GEN x, GEN y); | int cmp_pol(GEN x, GEN y); | ||
GEN | GEN | ||
get_bnfpol(GEN x, GEN *bnf, GEN *nf) | |||
{ | |||
*bnf = _checkbnf(x); | |||
*nf = _checknf(x); | |||
if (*nf) return (GEN)(*nf)[1]; | |||
if (typ(x) != t_POL) err(idealer1); | |||
return x; | |||
} | |||
GEN | |||
get_nfpol(GEN x, GEN *nf) | get_nfpol(GEN x, GEN *nf) | ||
{ | { | ||
if (typ(x) == t_POL) { *nf = NULL; return x; } | if (typ(x) == t_POL) { *nf = NULL; return x; } | ||
|
|
||
static GEN | static GEN | ||
nfiso0(GEN a, GEN b, long fliso) | nfiso0(GEN a, GEN b, long fliso) | ||
{ | { | ||
ulong av = avma; | gpmem_t av = avma; | ||
long n,m,i,vb,lx; | long n,m,i,vb,lx; | ||
GEN nfa,nfb,p1,y,la,lb; | GEN nfa,nfb,p1,y,la,lb; | ||
|
|
||
} | } | ||
else | else | ||
{ | { | ||
GEN da = nfa? (GEN)nfa[3]: discsr(a); | GEN da = nfa? (GEN)nfa[3]: ZX_disc(a); | ||
GEN db = nfb? (GEN)nfb[3]: discsr(b); | GEN db = nfb? (GEN)nfb[3]: ZX_disc(b); | ||
if (fliso) | if (fliso) | ||
{ | { | ||
p1=gdiv(da,db); | p1=gdiv(da,db); | ||
|
|
||
} | } | ||
y = gen_sort(y, 0, cmp_pol); | y = gen_sort(y, 0, cmp_pol); | ||
settyp(y, t_VEC); | settyp(y, t_VEC); | ||
if (vb == 0) delete_var(); | if (vb == 0) (void)delete_var(); | ||
} | } | ||
lx = lg(y); if (lx==1) { avma=av; return gzero; } | lx = lg(y); if (lx==1) { avma=av; return gzero; } | ||
for (i=1; i<lx; i++) | for (i=1; i<lx; i++) | ||
|
|
||
} | } | ||
GEN | GEN | ||
nfisisom(GEN a, GEN b) | nfisisom(GEN a, GEN b) { return nfiso0(a,b,1); } | ||
{ | |||
return nfiso0(a,b,1); | |||
} | |||
GEN | GEN | ||
nfisincl(GEN a, GEN b) | nfisincl(GEN a, GEN b) { return nfiso0(a,b,0); } | ||
{ | |||
return nfiso0(a,b,0); | |||
} | |||
/*************************************************************************/ | /*************************************************************************/ | ||
/** **/ | /** **/ | ||
/** INITALG **/ | /** INITALG **/ | ||
/** **/ | /** **/ | ||
/*************************************************************************/ | /*************************************************************************/ | ||
extern GEN LLL_nfbasis(GEN *x, GEN polr, GEN base, long prec); | |||
extern GEN eleval(GEN f,GEN h,GEN a); | |||
extern int canon_pol(GEN z); | |||
extern GEN mat_to_vecpol(GEN x, long v); | |||
extern GEN vecpol_to_mat(GEN v, long n); | |||
GEN | |||
get_roots(GEN x,long r1,long prec) | |||
{ | |||
GEN roo = (typ(x)==t_VEC)? dummycopy(x): roots(x,prec); | |||
long i, ru = (lg(roo)-1 + r1) >> 1; | |||
for (i=1; i<=r1; i++) roo[i]=lreal((GEN)roo[i]); | |||
for ( ; i<=ru; i++) roo[i]=roo[(i<<1)-r1]; | |||
roo[0]=evaltyp(t_VEC)|evallg(ru+1); return roo; | |||
} | |||
/* For internal use. compute trace(x mod pol), sym=polsym(pol,deg(pol)-1) */ | /* For internal use. compute trace(x mod pol), sym=polsym(pol,deg(pol)-1) */ | ||
GEN | GEN | ||
quicktrace(GEN x, GEN sym) | quicktrace(GEN x, GEN sym) | ||
|
|
||
static GEN | static GEN | ||
trace_col(GEN x, GEN T) | trace_col(GEN x, GEN T) | ||
{ | { | ||
ulong av = avma; | gpmem_t av = avma; | ||
GEN t = gzero; | GEN t = gzero; | ||
long i, l = lg(x); | long i, l = lg(x); | ||
|
|
||
return gerepileuptoint(av, t); | return gerepileuptoint(av, t); | ||
} | } | ||
/* Seek a new, simpler, polynomial pol defining the same number field as | |||
* x (assumed to be monic at this point). | |||
* *ptx receives pol | |||
* *ptdx receives disc(pol) | |||
* *ptrev expresses the new root in terms of the old one. | |||
* base updated in place. | |||
* prec = 0 <==> field is totally real | |||
*/ | |||
static void | |||
nfinit_reduce(long flag, GEN *px, GEN *pdx, GEN *prev, GEN *pbase, long prec) | |||
{ | |||
GEN chbas,a,phimax,dxn,s,sn,p1,p2,p3,polmax,rev,polr; | |||
GEN x = *px, dx = *pdx, base = *pbase; | |||
long i,j,nmax,numb,flc,v=varn(x), n=lg(base)-1; | |||
if (n == 1) | |||
{ | |||
*px = gsub(polx[v],gun); *pdx = gun; | |||
*prev = polymodrecip(gmodulcp(gun, x)); return; | |||
} | |||
polr = prec? roots(x, prec): NULL; | |||
if (polr) | |||
for (s=gzero,i=1; i<=n; i++) s = gadd(s,gnorm((GEN)polr[i])); | |||
else | |||
s = subii(sqri((GEN)x[n+1]), shifti((GEN)x[n],1)); | |||
chbas = LLL_nfbasis(&x,polr,base,prec); | |||
if (DEBUGLEVEL) msgtimer("LLL basis"); | |||
phimax=polx[v]; polmax=dummycopy(x); | |||
nmax=(flag & nf_PARTIAL)?min(n,3):n; | |||
for (numb=0,i=2; i<=nmax || (!numb && i<=n); i++) | |||
{ /* cf pols_for_polred */ | |||
if (DEBUGLEVEL>2) { fprintferr("i = %ld\n",i); flusherr(); } | |||
p1 = a = gmul(base,(GEN)chbas[i]); p3=content(p1); | |||
if (gcmp1(p3)) p3 = NULL; else p1 = gdiv(p1,p3); | |||
p1 = caract2(x,p1,v); | |||
if (p3) | |||
for (p2=gun, j=lgef(p1)-2; j>=2; j--) | |||
{ | |||
p2 = gmul(p2,p3); p1[j] = lmul((GEN)p1[j], p2); | |||
} | |||
p2 = modulargcd(derivpol(p1),p1); | |||
if (lgef(p2) > 3) continue; | |||
if (DEBUGLEVEL>3) outerr(p1); | |||
dxn=discsr(p1); flc=absi_cmp(dxn,dx); numb++; | |||
if (flc>0) continue; | |||
if (polr) | |||
for (sn=gzero,j=1; j<=n; j++) | |||
sn = gadd(sn,gnorm(poleval(a,(GEN)polr[j]))); | |||
else | |||
sn = subii(sqri((GEN)p1[n+1]), shifti((GEN)p1[n],1)); | |||
if (flc>=0) | |||
{ | |||
flc=gcmp(sn,s); | |||
if (flc>0 || (!flc && gpolcomp(p1,polmax) >= 0)) continue; | |||
} | |||
dx=dxn; s=sn; polmax=p1; phimax=a; | |||
} | |||
if (!numb) | |||
{ | |||
if (gisirreducible(x)!=gun) err(redpoler,"nfinit_reduce"); | |||
err(talker,"you have found a counter-example to a conjecture, " | |||
"please send us\nthe polynomial as soon as possible"); | |||
} | |||
if (phimax == polx[v]) /* no improvement */ | |||
rev = gmodulcp(phimax,x); | |||
else | |||
{ | |||
if (canon_pol(polmax) < 0) phimax = gneg_i(phimax); | |||
if (DEBUGLEVEL>1) fprintferr("polmax = %Z\n",polmax); | |||
p2 = gmodulcp(phimax,x); rev = polymodrecip(p2); | |||
a = base; base = cgetg(n+1,t_VEC); | |||
p1 = (GEN)rev[2]; | |||
for (i=1; i<=n; i++) | |||
base[i] = (long)eleval(polmax, (GEN)a[i], p1); | |||
p1 = vecpol_to_mat(base,n); p2=denom(p1); p1=gmul(p2,p1); | |||
p1 = gdiv(hnfmodid(p1,p2), p2); | |||
base = mat_to_vecpol(p1,v); | |||
} | |||
*px=polmax; *pdx=dx; *prev=rev; *pbase = base; | |||
} | |||
/* pol belonging to Z[x], return a monic polynomial generating the same field | /* pol belonging to Z[x], return a monic polynomial generating the same field | ||
* as pol (x-> ax+b)) set lead = NULL if pol was monic (after dividing | * as pol (x-> ax+b)) set lead = NULL if pol was monic (after dividing | ||
* by the content), and to to leading coeff otherwise. | * by the content), and to to leading coeff otherwise. | ||
|
|
||
pol_to_monic(GEN pol, GEN *lead) | pol_to_monic(GEN pol, GEN *lead) | ||
{ | { | ||
long n = lgef(pol)-1; | long n = lgef(pol)-1; | ||
GEN p1; | |||
if (n==1 || gcmp1((GEN)pol[n])) { *lead = NULL; return pol; } | if (n==1 || gcmp1((GEN)pol[n])) { *lead = NULL; return pol; } | ||
return primitive_pol_to_monic(primpart(pol), lead); | |||
p1=content(pol); if (!gcmp1(p1)) pol = gdiv(pol,p1); | |||
return primitive_pol_to_monic(pol,lead); | |||
} | } | ||
/* NB: TI integral, det(TI) = d_K, ideal/dideal = codifferent */ | /* Let T = (Tr(wi wj)), then T^(-1) Z^n = codifferent = coD/dcoD | ||
GEN | * NB: det(T) = dK, TI = dK T^(-1) integral */ | ||
make_MDI(GEN nf, GEN TI, GEN *ideal, GEN *dideal) | static GEN | ||
make_MDI(GEN nf, GEN TI, GEN *coD, GEN *dcoD) | |||
{ | { | ||
GEN c = content(TI); | GEN c, MDI, dK = (GEN)nf[3]; | ||
GEN p1; | |||
*dideal = divii((GEN)nf[3], c); | TI = Q_primitive_part(TI, &c); | ||
*ideal = p1 = hnfmodid(gdiv(TI,c), *dideal); | *dcoD = c? diviiexact(dK, c): dK; | ||
return gmul(c, ideal_two_elt(nf, p1)); | *coD = hnfmodid(TI, *dcoD); | ||
MDI = ideal_two_elt(nf, *coD); | |||
return c? gmul(c, MDI): MDI; | |||
} | } | ||
/* [bas[i]/den[i]]= integer basis. roo = real part of the roots */ | |||
GEN | |||
make_M(GEN basden,GEN roo) | |||
{ | |||
GEN p1,d,M, bas=(GEN)basden[1], den=(GEN)basden[2]; | |||
long i,j, ru = lg(roo), n = lg(bas); | |||
M = cgetg(n,t_MAT); | |||
for (j=1; j<n; j++) | |||
{ | |||
p1=cgetg(ru,t_COL); M[j]=(long)p1; | |||
for (i=1; i<ru; i++) | |||
p1[i]=(long)poleval((GEN)bas[j],(GEN)roo[i]); | |||
} | |||
if (den) | |||
{ | |||
long prec = precision((GEN)roo[1]); | |||
GEN invd,rd = cgetr(prec+1); | |||
for (j=1; j<n; j++) | |||
{ | |||
d = (GEN)den[j]; if (!d) continue; | |||
p1 = (GEN)M[j]; affir(d,rd); invd = ginv(rd); | |||
for (i=1; i<ru; i++) p1[i] = lmul((GEN)p1[i], invd); | |||
} | |||
} | |||
if (DEBUGLEVEL>4) msgtimer("matrix M"); | |||
return M; | |||
} | |||
GEN | |||
make_MC(long r1,GEN M) | |||
{ | |||
long i,j,av,tetpil, n = lg(M), ru = lg(M[1]); | |||
GEN p1,p2,MC=cgetg(ru,t_MAT); | |||
for (j=1; j<ru; j++) | |||
{ | |||
p1=cgetg(n,t_COL); MC[j]=(long)p1; | |||
for (i=1; i<n; i++) | |||
{ | |||
av=avma; p2=gconj(gcoeff(M,j,i)); tetpil=avma; | |||
p1[i] = (j<=r1)? (long)p2: lpile(av,tetpil,gmul2n(p2,1)); | |||
} | |||
} | |||
if (DEBUGLEVEL>4) msgtimer("matrix MC"); | |||
return MC; | |||
} | |||
GEN | |||
get_roots(GEN x,long r1,long ru,long prec) | |||
{ | |||
GEN roo = (typ(x)==t_VEC)? dummycopy(x): roots(x,prec); | |||
long i; | |||
for (i=1; i<=r1; i++) roo[i]=lreal((GEN)roo[i]); | |||
for ( ; i<=ru; i++) roo[i]=roo[(i<<1)-r1]; | |||
roo[0]=evaltyp(t_VEC)|evallg(ru+1); return roo; | |||
} | |||
/* assume lg(nf) > 3 && typ(nf) = container [hopefully a genuine nf] */ | /* assume lg(nf) > 3 && typ(nf) = container [hopefully a genuine nf] */ | ||
long | long | ||
nf_get_r1(GEN nf) | nf_get_r1(GEN nf) | ||
|
|
||
err(talker,"false nf in nf_get_r2"); | err(talker,"false nf in nf_get_r2"); | ||
return itos((GEN)x[2]); | return itos((GEN)x[2]); | ||
} | } | ||
void | |||
nf_get_sign(GEN nf, long *r1, long *r2) | |||
{ | |||
GEN x = (GEN)nf[2]; | |||
if (typ(x) != t_VEC || lg(x) != 3 | |||
|| typ(x[1]) != t_INT || typ(x[2]) != t_INT) | |||
err(talker,"false nf in nf_get_sign"); | |||
*r1 = itos((GEN)x[1]); | |||
*r2 = (degpol(nf[1]) - *r1) >> 1; | |||
} | |||
extern GEN mulmat_pol(GEN A, GEN x); | |||
static GEN | static GEN | ||
get_T(GEN mul, GEN x, GEN bas, GEN den) | get_Tr(GEN mul, GEN x, GEN basden) | ||
{ | { | ||
GEN tr,T,T1,sym, bas = (GEN)basden[1], den = (GEN)basden[2]; | |||
long i,j,n = lg(bas)-1; | long i,j,n = lg(bas)-1; | ||
GEN tr,T,T1,sym; | |||
T = cgetg(n+1,t_MAT); | T = cgetg(n+1,t_MAT); | ||
T1 = cgetg(n+1,t_COL); | T1 = cgetg(n+1,t_COL); | ||
sym = polsym(x, n-1); | sym = polsym(x, n-1); | ||
|
|
||
GEN | GEN | ||
get_bas_den(GEN bas) | get_bas_den(GEN bas) | ||
{ | { | ||
GEN z,d,den, dbas = dummycopy(bas); | GEN z,b,d,den, dbas = dummycopy(bas); | ||
long i, c = 0, l = lg(bas); | long i, l = lg(bas); | ||
int power = 1; | |||
den = cgetg(l,t_VEC); | den = cgetg(l,t_VEC); | ||
for (i=1; i<l; i++) | for (i=1; i<l; i++) | ||
{ | { | ||
d = denom(content((GEN)dbas[i])); | b = Q_remove_denom((GEN)bas[i], &d); | ||
if (is_pm1(d)) d = NULL; else { dbas[i] = lmul((GEN)dbas[i],d); c++; } | dbas[i]= (long)b; | ||
den[i] = (long)d; | den[i] = (long)d; if (d) power = 0; | ||
} | } | ||
if (!c) den = NULL; /* power basis */ | if (power) den = NULL; /* power basis */ | ||
z = cgetg(3,t_VEC); | z = cgetg(3,t_VEC); | ||
z[1] = (long)dbas; | z[1] = (long)dbas; | ||
z[2] = (long)den; return z; | z[2] = (long)den; return z; | ||
|
|
||
} | } | ||
GEN | GEN | ||
get_mul_table(GEN x,GEN basden,GEN invbas,GEN *T) | get_mul_table(GEN x,GEN basden,GEN invbas) | ||
{ | { | ||
long i,j, n = degpol(x); | long i,j, n = degpol(x); | ||
GEN z,d, mul = cgetg(n*n+1,t_MAT), bas=(GEN)basden[1], den=(GEN)basden[2]; | GEN z,d, mul = cgetg(n*n+1,t_MAT), bas=(GEN)basden[1], den=(GEN)basden[2]; | ||
for (j=1; j<=n*n; j++) mul[j]=lgetg(n+1,t_COL); | for (j=1; j<=n*n; j++) mul[j]=lgetg(n+1,t_COL); | ||
for (i=1; i<=n; i++) | for (i=1; i<=n; i++) | ||
for (j=i; j<=n; j++) | for (j=i; j<=n; j++) | ||
{ | { | ||
ulong av = avma; | gpmem_t av = avma; | ||
z = gres(gmul((GEN)bas[j],(GEN)bas[i]), x); | z = gres(gmul((GEN)bas[j],(GEN)bas[i]), x); | ||
z = mulmat_pol(invbas, z); /* integral column */ | z = mulmat_pol(invbas, z); /* integral column */ | ||
if (den) | if (den) | ||
|
|
||
} | } | ||
mul[j+(i-1)*n] = mul[i+(j-1)*n] = lpileupto(av,z); | mul[j+(i-1)*n] = mul[i+(j-1)*n] = lpileupto(av,z); | ||
} | } | ||
if (T) *T = get_T(mul,x,bas,den); | |||
return mul; | return mul; | ||
} | } | ||
/* fill mat = nf[5], as well as nf[8] and nf[9] | /* as get_Tr, mul_table not precomputed */ | ||
* If (small) only compute a subset (use dummy 0s for the rest) */ | static GEN | ||
make_Tr(GEN x, GEN w) | |||
{ | |||
long i,j, n = degpol(x); | |||
GEN p1,p2,t,d; | |||
GEN sym = cgetg(n+2,t_VEC); | |||
GEN den = cgetg(n+1,t_VEC); | |||
GEN T = cgetg(n+1,t_MAT); | |||
gpmem_t av; | |||
sym = polsym(x, n-1); | |||
p1 = get_bas_den(w); | |||
w = (GEN)p1[1]; | |||
den = (GEN)p1[2]; | |||
for (i=1; i<=n; i++) | |||
{ | |||
p1 = cgetg(n+1,t_COL); T[i] = (long)p1; | |||
for (j=1; j<i ; j++) p1[j] = coeff(T,i,j); | |||
for ( ; j<=n; j++) | |||
{ | |||
av = avma; | |||
p2 = gres(gmul((GEN)w[i],(GEN)w[j]), x); | |||
t = quicktrace(p2, sym); | |||
if (den) | |||
{ | |||
d = _mulii((GEN)den[i],(GEN)den[j]); | |||
if (d) t = diviiexact(t, d); | |||
} | |||
p1[j] = lpileuptoint(av, t); | |||
} | |||
} | |||
return T; | |||
} | |||
/* compute roots so that _absolute_ accuracy of M >= prec [also holds for G] */ | |||
static void | |||
get_roots_for_M(nffp_t *F) | |||
{ | |||
long n, eBD, er, PREC; | |||
if (F->extraprec < 0) | |||
{ /* not initialized yet */ | |||
n = degpol(F->x); | |||
eBD = 1 + gexpo((GEN)F->basden[1]); | |||
er = 1 + gexpo(F->ro? F->ro: cauchy_bound(F->x)); | |||
if (er < 0) er = 0; | |||
F->extraprec = ((n*er + eBD + (long)log2(n)) >> TWOPOTBITS_IN_LONG); | |||
} | |||
PREC = F->prec + F->extraprec; | |||
if (F->ro && gprecision((GEN)F->ro[1]) >= PREC) return; | |||
F->ro = get_roots(F->x, F->r1, PREC); | |||
} | |||
/* [bas[i]/den[i]]= integer basis. roo = real part of the roots */ | |||
static void | |||
make_M(nffp_t *F, int trunc) | |||
{ | |||
GEN bas = (GEN)F->basden[1], den = (GEN)F->basden[2], ro = F->ro; | |||
GEN m, d, M; | |||
long i, j, l = lg(ro), n = lg(bas); | |||
M = cgetg(n,t_MAT); | |||
m = cgetg(l, t_COL); M[1] = (long)m; | |||
for (i=1; i<l; i++) m[i] = un; /* bas[1] = 1 */ | |||
for (j=2; j<n; j++) | |||
{ | |||
m = cgetg(l,t_COL); M[j] = (long)m; | |||
for (i=1; i<l; i++) m[i] = (long)poleval((GEN)bas[j], (GEN)ro[i]); | |||
} | |||
if (den) | |||
{ | |||
GEN invd, rd = cgetr(F->prec + F->extraprec); | |||
for (j=2; j<n; j++) | |||
{ | |||
d = (GEN)den[j]; if (!d) continue; | |||
m = (GEN)M[j]; affir(d,rd); invd = ginv(rd); | |||
for (i=1; i<l; i++) m[i] = lmul((GEN)m[i], invd); | |||
} | |||
} | |||
if (trunc && gprecision(M) > F->prec) | |||
{ | |||
M = gprec_w(M, F->prec); | |||
F->ro = gprec_w(ro,F->prec); | |||
} | |||
if (DEBUGLEVEL>4) msgtimer("matrix M"); | |||
F->M = M; | |||
} | |||
/* return G real such that G~ * G = T_2 */ | |||
static void | |||
make_G(nffp_t *F) | |||
{ | |||
GEN G, m, g, r, M = F->M, sqrt2 = gsqrt(gdeux, F->prec + F->extraprec); | |||
long i, j, k, r1 = F->r1, l = lg(M); | |||
G = cgetg(l, t_MAT); | |||
for (j=1; j<l; j++) | |||
{ | |||
g = cgetg(l, t_COL); | |||
G[j] = (long)g; m = (GEN)M[j]; | |||
for (k=i=1; i<=r1; i++) g[k++] = m[i]; | |||
for ( ; k < l; i++) | |||
{ | |||
r = (GEN)m[i]; | |||
if (typ(r) == t_COMPLEX) | |||
{ | |||
g[k++] = lmpmul(sqrt2, (GEN)r[1]); | |||
g[k++] = lmpmul(sqrt2, (GEN)r[2]); | |||
} | |||
else | |||
{ | |||
g[k++] = lmpmul(sqrt2, r); | |||
g[k++] = zero; | |||
} | |||
} | |||
} | |||
F->G = G; | |||
} | |||
static void | |||
make_M_G(nffp_t *F, int trunc) | |||
{ | |||
get_roots_for_M(F); | |||
make_M(F, trunc); | |||
make_G(F); | |||
} | |||
void | void | ||
get_nf_matrices(GEN nf, long small) | remake_GM(GEN nf, nffp_t *F, long prec) | ||
{ | { | ||
GEN x=(GEN)nf[1],dK=(GEN)nf[3],index=(GEN)nf[4],ro=(GEN)nf[6],bas=(GEN)nf[7]; | F->x = (GEN)nf[1]; | ||
GEN basden,mul,invbas,M,MC,T,MDI,D,TI,A,dA,mat; | F->ro = NULL; | ||
long r1 = nf_get_r1(nf), n = lg(bas)-1; | F->r1 = nf_get_r1(nf); | ||
F->basden = get_bas_den((GEN)nf[7]); | |||
F->extraprec = -1; | |||
F->prec = prec; make_M_G(F, 1); | |||
} | |||
mat = cgetg(small? 4: 8,t_VEC); nf[5] = (long)mat; | static void | ||
basden = get_bas_den(bas); | nffp_init(nffp_t *F, nfbasic_t *T, GEN ro, long prec) | ||
M = make_M(basden,ro); | { | ||
MC = make_MC(r1,M); | F->x = T->x; | ||
mat[1]=(long)M; | F->ro = ro; | ||
mat[3]=(long)mulmat_real(MC,M); | F->r1 = T->r1; | ||
if (small) { nf[8]=nf[9]=mat[2]=zero; return; } | if (!T->basden) T->basden = get_bas_den(T->bas); | ||
F->basden = T->basden; | |||
F->extraprec = -1; | |||
F->prec = prec; | |||
} | |||
invbas = QM_inv(vecpol_to_mat(bas,n), gun); | static void | ||
mul = get_mul_table(x,basden,invbas,&T); | get_nf_fp_compo(nfbasic_t *T, nffp_t *F, GEN ro, long prec) | ||
{ | |||
nffp_init(F,T,ro,prec); | |||
make_M_G(F, 1); | |||
} | |||
static GEN | |||
get_sign(long r1, long n) | |||
{ | |||
GEN s = cgetg(3, t_VEC); | |||
s[1] = lstoi(r1); | |||
s[2] = lstoi((n - r1) >> 1); return s; | |||
} | |||
GEN | |||
nfbasic_to_nf(nfbasic_t *T, GEN ro, long prec) | |||
{ | |||
GEN nf = cgetg(10,t_VEC); | |||
GEN x = T->x; | |||
GEN invbas,Tr,D,TI,A,dA, mat = cgetg(8,t_VEC); | |||
nffp_t F; | |||
get_nf_fp_compo(T, &F, ro, prec); | |||
nf[1] = (long)T->x; | |||
nf[2] = (long)get_sign(T->r1, degpol(T->x)); | |||
nf[3] = (long)T->dK; | |||
nf[4] = (long)T->index; | |||
nf[6] = (long)F.ro; | |||
nf[5] = (long)mat; | |||
nf[7] = (long)T->bas; | |||
mat[1] = (long)F.M; | |||
mat[2] = (long)F.G; | |||
invbas = QM_inv(vecpol_to_mat(T->bas, lg(T->bas)-1), gun); | |||
nf[8] = (long)invbas; | |||
nf[9] = (long)get_mul_table(x, F.basden, invbas); | |||
if (DEBUGLEVEL) msgtimer("mult. table"); | if (DEBUGLEVEL) msgtimer("mult. table"); | ||
nf[8]=(long)invbas; | |||
nf[9]=(long)mul; | |||
TI = ZM_inv(T, dK); /* dK T^-1 */ | Tr = get_Tr((GEN)nf[9], x, F.basden); | ||
MDI = make_MDI(nf,TI, &A, &dA); | TI = ZM_inv(Tr, T->dK); /* dK T^-1 */ | ||
mat[6]=(long)TI; | mat[6] = (long)TI; | ||
mat[7]=(long)MDI; /* needed in idealinv below */ | mat[7] = (long)make_MDI(nf,TI, &A, &dA); /* needed in idealinv below */ | ||
if (is_pm1(index)) | if (is_pm1(T->index)) | ||
D = idealhermite_aux(nf, derivpol(x)); | D = idealhermite_aux(nf, derivpol(x)); | ||
else | else | ||
D = gmul(dA, idealinv(nf, A)); | D = gmul(dA, idealinv(nf, A)); | ||
mat[2]=(long)MC; | mat[3] = zero; /* FIXME: was gram matrix of current mat[2]. Useless */ | ||
mat[4]=(long)T; | mat[4] = (long)Tr; | ||
mat[5]=(long)D; | mat[5] = (long)D; | ||
if (DEBUGLEVEL) msgtimer("matrices"); | if (DEBUGLEVEL) msgtimer("matrices"); | ||
return nf; | |||
} | } | ||
/* Initialize the number field defined by the polynomial x (in variable v) | static GEN | ||
* flag & nf_REGULAR | hnffromLLL(GEN nf) | ||
* regular behaviour. | { | ||
* flag & nf_SMALL | GEN d, x; | ||
* compute only nf[1] (pol), nf[2] (signature), nf[5][3] (T2) and | x = vecpol_to_mat((GEN)nf[7], degpol(nf[1])); | ||
* nf[7] (integer basis), the other components are filled with gzero. | x = Q_remove_denom(x, &d); | ||
* flag & nf_REDUCE | if (!d) return x; /* power basis */ | ||
* try a polred first. | return gauss(hnfmodid(x, d), x); | ||
* flag & nf_PARTIAL | } | ||
* do a partial polred, not a polredabs | |||
* flag & nf_ORIG | |||
* do a polred and return [nfinit(x),Mod(a,red)], where | |||
* Mod(a,red)=Mod(v,x) (i.e return the base change). | |||
*/ | |||
/* here x can be a polynomial, an nf or a bnf */ | static GEN | ||
nfbasechange(GEN u, GEN x) | |||
{ | |||
long i,lx; | |||
GEN y; | |||
switch(typ(x)) | |||
{ | |||
case t_COL: /* nfelt */ | |||
return gmul(u, x); | |||
case t_MAT: /* ideal */ | |||
y = dummycopy(x); | |||
lx = lg(x); | |||
for (i=1; i<lx; i++) y[i] = lmul(u, (GEN)y[i]); | |||
break; | |||
case t_VEC: /* pr */ | |||
checkprimeid(x); | |||
y = dummycopy(x); | |||
y[2] = lmul(u, (GEN)y[2]); | |||
y[5] = lmul(u, (GEN)y[5]); | |||
break; | |||
default: y = x; | |||
} | |||
return y; | |||
} | |||
GEN | GEN | ||
initalgall0(GEN x, long flag, long prec) | nffromhnfbasis(GEN nf, GEN x) | ||
{ | { | ||
GEN lead = NULL,nf,ro,bas,mat,rev,dK,dx,index,res; | long tx = typ(x); | ||
long av=avma,n,i,r1,r2,ru,PRECREG; | gpmem_t av = avma; | ||
GEN u; | |||
if (!is_vec_t(tx)) return gcopy(x); | |||
nf = checknf(nf); | |||
u = hnffromLLL(nf); | |||
return gerepilecopy(av, nfbasechange(u, x)); | |||
} | |||
if (DEBUGLEVEL) timer2(); | GEN | ||
if (typ(x)==t_POL) | nftohnfbasis(GEN nf, GEN x) | ||
{ | |||
long tx = typ(x); | |||
gpmem_t av = avma; | |||
GEN u; | |||
if (!is_vec_t(tx)) return gcopy(x); | |||
nf = checknf(nf); | |||
u = ZM_inv(hnffromLLL(nf), gun); | |||
return gerepilecopy(av, nfbasechange(u, x)); | |||
} | |||
static GEN | |||
get_red_G(nfbasic_t *T, GEN *pro) | |||
{ | |||
GEN G, u, u0 = NULL; | |||
gpmem_t av; | |||
long i, prec, extraprec, n = degpol(T->x); | |||
nffp_t F; | |||
extraprec = (long) (0.5 * n * sizeof(long) / 8.); | |||
prec = DEFAULTPREC + extraprec; | |||
nffp_init(&F, T, *pro, prec); | |||
av = avma; | |||
for (i=1; ; i++) | |||
{ | { | ||
n=degpol(x); if (n<=0) err(constpoler,"initalgall0"); | F.prec = prec; make_M_G(&F, 0); G = F.G; | ||
check_pol_int(x,"initalg"); | if (u0) G = gmul(G, u0); | ||
if (gisirreducible(x) == gzero) err(redpoler,"nfinit"); | if (DEBUGLEVEL) | ||
if (!gcmp1((GEN)x[n+2])) | fprintferr("get_red_G: starting LLL, prec = %ld (%ld + %ld)\n", | ||
prec + F.extraprec, prec, F.extraprec); | |||
if ((u = lllfp_marked(1, G, 100, 2, prec, 0))) | |||
{ | { | ||
x = pol_to_monic(x,&lead); | if (typ(u) == t_MAT) break; | ||
if (!(flag & nf_SMALL)) | u = (GEN)u[1]; | ||
{ | if (u0) u0 = gerepileupto(av, gmul(u0,u)); | ||
if (!(flag & nf_REDUCE)) | else u0 = gerepilecopy(av, u); | ||
err(warner,"non-monic polynomial. Result of the form [nf,c]"); | |||
flag = flag | nf_REDUCE | nf_ORIG; | |||
} | |||
} | } | ||
bas = allbase4(x,0,&dK,NULL); | if (i == MAXITERPOL) err(accurer,"red_T2"); | ||
if (DEBUGLEVEL) msgtimer("round4"); | prec = (prec<<1)-2 + (gexpo(u0) >> TWOPOTBITS_IN_LONG); | ||
dx = discsr(x); r1 = sturm(x); | F.ro = NULL; | ||
if (DEBUGLEVEL) err(warnprec,"get_red_G", prec); | |||
} | } | ||
*pro = F.ro; return u0? gmul(u0,u): u; | |||
} | |||
/* Return the base change matrix giving an LLL-reduced basis for the | |||
* integer basis of nf(T). | |||
* *ro = roots of x, computed to precision prec [or NULL] */ | |||
static GEN | |||
get_LLL_basis(nfbasic_t *T, GEN *pro) | |||
{ | |||
GEN u; | |||
if (T->r1 != degpol(T->x)) u = get_red_G(T, pro); | |||
else | else | ||
{ | { | ||
i = lg(x); | u = lllint_marked(1, make_Tr(T->x, T->bas), 100, 1, &u,NULL,NULL); | ||
if (typ(x) == t_VEC && i<=4 && i>=3 && typ(x[1])==t_POL) | if (!u) u = idmat(1); | ||
{ /* polynomial + integer basis */ | |||
bas=(GEN)x[2]; x=(GEN)x[1]; n=degpol(x); | |||
if (typ(bas) == t_MAT) { mat = bas; bas = mat_to_vecpol(mat,varn(x)); } | |||
else mat = vecpol_to_mat(bas,n); | |||
dx = discsr(x); r1 = sturm(x); | |||
dK = gmul(dx, gsqr(det2(mat))); | |||
} | |||
else | |||
{ | |||
nf=checknf(x); | |||
bas=(GEN)nf[7]; x=(GEN)nf[1]; n=degpol(x); | |||
dK=(GEN)nf[3]; dx=mulii(dK, sqri((GEN)nf[4])); | |||
r1 = nf_get_r1(nf); | |||
} | |||
bas[1]=lpolun[varn(x)]; /* it may be gun => SEGV later */ | |||
} | } | ||
r2=(n-r1)>>1; ru=r1+r2; | return u; | ||
PRECREG = prec + (expi(dK)>>(TWOPOTBITS_IN_LONG+1)) | } | ||
+ (long)((sqrt((double)n)+3) / sizeof(long) * 4); | |||
rev = NULL; | static void | ||
if (flag & nf_REDUCE) | set_LLL_basis(nfbasic_t *T, GEN *pro) | ||
{ | |||
T->bas = gmul(T->bas, get_LLL_basis(T, pro)); | |||
T->basden = NULL; /* recompute */ | |||
} | |||
typedef struct { | |||
GEN xbest, dxbest; | |||
long ind, indmax, indbest; | |||
} ok_pol_t; | |||
/* is xn better than x ? */ | |||
static int | |||
better_pol(GEN xn, GEN dxn, GEN x, GEN dx) | |||
{ | |||
int fl; | |||
if (!x) return 1; | |||
fl = absi_cmp(dxn, dx); | |||
return (fl < 0 || (!fl && gpolcomp(xn, x) < 0)); | |||
} | |||
static GEN | |||
ok_pol(void *TT, GEN xn) | |||
{ | |||
ok_pol_t *T = (ok_pol_t*)TT; | |||
GEN dxn; | |||
if (++T->ind > T->indmax && T->xbest) return T->xbest; | |||
if (!ZX_is_squarefree(xn)) return (T->ind == T->indmax)? T->xbest: NULL; | |||
if (DEBUGLEVEL>3) outerr(xn); | |||
dxn = ZX_disc(xn); | |||
if (better_pol(xn, dxn, T->xbest, T->dxbest)) | |||
{ | { | ||
nfinit_reduce(flag, &x, &dx, &rev, &bas, r1==n? 0: prec); | T->dxbest = dxn; T->xbest = xn; T->indbest = T->ind; | ||
if (DEBUGLEVEL) msgtimer("polred"); | |||
} | } | ||
if (!carrecomplet(divii(dx,dK),&index)) | if (T->ind >= T->indmax) return T->xbest; | ||
err(bugparier,"nfinit (incorrect discriminant)"); | return NULL; | ||
} | |||
ro=get_roots(x,r1,ru,PRECREG); | /* z in Z[X] with positive leading coeff. Set z := z(-X) or z(X) so that | ||
if (DEBUGLEVEL) msgtimer("roots"); | * second coeff is > 0. In place. */ | ||
static int | |||
canon_pol(GEN z) | |||
{ | |||
long i,s; | |||
nf=cgetg(10,t_VEC); | for (i=lgef(z)-2; i>=2; i-=2) | ||
nf[1]=(long)x; | |||
nf[2]=lgetg(3,t_VEC); | |||
mael(nf,2,1)=lstoi(r1); | |||
mael(nf,2,2)=lstoi(r2); | |||
nf[3]=(long)dK; | |||
nf[4]=(long)index; | |||
nf[6]=(long)ro; | |||
nf[7]=(long)bas; | |||
get_nf_matrices(nf, flag & nf_SMALL); | |||
if (flag & nf_ORIG) | |||
{ | { | ||
if (!rev) err(talker,"bad flag in initalgall0"); | s = signe(z[i]); | ||
res = cgetg(3,t_VEC); | if (s > 0) | ||
res[1]=(long)nf; nf = res; | { | ||
res[2]=lead? ldiv(rev,lead): (long)rev; | for (; i>=2; i-=2) z[i]=lnegi((GEN)z[i]); | ||
return -1; | |||
} | |||
if (s) return 1; | |||
} | } | ||
return gerepilecopy(av, nf); | return 0; | ||
} | } | ||
static GEN _polred(GEN x, GEN a, GEN *pta, FP_chk_fun *CHECK); | |||
/* Seek a simpler, polynomial pol defining the same number field as | |||
* x (assumed to be monic at this point) */ | |||
static GEN | |||
nfpolred(int part, nfbasic_t *T) | |||
{ | |||
GEN x = T->x, dx = T->dx, a = T->bas; | |||
GEN phi, xbest, dxbest, mat, d, rev; | |||
long i, n = lg(a)-1, v = varn(x); | |||
ok_pol_t O; | |||
FP_chk_fun chk; | |||
if (degpol(x) == 1) { T->x = gsub(polx[v],gun); return gun; } | |||
if (!dx) dx = mulii(T->dK, sqri(T->index)); | |||
O.ind = 0; | |||
O.indmax = part? min(n,3): n; | |||
O.xbest = NULL; | |||
chk.f = &ok_pol; | |||
chk.data = (void*)&O; | |||
if (!_polred(x, a, NULL, &chk)) | |||
err(talker,"you found a counter-example to a conjecture, please report!"); | |||
xbest = O.xbest; dxbest = O.dxbest; | |||
if (!better_pol(xbest, dxbest, x, dx)) return NULL; /* no improvement */ | |||
/* update T */ | |||
phi = (GEN)a[O.indbest]; | |||
if (canon_pol(xbest) < 0) phi = gneg_i(phi); | |||
if (DEBUGLEVEL>1) fprintferr("xbest = %Z\n",xbest); | |||
rev = modreverse_i(phi, x); | |||
for (i=1; i<=n; i++) a[i] = (long)RX_RXQ_compo((GEN)a[i], rev, xbest); | |||
mat = vecpol_to_mat(Q_remove_denom(a, &d), n); | |||
if (!is_pm1(d)) mat = gdiv(hnfmodid(mat,d), d); else mat = idmat(n); | |||
(void)carrecomplet(diviiexact(dxbest,T->dK), &(T->index)); | |||
T->bas= mat_to_vecpol(mat, v); | |||
T->dx = dxbest; | |||
T->x = xbest; return rev; | |||
} | |||
GEN | GEN | ||
initalgred(GEN x, long prec) | get_nfindex(GEN bas) | ||
{ | { | ||
return initalgall0(x,nf_REDUCE,prec); | gpmem_t av = avma; | ||
long n = lg(bas)-1; | |||
GEN d, mat = vecpol_to_mat(Q_remove_denom(bas, &d), n); | |||
if (!d) { avma = av; return gun; } | |||
return gerepileuptoint(av, diviiexact(gpowgs(d, n), det(mat))); | |||
} | } | ||
void | |||
nfbasic_init(GEN x, long flag, GEN fa, nfbasic_t *T) | |||
{ | |||
GEN bas, dK, dx, index; | |||
long r1; | |||
T->basden = NULL; | |||
T->lead = NULL; | |||
if (DEBUGLEVEL) (void)timer2(); | |||
if (typ(x) == t_POL) | |||
{ | |||
check_pol_int(x, "nfinit"); | |||
if (gisirreducible(x) == gzero) err(redpoler, "nfinit"); | |||
x = pol_to_monic(x, &(T->lead)); | |||
bas = allbase(x, flag, &dx, &dK, &index, &fa); | |||
if (DEBUGLEVEL) msgtimer("round4"); | |||
r1 = sturm(x); | |||
} | |||
else if (typ(x) == t_VEC && lg(x)==3 && typ(x[1])==t_POL) | |||
{ /* monic integral polynomial + integer basis */ | |||
GEN mat; | |||
bas = (GEN)x[2]; x = (GEN)x[1]; | |||
if (typ(bas) == t_MAT) | |||
{ mat = bas; bas = mat_to_vecpol(mat,varn(x)); } | |||
else | |||
mat = vecpol_to_mat(bas, lg(bas)-1); | |||
index = get_nfindex(bas); | |||
dx = ZX_disc(x); | |||
dK = diviiexact(dx, sqri(index)); | |||
r1 = sturm(x); | |||
} | |||
else | |||
{ /* nf, bnf, bnr */ | |||
GEN nf = checknf(x); | |||
x = (GEN)nf[1]; | |||
dK = (GEN)nf[3]; | |||
index = (GEN)nf[4]; | |||
bas = (GEN)nf[7]; | |||
r1 = nf_get_r1(nf); dx = NULL; | |||
} | |||
T->x = x; | |||
T->r1 = r1; | |||
T->dx = dx; | |||
T->dK = dK; | |||
T->bas = bas; | |||
T->index = index; | |||
} | |||
/* Initialize the number field defined by the polynomial x (in variable v) | |||
* flag & nf_RED: try a polred first. | |||
* flag & nf_PARTRED: as nf_RED but check the first two elements only. | |||
* flag & nf_ORIG | |||
* do a polred and return [nfinit(x), Mod(a,red)], where | |||
* Mod(a,red) = Mod(v,x) (i.e return the base change). */ | |||
GEN | GEN | ||
initalgred2(GEN x, long prec) | _initalg(GEN x, long flag, long prec) | ||
{ | { | ||
return initalgall0(x,nf_REDUCE|nf_ORIG,prec); | const gpmem_t av = avma; | ||
GEN nf, rev = NULL, ro = NULL; | |||
nfbasic_t T; | |||
nfbasic_init(x, flag, NULL, &T); | |||
if (T.lead && !(flag & (nf_RED|nf_PARTRED))) | |||
{ | |||
err(warner,"non-monic polynomial. Result of the form [nf,c]"); | |||
flag |= nf_PARTRED | nf_ORIG; | |||
} | |||
if (flag & (nf_RED|nf_PARTRED)) | |||
{ | |||
set_LLL_basis(&T, &ro); | |||
rev = nfpolred(flag & nf_PARTRED, &T); | |||
if (rev) ro = NULL; /* changed T.x */ | |||
if (flag & nf_ORIG) | |||
{ | |||
if (!rev) rev = polx[varn(T.x)]; /* no improvement */ | |||
if (T.lead) rev = gdiv(rev, T.lead); | |||
rev = to_polmod(rev, T.x); | |||
} | |||
if (DEBUGLEVEL) msgtimer("polred"); | |||
} | |||
set_LLL_basis(&T, &ro); | |||
if (DEBUGLEVEL) msgtimer("LLL basis"); | |||
nf = nfbasic_to_nf(&T, ro, prec); | |||
if (flag & nf_ORIG) | |||
{ | |||
GEN res = cgetg(3,t_VEC); | |||
res[1] = (long)nf; | |||
res[2] = (long)rev; nf = res; | |||
} | |||
return gerepilecopy(av, nf); | |||
} | } | ||
GEN | GEN | ||
initalgred(GEN x, long prec) { return _initalg(x, nf_RED, prec); } | |||
GEN | |||
initalgred2(GEN x, long prec) { return _initalg(x, nf_RED|nf_ORIG, prec); } | |||
GEN | |||
initalg(GEN x, long prec) { return _initalg(x, 0, prec); } | |||
GEN | |||
nfinit0(GEN x, long flag,long prec) | nfinit0(GEN x, long flag,long prec) | ||
{ | { | ||
switch(flag) | switch(flag) | ||
{ | { | ||
case 0: | case 0: | ||
case 1: return initalgall0(x,nf_REGULAR,prec); | case 1: return _initalg(x,0,prec); | ||
case 2: return initalgall0(x,nf_REDUCE,prec); | case 2: return _initalg(x,nf_RED,prec); | ||
case 3: return initalgall0(x,nf_REDUCE|nf_ORIG,prec); | case 3: return _initalg(x,nf_RED|nf_ORIG,prec); | ||
case 4: return initalgall0(x,nf_REDUCE|nf_PARTIAL,prec); | case 4: return _initalg(x,nf_PARTRED,prec); | ||
case 5: return initalgall0(x,nf_REDUCE|nf_ORIG|nf_PARTIAL,prec); | case 5: return _initalg(x,nf_PARTRED|nf_ORIG,prec); | ||
case 6: return initalgall0(x,nf_SMALL,prec); | |||
default: err(flagerr,"nfinit"); | default: err(flagerr,"nfinit"); | ||
} | } | ||
return NULL; /* not reached */ | return NULL; /* not reached */ | ||
} | } | ||
GEN | |||
initalg(GEN x, long prec) | |||
{ | |||
return initalgall0(x,nf_REGULAR,prec); | |||
} | |||
/* assume x a bnr/bnf/nf */ | /* assume x a bnr/bnf/nf */ | ||
long | long | ||
nfgetprec(GEN x) | nfgetprec(GEN x) | ||
{ | { | ||
GEN nf = checknf(x), ro = (GEN)nf[6]; | GEN nf = checknf(x), ro = (GEN)nf[6]; | ||
return (typ(ro)==t_VEC)?precision((GEN)ro[1]): DEFAULTPREC; | return (typ(ro)==t_VEC)? precision((GEN)ro[1]): DEFAULTPREC; | ||
} | } | ||
GEN | GEN | ||
nfnewprec(GEN nf, long prec) | nfnewprec(GEN nf, long prec) | ||
{ | { | ||
long av=avma,r1,r2,ru,n; | gpmem_t av = avma; | ||
GEN y,pol,ro,basden,MC,mat,M; | GEN NF; | ||
nffp_t F; | |||
if (typ(nf) != t_VEC) err(talker,"incorrect nf in nfnewprec"); | if (typ(nf) != t_VEC) err(talker,"incorrect nf in nfnewprec"); | ||
if (lg(nf) == 11) return bnfnewprec(nf,prec); | if (lg(nf) == 11) return bnfnewprec(nf,prec); | ||
if (lg(nf) == 7) return bnrnewprec(nf,prec); | if (lg(nf) == 7) return bnrnewprec(nf,prec); | ||
(void)checknf(nf); | (void)checknf(nf); | ||
y = dummycopy(nf); | NF = dummycopy(nf); | ||
pol=(GEN)nf[1]; n=degpol(pol); | NF[5] = (long)dummycopy((GEN)NF[5]); | ||
r1 = nf_get_r1(nf); | remake_GM(NF, &F, prec); | ||
r2 = (n - r1) >> 1; ru = r1+r2; | NF[6] = (long)F.ro; | ||
mat=dummycopy((GEN)nf[5]); | mael(NF,5,1) = (long)F.M; | ||
ro=get_roots(pol,r1,ru,prec); | mael(NF,5,2) = (long)F.G; | ||
y[5]=(long)mat; | return gerepilecopy(av, NF); | ||
y[6]=(long)ro; | } | ||
basden = get_bas_den((GEN)nf[7]); | |||
M = make_M(basden,ro); | /********************************************************************/ | ||
MC = make_MC(r1,M); | /** **/ | ||
mat[1]=(long)M; | /** POLRED **/ | ||
if (typ(nf[8]) != t_INT) mat[2]=(long)MC; /* not a small nf */ | /** **/ | ||
mat[3]=(long)mulmat_real(MC,M); | /********************************************************************/ | ||
/* remove duplicate polynomials in y, updating a (same indexes), in place */ | |||
static void | |||
remove_duplicates(GEN y, GEN a) | |||
{ | |||
long k, i, l = lg(y); | |||
gpmem_t av = avma; | |||
GEN z; | |||
if (l < 2) return; | |||
z = new_chunk(3); | |||
z[1] = (long)y; | |||
z[2] = (long)a; (void)sort_factor(z, gcmp); | |||
for (k=1, i=2; i<l; i++) | |||
if (!gegal((GEN)y[i], (GEN)y[k])) | |||
{ | |||
k++; | |||
a[k] = a[i]; | |||
y[k] = y[i]; | |||
} | |||
l = k+1; setlg(a,l); setlg(y,l); | |||
avma = av; | |||
} | |||
/* if CHECK != NULL, return the first polynomial pol found such that | |||
* CHECK->f(CHECK->data, pol) != 0; return NULL if there are no such pol */ | |||
static GEN | |||
_polred(GEN x, GEN a, GEN *pta, FP_chk_fun *CHECK) | |||
{ | |||
long i, v = varn(x), l = lg(a); | |||
GEN ch, d, y = cgetg(l,t_VEC); | |||
for (i=1; i<l; i++) | |||
{ | |||
if (DEBUGLEVEL>2) { fprintferr("i = %ld\n",i); flusherr(); } | |||
ch = QX_caract(x, (GEN)a[i], v); | |||
if (CHECK) | |||
{ | |||
ch = CHECK->f(CHECK->data, ch); | |||
if (!ch) continue; | |||
return ch; | |||
} | |||
d = modulargcd(derivpol(ch), ch); | |||
if (degpol(d)) ch = gdivexact(ch,d); | |||
if (canon_pol(ch) < 0 && pta) a[i] = (long)gneg_i((GEN)a[i]); | |||
if (DEBUGLEVEL > 3) outerr(ch); | |||
y[i] = (long)ch; | |||
} | |||
if (CHECK) return NULL; /* no suitable polynomial found */ | |||
remove_duplicates(y,a); | |||
if (pta) *pta = a; | |||
return y; | |||
} | |||
static GEN | |||
allpolred(GEN x, long flag, GEN fa, GEN *pta, FP_chk_fun *CHECK) | |||
{ | |||
GEN ro = NULL; | |||
nfbasic_t T; nfbasic_init(x, flag, fa, &T); | |||
if (T.lead) err(impl,"polred for non-monic polynomial"); | |||
set_LLL_basis(&T, &ro); | |||
return _polred(T.x, T.bas, pta, CHECK); | |||
} | |||
/* FIXME: backward compatibility */ | |||
#define red_PARTIAL 1 | |||
#define red_ORIG 2 | |||
GEN | |||
polred0(GEN x, long flag, GEN fa) | |||
{ | |||
gpmem_t av = avma; | |||
GEN y, a; | |||
int fl = 0; | |||
if (fa && gcmp0(fa)) fa = NULL; /* compatibility */ | |||
if (flag & red_PARTIAL) fl |= nf_PARTIALFACT; | |||
if (flag & red_ORIG) fl |= nf_ORIG; | |||
y = allpolred(x, fl, fa, &a, NULL); | |||
if (fl & nf_ORIG) { | |||
GEN z = cgetg(3,t_MAT); | |||
z[1] = (long)a; | |||
z[2] = (long)y; y = z; | |||
} | |||
return gerepilecopy(av, y); | return gerepilecopy(av, y); | ||
} | } | ||
GEN | |||
ordred(GEN x) | |||
{ | |||
gpmem_t av = avma; | |||
GEN y; | |||
if (typ(x) != t_POL) err(typeer,"ordred"); | |||
if (!gcmp1(leading_term(x))) err(impl,"ordred"); | |||
if (!signe(x)) return gcopy(x); | |||
y = cgetg(3,t_VEC); | |||
y[1] = (long)x; | |||
y[2] = (long)idmat(degpol(x)); | |||
return gerepileupto(av, polred(y)); | |||
} | |||
GEN | |||
polredfirstpol(GEN x, long flag, FP_chk_fun *CHECK) | |||
{ return allpolred(x,flag,NULL,NULL,CHECK); } | |||
GEN | |||
polred(GEN x) | |||
{ return polred0(x, 0, NULL); } | |||
GEN | |||
smallpolred(GEN x) | |||
{ return polred0(x, red_PARTIAL, NULL); } | |||
GEN | |||
factoredpolred(GEN x, GEN fa) | |||
{ return polred0(x, 0, fa); } | |||
GEN | |||
polred2(GEN x) | |||
{ return polred0(x, red_ORIG, NULL); } | |||
GEN | |||
smallpolred2(GEN x) | |||
{ return polred0(x, red_PARTIAL|red_ORIG, NULL); } | |||
GEN | |||
factoredpolred2(GEN x, GEN fa) | |||
{ return polred0(x, red_PARTIAL, fa); } | |||
/********************************************************************/ | |||
/** **/ | |||
/** POLREDABS **/ | |||
/** **/ | |||
/********************************************************************/ | |||
GEN | |||
T2_from_embed_norm(GEN x, long r1) | |||
{ | |||
GEN p = sum(x, 1, r1); | |||
GEN q = sum(x, r1+1, lg(x)-1); | |||
if (q != gzero) p = gadd(p, gmul2n(q,1)); | |||
return p; | |||
} | |||
GEN | |||
T2_from_embed(GEN x, long r1) | |||
{ | |||
return T2_from_embed_norm(gnorm(x), r1); | |||
} | |||
typedef struct { | |||
long r1, v, prec; | |||
GEN ZKembed; /* embeddings of fincke-pohst-reduced Zk basis */ | |||
GEN u; /* matrix giving fincke-pohst-reduced Zk basis */ | |||
GEN M; /* embeddings of initial (LLL-reduced) Zk basis */ | |||
GEN bound; /* T2 norm of the polynomial defining nf */ | |||
} CG_data; | |||
/* characteristic pol of x */ | |||
static GEN | |||
get_polchar(CG_data *d, GEN x) | |||
{ | |||
GEN g = gmul(d->ZKembed,x); | |||
long e; | |||
g = grndtoi(roots_to_pol_r1r2(g, d->r1, d->v), &e); | |||
if (e > -5) err(precer, "get_polchar"); | |||
return g; | |||
} | |||
/* return a defining polynomial for Q(x) */ | |||
static GEN | |||
get_polmin(CG_data *d, GEN x) | |||
{ | |||
GEN g = get_polchar(d,x); | |||
GEN h = modulargcd(g,derivpol(g)); | |||
if (degpol(h)) g = gdivexact(g,h); | |||
return g; | |||
} | |||
/* does x generate the correct field ? */ | |||
static GEN | |||
chk_gen(void *data, GEN x) | |||
{ | |||
gpmem_t av = avma; | |||
GEN g = get_polchar((CG_data*)data,x); | |||
GEN h = modulargcd(g,derivpol(g)); | |||
if (degpol(h)) { avma = av; return NULL; } | |||
if (DEBUGLEVEL>3) fprintferr(" generator: %Z\n",g); | |||
return g; | |||
} | |||
/* mat = base change matrix, gram = LLL-reduced T2 matrix */ | |||
static GEN | |||
chk_gen_init(FP_chk_fun *chk, GEN gram, GEN mat) | |||
{ | |||
CG_data *d = (CG_data*)chk->data; | |||
GEN P,bound,prev,x,B; | |||
long l = lg(gram), N = l-1,i,prec,prec2; | |||
int skipfirst = 0; | |||
d->u = mat; | |||
d->ZKembed = gmul(d->M, mat); | |||
bound = d->bound; | |||
prev = NULL; | |||
x = zerocol(N); | |||
for (i=1; i<l; i++) | |||
{ | |||
B = gcoeff(gram,i,i); | |||
if (gcmp(B,bound) >= 0) continue; /* don't assume increasing norms */ | |||
x[i] = un; P = get_polmin(d,x); | |||
x[i] = zero; | |||
if (degpol(P) == N) | |||
{ | |||
if (gcmp(B,bound) < 0) bound = B ; | |||
continue; | |||
} | |||
if (DEBUGLEVEL>2) fprintferr("chk_gen_init: subfield %Z\n",P); | |||
if (skipfirst != i-1) continue; | |||
/* Q(w[1],...,w[i-1]) = Q[X]/(prev) is a strict subfield of nf */ | |||
if (prev && degpol(prev) < N && !gegal(prev,P)) | |||
{ | |||
if (degpol(prev) * degpol(P) > 64) continue; /* too expensive */ | |||
P = compositum(prev,P); | |||
P = (GEN)P[lg(P)-1]; | |||
if (degpol(P) == N) continue; | |||
if (DEBUGLEVEL>2 && degpol(P) != degpol(prev)) | |||
fprintferr("chk_gen_init: subfield %Z\n",P); | |||
} | |||
skipfirst++; prev = P; | |||
} | |||
/* x_1,...,x_skipfirst generate a strict subfield [unless N=skipfirst=1] */ | |||
chk->skipfirst = skipfirst; | |||
if (DEBUGLEVEL>2) fprintferr("chk_gen_init: skipfirst = %ld\n",skipfirst); | |||
/* should be gexpo( max_k C^n_k (bound/k)^(k/2) ) */ | |||
prec2 = (1 + (((gexpo(bound)*N)/2) >> TWOPOTBITS_IN_LONG)); | |||
if (prec2 < 0) prec2 = 0; | |||
prec = 3 + prec2; | |||
if (DEBUGLEVEL) | |||
fprintferr("chk_gen_init: new prec = %ld (initially %ld)\n", prec, d->prec); | |||
if (prec > d->prec) err(bugparier, "precision problem in polredabs"); | |||
if (prec < d->prec) d->ZKembed = gprec_w(d->ZKembed, prec); | |||
return bound; | |||
} | |||
/* store phi(beta mod z). Suitable for gerepileupto */ | |||
static GEN | |||
storeeval(GEN beta, GEN z, GEN lead) | |||
{ | |||
GEN y = cgetg(3,t_VEC); | |||
z = gcopy(z); | |||
beta = lead? gdiv(beta, lead): gcopy(beta); | |||
y[1] = (long)z; | |||
y[2] = (long)to_polmod(beta,z); | |||
return y; | |||
} | |||
static GEN | |||
storeraw(GEN beta, GEN z) | |||
{ | |||
GEN y = cgetg(3,t_VEC); | |||
y[1] = lcopy(z); | |||
y[2] = lcopy(beta); return y; | |||
} | |||
static void | |||
findmindisc(GEN *py, GEN *pa) | |||
{ | |||
GEN v, dmin, z, b, discs, y = *py, a = *pa; | |||
long i,k, l = lg(y); | |||
if (l == 2) { *py = (GEN)y[1]; *pa = (GEN)a[1]; return; } | |||
discs = cgetg(l,t_VEC); | |||
for (i=1; i<l; i++) discs[i] = labsi(ZX_disc((GEN)y[i])); | |||
v = sindexsort(discs); | |||
k = v[1]; | |||
dmin = (GEN)discs[k]; | |||
z = (GEN)y[k]; | |||
b = (GEN)a[k]; | |||
for (i=2; i<l; i++) | |||
{ | |||
k = v[i]; | |||
if (!egalii((GEN)discs[k], dmin)) break; | |||
if (gpolcomp((GEN)y[k],z) < 0) { z = (GEN)y[k]; b = (GEN)a[k]; } | |||
} | |||
*py = z; *pa = b; | |||
} | |||
static GEN | |||
storepol(GEN x, GEN z, GEN a, GEN lead, long flag) | |||
{ | |||
GEN y, b; | |||
if (flag & nf_RAW) | |||
y = storeraw(a, z); | |||
else if (flag & nf_ORIG) | |||
{ | |||
b = modreverse_i(a, x); | |||
y = storeeval(b,z, lead); | |||
} | |||
else | |||
y = gcopy(z); | |||
return y; | |||
} | |||
static GEN | |||
storeallpol(GEN x, GEN z, GEN a, GEN lead, long flag) | |||
{ | |||
GEN y, b; | |||
if (flag & nf_RAW) | |||
{ | |||
long i, c = lg(z); | |||
y = cgetg(c,t_VEC); | |||
for (i=1; i<c; i++) y[i] = (long)storeraw((GEN)a[i], (GEN)z[i]); | |||
} | |||
else if (flag & nf_ORIG) | |||
{ | |||
long i, c = lg(z); | |||
b = cgetg(c, t_VEC); | |||
for (i=1; i<c; i++) | |||
b[i] = (long)modreverse_i((GEN)a[i], x); | |||
y = cgetg(c,t_VEC); | |||
for (i=1; i<c; i++) y[i] = (long)storeeval((GEN)b[i], (GEN)z[i], lead); | |||
} | |||
else | |||
y = gcopy(z); | |||
return y; | |||
} | |||
GEN | |||
_polredabs(nfbasic_t *T, GEN *u) | |||
{ | |||
long i, prec, e, n = degpol(T->x); | |||
GEN v, ro = NULL; | |||
FP_chk_fun chk = { &chk_gen, &chk_gen_init, NULL, 0 }; | |||
nffp_t F; | |||
CG_data d; chk.data = (void*)&d; | |||
set_LLL_basis(T, &ro); | |||
if (DEBUGLEVEL) msgtimer("LLL basis"); | |||
/* || polchar ||_oo < 2^e */ | |||
e = n * (gexpo(gmulsg(n, cauchy_bound(T->x))) + 1); | |||
prec = DEFAULTPREC + (e >> TWOPOTBITS_IN_LONG); | |||
get_nf_fp_compo(T, &F, ro, prec); | |||
d.v = varn(T->x); | |||
d.r1= T->r1; | |||
d.bound = gmul(T2_from_embed(F.ro, d.r1), dbltor(1.00000001)); | |||
for (i=1; ; i++) | |||
{ | |||
GEN R = R_from_QR(F.G, prec); | |||
d.prec = prec; | |||
d.M = F.M; | |||
if (R) | |||
{ | |||
v = fincke_pohst(_vec(R),NULL,5000,1, 0, &chk); | |||
if (v) break; | |||
} | |||
if (i==MAXITERPOL) err(accurer,"polredabs0"); | |||
prec = (prec<<1)-2; | |||
get_nf_fp_compo(T, &F, NULL, prec); | |||
if (DEBUGLEVEL) err(warnprec,"polredabs0",prec); | |||
} | |||
*u = d.u; return v; | |||
} | |||
GEN | |||
polredabs0(GEN x, long flag) | |||
{ | |||
long i, l, vx; | |||
gpmem_t av = avma; | |||
GEN y, a, u; | |||
nfbasic_t T; | |||
nfbasic_init(x, flag & nf_PARTIALFACT, NULL, &T); | |||
x = T.x; vx = varn(x); | |||
if (degpol(x) == 1) | |||
{ | |||
u = NULL; | |||
y = _vec(polx[vx]); | |||
a = _vec(gsub((GEN)y[1], (GEN)x[2])); | |||
} | |||
else | |||
{ | |||
GEN v = _polredabs(&T, &u); | |||
y = (GEN)v[1]; | |||
a = (GEN)v[2]; | |||
} | |||
l = lg(a); | |||
for (i=1; i<l; i++) | |||
if (canon_pol((GEN)y[i]) < 0) a[i] = (long)gneg_i((GEN)a[i]); | |||
remove_duplicates(y,a); | |||
l = lg(a); | |||
if (l == 1) | |||
{ | |||
y = _vec(x); | |||
a = _vec(polx[vx]); | |||
} | |||
if (DEBUGLEVEL) fprintferr("%ld minimal vectors found.\n",l-1); | |||
if (flag & nf_ALL) | |||
{ | |||
if (u) for (i=1; i<l; i++) a[i] = lmul(T.bas, gmul(u, (GEN)a[i])); | |||
y = storeallpol(x,y,a,T.lead,flag); | |||
} | |||
else | |||
{ | |||
findmindisc(&y, &a); | |||
if (u) a = gmul(T.bas, gmul(u, a)); | |||
y = storepol(x,y,a,T.lead,flag); | |||
} | |||
return gerepileupto(av, y); | |||
} | |||
GEN | |||
polredabsall(GEN x, long flun) { return polredabs0(x, flun | nf_ALL); } | |||
GEN | |||
polredabs(GEN x) { return polredabs0(x,0); } | |||
GEN | |||
polredabs2(GEN x) { return polredabs0(x,nf_ORIG); } | |||
static long | static long | ||
nf_pm1(GEN y) | nf_pm1(GEN y) | ||
{ | { | ||
|
|
||
for (i=2; i<l; i++) | for (i=2; i<l; i++) | ||
if (signe(y[i])) return 0; | if (signe(y[i])) return 0; | ||
return signe(y[1]); | return signe(y[1]); | ||
} | } | ||
static GEN | static GEN | ||
|
|
||
return x; | return x; | ||
} | } | ||
/* only roots of 1 are +/- 1 */ | |||
static GEN | |||
trivroots(GEN nf) | |||
{ | |||
GEN y = cgetg(3, t_VEC); | |||
y[1] = deux; | |||
y[2] = (long)gscalcol_i(negi(gun), degpol(nf[1])); | |||
return y; | |||
} | |||
GEN | GEN | ||
rootsof1(GEN nf) | rootsof1(GEN nf) | ||
{ | { | ||
ulong av; | gpmem_t av = avma; | ||
long N,k,i,ws,prec; | long N,k,i,ws,prec; | ||
GEN algun,p1,y,R1,d,list,w; | GEN p1,y,d,list,w; | ||
y=cgetg(3,t_VEC); av=avma; nf=checknf(nf); | nf = checknf(nf); | ||
R1=gmael(nf,2,1); algun=gmael(nf,8,1); | if ( nf_get_r1(nf) ) return trivroots(nf); | ||
if (signe(R1)) | |||
{ | N = degpol(nf[1]); prec = nfgetprec(nf); | ||
y[1]=deux; | |||
y[2]=lneg(algun); return y; | |||
} | |||
N=degpol(nf[1]); prec=gprecision((GEN)nf[6]); | |||
#ifdef LONG_IS_32BIT | |||
if (prec < 10) prec = 10; | |||
#else | |||
if (prec < 6) prec = 6; | |||
#endif | |||
for (i=1; ; i++) | for (i=1; ; i++) | ||
{ | { | ||
p1 = fincke_pohst(nf,stoi(N),1000,1,prec,NULL); | GEN R = R_from_QR(gmael(nf,5,2), prec); | ||
if (p1) break; | if (R) | ||
{ | |||
p1 = fincke_pohst(_vec(R),stoi(N),1000,1, 0, NULL); | |||
if (p1) break; | |||
} | |||
if (i == MAXITERPOL) err(accurer,"rootsof1"); | if (i == MAXITERPOL) err(accurer,"rootsof1"); | ||
prec=(prec<<1)-2; | prec = (prec<<1)-2; | ||
if (DEBUGLEVEL) err(warnprec,"rootsof1",prec); | if (DEBUGLEVEL) err(warnprec,"rootsof1",prec); | ||
nf=nfnewprec(nf,prec); | nf = nfnewprec(nf,prec); | ||
} | } | ||
if (itos(ground((GEN)p1[2])) != N) err(bugparier,"rootsof1 (bug1)"); | if (itos(ground((GEN)p1[2])) != N) err(bugparier,"rootsof1 (bug1)"); | ||
w=(GEN)p1[1]; ws = itos(w); | w = (GEN)p1[1]; ws = itos(w); | ||
if (ws == 2) | if (ws == 2) { avma = av; return trivroots(nf); } | ||
{ | |||
y[1]=deux; avma=av; | |||
y[2]=lneg(algun); return y; | |||
} | |||
d = decomp(w); list = (GEN)p1[3]; k = lg(list); | d = decomp(w); list = (GEN)p1[3]; k = lg(list); | ||
for (i=1; i<k; i++) | for (i=1; i<k; i++) | ||
{ | { | ||
p1 = (GEN)list[i]; | p1 = (GEN)list[i]; | ||
p1 = is_primitive_root(nf,d,p1,ws); | p1 = is_primitive_root(nf,d,p1,ws); | ||
if (p1) | if (p1) break; | ||
{ | |||
y[2]=lpilecopy(av,p1); | |||
y[1]=lstoi(ws); return y; | |||
} | |||
} | } | ||
err(bugparier,"rootsof1"); | if (!p1) err(bugparier,"rootsof1"); | ||
return NULL; /* not reached */ | y = cgetg(3, t_VEC); | ||
y[1] = lstoi(ws); | |||
y[2] = (long)p1; return gerepilecopy(av, y); | |||
} | } | ||
/*******************************************************************/ | /*******************************************************************/ | ||
|
|
||
/* DEDEKIND ZETA FUNCTION */ | /* DEDEKIND ZETA FUNCTION */ | ||
/* */ | /* */ | ||
/*******************************************************************/ | /*******************************************************************/ | ||
ulong smulss(ulong x, ulong y, ulong *rem); | |||
static GEN | static GEN | ||
dirzetak0(GEN nf, long N0) | dirzetak0(GEN nf, long N0) | ||
{ | { | ||
GEN vect,p1,pol,disc,c,c2; | GEN vect,p1,pol,disc,c,c2; | ||
long av=avma,i,j,k,limk,lx; | gpmem_t av=avma; | ||
long i,j,k,limk,lx; | |||
ulong q,p,rem; | ulong q,p,rem; | ||
byteptr d=diffptr; | byteptr d=diffptr; | ||
long court[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),0}; | long court[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0}; | ||
pol=(GEN)nf[1]; disc=(GEN)nf[4]; | pol=(GEN)nf[1]; disc=(GEN)nf[4]; | ||
c = (GEN) gpmalloc((N0+1)*sizeof(long)); | c = (GEN) gpmalloc((N0+1)*sizeof(long)); | ||
|
|
||
while (court[2]<=N0) | while (court[2]<=N0) | ||
{ | { | ||
court[2] += *d++; if (! *d) err(primer1); | NEXT_PRIME_VIADIFF_CHECK(court[2], d); | ||
if (smodis(disc,court[2])) /* court does not divide index */ | if (smodis(disc,court[2])) /* court does not divide index */ | ||
{ vect = (GEN) simplefactmod(pol,court)[1]; lx=lg(vect); } | { vect = (GEN) simplefactmod(pol,court)[1]; lx=lg(vect); } | ||
else | else | ||
|
|
||
GEN nfz,nf,alpha,beta,mu,gr1,gr2,gru,p1,p2,cst,A0,c0,c1,c2,eps,coef; | GEN nfz,nf,alpha,beta,mu,gr1,gr2,gru,p1,p2,cst,A0,c0,c1,c2,eps,coef; | ||
GEN limx,bnf,resi,zet,C,coeflog,racpi,aij,tabj,colzero, *tabcstn, *tabcstni; | GEN limx,bnf,resi,zet,C,coeflog,racpi,aij,tabj,colzero, *tabcstn, *tabcstni; | ||
GEN c_even,ck_even,c_odd,ck_odd,serie_even,serie_odd,serie_exp,Pi; | GEN c_even,ck_even,c_odd,ck_odd,serie_even,serie_odd,serie_exp,Pi; | ||
long N0,imin,imax,r1,r2,ru,R,N,i,j,k,n, av,av2,tetpil; | long N0,imin,imax,r1,r2,ru,R,N,i,j,k,n; | ||
long court[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),0}; | gpmem_t av,av2,tetpil; | ||
long court[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0}; | |||
stackzone *zone, *zone0, *zone1; | stackzone *zone, *zone0, *zone1; | ||
/*************** Calcul du residu et des constantes ***************/ | /*************** Calcul du residu et des constantes ***************/ | ||
eps=gmul2n(gun,-bit_accuracy(prec)-6); p1=dbltor(0.5); | eps=gmul2n(gun,-bit_accuracy(prec)-6); p1=dbltor(0.5); | ||
nfz=cgetg(10,t_VEC); | nfz=cgetg(10,t_VEC); | ||
bnf=buchinit(pol,p1,p1,prec+1); prec=(prec<<1)-1; | bnf = bnfinit0(pol,2,NULL,prec+1); prec=(prec<<1)-1; | ||
bnf = checkbnf_discard(bnf); | bnf = checkbnf_discard(bnf); | ||
Pi = mppi(prec); racpi=gsqrt(Pi,prec); | Pi = mppi(prec); racpi=gsqrt(Pi,prec); | ||
|
|
||
if (is_bigint(p1) || N0 > 10000000) | if (is_bigint(p1) || N0 > 10000000) | ||
err(talker,"discriminant too large for initzeta, sorry"); | err(talker,"discriminant too large for initzeta, sorry"); | ||
if (DEBUGLEVEL>=2) | if (DEBUGLEVEL>=2) | ||
{ fprintferr("\ninitzeta:\nN0 = %ld\n",N0); flusherr(); timer2(); } | { fprintferr("\ninitzeta:\nN0 = %ld\n",N0); flusherr(); (void)timer2(); } | ||
/* Calcul de imax */ | /* Calcul de imax */ | ||
|
|
||
zone = switch_stack(NULL,i + 2*(N0+1) + 6*prec); | zone = switch_stack(NULL,i + 2*(N0+1) + 6*prec); | ||
zone1 = switch_stack(NULL,2*i); | zone1 = switch_stack(NULL,2*i); | ||
zone0 = switch_stack(NULL,2*i); | zone0 = switch_stack(NULL,2*i); | ||
switch_stack(zone,1); | (void)switch_stack(zone,1); | ||
tabcstn = (GEN*) cgetg(N0+1,t_VEC); | tabcstn = (GEN*) cgetg(N0+1,t_VEC); | ||
tabcstni = (GEN*) cgetg(N0+1,t_VEC); | tabcstni = (GEN*) cgetg(N0+1,t_VEC); | ||
p1 = ginv(cst); | p1 = ginv(cst); | ||
for (i=1; i<=N0; i++) { tabcstni[i] = gun; tabcstn[i] = mulsr(i,p1); } | for (i=1; i<=N0; i++) { tabcstni[i] = gun; tabcstn[i] = mulsr(i,p1); } | ||
switch_stack(zone,0); | (void)switch_stack(zone,0); | ||
/********** Calcul des coefficients a(i,j) independants de s **********/ | /********** Calcul des coefficients a(i,j) independants de s **********/ | ||
|
|
||
if (k&1) ck_odd[k]=lneg((GEN)ck_odd[k]); | if (k&1) ck_odd[k]=lneg((GEN)ck_odd[k]); | ||
ck_odd[k]=ladd(gru,(GEN)ck_odd[k]); | ck_odd[k]=ladd(gru,(GEN)ck_odd[k]); | ||
} | } | ||
ck_odd[1]=lsub((GEN)ck_odd[1],gmul(gr1,glog(gdeux,prec))); | ck_odd[1]=lsub((GEN)ck_odd[1],gmul(gr1,mplog2(prec))); | ||
serie_even =cgetg(ru+3,t_SER); serie_odd=cgetg(r2+3,t_SER); | serie_even =cgetg(ru+3,t_SER); serie_odd=cgetg(r2+3,t_SER); | ||
serie_even[1] = serie_odd[1] = evalsigne(1)+evalvalp(1); | serie_even[1] = serie_odd[1] = evalsigne(1)+evalvalp(1); | ||
i=0; | i=0; | ||
|
|
||
t[1] = lstoi(coef[i]); t[2] = lmul((GEN)t[1],p1); | t[1] = lstoi(coef[i]); t[2] = lmul((GEN)t[1],p1); | ||
for (j=2; j<=ru; j++) | for (j=2; j<=ru; j++) | ||
{ | { | ||
long av2 = avma; p2 = gmul((GEN)t[j],p1); | gpmem_t av2 = avma; | ||
p2 = gmul((GEN)t[j], p1); | |||
t[j+1] = (long)gerepileuptoleaf(av2, divrs(p2,j)); | t[j+1] = (long)gerepileuptoleaf(av2, divrs(p2,j)); | ||
} | } | ||
/* tabj[n,j]=coef(n)*ln(c/n)^(j-1)/(j-1)! */ | /* tabj[n,j]=coef(n)*ln(c/n)^(j-1)/(j-1)! */ | ||
|
|
||
} | } | ||
/* use a parallel stack */ | /* use a parallel stack */ | ||
z = i&1? zone1: zone0; | z = i&1? zone1: zone0; | ||
switch_stack(z, 1); | (void)switch_stack(z, 1); | ||
for (n=1; n<=N0; n++) | for (n=1; n<=N0; n++) | ||
if (coef[n]) tabcstni[n] = mpmul(tabcstni[n],tabcstn[n]); | if (coef[n]) tabcstni[n] = mpmul(tabcstni[n],tabcstn[n]); | ||
/* come back */ | /* come back */ | ||
switch_stack(z, 0); | (void)switch_stack(z, 0); | ||
} | } | ||
nfz[4] = (long) C; | nfz[4] = (long) C; | ||
if (DEBUGLEVEL>=2) msgtimer("Cik"); | if (DEBUGLEVEL>=2) msgtimer("Cik"); | ||
|
|
||
GEN resi,C,cst,cstlog,coeflog,cs,coef; | GEN resi,C,cst,cstlog,coeflog,cs,coef; | ||
GEN lambd,gammas,gammaunmoins,gammas2,gammaunmoins2,var1,var2; | GEN lambd,gammas,gammaunmoins,gammas2,gammaunmoins2,var1,var2; | ||
GEN p1,unmoins,gexpro,gar,val,valm,valk,valkm; | GEN p1,unmoins,gexpro,gar,val,valm,valk,valkm; | ||
long ts = typ(s), r1,r2,ru,imax,i,j,k,N0,sl,prec,bigprec, av = avma; | long ts = typ(s), r1,r2,ru,imax,i,j,k,N0,sl,prec,bigprec; | ||
gpmem_t av = avma; | |||
if (typ(nfz)!=t_VEC || lg(nfz)!=10 || typ(nfz[1]) != t_VEC) | if (typ(nfz)!=t_VEC || lg(nfz)!=10 || typ(nfz[1]) != t_VEC) | ||
err(talker,"not a zeta number field in zetakall"); | err(talker,"not a zeta number field in zetakall"); | ||
|
|
||
for (k=1; k<=ru; k++) | for (k=1; k<=ru; k++) | ||
{ | { | ||
p1 = gcoeff(C,i,k); | p1 = gcoeff(C,i,k); | ||
lambd = gsub(lambd,gdiv(p1,valk )); valk = mulii(val,valk); | lambd = gsub(lambd,gdiv(p1,valk )); valk = mulii(val,valk); | ||
lambd = gsub(lambd,gdiv(p1,valkm)); valkm = mulii(valm,valkm); | lambd = gsub(lambd,gdiv(p1,valkm)); valkm = mulii(valm,valkm); | ||
} | } | ||
val = addis(val,1); | val = addis(val,1); |