Annotation of OpenXM_contrib/pari-2.2/src/modules/subfield.c, Revision 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>