Annotation of OpenXM_contrib/pari-2.2/src/modules/subfield.c, Revision 1.2
1.2 ! noro 1: /* $Id: subfield.c,v 1.55 2002/09/07 19:15:13 karim Exp $
1.1 noro 2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10: ANY WARRANTY WHATSOEVER.
11:
12: Check the License for details. You should have received a copy of it, along
13: with the package; see the file 'COPYING'. If not, write to the Free Software
14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
16: /*******************************************************************/
17: /* */
18: /* SUBFIELDS OF A NUMBER FIELD */
19: /* J. Klueners and M. Pohst, J. Symb. Comp. (1996), vol. 11 */
20: /* */
21: /*******************************************************************/
22: #include "pari.h"
23: extern GEN matrixpow(long n, long m, GEN y, GEN P,GEN l);
24: extern GEN Fp_factor_irred(GEN P,GEN l, GEN Q);
25: extern GEN FpX_rand(long d1, long v, GEN p);
26: extern GEN hensel_lift_fact(GEN pol, GEN Q, GEN T, GEN p, GEN pe, long e);
1.2 ! noro 27: extern GEN ZX_disc_all(GEN x, ulong bound);
! 28: extern GEN indexpartial(GEN P, GEN DP);
! 29: extern GEN initgaloisborne(GEN T, GEN dn, long prec, GEN *ptL, GEN *ptprep, GEN *ptdis);
! 30: extern long ZX_get_prec(GEN x);
1.1 noro 31:
32: static GEN print_block_system(long N,GEN Y,long d,GEN vbs,long maxl);
33:
34: /* Computation of potential block systems of given size d associated to a
35: * rational prime p: give a row vector of row vectors containing the
36: * potential block systems of imprimitivity; a potential block system is a
37: * vector of row vectors (enumeration of the roots).
38: */
39: #define BIL 32 /* for 64bit machines also */
40: static GEN
41: calc_block(long N,GEN Z,long d,GEN Y,GEN vbs,ulong maxl)
42: {
43: long r,lK,i,j,t,tp,T,lpn,u,nn,lnon,lY;
44: GEN K,n,non,pn,pnon,e,Yp,Zp,Zpp;
45:
46: if (DEBUGLEVEL>3)
47: {
48: long l = vbs? lg(vbs): 0;
49: fprintferr("avma = %ld, lg(Z) = %ld, lg(Y) = %ld, lg(vbs) = %ld\n",
50: avma,lg(Z),lg(Y),l);
51: if (DEBUGLEVEL > 5)
52: {
53: fprintferr("Z = %Z\n",Z);
54: fprintferr("Y = %Z\n",Y);
55: if (DEBUGLEVEL > 7) fprintferr("vbs = %Z\n",vbs);
56: }
57: }
58: r=lg(Z); lnon = min(BIL, r);
59: e = new_chunk(BIL);
60: n = new_chunk(r);
61: non = new_chunk(lnon);
62: pnon = new_chunk(lnon);
63: pn = new_chunk(lnon);
64: Zp = cgetg(lnon,t_VEC);
65: Zpp = cgetg(lnon,t_VEC);
66: for (i=1; i<r; i++) n[i] = lg(Z[i])-1;
67:
68: K=divisors(stoi(n[1])); lK=lg(K);
69: for (i=1; i<lK; i++)
70: {
71: long ngcd = n[1], k = itos((GEN)K[i]), dk = d*k;
72: lpn=0;
73: for (j=2; j<r; j++)
74: if (n[j]%k == 0)
75: {
76: if (++lpn >= BIL) err(talker,"overflow in calc_block");
77: pn[lpn] = n[j]; pnon[lpn] = j;
78: ngcd = cgcd(ngcd, n[j]);
79: }
80: if (dk % ngcd) continue;
81: T = 1<<lpn; if (!lpn) lpn = 1;
82: for (t=0; t<T; t++)
83: {
84: for (nn=n[1],tp=t, u=1; u<=lpn; u++,tp>>=1)
85: {
86: if (tp&1) { nn += pn[u]; e[u]=1; } else e[u]=0;
87: }
88: if (dk == nn)
89: {
1.2 ! noro 90: gpmem_t av=avma;
1.1 noro 91: int Z_equal_Zp = 1;
92:
93: for (j=1; j<lnon; j++) non[j]=0;
94: Zp[1]=Z[1];
95: for (u=2,j=1; j<=lpn; j++)
96: if (e[j])
97: {
98: Zp[u]=Z[pnon[j]]; non[pnon[j]]=1;
99: if (Zp[u] != Z[u]) Z_equal_Zp = 0;
100: u++;
101: }
102: setlg(Zp, u);
103: lY=lg(Y); Yp = cgetg(lY+1,t_VEC);
104: for (j=1; j<lY; j++) Yp[j]=Y[j];
105: Yp[lY]=(long)Zp;
106: if (r == u && Z_equal_Zp)
107: vbs = print_block_system(N,Yp,d,vbs,maxl);
108: else
109: {
110: for (u=1,j=2; j<r; j++)
111: if (!non[j]) Zpp[u++] = Z[j];
112: setlg(Zpp, u);
113: vbs = calc_block(N,Zpp,d,Yp,vbs,maxl);
114: }
1.2 ! noro 115: if (vbs && (ulong)lg(vbs) > maxl) return vbs;
1.1 noro 116: avma=av;
117: }
118: }
119: }
120: return vbs;
121: }
122:
123: static GEN
124: potential_block_systems(long N, long d, GEN n, ulong maxl)
125: {
1.2 ! noro 126: long r, i, j, k;
! 127: gpmem_t av=avma;
1.1 noro 128: GEN p1,vbs,Z;
129:
130: r=lg(n); Z=cgetg(r,t_VEC);
131: for (k=0,i=1; i<r; i++)
132: {
133: p1=cgetg(n[i]+1,t_VECSMALL); Z[i]=(long)p1;
134: for (j=1; j<=n[i]; j++) p1[j]= ++k;
135: }
136: vbs=calc_block(N,Z,d, cgetg(1,t_VEC), NULL, maxl);
137: avma=av; return vbs;
138: }
139:
140: /* product of permutations. Put the result in perm1. */
141: static void
1.2 ! noro 142: perm_mul_i(GEN perm1,GEN perm2)
1.1 noro 143: {
1.2 ! noro 144: long i, N = lg(perm1);
! 145: gpmem_t av = avma;
1.1 noro 146: GEN perm=new_chunk(N);
147: for (i=1; i<N; i++) perm[i]=perm1[perm2[i]];
148: for (i=1; i<N; i++) perm1[i]=perm[i];
149: avma=av;
150: }
151:
152: /* cy is a cycle; compute cy^l as a permutation */
153: static GEN
154: cycle_power_to_perm(GEN perm,GEN cy,long l)
155: {
156: long lp,i,j,b, N = lg(perm), lcy = lg(cy)-1;
157:
158: lp = l % lcy;
159: for (i=1; i<N; i++) perm[i] = i;
160: if (lp)
161: {
1.2 ! noro 162: gpmem_t av = avma;
1.1 noro 163: GEN p1 = new_chunk(N);
164: b = cy[1];
165: for (i=1; i<lcy; i++) b = (perm[b] = cy[i+1]);
166: perm[b] = cy[1];
167: for (i=1; i<N; i++) p1[i] = perm[i];
168:
1.2 ! noro 169: for (j=2; j<=lp; j++) perm_mul_i(perm,p1);
1.1 noro 170: avma = av;
171: }
172: return perm;
173: }
174:
175: /* image du block system D par la permutation perm */
176: static GEN
177: im_block_by_perm(GEN D,GEN perm)
178: {
179: long i,j,lb,lcy;
180: GEN Dn,cy,p1;
181:
182: lb=lg(D); Dn=cgetg(lb,t_VEC);
183: for (i=1; i<lb; i++)
184: {
185: cy=(GEN)D[i]; lcy=lg(cy);
186: Dn[i]=lgetg(lcy,t_VECSMALL); p1=(GEN)Dn[i];
187: for (j=1; j<lcy; j++) p1[j] = perm[cy[j]];
188: }
189: return Dn;
190: }
191:
192: /* cy is a cycle; recturn cy(a) */
193: static long
194: im_by_cy(long a,GEN cy)
195: {
196: long k, l = lg(cy);
197:
198: k=1; while (k<l && cy[k] != a) k++;
199: if (k == l) return a;
200: k++; if (k == l) k = 1;
201: return cy[k];
202: }
203:
204: /* vbs[0] = current cardinality+1, vbs[1] = max number of elts */
205: static GEN
206: append_vbs(GEN vbs, GEN D)
207: {
208: long l,maxl,i,j,n, lD = lg(D);
209: GEN Dn, last;
210:
211: n = 0; for (i=1; i<lD; i++) n += lg(D[i]);
212: Dn = (GEN)gpmalloc((lD + n) * sizeof(long));
213: last = Dn + lD; Dn[0] = D[0];
214: for (i=1; i<lD; i++)
215: {
216: GEN cy = (GEN)D[i], cn = last;
217: for (j=0; j<lg(cy); j++) cn[j] = cy[j];
218: Dn[i] = (long)cn; last = cn + j;
219: }
220:
221: if (!vbs)
222: {
223: maxl = 1024;
224: vbs = (GEN)gpmalloc((2 + maxl)*sizeof(GEN));
225: *vbs = maxl; vbs++; setlg(vbs, 1);
226: }
227: l = lg(vbs); maxl = vbs[-1];
228: if (l == maxl)
229: {
1.2 ! noro 230: vbs = (GEN)gprealloc((void*)(vbs-1), (2 + (maxl<<1))*sizeof(GEN));
1.1 noro 231: *vbs = maxl<<1; vbs++; setlg(vbs, 1);
232: }
233: if (DEBUGLEVEL>5) fprintferr("appending D = %Z\n",D);
234: vbs[l] = (long)Dn; setlg(vbs, l+1); return vbs;
235: }
236:
237: GEN
238: myconcat(GEN D, long a)
239: {
240: long i,l = lg(D);
241: GEN x = cgetg(l+1,t_VECSMALL);
242: for (i=1; i<l; i++) x[i]=D[i];
243: x[l] = a; return x;
244: }
245:
246: void
247: myconcat2(GEN D, GEN a)
248: {
249: long i,l = lg(D), m = lg(a);
250: GEN x = D + (l-1);
251: for (i=1; i<m; i++) x[i]=a[i];
252: setlg(D, l+m-1);
253: }
254:
255: static GEN
256: print_block_system(long N,GEN Y,long d,GEN vbs,long maxl)
257: {
258: long a,i,j,l,ll,*k,*n,lp,**e,u,v,*t,ns, r = lg(Y);
259: GEN D,De,Z,cyperm,perm,p1,empty;
260:
261: if (DEBUGLEVEL>5) fprintferr("Y = %Z\n",Y);
262: empty = cgetg(1,t_VEC);
263: n = new_chunk(N+1);
1.2 ! noro 264: D = cget1(N+1, t_VEC);
1.1 noro 265: t=new_chunk(r+1); k=new_chunk(r+1); Z=cgetg(r+1,t_VEC);
266: for (ns=0,i=1; i<r; i++)
267: {
268: GEN Yi = (GEN)Y[i], cy;
269: long ki = 0, si = lg(Yi)-1;
270:
271: for (j=1; j<=si; j++) { n[j]=lg(Yi[j])-1; ki += n[j]; }
272: ki /= d;
273: De=cgetg(ki+1,t_VEC);
274: for (j=1; j<=ki; j++) De[j]=(long)empty;
275: for (j=1; j<=si; j++)
276: {
277: cy = (GEN)Yi[j]; a = cy[1];
278: for (l=1,lp=0; l<=n[j]; l++)
279: {
280: lp++; if (lp>ki) lp = 1;
281: a = im_by_cy(a, cy);
282: De[lp] = (long)myconcat((GEN)De[lp], a);
283: }
284: }
285: if (si>1 && ki>1)
286: {
287: ns++; t[ns]=si-1; k[ns]=ki;
288: Z[ns]=lgetg(si,t_VEC); p1=(GEN)Z[ns];
289: for (j=2; j<=si; j++) p1[j-1]=Yi[j];
290: }
291: myconcat2(D,De);
292: }
293: if (DEBUGLEVEL>2) { fprintferr("\nns = %ld\n",ns); flusherr(); }
294: if (!ns) return append_vbs(vbs,D);
295:
296: setlg(Z, ns+1);
297: e=(long**)new_chunk(ns+1);
298: for (i=1; i<=ns; i++)
299: {
300: e[i]=new_chunk(t[i]+1);
301: for (j=1; j<=t[i]; j++) e[i][j]=0;
302: }
303: cyperm = cgetg(N+1,t_VEC);
304: perm = cgetg(N+1,t_VEC); i=ns;
305: do
306: {
1.2 ! noro 307: gpmem_t av = avma;
1.1 noro 308: if (DEBUGLEVEL>5)
309: {
310: for (l=1; l<=ns; l++)
311: {
312: for (ll=1; ll<=t[l]; ll++)
313: fprintferr("e[%ld][%ld] = %ld, ",l,ll,e[l][ll]);
314: fprintferr("\n");
315: }
316: fprintferr("\n"); flusherr();
317: }
318: for (u=1; u<=N; u++) perm[u]=u;
319: for (u=1; u<=ns; u++)
320: for (v=1; v<=t[u]; v++)
1.2 ! noro 321: perm_mul_i(perm, cycle_power_to_perm(cyperm,gmael(Z,u,v),e[u][v]));
1.1 noro 322: vbs = append_vbs(vbs, im_block_by_perm(D,perm));
323: if (lg(vbs) > maxl) return vbs;
324: avma = av;
325:
326: e[ns][t[ns]]++;
327: if (e[ns][t[ns]] >= k[ns])
328: {
329: j=t[ns]-1;
330: while (j>=1 && e[ns][j] == k[ns]-1) j--;
331: if (j>=1) { e[ns][j]++; for (l=j+1; l<=t[ns]; l++) e[ns][l]=0; }
332: else
333: {
334: i=ns-1;
335: while (i>=1)
336: {
337: j=t[i];
338: while (j>=1 && e[i][j] == k[i]-1) j--;
339: if (j<1) i--;
340: else
341: {
342: e[i][j]++;
343: for (l=j+1; l<=t[i]; l++) e[i][l]=0;
344: for (ll=i+1; ll<=ns; ll++)
345: for (l=1; l<=t[ll]; l++) e[ll][l]=0;
346: break;
347: }
348: }
349: }
350: }
351: }
352: while (i>0);
353: return vbs;
354: }
355:
356: static GEN
357: polsimplify(GEN x)
358: {
359: long i,lx = lgef(x);
360: for (i=2; i<lx; i++)
361: if (typ(x[i]) == t_POL) x[i] = (long)constant_term(x[i]);
362: return x;
363: }
364:
365: /* return 0 if |g[i]| > M[i] for some i; 1 otherwise */
366: static long
367: ok_coeffs(GEN g,GEN M)
368: {
369: long i, lg = lgef(g)-1; /* g is monic, and cst term is ok */
370: for (i=3; i<lg; i++)
371: if (absi_cmp((GEN)g[i], (GEN)M[i]) > 0) return 0;
372: return 1;
373: }
374:
375: /* Return a polynomial g defining a potential subfield, or
376: * 0: if p | d(g)
377: * 1: if coeffs of g are too large
378: * 2: same, detected by cheap d-1 test */
379: static GEN
380: cand_for_subfields(GEN A,GEN DATA,GEN *ptlistdelta)
381: {
382: long N,m,i,j,d,lf;
383: GEN M,T,pe,p,pol,fhk,g;
1.2 ! noro 384: GEN d_1_term,delta,listdelta,whichdelta,firstroot;
1.1 noro 385:
386: pol=(GEN)DATA[1];
387: p = (GEN)DATA[2];
388: pe= (GEN)DATA[3];
389: T = (GEN)DATA[4];
390: fhk=(GEN)DATA[5];
391: M = (GEN)DATA[8]; N=degpol(pol); m=lg(A)-1; d=N/m; /* m | M */
392: firstroot = (GEN)DATA[13];
393:
394: delta = cgetg(m+1,t_VEC);
395: lf = lg(firstroot);
396: listdelta = cgetg(lf, t_VEC);
397: whichdelta = cgetg(N+1, t_VECSMALL);
1.2 ! noro 398: d_1_term = gzero;
1.1 noro 399: for (i=1; i<=m; i++)
400: {
401: GEN Ai = (GEN)A[i], p1 = gun;
402: for (j=1; j<=d; j++)
1.2 ! noro 403: p1 = Fq_mul(p1, (GEN)fhk[Ai[j]], T,pe);
1.1 noro 404: delta[i] = (long)p1;
405: if (DEBUGLEVEL>2) fprintferr("delta[%ld] = %Z\n",i,p1);
406: /* g = prod (X - delta[i])
407: * if g o h = 0 (pol), we'll have h(Ai[j]) = delta[i] for all j */
408: /* fk[k] belongs to block number whichdelta[k] */
409: for (j=1; j<=d; j++) whichdelta[Ai[j]] = i;
410: if (typ(p1) == t_POL) p1 = constant_term(p1);
1.2 ! noro 411: d_1_term = addii(d_1_term, p1);
1.1 noro 412: }
1.2 ! noro 413: d_1_term = centermod(d_1_term, pe); /* Tr(g) */
! 414: if (absi_cmp(d_1_term, (GEN)M[3]) > 0) return gdeux; /* d-1 test */
1.1 noro 415: g = FqV_roots_to_pol(delta, T, pe, 0);
416: g = centermod(polsimplify(g), pe); /* assume g in Z[X] */
417: if (DEBUGLEVEL>2) fprintferr("pol. found = %Z\n",g);
418: if (!ok_coeffs(g,M)) return gun;
419: if (!FpX_is_squarefree(g, p)) return gzero;
420:
421: for (i=1; i<lf; i++) listdelta[i] = delta[whichdelta[firstroot[i]]];
422: *ptlistdelta = listdelta; return g;
423: }
424:
425: /* return U list of polynomials s.t U[i] = 1 mod fk[i] and 0 mod fk[j] for all
426: * other j */
427: static GEN
428: get_bezout(GEN pol, GEN fk, GEN p)
429: {
430: GEN A,B,d,u,v,bezout_coeff;
431: long i, l = lg(fk);
432:
433: pol = FpX_red(pol, p);
434: bezout_coeff = cgetg(l, t_VEC);
435: for (i=1; i<l; i++)
436: {
437: A = (GEN)fk[i];
438: B = FpX_div(pol, A, p);
439: d = FpX_extgcd(A,B,p, &u, &v);
440: if (degpol(d) > 0) err(talker, "relatively prime polynomials expected");
441: d = constant_term(d);
442: if (!gcmp1(d))
443: {
444: d = mpinvmod(d, p);
445: v = FpX_Fp_mul(v,d, p);
446: }
447: bezout_coeff[i] = (long)FpX_mul(B,v, p);
448: }
449: return bezout_coeff;
450: }
451:
452: /* assume x in Fq, return Tr_{Fq/Fp}(x) as a t_INT */
453: static GEN
454: trace(GEN x, GEN Trq, GEN p)
455: {
456: long i, l;
457: GEN s;
458: if (typ(x) == t_INT) return modii(mulii(x, (GEN)Trq[1]), p);
459: l = lgef(x)-1; if (l == 1) return gzero;
460: x++; s = mulii((GEN)x[1], (GEN)Trq[1]);
461: for (i=2; i<l; i++)
462: s = addii(s, mulii((GEN)x[i], (GEN)Trq[i]));
463: return modii(s, p);
464: }
465:
466: /* assume x in Fq[X], return Tr_{Fq[X]/Fp[X]}(x), varn(X) = 0 */
467: static GEN
468: poltrace(GEN x, GEN Trq, GEN p)
469: {
470: long i,l;
471: GEN y;
472: if (typ(x) == t_INT || varn(x) != 0) return trace(x, Trq, p);
473: l = lgef(x); y = cgetg(l,t_POL); y[1]=x[1];
474: for (i=2; i<l; i++) y[i] = (long)trace((GEN)x[i],Trq,p);
475: return y;
476: }
477:
478: /* Find h in Fp[X] such that h(a[i]) = listdelta[i] for all modular factors
479: * ff[i], where a[i] is a fixed root of ff[i] in Fq = Z[Y]/(p,T) [namely the
480: * first one in Fp_factor_irred output]. Let f = ff[i], A the given root, then
481: * h mod f is Tr_Fq/Fp ( h(A) f(X)/(X-A)f'(A) ), most of the expression being
482: * precomputed. The complete h is recovered via chinese remaindering */
483: static GEN
484: chinese_retrieve_pol(GEN DATA, GEN listdelta)
485: {
486: GEN interp,Trk,bezoutC,T,p, S,p1;
487: long i,l;
488: p = (GEN)DATA[2];
489: T = (GEN)DATA[4];
490: interp = (GEN)DATA[10];
491: Trk = (GEN)DATA[11];
492: bezoutC = (GEN)DATA[12];
493:
494: S = NULL; l = lg(interp);
495: for (i=1; i<l; i++)
496: { /* h(firstroot[i]) = listdelta[i] */
497: p1 = FpXQX_FpXQ_mul((GEN)interp[i], (GEN)listdelta[i], T,p);
498: p1 = poltrace(p1, (GEN)Trk[i], p);
499: p1 = gmul(p1, (GEN)bezoutC[i]);
500: S = S? gadd(S,p1): p1;
501: }
502: return FpX_res(FpX_red(S, p), FpX_red((GEN)DATA[1],p), p);
503: }
504:
1.2 ! noro 505: /* return P(X + c) using destructive Horner */
! 506: static GEN
! 507: TR_pol(GEN P, GEN c)
! 508: {
! 509: GEN Q = dummycopy(P), *R;
! 510: long i,k,n;
! 511:
! 512: if (!signe(P)) return Q;
! 513: R = (GEN*)(Q+2); n = degpol(P);
! 514: if (gcmp1(c))
! 515: {
! 516: for (i=1; i<=n; i++)
! 517: for (k=n-i; k<n; k++)
! 518: R[k] = gadd(R[k], R[k+1]);
! 519: }
! 520: else if (gcmp_1(c))
! 521: {
! 522: for (i=1; i<=n; i++)
! 523: for (k=n-i; k<n; k++)
! 524: R[k] = gsub(R[k], R[k+1]);
! 525: }
! 526: else
! 527: {
! 528: for (i=1; i<=n; i++)
! 529: for (k=n-i; k<n; k++)
! 530: R[k] = gadd(R[k], gmul(c, R[k+1]));
! 531: }
! 532: return Q;
! 533: }
! 534:
1.1 noro 535: /* g in Z[X] potentially defines a subfield of Q[X]/f. It is a subfield iff A
536: * (cf cand_for_subfields) was a block system; then there
537: * exists h in Q[X] such that f | g o h. listdelta determines h s.t f | g o h
538: * in Fp[X] (cf chinese_retrieve_pol). We try to lift it. */
539: static GEN
540: embedding_of_potential_subfields(GEN g,GEN DATA,GEN listdelta)
541: {
542: GEN TR,w0_Q,w0,w1_Q,w1,wpow,h0,gp,T,q2,q,p,ind,maxp,a;
543: long rt;
1.2 ! noro 544: gpmem_t av;
1.1 noro 545:
546: T = (GEN)DATA[1]; rt = brent_kung_optpow(degpol(T), 2);
547: p = (GEN)DATA[2];
548: maxp = (GEN)DATA[7];
549: ind = (GEN)DATA[9];
550: gp = derivpol(g); av = avma;
551: w0 = chinese_retrieve_pol(DATA,listdelta);
552: w0_Q = centermod(gmul(w0,ind), p);
553: h0 = FpXQ_inv(FpX_FpXQ_compo(gp,w0, T,p), T,p); /* = 1/g'(w0) mod (T,p) */
554: wpow = NULL; q = sqri(p);
555: for(;;)
556: {/* Given g,w0,h0 in Z[x], s.t. h0.g'(w0) = 1 and g(w0) = 0 mod (T,p), find
557: * [w1,h1] satisfying the same conditions mod p^2, [w1,h1] = [w0,h0] (mod p)
558: * (cf. Dixon: J. Austral. Math. Soc., Series A, vol.49, 1990, p.445) */
559: if (DEBUGLEVEL>1) fprintferr("lifting embedding mod p = %Z\n",q);
560:
561: /* w1 := w0 - h0 g(w0) mod (T,q) */
562: if (wpow)
563: a = FpX_FpXQV_compo(g,wpow, T,q);
564: else
565: a = FpX_FpXQ_compo(g,w0, T,q); /* first time */
566: /* now, a = 0 (p) */
567: a = gmul(gneg(h0), gdivexact(a, p));
568: w1 = gadd(w0, gmul(p, FpX_res(a, T,p)));
569:
570: w1_Q = centermod(gmul(w1, resii(ind,q)), q);
571: if (gegal(w1_Q, w0_Q) || cmpii(q,maxp) > 0)
572: {
1.2 ! noro 573: GEN G = gcmp1(ind)? g: rescale_pol(g,ind);
1.1 noro 574: if (gcmp0(poleval(G, gmodulcp(w1_Q,T)))) break;
575: }
576: if (cmpii(q, maxp) > 0)
577: {
578: if (DEBUGLEVEL) fprintferr("coeff too big for embedding\n");
579: return NULL;
580: }
581: {
582: GEN *gptr[5]; gptr[0]=&w1; gptr[1]=&h0; gptr[2]=&w1_Q;
583: gptr[3]=&q; gptr[4]=&p;
584: gerepilemany(av,gptr,5);
585: }
586:
587: q2 = sqri(q);
588: wpow = FpXQ_powers(w1, rt, T, q2);
589: /* h0 := h0 * (2 - h0 g'(w1)) mod (T,q)
590: * = h0 + h0 * (1 - h0 g'(w1)) */
591: a = gmul(gneg(h0), FpX_FpXQV_compo(gp, FpXV_red(wpow,q),T,q));
592: a = ZX_s_add(FpX_res(a, T,q), 1); /* 1 - h0 g'(w1) = 0 (p) */
593: a = gmul(h0, gdivexact(a, p));
594: h0 = gadd(h0, gmul(p, FpX_res(a, T,p)));
595: w0 = w1; w0_Q = w1_Q; p = q; q = q2;
596: }
597: TR = (GEN)DATA[14];
1.2 ! noro 598: if (!gcmp0(TR)) w1_Q = TR_pol(w1_Q, TR);
1.1 noro 599: return gdiv(w1_Q,ind);
600: }
601:
602: static GEN
603: choose_prime(GEN pol,GEN dpol,long d,GEN *ptff,GEN *ptlistpotbl, long *ptlcm)
604: {
1.2 ! noro 605: gpmem_t av;
! 606: byteptr di=diffptr;
! 607: long j,k,oldllist,llist,r,lcm,oldlcm,pp,minp, N = degpol(pol), m = N/d;
1.1 noro 608: GEN p,listpotbl,oldlistpotbl,ff,oldff,n,oldn;
1.2 ! noro 609:
! 610: minp = N*(m-1)/2;
! 611: if (DEBUGLEVEL) (void)timer2();
! 612: di++; p = stoi(2);
! 613: while (p[2]<=minp)
! 614: NEXT_PRIME_VIADIFF(p[2], di);
1.1 noro 615: oldllist = 100000;
616: oldlcm = 0;
617: oldlistpotbl = oldff = oldn = NULL; pp = 0; /* gcc -Wall */
618: av = avma;
1.2 ! noro 619: for(k=1; k<11 || !oldn; k++,avma = av)
1.1 noro 620: {
1.2 ! noro 621: while (!smodis(dpol,p[2]))
! 622: NEXT_PRIME_VIADIFF(p[2], di);
1.1 noro 623: if (k > 50) err(talker,"sorry, too many block systems in nfsubfields");
1.2 ! noro 624: ff=(GEN)factmod0(pol,p)[1]; r=lg(ff)-1;
! 625: if (r == 1 || r == N) goto repeat;
1.1 noro 626:
627: n = cgetg(r+1, t_VECSMALL);
628: lcm = n[1] = degpol(ff[1]);
629: for (j=2; j<=r; j++) { n[j] = degpol(ff[j]); lcm = clcm(lcm,n[j]); }
1.2 ! noro 630: if (lcm < oldlcm) goto repeat; /* false when oldlcm = 0 */
! 631: if (r >= BIL) { err(warner,"subfields: overflow in calc_block"); goto repeat; }
1.1 noro 632: if (DEBUGLEVEL) fprintferr("p = %ld,\tlcm = %ld,\torbits: %Z\n",p[2],lcm,n);
1.2 ! noro 633: if (oldn && gegal(n,oldn)) goto repeat;
1.1 noro 634:
635: listpotbl = potential_block_systems(N,d,n, oldllist);
636: if (!listpotbl) { oldlistpotbl = NULL; pp = p[2]; break; }
637: llist = lg(listpotbl)-1;
638: if (llist >= oldllist)
639: {
640: if (DEBUGLEVEL) msgtimer("#pbs >= %ld [aborted]",oldllist);
641: for (j=1; j<llist; j++) free((void*)listpotbl[j]);
1.2 ! noro 642: free((void*)(listpotbl-1)); goto repeat;
1.1 noro 643: }
644: oldllist = llist; oldlistpotbl = listpotbl;
645: pp = p[2]; oldff = ff; oldlcm = lcm; oldn = n;
646: if (DEBUGLEVEL) msgtimer("#pbs = %ld",llist);
647: av = avma;
1.2 ! noro 648: repeat:
! 649: NEXT_PRIME_VIADIFF(p[2], di);
1.1 noro 650: }
1.2 ! noro 651: if (DEBUGLEVEL) fprintferr("Chosen prime: p = %ld\n",pp);
1.1 noro 652: *ptlistpotbl=oldlistpotbl; *ptff=oldff; *ptlcm=oldlcm; return stoi(pp);
653: }
654:
655: static GEN
656: bound_for_coeff(long m,GEN rr, GEN *maxroot)
657: {
658: long i,r1, lrr=lg(rr);
659: GEN p1,b1,b2,B,M, C = matpascal(m-1);
660:
1.2 ! noro 661: for (r1=1; r1 < lrr; r1++)
! 662: if (typ(rr[r1]) != t_REAL) break;
! 663: r1--;
1.1 noro 664:
665: rr = gabs(rr,0); *maxroot = vecmax(rr);
666: for (i=1; i<lrr; i++)
667: if (gcmp((GEN)rr[i], gun) < 0) rr[i] = un;
668: for (b1=gun,i=1; i<=r1; i++) b1 = gmul(b1, (GEN)rr[i]);
669: for (b2=gun ; i<lrr; i++) b2 = gmul(b2, (GEN)rr[i]);
670: B = gmul(b1, gsqr(b2)); /* Mahler measure */
671: M = cgetg(m+2, t_VEC); M[1]=M[2]=zero; /* unused */
672: for (i=1; i<m; i++)
673: {
674: p1 = gadd(gmul(gcoeff(C, m, i+1), B),/* binom(m-1, i) */
675: gcoeff(C, m, i)); /* binom(m-1, i-1) */
676: M[i+2] = (long)ceil_safe(p1);
677: }
678: return M;
679: }
680:
681: static GEN
682: init_traces(GEN ff, GEN T, GEN p)
683: {
684: long N = degpol(T),i,j,k, r = lg(ff);
685: GEN Frob = matrixpow(N,N, FpXQ_pow(polx[varn(T)],p, T,p), T,p);
686: GEN y,p1,p2,Trk,pow,pow1;
687:
688: k = degpol(ff[r-1]); /* largest degree in modular factorization */
689: pow = cgetg(k+1, t_VEC);
690: pow[1] = (long)zero; /* dummy */
691: pow[2] = (long)Frob;
692: pow1= cgetg(k+1, t_VEC); /* 1st line */
693: for (i=3; i<=k; i++)
694: pow[i] = (long)FpM_mul((GEN)pow[i-1], Frob, p);
695: pow1[1] = (long)zero; /* dummy */
696: for (i=2; i<=k; i++)
697: {
698: p1 = cgetg(N+1, t_VEC);
699: pow1[i] = (long)p1; p2 = (GEN)pow[i];
700: for (j=1; j<=N; j++) p1[j] = coeff(p2,1,j);
701: }
702: p1 = cgetg(N+1,t_VEC); p1[1] = un;
703: for (i=2; i<=N; i++) p1[i] = zero;
704: /* Trk[i] = line 1 of x -> x + x^p + ... + x^{p^(i-1)} */
705: Trk = pow; /* re-use (destroy) pow */
706: Trk[1] = (long)p1;
707: for (i=2; i<=k; i++)
708: Trk[i] = ladd((GEN)Trk[i-1], (GEN)pow1[i]);
709: y = cgetg(r, t_VEC);
710: for (i=1; i<r; i++) y[i] = Trk[degpol(ff[i])];
711: return y;
712: }
713:
714: /* return C in Z[X]/(p,T), C[ D[1] ] = 1, C[ D[i] ] = 0 otherwise. H is the
715: * list of degree 1 polynomials X - D[i] (come directly from factorization) */
716: static GEN
717: interpol(GEN H, GEN T, GEN p)
718: {
719: long i, m = lg(H);
720: GEN X = polx[0],d,p1,p2,a;
721:
722: p1=polun[0]; p2=gun; a = gneg(constant_term(H[1])); /* = D[1] */
723: for (i=2; i<m; i++)
724: {
725: d = constant_term(H[i]); /* -D[i] */
726: p1 = FpXQX_mul(p1,gadd(X,d), T,p);
1.2 ! noro 727: p2 = Fq_mul(p2, gadd(a,d), T,p);
1.1 noro 728: }
729: p2 = FpXQ_inv(p2, T,p);
730: return FpXQX_FpXQ_mul(p1,p2, T,p);
731: }
732:
733: static GEN
734: roots_from_deg1(GEN x)
735: {
736: long i,l = lg(x);
737: GEN r = cgetg(l,t_VEC);
738: for (i=1; i<l; i++) r[i] = lneg(constant_term(x[i]));
739: return r;
740: }
741:
742: struct poldata
743: {
744: GEN pol;
745: GEN dis; /* |disc(pol)| */
746: GEN roo; /* roots(pol) */
747: GEN den; /* multiple of index(pol) */
748: };
749:
750: /* ff = factmod(nf[1], p), d = requested degree for subfield. Return DATA,
751: * valid for given nf, p and d
752: * 1: polynomial nf[1],
753: * 2: prime p,
754: * 3: exponent e (for Hensel lifts) such that p^e > max(M),
755: * 4: polynomial defining the field F_(p^nn),
756: * 5: Hensel lift to precision p^e of DATA[6]
757: * 6: roots of nf[1] in F_(p^nn),
758: * 7: Hadamard bound for coefficients of h(x) such that g o h = 0 mod nf[1].
759: * 8: bound M for polynomials defining subfields x DATA[9]
760: * 9: multiple of index of nf
761: *10: *[i] = interpolation polynomial for ff[i] [= 1 on the first root
762: firstroot[i], 0 on the others]
763: *11: Trk used to compute traces (cf poltrace)
764: *12: Bezout coefficients associated to the ff[i]
765: *13: *[i] = index of first root of ff[i] (in DATA[6])
766: *14: number of polynomial changes (translations) */
767: static GEN
768: compute_data(GEN DATA, struct poldata PD, long d, GEN ff, GEN T,GEN p)
769: {
1.2 ! noro 770: long i,j,l,e,N, lff = lg(ff);
! 771: GEN ffL,den,roo,pe,p1,p2,fk,fhk,MM,maxroot,pol,interp,bezoutC;
1.1 noro 772:
773: if (DEBUGLEVEL>1) { fprintferr("Entering compute_data()\n\n"); flusherr(); }
774: pol = PD.pol; N = degpol(pol);
775: roo = PD.roo;
776: den = PD.den;
777: if (DATA) /* update (translate) an existing DATA */
778: {
779: GEN Xm1 = gsub(polx[varn(pol)], gun);
780: GEN TR = addis((GEN)DATA[14],1);
1.2 ! noro 781: GEN mTR = negi(TR), mun = negi(gun);
1.1 noro 782: DATA[14] = (long)TR;
1.2 ! noro 783: pol = TR_pol(pol, mTR);
1.1 noro 784: p1 = dummycopy(roo); l = lg(p1);
785: for (i=1; i<l; i++) p1[i] = ladd(TR, (GEN)p1[i]);
786: roo = p1;
787:
788: fk = (GEN)DATA[6]; l=lg(fk);
789: for (i=1; i<l; i++) fk[i] = lsub(Xm1, (GEN)fk[i]);
790:
791: interp = (GEN)DATA[10];
792: bezoutC = (GEN)DATA[12]; l = lg(interp);
793: for (i=1; i<l; i++)
794: {
795: if (degpol(interp[i]) > 0) /* do not turn polun[0] into gun */
796: {
1.2 ! noro 797: p1 = TR_pol((GEN)interp[i], mun);
1.1 noro 798: interp[i] = (long)FpXQX_red(p1, NULL,p);
799: }
800: if (degpol(bezoutC[i]) > 0)
801: {
1.2 ! noro 802: p1 = TR_pol((GEN)bezoutC[i], mun);
1.1 noro 803: bezoutC[i] = (long)FpXQX_red(p1, NULL,p);
804: }
805: }
1.2 ! noro 806: p1 = cgetg(lff, t_VEC);
! 807: for (i=1; i<lff; i++)
! 808: p1[i] = (long)FpX_red(TR_pol((GEN)ff[i], mTR), p);
! 809: ff = p1;
1.1 noro 810: }
811: else
812: {
1.2 ! noro 813: GEN firstroot = cgetg(lff, t_VECSMALL);
1.1 noro 814: DATA = cgetg(15,t_VEC);
815: DATA[2] = (long)p;
816: DATA[4] = (long)T;
1.2 ! noro 817: interp = cgetg(lff, t_VEC);
1.1 noro 818: fk = cgetg(N+1,t_VEC);
1.2 ! noro 819: for (l=1,j=1; j<lff; j++)
1.1 noro 820: { /* compute roots and fix ordering (Frobenius cycles) */
1.2 ! noro 821: p1 = Fp_factor_irred((GEN)ff[j], p,T);
1.1 noro 822: interp[j] = (long)interpol(p1,T,p);
823: firstroot[j] = l;
824: for (i=1; i<lg(p1); i++,l++) fk[l] = p1[i];
825: }
826: DATA[9] = (long)PD.den;
827: DATA[10]= (long)interp;
828: DATA[11]= (long)init_traces(ff, T,p);
829: DATA[12]= (long)get_bezout(pol,ff,p);
830: DATA[13]= (long)firstroot;
831: DATA[14]= zero;
832: }
833: DATA[1] = (long)pol;
834: MM = bound_for_coeff(d, roo, &maxroot);
835: MM = gmul2n(MM,1);
836: DATA[8] = (long)MM;
837: e = logint(shifti(vecmax(MM),20), p, &pe); /* overlift 2^20 [for d-1 test] */
838: DATA[3] = (long)pe;
839: DATA[6] = (long)roots_from_deg1(fk);
1.2 ! noro 840:
! 841: #if 0
1.1 noro 842: fhk = hensel_lift_fact(pol,fk,(GEN)DATA[4],p,pe,e);
1.2 ! noro 843: #else
! 844: /* first lift in Zp to precision p^e */
! 845: ffL = hensel_lift_fact(pol,ff, NULL,p,pe,e);
! 846: fhk = NULL;
! 847: for (l=i=1; i<lff; i++)
! 848: { /* lift factorization of ff[i] in Qp[X] / T */
! 849: GEN L = (GEN)ffL[i];
! 850: long di = degpol(L);
! 851: p2 = cgetg(di+1, t_VEC);
! 852: for (j=1; j<=di; j++) p2[j] = fk[l++];
! 853: p1 = hensel_lift_fact(L,p2,T,p,pe,e);
! 854: fhk = fhk? concatsp(fhk, p1): p1;
! 855: }
! 856: #endif
1.1 noro 857: DATA[5] = (long)roots_from_deg1(fhk);
858:
859: p1 = gmul(stoi(N), gsqrt(gpuigs(stoi(N-1),N-1),DEFAULTPREC));
860: p2 = gpuigs(maxroot, N/d + N*(N-1)/2);
861: p1 = gdiv(gmul(p1,p2), gsqrt(PD.dis,DEFAULTPREC));
862: p1 = shifti(ceil_safe(p1), 1);
863: DATA[7] = lmulii(p1, PD.den);
864:
865: if (DEBUGLEVEL>1) {
866: fprintferr("f = %Z\n",DATA[1]);
867: fprintferr("p = %Z, lift to p^%ld\n",DATA[2], e);
868: fprintferr("Fq defined by %Z\n",DATA[4]);
869: fprintferr("2 * Hadamard bound * ind = %Z\n",DATA[7]);
870: fprintferr("2 * M = %Z\n",DATA[8]);
871: }
872: return DATA;
873: }
874:
875: /* g = polynomial, h = embedding. Return [[g,h]] */
876: static GEN
877: _subfield(GEN g, GEN h)
878: {
879: GEN x = cgetg(3,t_VEC);
880: x[1] = (long)g;
881: x[2] = (long)h; return _vec(x);
882: }
883:
884: /* subfields of degree d */
885: static GEN
886: subfields_of_given_degree(struct poldata PD,long d)
887: {
1.2 ! noro 888: gpmem_t av, av2;
! 889: long llist,i,nn,v;
1.1 noro 890: GEN listpotbl,ff,A,CSF,ESF,LSB,p,T,DATA,listdelta;
891: GEN pol = PD.pol, dpol = PD.dis;
892:
893: if (DEBUGLEVEL) fprintferr("\n*** Look for subfields of degree %ld\n\n", d);
894: av = avma;
895: p = choose_prime(pol,dpol,degpol(pol)/d,&ff,&listpotbl,&nn);
896: if (!listpotbl) { avma=av; return cgetg(1,t_VEC); }
1.2 ! noro 897: v = fetch_var(); name_var(v,"y");
! 898: T = lift_intern(ffinit(p,nn, v));
1.1 noro 899: DATA = NULL; LSB = cgetg(1,t_VEC);
900: i = 1; llist = lg(listpotbl);
901: CHANGE:
902: DATA = compute_data(DATA,PD,d, ff,T,p);
903: for (; i<llist; i++)
904: {
905: av2 = avma; A = (GEN)listpotbl[i];
906: if (DEBUGLEVEL > 1) fprintferr("\n* Potential block # %ld: %Z\n",i,A);
907: CSF = cand_for_subfields(A,DATA,&listdelta);
908: if (typ(CSF)==t_INT)
909: {
910: avma = av2;
911: if (DEBUGLEVEL > 1) switch(itos(CSF))
912: {
913: case 0: fprintferr("changing f(x): p divides disc(g(x))\n"); break;
914: case 1: fprintferr("coeff too big for pol g(x)\n"); break;
915: case 2: fprintferr("d-1 test failed\n"); break;
916: }
917: if (CSF==gzero) goto CHANGE;
918: }
919: else
920: {
921: if (DEBUGLEVEL) fprintferr("candidate = %Z\n",CSF);
922: ESF = embedding_of_potential_subfields(CSF,DATA,listdelta);
923: if (!ESF) { avma = av2; continue; }
924:
925: if (DEBUGLEVEL) fprintferr("embedding = %Z\n",ESF);
926: LSB = gerepileupto(av2, concat(LSB, _subfield(CSF,ESF)));
927: }
928: }
1.2 ! noro 929: (void)delete_var();
1.1 noro 930: for (i=1; i<llist; i++) free((void*)listpotbl[i]);
931: free((void*)(listpotbl-1));
932: if (DEBUGLEVEL) fprintferr("\nSubfields of degree %ld: %Z\n",d, LSB);
933: return gerepilecopy(av, LSB);
934: }
935:
936: static GEN
937: fix_var(GEN x, long v)
938: {
939: long i, l = lg(x);
940: x = gcopy(x); if (!v) return x;
941: for (i=1; i<l; i++) { GEN t = (GEN)x[i]; setvarn(t[1],v); setvarn(t[2],v); }
942: return x;
943: }
944:
945: void
946: subfields_poldata(GEN T, struct poldata *PD)
947: {
1.2 ! noro 948: GEN nf,L,dis;
1.1 noro 949:
950: T = get_nfpol(T, &nf);
951: PD->pol = dummycopy(T); /* may change variable number later */
952: if (nf)
953: {
1.2 ! noro 954: PD->den = Q_denom((GEN)nf[7]);
1.1 noro 955: PD->roo = (GEN)nf[6];
956: PD->dis = mulii(absi((GEN)nf[3]), sqri((GEN)nf[4]));
957: }
1.2 ! noro 958: else
! 959: {
! 960: PD->den = initgaloisborne(T,NULL,ZX_get_prec(T), &L,NULL,&dis);
! 961: PD->roo = L;
! 962: PD->dis = absi(dis);
! 963: }
1.1 noro 964: }
965:
966: GEN
967: subfields(GEN nf,GEN d)
968: {
1.2 ! noro 969: gpmem_t av = avma;
! 970: long di, N, v0;
! 971: GEN LSB, pol, G;
1.1 noro 972: struct poldata PD;
973:
974: pol = get_nfpol(nf, &nf); /* in order to treat trivial cases */
975: v0=varn(pol); N=degpol(pol); di=itos(d);
976: if (di==N) return gerepilecopy(av, _subfield(pol, polx[v0]));
977: if (di==1) return gerepilecopy(av, _subfield(polx[v0], pol));
978: if (di < 1 || di > N || N % di) return cgetg(1,t_VEC);
979:
1.2 ! noro 980: /* much easier if nf is Galois (WSS) */
! 981: G = galoisconj4(nf? nf: pol, NULL, 1, 0);
! 982: if (typ(G) != t_INT)
! 983: { /* Bingo */
! 984: GEN L = galoissubgroups(G), F;
! 985: long k,i, l = lg(L), o = N/di;
! 986: F = cgetg(l, t_VEC);
! 987: k = 1;
! 988: for (i=1; i<l; i++)
! 989: {
! 990: GEN H = (GEN)L[i];
! 991: if (group_order(H) == o)
! 992: F[k++] = (long) lift_intern(galoisfixedfield(G, (GEN)H[1], 0, v0));
! 993: }
! 994: setlg(F, k);
! 995: return gerepilecopy(av, F);
! 996: }
! 997:
! 998: subfields_poldata(nf? nf:pol, &PD);
1.1 noro 999: pol = PD.pol;
1000: setvarn(pol, 0);
1001: LSB = subfields_of_given_degree(PD, di);
1002: return gerepileupto(av, fix_var(LSB,v0));
1003: }
1004:
1005: static GEN
1006: subfieldsall(GEN nf)
1007: {
1.2 ! noro 1008: gpmem_t av = avma;
! 1009: long N, ld, i, v0;
! 1010: GEN G, pol, dg, LSB, NLSB;
1.1 noro 1011: struct poldata PD;
1012:
1.2 ! noro 1013: /* much easier if nf is Galois (WSS) */
! 1014: G = galoisconj4(nf, NULL, 1, 0);
! 1015: if (typ(G) != t_INT)
! 1016: {
! 1017: pol = get_nfpol(nf, &nf);
! 1018: return gerepilecopy(av, lift_intern( galoissubfields(G, 0, varn(pol)) ));
! 1019: }
! 1020:
1.1 noro 1021: subfields_poldata(nf, &PD);
1022: pol = PD.pol;
1023:
1024: v0 = varn(pol); N = degpol(pol);
1025: dg = divisors(stoi(N)); ld = lg(dg)-1;
1026: if (DEBUGLEVEL) fprintferr("\n***** Entering subfields\n\npol = %Z\n",pol);
1027:
1028: LSB = _subfield(pol,polx[0]);
1029: if (ld > 2)
1030: {
1031: setvarn(pol, 0);
1032: for (i=2; i<ld; i++)
1033: {
1.2 ! noro 1034: gpmem_t av1 = avma;
1.1 noro 1035: NLSB = subfields_of_given_degree(PD, N / itos((GEN)dg[i]));
1036: if (lg(NLSB) > 1) LSB = concatsp(LSB,NLSB); else avma = av1;
1037: }
1038: }
1039: LSB = concatsp(LSB, _subfield(polx[0],pol));
1040: if (DEBUGLEVEL) fprintferr("\n***** Leaving subfields\n\n");
1041: return gerepileupto(av, fix_var(LSB,v0));
1042: }
1043:
1044: GEN
1045: subfields0(GEN nf,GEN d)
1046: {
1047: return d? subfields(nf,d): subfieldsall(nf);
1048: }
1049:
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-2.2 / src / modules