Annotation of OpenXM_contrib/pari/src/modules/subfield.c, Revision 1.1.1.1
1.1 maekawa 1: /*******************************************************************/
2: /* */
3: /* SUBFIELDS OF A NUMBER FIELD */
4: /* */
5: /* J. Klueners and M. Pohst, J. Symb. Comp. (1996), vol. 11 */
6: /* */
7: /*******************************************************************/
8: /* $Id: subfield.c,v 1.1.1.1 1999/09/16 13:48:20 karim Exp $ */
9: #include "pari.h"
10: #ifdef __WIN32
11: # include <io.h> /* for open, read, close */
12: #endif
13: GEN roots_to_pol(GEN a, long v);
14:
15: static long TR; /* nombre de changements de polynomes (degre fixe) */
16: static GEN FACTORDL; /* factorisation of |disc(L)| */
17:
18: static GEN print_block_system(long N,GEN Y,long d, GEN vbs);
19: static GEN compute_data(GEN nf,GEN ff,GEN p,long m,long nn);
20:
21: /* Computation of potential block systems of given size d associated to a
22: * rational prime p: give a row vector of row vectors containing the
23: * potential block systems of imprimitivity; a potential block system is a
24: * vector of row vectors (enumeration of the roots).
25: */
26: #define BIL 32 /* for 64bit machines also */
27: static GEN
28: calc_block(long N,GEN Z,long d,GEN Y,GEN vbs)
29: {
30: long r,lK,i,j,k,t,tp,T,lpn,u,nn,lnon,lY;
31: GEN K,n,non,pn,pnon,e,Yp,Zp,Zpp;
32:
33: if (DEBUGLEVEL>3)
34: {
35: long l = vbs? lg(vbs): 0;
36: fprintferr("avma = %ld, lg(Z) = %ld, lg(Y) = %ld, lg(vbs) = %ld\n",
37: avma,lg(Z),lg(Y),l);
38: if (DEBUGLEVEL > 5)
39: {
40: fprintferr("Z = %Z\n",Z);
41: fprintferr("Y = %Z\n",Y);
42: if (DEBUGLEVEL > 7) fprintferr("vbs = %Z\n",vbs);
43: }
44: }
45: r=lg(Z); lnon = min(BIL, r);
46: e = new_chunk(BIL);
47: n = new_chunk(r);
48: non = new_chunk(lnon);
49: pnon = new_chunk(lnon);
50: pn = new_chunk(lnon);
51: Zp = cgetg(lnon,t_VEC);
52: Zpp = cgetg(lnon,t_VEC);
53: for (i=1; i<r; i++) n[i] = lg(Z[i])-1;
54:
55: K=divisors(stoi(n[1])); lK=lg(K);
56: for (i=1; i<lK; i++)
57: {
58: lpn=0; k = itos((GEN)K[i]);
59: for (j=2; j<r; j++)
60: if (n[j]%k == 0)
61: {
62: if (++lpn >= BIL) err(talker,"overflow in calc_block");
63: pn[lpn]=n[j]; pnon[lpn]=j;
64: }
65: if (!lpn)
66: {
67: if (d*k != n[1]) continue;
68: T=1; lpn=1;
69: }
70: else
71: T = 1<<lpn;
72: for (t=0; t<T; t++)
73: {
74: for (nn=n[1],tp=t, u=1; u<=lpn; u++,tp>>=1)
75: {
76: if (tp&1) { nn += pn[u]; e[u]=1; } else e[u]=0;
77: }
78: if (d*k == nn)
79: {
80: long av=avma;
81: int Z_equal_Zp = 1;
82:
83: for (j=1; j<lnon; j++) non[j]=0;
84: Zp[1]=Z[1];
85: for (u=2,j=1; j<=lpn; j++)
86: if (e[j])
87: {
88: Zp[u]=Z[pnon[j]]; non[pnon[j]]=1;
89: if (Zp[u] != Z[u]) Z_equal_Zp = 0;
90: u++;
91: }
92: setlg(Zp, u);
93: lY=lg(Y); Yp = cgetg(lY+1,t_VEC);
94: for (j=1; j<lY; j++) Yp[j]=Y[j];
95: Yp[lY]=(long)Zp;
96: if (r == u && Z_equal_Zp)
97: vbs = print_block_system(N,Yp,d,vbs);
98: else
99: {
100: for (u=1,j=2; j<r; j++)
101: if (!non[j]) Zpp[u++] = Z[j];
102: setlg(Zpp, u);
103: vbs = calc_block(N,Zpp,d,Yp,vbs);
104: }
105: avma=av;
106: }
107: }
108: }
109: return vbs;
110: }
111:
112: static GEN
113: potential_block_systems(long N, long d,GEN ff,long *n)
114: {
115: long av=avma,r,i,j,k;
116: GEN p1,vbs,Z;
117:
118: r=lg(ff); Z=cgetg(r,t_VEC);
119: for (k=0,i=1; i<r; i++)
120: {
121: p1=cgetg(n[i]+1,t_VECSMALL); Z[i]=(long)p1;
122: for (j=1; j<=n[i]; j++) p1[j]= ++k;
123: }
124: vbs=calc_block(N,Z,d, cgetg(1,t_VEC), NULL);
125: avma=av; return vbs;
126: }
127:
128: /* product of permutations. Put the result in perm1. */
129: static void
130: perm_mul(GEN perm1,GEN perm2)
131: {
132: long av = avma,i, N = lg(perm1);
133: GEN perm=new_chunk(N);
134: for (i=1; i<N; i++) perm[i]=perm1[perm2[i]];
135: for (i=1; i<N; i++) perm1[i]=perm[i];
136: avma=av;
137: }
138:
139: /* cy is a cycle; compute cy^l as a permutation */
140: static GEN
141: cycle_power_to_perm(GEN perm,GEN cy,long l)
142: {
143: long lp,i,j,b, N = lg(perm), lcy = lg(cy)-1;
144:
145: lp = l % lcy;
146: for (i=1; i<N; i++) perm[i] = i;
147: if (lp)
148: {
149: long av = avma;
150: GEN p1 = new_chunk(N);
151: b = cy[1];
152: for (i=1; i<lcy; i++) b = (perm[b] = cy[i+1]);
153: perm[b] = cy[1];
154: for (i=1; i<N; i++) p1[i] = perm[i];
155:
156: for (j=2; j<=lp; j++) perm_mul(perm,p1);
157: avma = av;
158: }
159: return perm;
160: }
161:
162: /* image du block system D par la permutation perm */
163: static GEN
164: im_block_by_perm(GEN D,GEN perm)
165: {
166: long i,j,lb,lcy;
167: GEN Dn,cy,p1;
168:
169: lb=lg(D); Dn=cgetg(lb,t_VEC);
170: for (i=1; i<lb; i++)
171: {
172: cy=(GEN)D[i]; lcy=lg(cy);
173: Dn[i]=lgetg(lcy,t_VECSMALL); p1=(GEN)Dn[i];
174: for (j=1; j<lcy; j++) p1[j] = perm[cy[j]];
175: }
176: return Dn;
177: }
178:
179: /* cy is a cycle; recturn cy(a) */
180: static long
181: im_by_cy(long a,GEN cy)
182: {
183: long k, l = lg(cy);
184:
185: k=1; while (k<l && cy[k] != a) k++;
186: if (k == l) return a;
187: k++; if (k == l) k = 1;
188: return cy[k];
189: }
190:
191: /* renvoie 0 si l'un des coefficients de g[i] est de module > M[i]; 1 sinon */
192: static long
193: ok_coeffs(GEN g,GEN M)
194: {
195: long i, lg = lgef(g)-1; /* g is monic, and cst term is ok */
196: for (i=3; i<lg; i++)
197: if (absi_cmp((GEN)g[i], (GEN)M[i]) > 0) return 0;
198: return 1;
199: }
200:
201: /* vbs[0] = current cardinality+1, vbs[1] = max number of elts */
202: static GEN
203: append_vbs(GEN vbs, GEN D)
204: {
205: long l,maxl,i,j,n, lD = lg(D);
206: GEN Dn, last;
207:
208: n = 0; for (i=1; i<lD; i++) n += lg(D[i]);
209: Dn = (GEN)gpmalloc((lD + n) * sizeof(long));
210: last = Dn + lD; Dn[0] = D[0];
211: for (i=1; i<lD; i++)
212: {
213: GEN cy = (GEN)D[i], cn = last;
214: for (j=0; j<lg(cy); j++) cn[j] = cy[j];
215: Dn[i] = (long)cn; last = cn + j;
216: }
217:
218: if (!vbs)
219: {
220: maxl = 1024;
221: vbs = (GEN)gpmalloc((2 + maxl)*sizeof(GEN));
222: *vbs = maxl; vbs++; setlg(vbs, 1);
223: }
224: l = lg(vbs); maxl = vbs[-1];
225: if (l == maxl)
226: {
227: vbs = (GEN)gprealloc((void*)(vbs-1), (2 + (maxl<<1))*sizeof(GEN),
228: (2 + maxl)*sizeof(GEN));
229: *vbs = maxl<<1; vbs++; setlg(vbs, 1);
230: }
231: if (DEBUGLEVEL>5) fprintferr("appending D = %Z\n",D);
232: vbs[l] = (long)Dn; setlg(vbs, l+1); return vbs;
233: }
234:
235: GEN
236: myconcat(GEN D, long a)
237: {
238: long i,l = lg(D);
239: GEN x = cgetg(l+1,t_VECSMALL);
240: for (i=1; i<l; i++) x[i]=D[i];
241: x[l] = a; return x;
242: }
243:
244: void
245: myconcat2(GEN D, GEN a)
246: {
247: long i,l = lg(D), m = lg(a);
248: GEN x = D + (l-1);
249: for (i=1; i<m; i++) x[i]=a[i];
250: setlg(D, l+m-1);
251: }
252:
253: static GEN
254: print_block_system(long N,GEN Y,long d, GEN vbs)
255: {
256: long a,i,j,l,ll,*k,*n,lp,**e,u,v,*t,ns, r = lg(Y);
257: GEN D,De,Z,cyperm,perm,p1,empty;
258:
259: if (DEBUGLEVEL>5) fprintferr("Y = %Z\n",Y);
260: empty = cgetg(1,t_VEC);
261: n = new_chunk(N+1);
262: D = cgetg(N+1, t_VEC); setlg(D,1);
263: t=new_chunk(r+1); k=new_chunk(r+1); Z=cgetg(r+1,t_VEC);
264: for (ns=0,i=1; i<r; i++)
265: {
266: GEN Yi = (GEN)Y[i], cy;
267: long ki = 0, si = lg(Yi)-1;
268:
269: for (j=1; j<=si; j++) { n[j]=lg(Yi[j])-1; ki += n[j]; }
270: ki /= d;
271: De=cgetg(ki+1,t_VEC);
272: for (j=1; j<=ki; j++) De[j]=(long)empty;
273: for (j=1; j<=si; j++)
274: {
275: a = mael(Yi,j,1); cy = (GEN)Yi[j];
276: for (l=1,lp=0; l<=n[j]; l++)
277: {
278: lp++; if (lp>ki) lp = 1;
279: a = im_by_cy(a, cy);
280: De[lp] = (long)myconcat((GEN)De[lp], a);
281: }
282: }
283: if (si>1 && ki>1)
284: {
285: ns++; t[ns]=si-1; k[ns]=ki;
286: Z[ns]=lgetg(si,t_VEC); p1=(GEN)Z[ns];
287: for (j=2; j<=si; j++) p1[j-1]=Yi[j];
288: }
289: myconcat2(D,De);
290: }
291: if (DEBUGLEVEL>2) { fprintferr("\nns = %ld\n",ns); flusherr(); }
292: if (!ns) return append_vbs(vbs,D);
293:
294: setlg(Z, ns+1);
295: e=(long**)new_chunk(ns+1);
296: for (i=1; i<=ns; i++)
297: {
298: e[i]=new_chunk(t[i]+1);
299: for (j=1; j<=t[i]; j++) e[i][j]=0;
300: }
301: cyperm = cgetg(N+1,t_VEC);
302: perm = cgetg(N+1,t_VEC); i=ns;
303: do
304: {
305: long av = avma;
306: if (DEBUGLEVEL>5)
307: {
308: for (l=1; l<=ns; l++)
309: {
310: for (ll=1; ll<=t[l]; ll++)
311: fprintferr("e[%ld][%ld] = %ld, ",l,ll,e[l][ll]);
312: fprintferr("\n");
313: }
314: fprintferr("\n"); flusherr();
315: }
316: for (u=1; u<=N; u++) perm[u]=u;
317: for (u=1; u<=ns; u++)
318: for (v=1; v<=t[u]; v++)
319: perm_mul(perm, cycle_power_to_perm(cyperm,gmael(Z,u,v),e[u][v]));
320: vbs = append_vbs(vbs, im_block_by_perm(D,perm));
321: avma = av;
322:
323: e[ns][t[ns]]++;
324: if (e[ns][t[ns]] >= k[ns])
325: {
326: j=t[ns]-1;
327: while (j>=1 && e[ns][j] == k[ns]-1) j--;
328: if (j>=1) { e[ns][j]++; for (l=j+1; l<=t[ns]; l++) e[ns][l]=0; }
329: else
330: {
331: i=ns-1;
332: while (i>=1)
333: {
334: j=t[i];
335: while (j>=1 && e[i][j] == k[i]-1) j--;
336: if (j<1) i--;
337: else
338: {
339: e[i][j]++;
340: for (l=j+1; l<=t[i]; l++) e[i][l]=0;
341: for (ll=i+1; ll<=ns; ll++)
342: for (l=1; l<=t[ll]; l++) e[ll][l]=0;
343: break;
344: }
345: }
346: }
347: }
348: }
349: while (i>0);
350: return vbs;
351: }
352:
353: /* rend le numero du cycle (a1,...,an) dans le support duquel se trouve a */
354: /* met dans *pt l'indice i tq ai = a */
355: static long
356: in_what_cycle(long a,GEN cys,long *pt)
357: {
358: long i,k,nk, lcys=lg(cys);
359:
360: for (k=1; k<lcys; k++)
361: {
362: GEN c = (GEN)cys[k]; nk = lg(c);
363: for (i=1; i<nk; i++)
364: if (a == c[i]) { *pt = i; return k; }
365: }
366: err(talker,"impossible to find %d in in_what_cycle",a);
367: return 0; /* not reached */
368: }
369:
370: /* Return common factors to A and B + s the prime to A part of B */
371: static GEN
372: commonfactor(GEN A, GEN B)
373: {
374: GEN f,p1,p2,s, y = cgetg(3,t_MAT);
375: long lf,i;
376:
377: s = absi(B); f=(GEN)A[1]; lf=lg(f);
378: p1=cgetg(lf+1,t_COL); y[1]=(long)p1;
379: p2=cgetg(lf+1,t_COL); y[2]=(long)p2;
380: for (i=1; i<lf; i++)
381: {
382: p1[i] = f[i];
383: p2[i] = lstoi(pvaluation(s,(GEN)f[i], &s));
384: }
385: p1[i] = (long)s;
386: p2[i] = un; return y;
387: }
388:
389: static void
390: polsimplify(GEN x)
391: {
392: long i,lx = lgef(x);
393: for (i=2; i<lx; i++)
394: if (typ(x[i]) == t_POL) x[i] = signe(x[i])? mael(x,i,2): zero;
395: }
396:
397: /* Renvoie un polynome g definissant un sous-corps potentiel, ou
398: * 0: si le polynome trouve n'est pas separable,
399: * 1: si les coefficients du polynome trouve sont plus grands que la borne M,
400: * 2: si p divise le discriminant de g,
401: * 3: si le discriminant de g est nul,
402: * 4: si la partie s de d(g) premiere avec d(L) n'est pas un carre,
403: * 5: si s est un carre et si un des facteurs premiers communs a d(g) et d(L)
404: * a un exposant impair dans d(g) et un exposant plus petit que d dans d(L),
405: * 6: si le discriminant du corps defini par g a la puissance d ne divise pas
406: * le discriminant du corps nf (soit L).
407: */
408: static GEN
409: cand_for_subfields(GEN A,GEN DATA,GEN *ptdelta,GEN *ptrootsA)
410: {
411: long av=avma,N,m,i,j,d,lf;
412: GEN P,pe,p,pol,cys,tabroots,delta,g,dg,unmodpe,tabrA;
413: GEN factcommon,ff1,ff2,p1;
414: GEN *gptr[3];
415:
416: pol=(GEN)DATA[1]; N=lgef(pol)-3; m=lg(A)-1; d=N/m;
417: if (N%m) err(talker,"incompatible block system in cand_for_subfields");
418: p = (GEN)DATA[2];
419: cys=(GEN)DATA[5];
420: tabroots=(GEN)DATA[10];
421: pe = gclone((GEN)DATA[9]);
422: unmodpe = cgetg(3,t_INTMOD); unmodpe[1]=(long)pe; unmodpe[2]=un;
423:
424: delta = cgetg(m+1,t_VEC);
425: tabrA = cgetg(m+1,t_VEC);
426: for (i=1; i<=m; i++)
427: {
428: GEN Ai=(GEN)A[i], col = cgetg(d+1,t_VEC);
429: long l,k;
430:
431: tabrA[i]=(long)col; p1 = unmodpe;
432: for (j=1; j<=d; j++)
433: {
434: l=in_what_cycle(Ai[j],cys,&k);
435: col[j] = mael(tabroots, l, k);
436: p1 = gmul(p1, (GEN)col[j]);
437: }
438: p1 = lift_intern((GEN)p1[2]);
439: for (j=1; j<i; j++)
440: if (gegal(p1,(GEN)delta[j])) { avma=av; return gzero; }
441: if (DEBUGLEVEL>2) fprintferr("delta[%ld] = %Z\n",i,p1);
442: delta[i] = (long)p1;
443: }
444: P = gmael3(tabroots,1,1,1);
445: for (i=1; i<=m; i++)
446: {
447: p1 = cgetg(3,t_POLMOD); p1[1]=(long)P; p1[2]=delta[i];
448: delta[i] = (long)p1;
449: }
450: g = roots_to_pol(gmul(unmodpe,delta),0);
451: g=centerlift(lift_intern(g)); polsimplify(g);
452: if (DEBUGLEVEL>2) fprintferr("pol. found = %Z\n",g);
453: if (!ok_coeffs(g,(GEN)DATA[8])) return gun;
454: dg=discsr(g);
455: if (!signe(dg)) return stoi(3);
456: if (!signe(resii(dg,p))) return gdeux;
457: factcommon=commonfactor(FACTORDL,dg);
458: ff1=(GEN)factcommon[1]; lf=lg(ff1)-1;
459: if (!carreparfait((GEN)ff1[lf])) return stoi(4);
460: ff2=(GEN)factcommon[2];
461: for (i=1; i<lf; i++)
462: if (mod2((GEN)ff2[i]) && itos(gmael(FACTORDL,2,i)) < d) return stoi(5);
463: gunclone(pe);
464:
465: *ptdelta=delta; *ptrootsA=tabrA;
466: gptr[0]=&g; gptr[1]=ptdelta; gptr[2]=ptrootsA;
467: gerepilemany(av,gptr,3); return g;
468: }
469:
470: /* a partir d'un polynome h(x) dont les coefficients sont definis mod p^k,
471: * on construit un polynome a coefficients dans Q dont les coefficients ont
472: * pour approximation p-adique les coefficients de h */
473: static GEN
474: retrieve_p_adique_polynomial_in_Q(GEN ind,GEN h)
475: {
476: return gdiv(centerlift(gmul(h,ind)), ind);
477: }
478:
479: /* interpolation polynomial P(x) s.t P(T[j][i]) = delta[i] mod p */
480: static GEN
481: interpolation_polynomial(GEN T, GEN delta)
482: {
483: long i,j,i1,j1, m = lg(T), d = lg(T[1]);
484: GEN P = NULL, x0 = gneg(polx[0]);
485:
486: for (j=1; j<m; j++)
487: {
488: GEN p3 = NULL;
489: for (i=1; i<d; i++)
490: {
491: GEN p1=gun, p2=gun, a = gneg(gmael(T,j,i));
492: for (j1=1; j1<m; j1++)
493: for (i1=1; i1<d; i1++)
494: if (i1 != i || j1 != j)
495: {
496: p1 = gmul(p1,gadd(gmael(T,j1,i1), x0));
497: p2 = gmul(p2,gadd(gmael(T,j1,i1), a));
498: }
499: p1 = gdiv(p1,p2);
500: p3 = p3? gadd(p3, p1): p1;
501: }
502: p3 = gmul((GEN)delta[j],p3);
503: P = P? gadd(P,p3): p3;
504: }
505: return P;
506: }
507:
508: /* nf est le corps de nombres, g un polynome de Z[x] candidat
509: * pour definir un sous-corps, p le nombre premier ayant servi a definir le
510: * potential block system rootsA donne par les racines avec une approximation
511: * convenable, e est la precision p-adique des elements de rootsA et delta la
512: * liste des racines de g dans une extension convenable en precision p^e.
513: * Renvoie un polynome h de Q[x] tel que f divise g o h et donc tel que le
514: * couple (g,h) definisse un sous-corps, ou bien gzero si rootsA n'est pas un
515: * block system
516: */
517: static GEN
518: embedding_of_potential_subfields(GEN nf,GEN g,GEN DATA,GEN rootsA,GEN delta)
519: {
520: GEN w0_inQ,w0,w1,h0,gp,p2,f,unmodp,p,ind, maxp;
521: long av = avma, av1;
522:
523: f=(GEN)nf[1]; ind=(GEN)nf[4]; p=(GEN)DATA[2];
524: maxp=mulii((GEN)DATA[11],ind);
525: gp=deriv(g,varn(g)); unmodp=gmodulsg(1,p);
526: av1 = avma;
527: w0 = interpolation_polynomial(gmul(rootsA,unmodp), delta);
528: w0 = lift_intern(w0); /* in Fp[x] */
529: polsimplify(w0);
530: w0_inQ = retrieve_p_adique_polynomial_in_Q(ind,w0);
531: (void)gbezout(poleval(gp,w0), gmul(unmodp,f), &h0, &p2);
532: w0 = lift_intern(w0); /* in Z[x] */
533: h0 = lift_intern(lift_intern(h0));
534: for(;;)
535: {
536: GEN p1;
537: /* Given g in Z[x], gp its derivative, p a prime, [w0,h0] in Z[x] s.t.
538: * h0(x).gp(w0(x)) = 1 and g(w0(x)) = 0 (mod f,mod p), return
539: * [w1,h1] satisfying the same condition mod p^2. Moreover,
540: * [w1,h1] = [w0,h0] (mod p)
541: * (cf. Dixon: J. Austral. Math. Soc., Series A, vol.49, 1990, p.445) */
542: if (DEBUGLEVEL>2)
543: {
544: fprintferr("w = "); outerr(w0);
545: fprintferr("h = "); outerr(h0);
546: }
547: p = sqri(p); unmodp[1] = (long)p;
548: p1 = gneg(gmul(h0, poleval(g,w0)));
549: w1 = gres(gmul(unmodp,gadd(w0,p1)), f);
550: p2 = retrieve_p_adique_polynomial_in_Q(ind,w1);
551: if ((gegal(p2, w0_inQ) || cmpii(p,maxp)) && gdivise(poleval(g,p2), f))
552: return gerepileupto(av, poleval(p2, gadd(polx[0],stoi(TR))));
553: if (DEBUGLEVEL>2)
554: {
555: fprintferr("Old Q-polynomial: "); outerr(w0_inQ);
556: fprintferr("New Q-polynomial: "); outerr(p2);
557: }
558: if (cmpii(p, maxp) > 0)
559: {
560: if (DEBUGLEVEL) fprintferr("coeff too big for embedding\n");
561: avma=av; return gzero;
562: }
563:
564: w1 = lift_intern(w1);
565: p1 = gneg(gmul(h0, poleval(gp,w1)));
566: p1 = gmul(h0, gadd(gdeux,p1));
567: h0 = lift_intern(gres(gmul(unmodp,p1), f));
568: w0 = w1; w0_inQ = p2;
569: {
570: GEN *gptr[4]; gptr[0]=&w0; gptr[1]=&h0; gptr[2]=&w0_inQ; gptr[3]=&p;
571: gerepilemany(av1,gptr,4);
572: }
573: }
574: }
575:
576: static long
577: choose_prime(GEN pol,GEN dpol,long d,GEN *ptff,GEN *ptlistpotbl, long *ptnn)
578: {
579: long j,k,oldllist,llist,r,nn,oldnn,*n,N,pp;
580: GEN p,listpotbl,oldlistpotbl,ff,oldff,p3;
581: byteptr di=diffptr;
582:
583: if (DEBUGLEVEL) timer2();
584: di++; p = stoi(2); N = lgef(pol)-3;
585: while (p[2]<=N) p[2] += *di++;
586: oldllist = oldnn = BIGINT;
587: n = new_chunk(N+1);
588: for(k=1; k<11 || oldnn == BIGINT; k++,p[2]+= *di++)
589: {
590: long av=avma;
591: while (!smodis(dpol,p[2])) p[2] += *di++;
592: ff=(GEN)factmod(pol,p)[1]; r=lg(ff)-1;
593: if (r>1 && r<N)
594: {
595: for (j=1; j<=r; j++) n[j]=lgef(ff[j])-3;
596: p3 = stoi(n[1]);
597: for (j=2; j<=r; j++) p3 = glcm(p3,stoi(n[j]));
598: nn=itos(p3);
599: if (nn > oldnn)
600: {
601: if (DEBUGLEVEL)
602: {
603: fprintferr("p = %ld,\tr = %ld,\tnn = %ld,\t#pbs = skipped\n",
604: p[2],r,nn);
605: }
606: continue;
607: }
608: listpotbl=potential_block_systems(N,d,ff,n);
609: if (!listpotbl) { oldlistpotbl = NULL; pp = p[2]; break; }
610: llist=lg(listpotbl)-1;
611: if (DEBUGLEVEL)
612: {
613: fprintferr("Time: %ldms,\tp = %ld,\tr = %ld,\tnn = %ld,\t#pbs = %ld\n",
614: timer2(),p[2],r,nn,llist);
615: flusherr();
616: }
617: if (nn<oldnn || llist<oldllist)
618: {
619: oldllist=llist; oldlistpotbl=listpotbl;
620: pp=p[2]; oldff=ff; oldnn=nn; continue;
621: }
622: for (j=1; j<llist; j++) free((void*)listpotbl[j]);
623: free((void*)(listpotbl-1));
624: }
625: avma = av;
626: }
627: if (DEBUGLEVEL)
628: {
629: fprintferr("Chosen prime: p = %ld\n",pp);
630: if (DEBUGLEVEL>2)
631: fprintferr("List of potential block systems of size %ld: %Z\n",
632: d,oldlistpotbl);
633: flusherr();
634: }
635: *ptlistpotbl=oldlistpotbl; *ptff=oldff; *ptnn=oldnn; return pp;
636: }
637:
638: static GEN
639: change_pol(GEN nf, GEN ff)
640: {
641: long i,l;
642: GEN pol = (GEN)nf[1], p1 = gsub(polx[0],gun);
643:
644: TR++; pol=poleval(pol, p1);
645: nf = dummycopy(nf);
646: nf[6] = (long)dummycopy((GEN)nf[6]);
647: l=lg(ff);
648: for (i=1; i<l; i++) ff[i]=(long)poleval((GEN)ff[i], p1);
649: l=lg(nf[6]); p1=(GEN)nf[6];
650: for (i=1; i<l; i++) p1[i]=ladd(gun,(GEN)p1[i]);
651: nf[1]=(long)pol; return nf;
652: }
653:
654: static GEN
655: bound_for_coeff(long m,GEN rr,long r1, GEN *maxroot)
656: {
657: long i, lrr=lg(rr);
658: GEN p1,b1,b2,B,M, C = matpascal(m-1);
659:
660: rr = gabs(rr,DEFAULTPREC); *maxroot = vecmax(rr);
661: for (i=1; i<lrr; i++)
662: if (gcmp((GEN)rr[i], gun) < 0) rr[i] = un;
663: for (b1=gun,i=1; i<=r1; i++) b1 = gmul(b1, (GEN)rr[i]);
664: for (b2=gun ; i<lrr; i++) b2 = gmul(b2, (GEN)rr[i]);
665: B = gmul(b1, gsqr(b2));
666: M = cgetg(m+2, t_VEC); M[1]=M[2]=zero; /* unused */
667: for (i=1; i<m; i++)
668: {
669: p1 = gadd(gmul(gcoeff(C, m, i), B),
670: gcoeff(C, m, i+1));
671: M[i+2] = lceil(p1);
672: }
673: return M;
674: }
675:
676: /* liste des sous corps de degre d du corps de nombres nf */
677: static GEN
678: subfields_of_given_degree(GEN nf,GEN dpol,long d)
679: {
680: long av,av1,av2,tetpil,pp,llist,i,nn,N;
681: GEN listpotbl,ff,A,delta,rootsA,CSF,ESF,p1,p2,LSB;
682: GEN DATA, pol = (GEN)nf[1];
683:
684: av=avma;
685: N = lgef(pol)-3;
686: pp=choose_prime(pol,dpol,N/d,&ff,&listpotbl,&nn);
687: if (!listpotbl) { avma=av; return cgetg(1,t_VEC); }
688: llist=lg(listpotbl);
689: LAB0:
690: av1=avma; LSB=cgetg(1,t_VEC);
691: DATA=compute_data(nf,ff,stoi(pp),d,nn);
692: for (i=1; i<llist; i++)
693: {
694: av2=avma; A=(GEN)listpotbl[i];
695: if (DEBUGLEVEL > 1)
696: fprintferr("\n* Potential block # %ld: %Z\n",i,A);
697: CSF=cand_for_subfields(A,DATA,&delta,&rootsA);
698: if (typ(CSF)==t_INT)
699: {
700: if (DEBUGLEVEL > 1) switch(itos(CSF))
701: {
702: case 0: fprintferr("changing f(x): non separable g(x)\n"); break;
703: case 1: fprintferr("coeff too big for pol g(x)\n"); break;
704: case 2: fprintferr("changing f(x): p divides disc(g(x))\n"); break;
705: case 3: fprintferr("non irreducible polynomial g(x)\n"); break;
706: case 4: fprintferr("prime to d(L) part of d(g) not a square\n"); break;
707: case 5: fprintferr("too small exponent of a prime factor in d(L)\n"); break;
708: case 6: fprintferr("the d-th power of d(K) does not divide d(L)\n");
709: }
710: switch(itos(CSF))
711: {
712: case 0: case 2:
713: avma=av1; nf = change_pol(nf,ff); pol = (GEN)nf[1];
714: if (DEBUGLEVEL) fprintferr("new f = %Z\n",pol);
715: goto LAB0;
716: }
717: avma=av2;
718: }
719: else
720: {
721: if (DEBUGLEVEL) fprintferr("candidate = %Z\n",CSF);
722: ESF=embedding_of_potential_subfields(nf,CSF,DATA,rootsA,delta);
723: if (ESF == gzero) avma=av2;
724: else
725: {
726: if (DEBUGLEVEL) fprintferr("embedding = %Z\n",ESF);
727: p1=cgetg(3,t_VEC); p2=cgetg(2,t_VEC); p2[1]=(long)p1;
728: p1[1]=(long)CSF;
729: p1[2]=(long)ESF; tetpil=avma;
730: LSB=gerepile(av2,tetpil, concat(LSB,p2));
731: }
732: }
733: }
734: for (i=1; i<llist; i++) free((void*)listpotbl[i]);
735: free((void*)(listpotbl-1)); tetpil=avma;
736: return gerepile(av,tetpil,gcopy(LSB));
737: }
738:
739: GEN
740: subfields(GEN nf,GEN d)
741: {
742: long av=avma,di,N,v0,lp1,i;
743: GEN dpol,p1,LSB,p2,pol;
744:
745: nf=checknf(nf); pol = (GEN)nf[1];
746: v0=varn(pol); N=lgef(pol)-3; di=itos(d);
747: if (di==N)
748: {
749: LSB=cgetg(2,t_VEC); p1=cgetg(3,t_VEC); LSB[1]=(long)p1;
750: p1[1]=lcopy(pol); p1[2]=lpolx[v0]; return LSB;
751: }
752: if (di==1)
753: {
754: LSB=cgetg(2,t_VEC); p1=cgetg(3,t_VEC); LSB[1]=(long)p1;
755: p1[1]=lpolx[v0]; p1[2]=lcopy(pol); return LSB;
756: }
757: if (di<=0 || di>N || N%di) return cgetg(1,t_VEC);
758:
759: TR=0; dpol=mulii((GEN)nf[3],sqri((GEN)nf[4]));
760: if (v0) nf=gsubst(nf,v0,polx[0]);
761: FACTORDL=factor(absi((GEN)nf[3]));
762: p1=subfields_of_given_degree(nf,dpol,di); lp1=lg(p1)-1;
763: if (v0)
764: for (i=1; i<=lp1; i++)
765: { p2=(GEN)p1[i]; setvarn(p2[1],v0); setvarn(p2[2],v0); }
766: return gerepileupto(av,p1);
767: }
768:
769: static GEN
770: subfieldsall(GEN nf)
771: {
772: long av=avma,av1,N,ld,d,i,j,lNLSB,v0,lp1;
773: GEN pol,dpol,dg,LSB,NLSB,p1,p2;
774:
775: nf=checknf(nf); pol = (GEN)nf[1];
776: v0=varn(pol); N=lgef(pol)-3;
777: if (isprime(stoi(N)))
778: {
779: avma=av; LSB=cgetg(3,t_VEC);
780: LSB[1]=lgetg(3,t_VEC); LSB[2]=lgetg(3,t_VEC);
781: p1=(GEN)LSB[1]; p1[1]=lcopy(pol); p1[2]=lpolx[v0];
782: p2=(GEN)LSB[2]; p2[1]=p1[2]; p2[2]=p1[1];
783: return LSB;
784: }
785: FACTORDL=factor(absi((GEN)nf[3])); dg=divisors(stoi(N));
786: dpol=mulii(sqri((GEN)nf[4]),(GEN)nf[3]);
787: if (DEBUGLEVEL>0)
788: {
789: fprintferr("\n***** Entering subfields\n\n");
790: fprintferr("pol = "); outerr(pol);
791: fprintferr("dpol = "); outerr(dpol);
792: fprintferr("divisors = "); outerr(dg);
793: }
794: ld=lg(dg)-1; LSB=cgetg(2,t_VEC); LSB[1]=lgetg(3,t_VEC);
795: p1=(GEN)LSB[1]; p1[1]=(long)pol; p1[2]=(long)polx[0];
796: if (v0) nf=gsubst(nf,v0,polx[0]);
797: for (i=2; i<ld; i++)
798: {
799: TR=0; av1=avma; d=itos((GEN)dg[i]);
800: if (DEBUGLEVEL>0)
801: {
802: fprintferr("\n*** Looking for subfields of degree %ld\n\n",N/d);
803: flusherr();
804: }
805: NLSB=subfields_of_given_degree(nf,dpol,N/d);
806: if (DEBUGLEVEL)
807: {
808: fprintferr("\nSubfields of degree %ld:\n",N/d);
809: lNLSB=lg(NLSB)-1; for (j=1; j<=lNLSB; j++) outerr((GEN)NLSB[j]);
810: }
811: if (lg(NLSB)>1) LSB = concatsp(LSB,NLSB); else avma=av1;
812: }
813: p1=cgetg(2,t_VEC); p1[1]=lgetg(3,t_VEC);p2=(GEN)p1[1];
814: p2[1]=(long)polx[0]; p2[2]=(long)pol;
815: LSB=concatsp(LSB,p1); lp1=lg(LSB)-1;
816: LSB = gerepileupto(av, gcopy(LSB));
817: if (v0)
818: for (i=1; i<=lp1; i++)
819: { p2=(GEN)LSB[i]; setvarn(p2[1],v0); setvarn(p2[2],v0); }
820: if (DEBUGLEVEL>0) fprintferr("\n***** Leaving subfields\n\n");
821: return LSB;
822: }
823:
824: GEN
825: subfields0(GEN nf,GEN d)
826: {
827: return d? subfields(nf,d): subfieldsall(nf);
828: }
829:
830: /* irreducible (unitary) polynomial of degree n over Fp[v] */
831: GEN
832: ffinit(GEN p,long n,long v)
833: {
834: long av,av1,tetpil,i,*a,j,l,pp;
835: GEN pol,fpol;
836:
837: if (n<=0) err(talker,"non positive degree in ffinit");
838: if (is_bigint(p)) err(talker,"prime field too big in ffinit");
839: if (v<0) v = 0;
840: av=avma; pp=itos(p); pol = cgetg(n+3,t_POL);
841: pol[1] = evalsigne(1)|evalvarn(v)|evallgef(n+3);
842: a=new_chunk(n+2);
843: a[1]=1; for (i=2; i<=n+1; i++) a[i]=0;
844: pol[n+2]=un; av1=avma;
845: for(;;)
846: {
847: a[n+1]++;
848: if (a[n+1]>=pp)
849: {
850: j=n; while (j>=2 && a[j]==pp-1) j--;
851: if (j>=2) { a[j]++; for (l=j+1; l<=n+1; l++) a[l]=0; }
852: }
853: for (i=2; i<=n+1; i++) pol[i]=lstoi(a[n+3-i]);
854: fpol=simplefactmod(pol,p);
855: if (lg(fpol[1])==2 && gcmp1(gmael(fpol,2,1))) break;
856: avma=av1;
857: }
858: tetpil=avma; return gerepile(av,tetpil,Fp_pol(pol,p));
859: }
860:
861: static GEN
862: lift_coeff(GEN x, GEN fq)
863: {
864: GEN r;
865: if (typ(x) == t_POLMOD) { r = x; x = (GEN)x[2]; }
866: else r = cgetg(3,t_POLMOD);
867: r[1]=(long)fq; r[2]=(long)lift_intern(x); return r;
868: }
869:
870: /* a is a polynomial whose coeffs are in Fq (= (Z/p)[y] / (fqbar), where
871: * fqbar is the reduction of fq mod p).
872: * Lift _in place_ the coeffs so that they belong to Z[y] / (fq)
873: */
874: static GEN
875: special_lift(GEN a,GEN fq)
876: {
877: long la,i;
878: GEN c;
879:
880: if (typ(a)==t_POL)
881: {
882: la=lgef(a); c=cgetg(la,t_POL); c[1]=a[1];
883: for (i=2; i<la; i++) c[i]=(long)lift_coeff((GEN)a[i],fq);
884: return c;
885: }
886: return lift_coeff(a,fq);
887: }
888:
889: /* Hensel lift: fk = vector of factors of pol (unramified) in finite field
890: * Fp / fkk. Lift it to the precision p^e. This is equivalent to working
891: * in precision pi^e in the unramified extension of Qp given by fkk.
892: */
893: GEN
894: hensel_lift(GEN pol,GEN fk,GEN fkk,GEN p,long e)
895: {
896: long av = avma, i, r = lg(fk)-1;
897: GEN p1,A,B,C,R,U,V,fklift,fklift2,fk2;
898: GEN unmodp = gmodulsg(1,p), fq = lift(fkk);
899:
900: fk2=cgetg(r+1,t_VEC);
901: fklift=cgetg(r+1,t_VEC);
902: fklift2=cgetg(r+1,t_VEC);
903: fk2[r] = fklift2[r] = un;
904: for (i=r; i>1; i--)
905: {
906: fk2[i-1] = lmul((GEN)fk2[i],(GEN)fk[i]);
907: fklift[i] = (long)special_lift(gcopy((GEN)fk[i]),fq);
908: fklift2[i-1] = lmul((GEN)fklift2[i],(GEN)fklift[i]);
909: }
910: fklift[1] = (long)special_lift(gcopy((GEN)fk[1]),fq);
911: R=cgetg(r+1,t_VEC); C=pol;
912: for (i=1; i<r; i++)
913: { /* treat factors two by two: fk[i] and fk2[i] = product fk[i+1..] */
914: long av1 = avma,tetpil1, ex = 1;
915: GEN pp;
916:
917: (void)gbezout((GEN)fk[i],(GEN)fk2[i],&U,&V);
918: A = (GEN)fklift[i]; U = special_lift(U,fq);
919: B = (GEN)fklift2[i]; V = special_lift(V,fq);
920: for (pp=p;; pp=sqri(pp))
921: { /* Algorithm 3.5.[5,6] H. Cohen page 137 (1995) */
922: GEN f,t,A0,B0,U0,V0;
923:
924: unmodp[1] = (long)pp;
925: p1 = gneg_i(gmul(A,B));
926: p1=gdiv(gadd(C,p1),pp);
927: f=gmul(p1,unmodp);
928: t=poldivres(gmul(V,f),A, &A0);
929: A0=special_lift(A0,fq);
930: B0=special_lift(gadd(gmul(U,f),gmul(B,t)),fq);
931: A0 = gmul(A0,pp);
932: B0 = gmul(B0,pp); tetpil1 = avma;
933: A = gadd(A, A0);
934: B = gadd(B, B0); ex <<= 1;
935: if (ex>=e)
936: {
937: GEN *gptr[2]; gptr[0]=&A; gptr[1]=&B;
938: gerepilemanysp(av1,tetpil1,gptr,2);
939: C = B; R[i] = (long)A; break;
940: }
941: p1 = gneg_i(gadd(gmul(U,A),gmul(V,B)));
942: p1=gdiv(gadd(gun,p1),pp);
943: f=gmul(p1,unmodp);
944: t=poldivres(gmul(V,f),A, &V0);
945: U0=special_lift(gadd(gmul(U,f),gmul(B,t)),fq);
946: V0=special_lift(V0,fq);
947: U = gadd(U, gmul(U0,pp));
948: V = gadd(V, gmul(V0,pp));
949: }
950: }
951: if (r==1) C = gcopy(C);
952: R[r] = (long)C; return gerepileupto(av,R);
953: }
954:
955: /* etant donne nf et p et la factorisation de nf[1] mod p, et le degre m des
956: * sous corps cherches, cree un vecteur ligne a 13 composantes:
957: * 1 : le polynome nf[1],
958: * 2 : le premier p,
959: * 3 : la factorisation ff,
960: * 4 : la longeur des cycles associes (n_1,...,n_r),
961: * 5 : les cycles associes,
962: * 6 : le corps F_(p^q),
963: * 7 : les racines de f dans F_(p^q) par facteur de ff,
964: * 8 : la borne M pour les sous-corps,
965: * 9 : l'exposant e telle que la precision des lifts soit p^e>2.M,
966: * 10: le lift de Hensel a la precision p^e de la factorisation en facteurs
967: * lineaires de nf[1] dans F_(p^q),
968: * 11: la borne de Hadamard pour les coefficients de h(x) tel que g o h = 0
969: * mod nf[1].
970: * ces donnees sont valides pour nf, p et m (d) donnes...
971: */
972: static GEN
973: compute_data(GEN nf, GEN ff, GEN p, long m, long nn)
974: {
975: long i,j,l,r,*n,e,N,pp,d,r1;
976: GEN DATA,p1,p2,cys,fhk,tabroots,MM,fk,dpol,maxroot,maxMM,pol;
977:
978: if (DEBUGLEVEL>1) { fprintferr("Entering compute_data()\n\n"); flusherr(); }
979: pol = (GEN)nf[1]; N = lgef(pol)-3;
980: DATA=cgetg(14,t_VEC);
981: DATA[1]=(long)pol;
982: DATA[2]=(long)p; r=lg(ff)-1;
983: DATA[3]=(long)ff;
984: n = cgetg(r+1, t_VECSMALL);
985: DATA[4]= (long)n;
986: for (j=1; j<=r; j++) n[j]=lgef(ff[j])-3;
987: cys=cgetg(r+1,t_VEC); l=0;
988: for (i=1; i<=r; i++)
989: {
990: p1 = cgetg(n[i]+1, t_VECSMALL);
991: cys[i] = (long)p1; for (j=1; j<=n[i]; j++) p1[j]=++l;
992: }
993: DATA[5]=(long)cys;
994: DATA[6]=(long)ffinit(p,nn,MAXVARN);
995: tabroots=cgetg(r+1,t_VEC);
996: for (j=1; j<=r; j++)
997: {
998: p1=(GEN)factmod9((GEN)ff[j],p,(GEN)DATA[6])[1];
999: p2=cgetg(n[j]+1,t_VEC); tabroots[j]=(long)p2;
1000: p2[1]=lneg(gmael(p1,1,2));
1001: for (i=2; i<=n[j]; i++) p2[i]=(long)powgi((GEN)p2[i-1],p);
1002: }
1003: DATA[7]=(long)tabroots;
1004: r1=itos(gmael(nf,2,1));
1005: MM = bound_for_coeff(m, (GEN)nf[6], r1, &maxroot);
1006: MM = gmul2n(MM,1);
1007: DATA[8]=(long)MM;
1008: pp=itos(p); maxMM = vecmax(MM);
1009: for (e=1,p1=p; cmpii(p1, maxMM) < 0; ) { p1 = mulis(p1,pp); e++; }
1010: DATA[9]=lpuigs(p,e); fk=cgetg(N+1,t_VEC);
1011: for (l=1,j=1; j<=r; j++)
1012: for (i=1; i<=n[j]; i++)
1013: fk[l++] = lsub(polx[0],gmael(tabroots,j,i));
1014: fhk = hensel_lift(pol,fk,(GEN)DATA[6],p,e);
1015: tabroots=cgetg(r+1,t_VEC);
1016: for (l=1,j=1; j<=r; j++)
1017: {
1018: p1 = cgetg(n[j]+1,t_VEC); tabroots[j]=(long)p1;
1019: for (i=1; i<=n[j]; i++,l++) p1[i] = lneg(gmael(fhk,l,2));
1020: }
1021: DATA[10]=(long)tabroots;
1022:
1023: d=N/m; p1=gmul(stoi(N), gsqrt(gpuigs(stoi(N-1),N-1),DEFAULTPREC));
1024: p2 = gpuigs(maxroot, d + N*(N-1)/2);
1025: dpol=mulii(sqri((GEN)nf[4]),(GEN)nf[3]);
1026: p1 = gdiv(gmul(p1,p2), gsqrt(absi(dpol),DEFAULTPREC));
1027: p1 = grndtoi(p1, &e);
1028: if (e>=0) p1 = addii(p1, shifti(gun, e));
1029: p1 = shifti(p1, 1);
1030: DATA[11]=(long)p1;
1031:
1032: if (DEBUGLEVEL>1)
1033: {
1034: fprintferr("DATA =\n");
1035: fprintferr("f = "); outerr((GEN)DATA[1]);
1036: fprintferr("p = "); outerr((GEN)DATA[2]);
1037: fprintferr("ff = "); outerr((GEN)DATA[3]);
1038: fprintferr("lcy = "); outerr((GEN)DATA[4]);
1039: fprintferr("cys = "); outerr((GEN)DATA[5]);
1040: fprintferr("bigfq = "); outerr((GEN)DATA[6]);
1041: fprintferr("roots = "); outerr((GEN)DATA[7]);
1042: fprintferr("2 * M = "); outerr((GEN)DATA[8]);
1043: fprintferr("p^e = "); outerr((GEN)DATA[9]);
1044: fprintferr("lifted roots = "); outerr((GEN)DATA[10]);
1045: fprintferr("2 * Hadamard bound = "); outerr((GEN)DATA[11]);
1046: }
1047: return DATA;
1048: }
1049:
1050: /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1051: /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1052: /* */
1053: /* AUTOMORPHISMS OF AN ABELIAN NUMBER FIELD */
1054: /* */
1055: /* V. Acciaro and J. Klueners (1996) */
1056: /* */
1057: /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1058: /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
1059:
1060: /* calcul du frobenius en p pour le corps abelien defini par le polynome pol,
1061: * par relevement de hensel du frobenius frobp de l'extension des corps
1062: * residuels (frobp est un polynome mod pol a coefficients dans F_p)
1063: */
1064: static GEN
1065: frobenius(GEN pol,GEN frobp,GEN p,GEN B,GEN d)
1066: {
1067: long av=avma,v0,deg,i,depas;
1068: GEN b0,b1,pold,polp,poldp,w0,w1,g0,g1,unmodp,polpp,v,pp,unmodpp,poldpp,bl0,bl1;
1069:
1070: v0=varn(pol); unmodp=gmodulsg(1,p); pold=deriv(pol,v0);
1071: b0=frobp; polp=gmul(unmodp,pol);
1072: poldp=gsubst(deriv(polp,v0),v0,frobp);
1073: w0=ginv(poldp);
1074: bl0=lift(b0); deg=lgef(bl0)-3;
1075: v=cgetg(deg+2,t_VEC);
1076: for (i=1; i<=deg+1; i++)
1077: v[i]=ldiv(centerlift(gmul(d,compo(bl0,deg+2-i))),d);
1078: g0=gtopoly(v,v0);
1079: if (DEBUGLEVEL>2)
1080: {
1081: fprintferr("val. initiales:\n");
1082: fprintferr("b0 = "); outerr(b0);
1083: fprintferr("w0 = "); outerr(w0);
1084: fprintferr("g0 = "); outerr(g0);
1085: }
1086: depas=1; pp=gsqr(p);
1087: for(;;)
1088: {
1089: if (gcmp(pp,B)>0) depas=0;
1090: unmodpp=gmodulsg(1,pp);
1091: polpp=gmul(unmodpp,pol); poldpp=gmul(unmodpp,pold);
1092: b0=gmodulcp(gmul(unmodpp,lift_intern(lift_intern(b0))),polpp);
1093: w0=gmodulcp(gmul(unmodpp,lift_intern(lift_intern(w0))),polpp);
1094: b1=gsub(b0,gmul(w0,gsubst(polpp,v0,b0)));
1095: w1=gmul(w0,gsub(gdeux,gmul(w0,gsubst(poldpp,v0,b1))));
1096: bl1=lift(b1); deg=lgef(bl1)-3;
1097: v=cgetg(deg+2,t_VEC);
1098: for (i=1; i<=deg+1; i++)
1099: v[i]=ldiv(centerlift(gmul(d,compo(bl1,deg+2-i))),d);
1100: g1=gtopoly(v,v0);
1101: if (DEBUGLEVEL>2)
1102: {
1103: fprintferr("pp = "); outerr(pp);
1104: fprintferr("b1 = "); outerr(b1);
1105: fprintferr("w1 = "); outerr(w1);
1106: fprintferr("g1 = "); outerr(g1);
1107: }
1108: if (gegal(g0,g1)) return gerepileupto(av,g1);
1109: pp=gsqr(pp); b0=b1; w0=w1; g0=g1;
1110: if (!depas) err(talker,"the number field is not an Abelian number field");
1111: }
1112: }
1113:
1114: static GEN
1115: compute_denom(GEN dpol)
1116: {
1117: long av=avma,lf,i,a;
1118: GEN d,f1,f2, f = decomp(dpol);
1119:
1120: f1=(GEN)f[1];
1121: f2=(GEN)f[2]; lf=lg(f1);
1122: for (d=gun,i=1; i<lf; i++)
1123: {
1124: a = itos((GEN)f2[i]) >> 1;
1125: d = mulii(d, gpuigs((GEN)f1[i],a));
1126: }
1127: return gerepileupto(av,d);
1128: }
1129:
1130: static GEN
1131: compute_bound_for_lift(GEN pol,GEN dpol,GEN d)
1132: {
1133: long av=avma,n,i;
1134: GEN p1,p2,p3;
1135:
1136: n=lgef(pol)-3;
1137: p1=gdiv(gmul(stoi(n),gpui(stoi(n-1),gdivgs(stoi(n-1),2),DEFAULTPREC)),
1138: gsqrt(dpol,DEFAULTPREC));
1139: p2=gzero;
1140: for (i=2; i<=n+2; i++) p2=gadd(p2,gsqr((GEN)pol[i]));
1141: p2=gpuigs(gsqrt(p2,DEFAULTPREC),n-1);
1142: p1=gmul(p1,p2); p2=gzero;
1143: for (i=2; i<=n+2; i++)
1144: {
1145: p3 = gabs((GEN)pol[i],DEFAULTPREC);
1146: if (gcmp(p3,p2)>0) p2 = p3;
1147: }
1148: p2=gmul(d,gadd(gun,p2));
1149: return gerepileupto(av, gmul2n(gsqr(gmul(p1,p2)),1));
1150:
1151: /* Borne heuristique de P. S. Wang, Math. Comp. 30, 1976, p. 332
1152: p2=gzero; for (i=2; i<=n+2; i++) p2=gadd(p2,gsqr((GEN)pol[i]));
1153: p1=gzero;
1154: for (i=2; i<=n+2; i++){ if (gcmp(gabs((GEN)pol[i],4),p1)>0) p1=gabs((GEN)pol[i],4); }
1155: if (gcmp(p2,p1)>0) p1=p2;
1156: p2=gmul(gdiv(mpfactr(n,4),gsqr(mpfactr(n/2,4))),d);
1157: B=gmul(p1,p2);
1158: tetpil=avma; return gerepile(av,tetpil,gcopy(B));
1159: */
1160: }
1161:
1162: static long
1163: isinlist(GEN T,long longT,GEN x)
1164: {
1165: long i;
1166: for (i=1; i<=longT; i++)
1167: if (gegal(x,(GEN)T[i])) return i;
1168: return 0;
1169: }
1170:
1171: /* renvoie 0 si frobp n'est pas dans la liste T; sinon le no de frobp dans T */
1172: static long
1173: isinlistmodp(GEN T,long longT,GEN frobp,GEN p)
1174: {
1175: long av=avma,i;
1176: GEN p1,p2,unmodp;
1177:
1178: p1=lift_intern(lift_intern(frobp)); unmodp=gmodulsg(1,p);
1179: for (i=1; i<=longT; i++)
1180: {
1181: p2=lift_intern(gmul(unmodp,(GEN)T[i]));
1182: if (gegal(p2,p1)) { avma=av; return i; }
1183: }
1184: avma=av; return 0;
1185: }
1186:
1187: /* renvoie le plus petit f tel que frobp^f est dans la liste T */
1188: static long
1189: minimalexponent(GEN T,long longT,GEN frobp,GEN p,long N)
1190: {
1191: long av=avma,i;
1192: GEN p1 = frobp;
1193: for (i=1; i<=N; i++)
1194: {
1195: if (isinlistmodp(T,longT,p1,p)) {avma=av; return i;}
1196: p1 = gpui(p1,p,DEFAULTPREC);
1197: }
1198: err(talker,"missing frobenius (field not abelian ?)");
1199: return 0; /* not reached */
1200: }
1201:
1202:
1203: /* Computation of all the automorphisms of the abelian number field
1204: defined by the monic irreducible polynomial pol with integral coefficients */
1205: GEN
1206: conjugates(GEN pol)
1207: {
1208: long av,tetpil,N,i,j,pp,bound_primes,nbprimes,longT,v0,flL,f,longTnew,*tab,nop,flnf;
1209: GEN T,S,p1,p2,p,dpol,modunp,polp,xbar,frobp,frob,d,B,nf;
1210: byteptr di;
1211:
1212: if (DEBUGLEVEL>2){ fprintferr("** Entree dans conjugates\n"); flusherr(); }
1213: flnf=0; if (typ(pol)!=t_POL){ nf=checknf(pol); flnf=1; pol=(GEN)nf[1]; }
1214: av=avma; N=lgef(pol)-3; v0=varn(pol);
1215: if (N==1) { S=cgetg(2,t_VEC); S[1]=(long)polx[v0]; return S; }
1216: if (N==2)
1217: {
1218: S=cgetg(3,t_VEC); S[1]=(long)polx[v0];
1219: S[2]=lsub(gneg(polx[v0]),(GEN)pol[3]);
1220: tetpil=avma; return gerepile(av,tetpil,gcopy(S));
1221: }
1222: dpol=absi(discsr(pol));
1223: if (DEBUGLEVEL>2)
1224: { fprintferr("discriminant du polynome: "); outerr(dpol); }
1225: d = flnf? (GEN)nf[4]: compute_denom(dpol);
1226: if (DEBUGLEVEL>2)
1227: { fprintferr("facteur carre du discriminant: "); outerr(d); }
1228: B=compute_bound_for_lift(pol,dpol,d);
1229: if (DEBUGLEVEL>2) { fprintferr("borne pour les lifts: "); outerr(B); }
1230: /* sous GRH il faut en fait 3.47*log(dpol) */
1231: p1=gfloor(glog(dpol,DEFAULTPREC));
1232: bound_primes=itos(p1);
1233: if (DEBUGLEVEL>2)
1234: { fprintferr("borne pour les premiers: %ld\n",bound_primes); flusherr(); }
1235: nbprimes=itos(gfloor(gmul(dbltor(1.25506),
1236: gdiv(p1,glog(p1,DEFAULTPREC)))));
1237: if (DEBUGLEVEL>2)
1238: { fprintferr("borne pour le nombre de premiers: %ld\n",nbprimes); flusherr(); }
1239: S=cgetg(nbprimes+1,t_VEC);
1240: di=diffptr; pp=*di; i=0;
1241: while (pp<=bound_primes)
1242: {
1243: if (smodis(dpol,pp)) { i++; S[i]=lstoi(pp); }
1244: pp = pp + (*(++di));
1245: }
1246: for (j=i+1; j<=nbprimes; j++) S[j]=zero;
1247: nbprimes=i; tab=new_chunk(nbprimes+1);
1248: for (i=1; i<=nbprimes; i++) tab[i]=0;
1249: if (DEBUGLEVEL>2)
1250: {
1251: fprintferr("nombre de premiers: %ld\n",nbprimes);
1252: fprintferr("table des premiers: "); outerr(S);
1253: }
1254: T=cgetg(N+1,t_VEC); T[1]=(long)polx[v0];
1255: for (i=2; i<=N; i++) T[i]=zero; longT=1;
1256: if (DEBUGLEVEL>2) { fprintferr("table initiale: "); outerr(T); }
1257: for(;;)
1258: {
1259: do
1260: {
1261: do
1262: {
1263: nop = 1+itos(shifti(mulss(mymyrand(),nbprimes),-(BITS_IN_RANDOM-1)));
1264: }
1265: while (tab[nop]);
1266: tab[nop]=1; p=(GEN)S[nop];
1267: if (DEBUGLEVEL>2) { fprintferr("\nnombre premier: "); outerr(p); }
1268: modunp=gmodulsg(1,p);
1269: polp=gmul(modunp,pol);
1270: xbar=gmodulcp(gmul(polx[v0],modunp),polp);
1271: frobp=gpui(xbar,p,4);
1272: if (DEBUGLEVEL>2) { fprintferr("frobenius mod p: "); outerr(frobp); }
1273: flL=isinlistmodp(T,longT,frobp,p);
1274: if (DEBUGLEVEL>2){ fprintferr("flL: %ld\n",flL); flusherr(); }
1275: }
1276: while (flL);
1277: f=minimalexponent(T,longT,frobp,p,N);
1278: if (DEBUGLEVEL>2){ fprintferr("exposant minimum: %ld\n",f); flusherr(); }
1279: frob=frobenius(pol,frobp,p,B,d);
1280: if (DEBUGLEVEL>2) { fprintferr("frobenius: "); outerr(frob); }
1281: /* Ce passage n'est vrai que si le corps est abelien !! */
1282: longTnew=longT;
1283: p2=gmodulcp(frob,pol);
1284: for (i=1; i<=longTnew; i++)
1285: for (j=1; j<f; j++)
1286: {
1287: p1=lift(gsubst((GEN)T[i],v0,gpuigs(p2,j)));
1288: if (DEBUGLEVEL>2)
1289: {
1290: fprintferr("test de la puissance (%ld,%ld): ",i,j); outerr(p1);
1291: }
1292: if (!isinlist(T,longTnew,p1))
1293: {
1294: longT++; T[longT]=(long)p1;
1295: if (longT==N)
1296: {
1297: if (DEBUGLEVEL>2)
1298: { fprintferr("** Sortie de conjugates\n"); flusherr(); }
1299: tetpil=avma; return gerepile(av,tetpil,gcopy(T));
1300: }
1301: }
1302: }
1303: if (DEBUGLEVEL>2) { fprintferr("nouvelle table: "); outerr(T); }
1304: }
1305: }
1306:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to subfield.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari / src / modules