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