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