Annotation of OpenXM_contrib/pari-2.2/src/basemath/galconj.c, Revision 1.1
1.1 ! noro 1: /* $Id: galconj.c,v 1.67 2001/10/01 16:41:42 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: /** GALOIS CONJUGATES **/
! 19: /** **/
! 20: /*************************************************************************/
! 21: #include "pari.h"
! 22: #define mylogint(x,y) logint((x),(y),NULL)
! 23: GEN
! 24: galoisconj(GEN nf)
! 25: {
! 26: GEN x, y, z;
! 27: long i, lz, v, av = avma;
! 28: nf = checknf(nf);
! 29: x = (GEN) nf[1];
! 30: v = varn(x);
! 31: if (v == 0)
! 32: nf = gsubst(nf, 0, polx[MAXVARN]);
! 33: else
! 34: {
! 35: x = dummycopy(x);
! 36: setvarn(x, 0);
! 37: }
! 38: z = nfroots(nf, x);
! 39: lz = lg(z);
! 40: y = cgetg(lz, t_COL);
! 41: for (i = 1; i < lz; i++)
! 42: {
! 43: GEN p1 = lift((GEN) z[i]);
! 44: setvarn(p1, v);
! 45: y[i] = (long) p1;
! 46: }
! 47: return gerepileupto(av, y);
! 48: }
! 49:
! 50: /* nbmax: maximum number of possible conjugates */
! 51: GEN
! 52: galoisconj2pol(GEN x, long nbmax, long prec)
! 53: {
! 54: long av = avma, i, n, v, nbauto;
! 55: GEN y, w, polr, p1, p2;
! 56: n = degpol(x);
! 57: if (n <= 0)
! 58: return cgetg(1, t_VEC);
! 59: if (gisirreducible(x) == gzero)
! 60: err(redpoler, "galoisconj2pol");
! 61: polr = roots(x, prec);
! 62: p1 = (GEN) polr[1];
! 63: nbauto = 1;
! 64: prec = (long) (bit_accuracy(prec) * L2SL10 * 0.75);
! 65: w = cgetg(n + 2, t_VEC);
! 66: w[1] = un;
! 67: for (i = 2; i <= n; i++)
! 68: w[i] = lmul(p1, (GEN) w[i - 1]);
! 69: v = varn(x);
! 70: y = cgetg(nbmax + 1, t_COL);
! 71: y[1] = (long) polx[v];
! 72: for (i = 2; i <= n && nbauto < nbmax; i++)
! 73: {
! 74: w[n + 1] = polr[i];
! 75: p1 = lindep2(w, prec);
! 76: if (signe(p1[n + 1]))
! 77: {
! 78: setlg(p1, n + 1);
! 79: p2 = gdiv(gtopolyrev(p1, v), negi((GEN) p1[n + 1]));
! 80: if (gdivise(poleval(x, p2), x))
! 81: {
! 82: y[++nbauto] = (long) p2;
! 83: if (DEBUGLEVEL > 1)
! 84: fprintferr("conjugate %ld: %Z\n", i, y[nbauto]);
! 85: }
! 86: }
! 87: }
! 88: setlg(y, 1 + nbauto);
! 89: return gerepileupto(av, gen_sort(y, 0, cmp_pol));
! 90: }
! 91:
! 92: GEN
! 93: galoisconj2(GEN nf, long nbmax, long prec)
! 94: {
! 95: long av = avma, i, j, n, r1, ru, nbauto;
! 96: GEN x, y, w, polr, p1, p2;
! 97: if (typ(nf) == t_POL)
! 98: return galoisconj2pol(nf, nbmax, prec);
! 99: nf = checknf(nf);
! 100: x = (GEN) nf[1];
! 101: n = degpol(x);
! 102: if (n <= 0)
! 103: return cgetg(1, t_VEC);
! 104: r1 = nf_get_r1(nf);
! 105: p1 = (GEN) nf[6];
! 106: prec = precision((GEN) p1[1]);
! 107: /* accuracy in decimal digits */
! 108: prec = (long) (bit_accuracy(prec) * L2SL10 * 0.75);
! 109: ru = (n + r1) >> 1;
! 110: nbauto = 1;
! 111: polr = cgetg(n + 1, t_VEC);
! 112: for (i = 1; i <= r1; i++)
! 113: polr[i] = p1[i];
! 114: for (j = i; i <= ru; i++)
! 115: {
! 116: polr[j++] = p1[i];
! 117: polr[j++] = lconj((GEN) p1[i]);
! 118: }
! 119: p2 = gmael(nf, 5, 1);
! 120: w = cgetg(n + 2, t_VEC);
! 121: for (i = 1; i <= n; i++)
! 122: w[i] = coeff(p2, 1, i);
! 123: y = cgetg(nbmax + 1, t_COL);
! 124: y[1] = (long) polx[varn(x)];
! 125: for (i = 2; i <= n && nbauto < nbmax; i++)
! 126: {
! 127: w[n + 1] = polr[i];
! 128: p1 = lindep2(w, prec);
! 129: if (signe(p1[n + 1]))
! 130: {
! 131: setlg(p1, n + 1);
! 132: settyp(p1, t_COL);
! 133: p2 = gdiv(gmul((GEN) nf[7], p1), negi((GEN) p1[n + 1]));
! 134: if (gdivise(poleval(x, p2), x))
! 135: {
! 136: y[++nbauto] = (long) p2;
! 137: if (DEBUGLEVEL > 1)
! 138: fprintferr("conjugate %ld: %Z\n", i, y[nbauto]);
! 139: }
! 140: }
! 141: }
! 142: setlg(y, 1 + nbauto);
! 143: return gerepileupto(av, gen_sort(y, 0, cmp_pol));
! 144: }
! 145: /*************************************************************************/
! 146: /** **/
! 147: /** GALOISCONJ4 **/
! 148: /** **/
! 149: /** **/
! 150: /*************************************************************************/
! 151: /*DEBUGLEVEL:
! 152: 1: timing
! 153: 2: outline
! 154: 4: complete outline
! 155: 6: detail
! 156: 7: memory
! 157: 9: complete detail
! 158: */
! 159: /* retourne la permutation identite */
! 160: GEN
! 161: permidentity(long l)
! 162: {
! 163: GEN perm;
! 164: int i;
! 165: perm = cgetg(l + 1, t_VECSMALL);
! 166: for (i = 1; i <= l; i++)
! 167: perm[i] = i;
! 168: return perm;
! 169: }
! 170:
! 171: GEN vandermondeinverseprepold(GEN T, GEN L)
! 172: {
! 173: int i, n = lg(L);
! 174: GEN V, dT;
! 175: dT = derivpol(T);
! 176: V = cgetg(n, t_VEC);
! 177: for (i = 1; i < n; i++)
! 178: V[i] = (long) poleval(dT, (GEN) L[i]);
! 179: return V;
! 180: }
! 181:
! 182: GEN vandermondeinverseprep(GEN T, GEN L)
! 183: {
! 184: int i, j, n = lg(L);
! 185: GEN V;
! 186: V = cgetg(n, t_VEC);
! 187: for (i = 1; i < n; i++)
! 188: {
! 189: ulong ltop=avma;
! 190: GEN W=cgetg(n,t_VEC);
! 191: for (j = 1; j < n; j++)
! 192: if (i==j)
! 193: W[j]=un;
! 194: else
! 195: W[j]=lsub((GEN)L[i],(GEN)L[j]);
! 196: V[i]=lpileupto(ltop,divide_conquer_prod(W,&gmul));
! 197: }
! 198: return V;
! 199: }
! 200:
! 201: /* Calcule l'inverse de la matrice de van der Monde de T multiplie par den */
! 202: GEN
! 203: vandermondeinverse(GEN L, GEN T, GEN den, GEN prep)
! 204: {
! 205: ulong ltop = avma, lbot;
! 206: int i, j, n = lg(L);
! 207: long x = varn(T);
! 208: GEN M, P;
! 209: if (!prep)
! 210: prep=vandermondeinverseprep(T,L);
! 211: M = cgetg(n, t_MAT);
! 212: for (i = 1; i < n; i++)
! 213: {
! 214: M[i] = lgetg(n, t_COL);
! 215: P = gdiv(gdeuc(T, gsub(polx[x], (GEN) L[i])), (GEN) prep[i]);
! 216: for (j = 1; j < n; j++)
! 217: ((GEN *) M)[i][j] = P[1 + j];
! 218: }
! 219: lbot = avma;
! 220: M = gmul(den, M);
! 221: return gerepile(ltop, lbot, M);
! 222: }
! 223:
! 224: /* Calcule les bornes sur les coefficients a chercher */
! 225: struct galois_borne
! 226: {
! 227: GEN l;
! 228: long valsol;
! 229: long valabs;
! 230: GEN bornesol;
! 231: GEN ladicsol;
! 232: GEN ladicabs;
! 233: GEN lbornesol;
! 234: };
! 235:
! 236:
! 237: GEN indexpartial(GEN,GEN);
! 238: GEN ZX_disc_all(GEN,long);
! 239:
! 240: GEN
! 241: galoisborne(GEN T, GEN dn, struct galois_borne *gb, long ppp)
! 242: {
! 243: ulong ltop = avma, av2;
! 244: GEN borne, borneroots, borneabs;
! 245: int i, j;
! 246: int n;
! 247: GEN L, M, z, prep, den;
! 248: long prec;
! 249: n = degpol(T);
! 250: prec = 1;
! 251: for (i = 2; i < lgef(T); i++)
! 252: if (lg(T[i]) > prec)
! 253: prec = lg(T[i]);
! 254: /*prec++;*/
! 255: if (DEBUGLEVEL>=4) gentimer(3);
! 256: L = roots(T, prec);
! 257: if (DEBUGLEVEL>=4) genmsgtimer(3,"roots");
! 258: for (i = 1; i <= n; i++)
! 259: {
! 260: z = (GEN) L[i];
! 261: if (signe(z[2]))
! 262: break;
! 263: L[i] = z[1];
! 264: }
! 265: if (DEBUGLEVEL>=4) gentimer(3);
! 266: prep=vandermondeinverseprep(T, L);
! 267: if (!dn)
! 268: {
! 269: GEN res = divide_conquer_prod(gabs(prep,prec), mpmul);
! 270: GEN dis;
! 271: disable_dbg(0);
! 272: dis = ZX_disc_all(T, 1+mylogint(res,gdeux));
! 273: disable_dbg(-1);
! 274: den = gclone(indexpartial(T,dis));
! 275: }
! 276: else
! 277: den=dn;
! 278: M = vandermondeinverse(L, gmul(T, realun(prec)), den, prep);
! 279: if (DEBUGLEVEL>=4) genmsgtimer(3,"vandermondeinverse");
! 280: borne = realzero(prec);
! 281: for (i = 1; i <= n; i++)
! 282: {
! 283: z = gzero;
! 284: for (j = 1; j <= n; j++)
! 285: z = gadd(z, gabs(gcoeff(M,i,j), prec));
! 286: if (gcmp(z, borne) > 0)
! 287: borne = z;
! 288: }
! 289: borneroots = realzero(prec);
! 290: for (i = 1; i <= n; i++)
! 291: {
! 292: z = gabs((GEN) L[i], prec);
! 293: if (gcmp(z, borneroots) > 0)
! 294: borneroots = z;
! 295: }
! 296: borneabs = addsr(1, gmulsg(n, gpowgs(borneroots, n/ppp)));
! 297: /*if (ppp == 1)
! 298: borneabs = addsr(1, gmulsg(n, gpowgs(borneabs, 2)));*/
! 299: borneroots = addsr(1, gmul(borne, borneroots));
! 300: av2 = avma;
! 301: /*We use d-1 test, so we must overlift to 2^BITS_IN_LONG*/
! 302: gb->valsol = mylogint(gmul2n(borneroots,2+BITS_IN_LONG), gb->l);
! 303: gb->valabs = mylogint(gmul2n(borneabs,2), gb->l);
! 304: gb->valabs = max(gb->valsol,gb->valabs);
! 305: if (DEBUGLEVEL >= 4)
! 306: fprintferr("GaloisConj:val1=%ld val2=%ld\n", gb->valsol, gb->valabs);
! 307: avma = av2;
! 308: gb->bornesol = gerepileupto(ltop, ceil_safe(borneroots));
! 309: if (DEBUGLEVEL >= 9)
! 310: fprintferr("GaloisConj: Bound %Z\n",borneroots);
! 311: gb->ladicsol = gpowgs(gb->l, gb->valsol);
! 312: gb->ladicabs = gpowgs(gb->l, gb->valabs);
! 313: gb->lbornesol = subii(gb->ladicsol,gb->bornesol);
! 314: if (!dn)
! 315: {
! 316: dn=forcecopy(den);
! 317: gunclone(den);
! 318: }
! 319: return dn;
! 320: }
! 321:
! 322: struct galois_lift
! 323: {
! 324: GEN T;
! 325: GEN den;
! 326: GEN p;
! 327: GEN L;
! 328: GEN Lden;
! 329: long e;
! 330: GEN Q;
! 331: GEN TQ;
! 332: struct galois_borne *gb;
! 333: };
! 334: /* Initialise la structure galois_lift */
! 335: GEN makeLden(GEN L,GEN den, struct galois_borne *gb)
! 336: {
! 337: ulong ltop=avma;
! 338: long i,l=lg(L);
! 339: GEN Lden;
! 340: Lden=cgetg(l,t_VEC);
! 341: for (i=1;i<l;i++)
! 342: Lden[i]=lmulii((GEN)L[i],den);
! 343: for (i=1;i<l;i++)
! 344: Lden[i]=lmodii((GEN)Lden[i],gb->ladicsol);
! 345: return gerepileupto(ltop,Lden);
! 346: }
! 347: void
! 348: initlift(GEN T, GEN den, GEN p, GEN L, GEN Lden, struct galois_borne *gb, struct galois_lift *gl)
! 349: {
! 350: ulong ltop, lbot;
! 351: gl->gb=gb;
! 352: gl->T = T;
! 353: gl->den = den;
! 354: gl->p = p;
! 355: gl->L = L;
! 356: gl->Lden = Lden;
! 357: ltop = avma;
! 358: gl->e = mylogint(gmul2n(gb->bornesol, 2+BITS_IN_LONG),p);
! 359: gl->e = max(2,gl->e);
! 360: lbot = avma;
! 361: gl->Q = gpowgs(p, gl->e);
! 362: gl->Q = gerepile(ltop, lbot, gl->Q);
! 363: gl->TQ = FpX_red(T,gl->Q);
! 364: }
! 365: #if 0
! 366: /*
! 367: * T est le polynome \sum t_i*X^i S est un Polmod T
! 368: * Retourne un vecteur Spow
! 369: * verifiant la condition: i>=1 et t_i!=0 ==> Spow[i]=S^(i-1)*t_i
! 370: * Essaie d'etre efficace sur les polynomes lacunaires
! 371: */
! 372: GEN
! 373: compoTS(GEN T, GEN S, GEN Q, GEN mod)
! 374: {
! 375: GEN Spow;
! 376: int i;
! 377: Spow = cgetg(lgef(T) - 2, t_VEC);
! 378: for (i = 3; i < lg(Spow); i++)
! 379: Spow[i] = 0;
! 380: Spow[1] = (long) polun[varn(S)];
! 381: Spow[2] = (long) S;
! 382: for (i = 2; i < lg(Spow) - 1; i++)
! 383: {
! 384: if (signe((GEN) T[i + 3]))
! 385: for (;;)
! 386: {
! 387: int k, l;
! 388: for (k = 1; k <= (i >> 1); k++)
! 389: if (Spow[k + 1] && Spow[i - k + 1])
! 390: break;
! 391: if ((k << 1) < i)
! 392: {
! 393: Spow[i + 1] = (long) FpXQ_mul((GEN) Spow[k + 1], (GEN) Spow[i - k + 1],Q,mod);
! 394: break;
! 395: }
! 396: else if ((k << 1) == i)
! 397: {
! 398: Spow[i + 1] = (long) FpXQ_sqr((GEN) Spow[k + 1],Q,mod);
! 399: break;
! 400: }
! 401: for (k = i - 1; k > 0; k--)
! 402: if (Spow[k + 1])
! 403: break;
! 404: if ((k << 1) < i)
! 405: {
! 406: Spow[(k << 1) + 1] = (long) FpXQ_sqr((GEN) Spow[k + 1],Q,mod);
! 407: continue;
! 408: }
! 409: for (l = i - k; l > 0; l--)
! 410: if (Spow[l + 1])
! 411: break;
! 412: if (Spow[i - l - k + 1])
! 413: Spow[i - k + 1] = (long) FpXQ_mul((GEN) Spow[i - l - k + 1], (GEN) Spow[l + 1],Q,mod);
! 414: else
! 415: Spow[l + k + 1] = (long) FpXQ_mul((GEN) Spow[k + 1], (GEN) Spow[l + 1],Q,mod);
! 416: }
! 417: }
! 418: for (i = 1; i < lg(Spow); i++)
! 419: if (signe((GEN) T[i + 2]))
! 420: Spow[i] = (long) FpX_Fp_mul((GEN) Spow[i], (GEN) T[i + 2],mod);
! 421: return Spow;
! 422: }
! 423:
! 424: /*
! 425: * Calcule T(S) a l'aide du vecteur Spow
! 426: */
! 427: static GEN
! 428: calcTS(GEN Spow, GEN P, GEN S, GEN Q, GEN mod)
! 429: {
! 430: GEN res = gzero;
! 431: int i;
! 432: for (i = 1; i < lg(Spow); i++)
! 433: if (signe((GEN) P[i + 2]))
! 434: res = FpX_add(res, (GEN) Spow[i],NULL);
! 435: res = FpXQ_mul(res,S,Q,mod);
! 436: res=FpX_Fp_add(res,(GEN)P[2],mod);
! 437: return res;
! 438: }
! 439:
! 440: /*
! 441: * Calcule T'(S) a l'aide du vecteur Spow
! 442: */
! 443: static GEN
! 444: calcderivTS(GEN Spow, GEN P, GEN mod)
! 445: {
! 446: GEN res = gzero;
! 447: int i;
! 448: for (i = 1; i < lg(Spow); i++)
! 449: if (signe((GEN) P[i + 2]))
! 450: res = FpX_add(res, FpX_Fp_mul((GEN) Spow[i], stoi(i),mod),NULL);
! 451: return FpX_red(res,mod);
! 452: }
! 453: #endif
! 454: /*
! 455: * Verifie que f est une solution presque surement et calcule sa permutation
! 456: */
! 457: static int
! 458: poltopermtest(GEN f, struct galois_lift *gl, GEN pf)
! 459: {
! 460: ulong ltop;
! 461: GEN fx, fp;
! 462: long i, j,ll;
! 463: for (i = 2; i< lgef(f); i++)
! 464: if (cmpii((GEN)f[i],gl->gb->bornesol)>0
! 465: && cmpii((GEN)f[i],gl->gb->lbornesol)<0)
! 466: {
! 467: if (DEBUGLEVEL>=4)
! 468: fprintferr("GaloisConj: Solution too large, discard it.\n");
! 469: return 0;
! 470: }
! 471: ll=lg(gl->L);
! 472: fp = cgetg(ll, t_VECSMALL);
! 473: ltop = avma;
! 474: for (i = 1; i < ll; i++)
! 475: fp[i] = 1;
! 476: for (i = 1; i < ll; i++)
! 477: {
! 478: fx = FpX_eval(f, (GEN) gl->L[i],gl->gb->ladicsol);
! 479: for (j = 1; j < ll; j++)
! 480: {
! 481: if (fp[j] && egalii(fx, (GEN) gl->Lden[j]))
! 482: {
! 483: pf[i] = j;
! 484: fp[j] = 0;
! 485: break;
! 486: }
! 487: }
! 488: if (j == ll)
! 489: return 0;
! 490: avma = ltop;
! 491: }
! 492: return 1;
! 493: }
! 494:
! 495: GEN polratlift(GEN P, GEN mod, GEN amax, GEN bmax, GEN denom);
! 496:
! 497: /*
! 498: * Soit P un polynome de \ZZ[X] , p un nombre premier , S\in\FF_p[X]/(Q) tel
! 499: * que P(S)=0 [p,Q] Relever S en S_0 tel que P(S_0)=0 [p^e,Q]
! 500: * Unclean stack.
! 501: */
! 502: static long monoratlift(GEN S, GEN q, GEN qm1old,struct galois_lift *gl, GEN frob)
! 503: {
! 504: ulong ltop=avma;
! 505: GEN tlift = polratlift(S,q,qm1old,qm1old,gl->den);
! 506: if (tlift)
! 507: {
! 508: if(DEBUGLEVEL>=4)
! 509: fprintferr("MonomorphismLift: trying early solution %Z\n",tlift);
! 510: /*Rationals coefficients*/
! 511: tlift=lift(gmul(tlift,gmodulcp(gl->den,gl->gb->ladicsol)));
! 512: if (poltopermtest(tlift, gl, frob))
! 513: {
! 514: if(DEBUGLEVEL>=4)
! 515: fprintferr("MonomorphismLift: true early solution.\n");
! 516: avma=ltop;
! 517: return 1;
! 518: }
! 519: if(DEBUGLEVEL>=4)
! 520: fprintferr("MonomorphismLift: false early solution.\n");
! 521: }
! 522: avma=ltop;
! 523: return 0;
! 524: }
! 525: GEN
! 526: monomorphismratlift(GEN P, GEN S, struct galois_lift *gl, GEN frob)
! 527: {
! 528: ulong ltop, lbot;
! 529: long rt;
! 530: GEN Q=gl->T, p=gl->p;
! 531: long e=gl->e, level=1;
! 532: long x;
! 533: GEN q, qold, qm1, qm1old;
! 534: GEN W, Pr, Qr, Sr, Wr = gzero, Prold, Qrold, Spow;
! 535: long i,nb,mask;
! 536: GEN *gptr[2];
! 537: if (DEBUGLEVEL == 1)
! 538: timer2();
! 539: x = varn(P);
! 540: rt = brent_kung_optpow(degpol(Q),1);
! 541: q = p; qm1 = gun; /*during the run, we have p*qm1=q*/
! 542: nb=hensel_lift_accel(e, &mask);
! 543: Pr = FpX_red(P,q);
! 544: Qr = (P==Q)?Pr:FpX_red(Q, q);/*A little speed up for automorphismlift*/
! 545: W=FpX_FpXQ_compo(deriv(Pr, x),S,Qr,q);
! 546: W=FpXQ_inv(W,Qr,q);
! 547: qold = p; qm1old=gun;
! 548: Prold = Pr;
! 549: Qrold = Qr;
! 550: gptr[0] = &S;
! 551: gptr[1] = &Wr;
! 552: for (i=0; i<nb;i++)
! 553: {
! 554: if (DEBUGLEVEL>=2)
! 555: {
! 556: level=(level<<1)-((mask>>i)&1);
! 557: timer2();
! 558: }
! 559: qm1 = (mask&(1<<i))?sqri(qm1):mulii(qm1, q);
! 560: q = mulii(qm1, p);
! 561: Pr = FpX_red(P, q);
! 562: Qr = (P==Q)?Pr:FpX_red(Q, q);/*A little speed up for automorphismlift*/
! 563: ltop = avma;
! 564: Sr = S;
! 565: /*Spow = compoTS(Pr, Sr, Qr, q);*/
! 566: Spow = FpXQ_powers(Sr, rt, Qr, q);
! 567:
! 568: if (i)
! 569: {
! 570: /*W = FpXQ_mul(Wr, calcderivTS(Spow, Prold,qold), Qrold, qold);*/
! 571: W = FpXQ_mul(Wr, FpX_FpXQV_compo(deriv(Pr,-1),FpXV_red(Spow,qold),Qrold,qold), Qrold, qold);
! 572: W = FpX_neg(W, qold);
! 573: W = FpX_Fp_add(W, gdeux, qold);
! 574: W = FpXQ_mul(Wr, W, Qrold, qold);
! 575: }
! 576: Wr = W;
! 577: /*S = FpXQ_mul(Wr, calcTS(Spow, Pr, Sr, Qr, q),Qr,q);*/
! 578: S = FpXQ_mul(Wr, FpX_FpXQV_compo(Pr, Spow, Qr, q),Qr,q);
! 579: S = FpX_sub(Sr, S, NULL);
! 580: lbot = avma;
! 581: Wr = gcopy(Wr);
! 582: S = FpX_red(S, q);
! 583: gerepilemanysp(ltop, lbot, gptr, 2);
! 584: if (i && i<nb-1 && frob && monoratlift(S,q,qm1old,gl,frob))
! 585: return NULL;
! 586: qold = q; qm1old=qm1; Prold = Pr; Qrold = Qr;
! 587: if (DEBUGLEVEL >= 2)
! 588: msgtimer("MonomorphismLift: lift to prec %d",level);
! 589: }
! 590: if (DEBUGLEVEL == 1)
! 591: msgtimer("monomorphismlift()");
! 592: return S;
! 593: }
! 594: /*
! 595: * Soit T un polynome de \ZZ[X] , p un nombre premier , S\in\FF_p[X]/(T) tel
! 596: * que T(S)=0 [p,T] Relever S en S_0 tel que T(S_0)=0 [T,p^e]
! 597: * Unclean stack.
! 598: */
! 599: GEN
! 600: automorphismlift(GEN S, struct galois_lift *gl, GEN frob)
! 601: {
! 602: return monomorphismratlift(gl->T, S, gl, frob);
! 603: }
! 604: GEN
! 605: monomorphismlift(GEN P, GEN S, GEN Q, GEN p, long e)
! 606: {
! 607: struct galois_lift gl;
! 608: gl.T=Q;
! 609: gl.p=p;
! 610: gl.e=e;
! 611: return monomorphismratlift(P,S,&gl,NULL);
! 612: }
! 613:
! 614: struct galois_testlift
! 615: {
! 616: long n;
! 617: long f;
! 618: long g;
! 619: GEN bezoutcoeff;
! 620: GEN pauto;
! 621: GEN C;
! 622: GEN Cd;
! 623: };
! 624: GEN bezout_lift_fact(GEN T, GEN Tmod, GEN p, long e);
! 625: static GEN
! 626: galoisdolift(struct galois_lift *gl, GEN frob)
! 627: {
! 628: ulong ltop=avma;
! 629: long v = varn(gl->T);
! 630: GEN Tp = FpX_red(gl->T, gl->p);
! 631: GEN S = FpXQ_pow(polx[v],gl->p, Tp,gl->p);
! 632: GEN plift = automorphismlift(S, gl, frob);
! 633: return gerepileupto(ltop,plift);
! 634: }
! 635:
! 636: static void
! 637: inittestlift( GEN plift, GEN Tmod, struct galois_lift *gl,
! 638: struct galois_testlift *gt)
! 639: {
! 640: long v = varn(gl->T);
! 641: gt->n = lg(gl->L) - 1;
! 642: gt->g = lg(Tmod) - 1;
! 643: gt->f = gt->n / gt->g;
! 644: gt->bezoutcoeff = bezout_lift_fact(gl->T, Tmod, gl->p, gl->e);
! 645: gt->pauto = cgetg(gt->f + 1, t_VEC);
! 646: gt->pauto[1] = (long) polx[v];
! 647: gt->pauto[2] = (long) gcopy(plift);
! 648: if (gt->f > 2)
! 649: {
! 650: ulong ltop=avma;
! 651: int i;
! 652: long nautpow=brent_kung_optpow(gt->n-1,gt->f-2);
! 653: GEN autpow;
! 654: if (DEBUGLEVEL >= 1) timer2();
! 655: autpow = FpXQ_powers(plift,nautpow,gl->TQ,gl->Q);
! 656: for (i = 3; i <= gt->f; i++)
! 657: gt->pauto[i] = (long) FpX_FpXQV_compo((GEN)gt->pauto[i-1],autpow,gl->TQ,gl->Q);
! 658: /*Somewhat paranoid with memory, but this function use a lot of stack.*/
! 659: gt->pauto=gerepileupto(ltop, gt->pauto);
! 660: if (DEBUGLEVEL >= 1) msgtimer("frobenius power");
! 661: }
! 662: }
! 663:
! 664: long intheadlong(GEN x, GEN mod)
! 665: {
! 666: ulong ltop=avma;
! 667: GEN r;
! 668: int s;
! 669: long res;
! 670: r=divii(shifti(x,BITS_IN_LONG),mod);
! 671: s=signe(r);
! 672: res=s?s>0?r[2]:-r[2]:0;
! 673: avma=ltop;
! 674: return res;
! 675: }
! 676:
! 677: GEN matheadlong(GEN W, GEN mod)
! 678: {
! 679: long i,j;
! 680: GEN V=cgetg(lg(W),t_VEC);
! 681: for(i=1;i<lg(W);i++)
! 682: {
! 683: V[i]=lgetg(lg(W[i]),t_VECSMALL);
! 684: for(j=1;j<lg(W[i]);j++)
! 685: mael(V,i,j)=intheadlong(gmael(W,i,j),mod);
! 686: }
! 687: return V;
! 688: }
! 689:
! 690: long polheadlong(GEN P, long n, GEN mod)
! 691: {
! 692: return (lgef(P)>n+2)?intheadlong((GEN)P[n+2],mod):0;
! 693: }
! 694: /*
! 695: *
! 696: */
! 697: long
! 698: frobeniusliftall(GEN sg, long el, GEN *psi, struct galois_lift *gl,
! 699: struct galois_testlift *gt, GEN frob)
! 700: {
! 701: ulong ltop = avma, av, ltop2;
! 702: long d, z, m, c, n, ord;
! 703: int i, j, k;
! 704: GEN pf, u, v;
! 705: GEN C, Cd;
! 706: long Z, Zv;
! 707: long stop=0;
! 708: GEN NN,NQ,NR;
! 709: long N1,N2,R1,Ni;
! 710: m = gt->g;
! 711: ord = gt->f;
! 712: n = lg(gl->L) - 1;
! 713: c = lg(sg) - 1;
! 714: d = m / c;
! 715: pf = cgetg(m, t_VECSMALL);
! 716: *psi = pf;
! 717: ltop2 = avma;
! 718: NN = gdiv(mpfact(m), gmul(stoi(c), gpowgs(mpfact(d), c)));
! 719: if (DEBUGLEVEL >= 4)
! 720: fprintferr("GaloisConj:I will try %Z permutations\n", NN);
! 721: N1=10000000;
! 722: NQ=dvmdis(NN,N1,&NR);
! 723: if (cmpis(NQ,1000000000)>0)
! 724: {
! 725: err(warner,"Combinatorics too hard : would need %Z tests!\n
! 726: I will skip it,but it may induce galoisinit to loop",NN);
! 727: avma = ltop;
! 728: *psi = NULL;
! 729: return 0;
! 730: }
! 731: N2=itos(NQ); R1=itos(NR); if(!N2) N1=R1;
! 732: if (DEBUGLEVEL>=4)
! 733: {
! 734: stop=N1/20;
! 735: timer2();
! 736: }
! 737: avma = ltop2;
! 738: C=gt->C;
! 739: Cd=gt->Cd;
! 740: v = FpX_Fp_mul(FpXQ_mul((GEN)gt->pauto[1+el%ord], (GEN)
! 741: gt->bezoutcoeff[m],gl->TQ,gl->Q),gl->den,gl->Q);
! 742: Zv=polheadlong(v,0,gl->Q);
! 743: for (i = 1; i < m; i++)
! 744: pf[i] = 1 + i / d;
! 745: av = avma;
! 746: for (Ni = 0, i = 0; ;i++)
! 747: {
! 748: Z=Zv;
! 749: for (j = 1; j < m; j++)
! 750: {
! 751: long h;
! 752: h=(el*sg[pf[j]])%ord + 1;
! 753: if (!mael(C,h,j))
! 754: {
! 755: ulong av3=avma;
! 756: GEN r;
! 757: r=FpX_Fp_mul(FpXQ_mul((GEN) gt->pauto[h],
! 758: (GEN) gt->bezoutcoeff[j],gl->TQ,gl->Q),gl->den,gl->Q);
! 759: mael(C,h,j) = lclone(r);
! 760: mael(Cd,h,j) = polheadlong(r,0,gl->Q);
! 761: avma=av3;
! 762: }
! 763: Z += mael(Cd,h,j);
! 764: }
! 765: if (labs(Z)<=n )
! 766: {
! 767: Z=polheadlong(v,1,gl->Q);
! 768: for (j = 1; j < m; j++)
! 769: Z += polheadlong(gmael(C,(el*sg[pf[j]])%ord + 1,j),1,gl->Q);
! 770: if (labs(Z)<=n )
! 771: {
! 772: u = v;
! 773: for (j = 1; j < m; j++)
! 774: u = FpX_add(u, gmael(C,1+(el*sg[pf[j]])%ord,j),NULL);
! 775: u = FpX_center(FpX_red(u, gl->Q), gl->Q);
! 776: if (poltopermtest(u, gl, frob))
! 777: {
! 778: if (DEBUGLEVEL >= 4 )
! 779: msgtimer("");
! 780: avma = ltop2;
! 781: return 1;
! 782: }
! 783: }
! 784: }
! 785: if (DEBUGLEVEL >= 4 && i==stop)
! 786: {
! 787: stop+=N1/20;
! 788: msgtimer("GaloisConj:Testing %Z", addis(mulss(Ni,N1),i));
! 789: }
! 790: avma = av;
! 791: if (i == N1 - 1)
! 792: {
! 793: if (Ni==N2-1)
! 794: N1=R1;
! 795: if (Ni==N2)
! 796: break;
! 797: Ni++;
! 798: i=0;
! 799: if (DEBUGLEVEL>=4)
! 800: {
! 801: stop=N1/20;
! 802: timer2();
! 803: }
! 804: }
! 805: for (j = 2; j < m && pf[j - 1] >= pf[j]; j++);
! 806: for (k = 1; k < j - k && pf[k] != pf[j - k]; k++)
! 807: {
! 808: z = pf[k];
! 809: pf[k] = pf[j - k];
! 810: pf[j - k] = z;
! 811: }
! 812: for (k = j - 1; pf[k] >= pf[j]; k--);
! 813: z = pf[j];
! 814: pf[j] = pf[k];
! 815: pf[k] = z;
! 816: }
! 817: avma = ltop;
! 818: *psi = NULL;
! 819: return 0;
! 820: }
! 821: /*
! 822: * applique une permutation t a un vecteur s sans duplication
! 823: * Propre si s est un VECSMALL
! 824: */
! 825: static GEN
! 826: permapply(GEN s, GEN t)
! 827: {
! 828: GEN u;
! 829: int i;
! 830: if (lg(s) < lg(t))
! 831: err(talker, "First permutation shorter than second in permapply");
! 832: u = cgetg(lg(s), typ(s));
! 833: for (i = 1; i < lg(s); i++)
! 834: u[i] = s[t[i]];
! 835: return u;
! 836: }
! 837:
! 838: /*
! 839: * alloue une ardoise pour n entiers de longueurs pour les test de
! 840: * permutation
! 841: */
! 842: GEN
! 843: alloue_ardoise(long n, long s)
! 844: {
! 845: GEN ar;
! 846: long i;
! 847: ar = cgetg(n + 1, t_VEC);
! 848: for (i = 1; i <= n; i++)
! 849: ar[i] = lgetg(s, t_INT);
! 850: return ar;
! 851: }
! 852:
! 853: /*
! 854: * structure contenant toutes les données pour le tests des permutations:
! 855: *
! 856: * ordre :ordre des tests pour verifie_test ordre[lg(ordre)]: numero du test
! 857: * principal borne : borne sur les coefficients a trouver ladic: modulo
! 858: * l-adique des racines lborne:ladic-borne TM:vecteur des ligne de M
! 859: * PV:vecteur des clones des matrices de test (Vmatrix) (ou NULL si non
! 860: * calcule) L,M comme d'habitude (voir plus bas)
! 861: */
! 862: struct galois_test
! 863: {
! 864: GEN ordre;
! 865: GEN borne, lborne, ladic;
! 866: GEN PV, TM;
! 867: GEN L, M;
! 868: };
! 869: /* Calcule la matrice de tests correspondant a la n-ieme ligne de V */
! 870: static GEN
! 871: Vmatrix(long n, struct galois_test *td)
! 872: {
! 873: ulong ltop = avma, lbot;
! 874: GEN V;
! 875: long i;
! 876: V = cgetg(lg(td->L), t_VEC);
! 877: for (i = 1; i < lg(V); i++)
! 878: V[i] = mael(td->M,i,n);
! 879: V = gmul(td->L, V);
! 880: lbot = avma;
! 881: V = gmod(V, td->ladic);
! 882: return gerepile(ltop, lbot, V);
! 883: }
! 884:
! 885: /*
! 886: * Initialise la structure galois_test
! 887: */
! 888: static void
! 889: inittest(GEN L, GEN M, GEN borne, GEN ladic, struct galois_test *td)
! 890: {
! 891: ulong ltop;
! 892: long i;
! 893: int n = lg(L) - 1;
! 894: if (DEBUGLEVEL >= 8)
! 895: fprintferr("GaloisConj:Entree Init Test\n");
! 896: td->ordre = cgetg(n + 1, t_VECSMALL);
! 897: for (i = 1; i <= n - 2; i++)
! 898: td->ordre[i] = i + 2;
! 899: for (; i <= n; i++)
! 900: td->ordre[i] = i - n + 2;
! 901: td->borne = borne;ltop = avma;
! 902: td->lborne = subii(ladic, borne);
! 903: td->ladic = ladic;
! 904: td->L = L;
! 905: td->M = M;
! 906: td->PV = cgetg(n + 1, t_VECSMALL); /* peut-etre t_VEC est correct ici */
! 907: for (i = 1; i <= n; i++)
! 908: td->PV[i] = 0;
! 909: ltop = avma;
! 910: td->PV[td->ordre[n]] = (long) gclone(Vmatrix(td->ordre[n], td));
! 911: avma = ltop;
! 912: td->TM = gtrans(M);
! 913: settyp(td->TM, t_VEC);
! 914: for (i = 1; i < lg(td->TM); i++)
! 915: settyp(td->TM[i], t_VEC);
! 916: if (DEBUGLEVEL >= 8)
! 917: fprintferr("GaloisConj:Sortie Init Test\n");
! 918: }
! 919:
! 920: /*
! 921: * liberer les clones de la structure galois_test
! 922: *
! 923: * Reservé a l'accreditation ultra-violet:Liberez les clones! Paranoia(tm)
! 924: */
! 925: static void
! 926: freetest(struct galois_test *td)
! 927: {
! 928: int i;
! 929: for (i = 1; i < lg(td->PV); i++)
! 930: if (td->PV[i])
! 931: {
! 932: gunclone((GEN) td->PV[i]);
! 933: td->PV[i] = 0;
! 934: }
! 935: }
! 936:
! 937: /*
! 938: * Test si le nombre padique P est proche d'un entier inferieur a td->borne
! 939: * en valeur absolue.
! 940: */
! 941: static long
! 942: padicisint(GEN P, struct galois_test *td)
! 943: {
! 944: long r;
! 945: ulong ltop = avma;
! 946: GEN U;
! 947: /*if (typ(P) != t_INT) err(typeer, "padicisint");*/
! 948: U = modii(P, td->ladic);
! 949: r = gcmp(U, (GEN) td->borne) <= 0 || gcmp(U, (GEN) td->lborne) >= 0;
! 950: avma = ltop;
! 951: return r;
! 952: }
! 953:
! 954: /*
! 955: * Verifie si pf est une vrai solution et non pas un "hop"
! 956: */
! 957: static long
! 958: verifietest(GEN pf, struct galois_test *td)
! 959: {
! 960: ulong av = avma;
! 961: GEN P, V;
! 962: int i, j;
! 963: int n = lg(td->L) - 1;
! 964: if (DEBUGLEVEL >= 8)
! 965: fprintferr("GaloisConj:Entree Verifie Test\n");
! 966: P = permapply(td->L, pf);
! 967: for (i = 1; i < n; i++)
! 968: {
! 969: long ord;
! 970: GEN PW;
! 971: ord = td->ordre[i];
! 972: PW = (GEN) td->PV[ord];
! 973: if (PW)
! 974: {
! 975: V = ((GEN **) PW)[1][pf[1]];
! 976: for (j = 2; j <= n; j++)
! 977: V = addii(V, ((GEN **) PW)[j][pf[j]]);
! 978: }
! 979: else
! 980: V = centermod(gmul((GEN) td->TM[ord], P),td->ladic);
! 981: if (!padicisint(V, td))
! 982: break;
! 983: }
! 984: if (i == n)
! 985: {
! 986: if (DEBUGLEVEL >= 8)
! 987: fprintferr("GaloisConj:Sortie Verifie Test:1\n");
! 988: avma = av;
! 989: return 1;
! 990: }
! 991: if (!td->PV[td->ordre[i]])
! 992: {
! 993: td->PV[td->ordre[i]] = (long) gclone(Vmatrix(td->ordre[i], td));
! 994: if (DEBUGLEVEL >= 4)
! 995: fprintferr("M");
! 996: }
! 997: if (DEBUGLEVEL >= 4)
! 998: fprintferr("%d.", i);
! 999: if (i > 1)
! 1000: {
! 1001: long z;
! 1002: z = td->ordre[i];
! 1003: for (j = i; j > 1; j--)
! 1004: td->ordre[j] = td->ordre[j - 1];
! 1005: td->ordre[1] = z;
! 1006: if (DEBUGLEVEL >= 8)
! 1007: fprintferr("%Z", td->ordre);
! 1008: }
! 1009: if (DEBUGLEVEL >= 8)
! 1010: fprintferr("GaloisConj:Sortie Verifie Test:0\n");
! 1011: avma = av;
! 1012: return 0;
! 1013: }
! 1014: /*Compute a*b/c when a*b will overflow*/
! 1015: static long muldiv(long a,long b,long c)
! 1016: {
! 1017: return (long)((double)a*(double)b/c);
! 1018: }
! 1019:
! 1020: /*
! 1021: * F est la decomposition en cycle de sigma, B est la decomposition en cycle
! 1022: * de cl(tau) Teste toutes les permutations pf pouvant correspondre a tau tel
! 1023: * que : tau*sigma*tau^-1=sigma^s et tau^d=sigma^t ou d est l'ordre de
! 1024: * cl(tau) --------- x: vecteur des choix y: vecteur des mises a jour G:
! 1025: * vecteur d'acces a F linéaire
! 1026: */
! 1027: GEN
! 1028: testpermutation(GEN F, GEN B, long s, long t, long cut,
! 1029: struct galois_test *td)
! 1030: {
! 1031: ulong av, avm = avma;
! 1032: int a, b, c, d, n;
! 1033: GEN pf, x, ar, y, *G;
! 1034: int p1, p2, p3, p4, p5, p6;
! 1035: long l1, l2, N1, N2, R1;
! 1036: long i, j, cx, hop = 0, start = 0;
! 1037: GEN W, NN, NQ, NR;
! 1038: long V;
! 1039: if (DEBUGLEVEL >= 1)
! 1040: timer2();
! 1041: a = lg(F) - 1;
! 1042: b = lg(F[1]) - 1;
! 1043: c = lg(B) - 1;
! 1044: d = lg(B[1]) - 1;
! 1045: n = a * b;
! 1046: s = (b + s) % b;
! 1047: pf = cgetg(n + 1, t_VECSMALL);
! 1048: av = avma;
! 1049: ar = cgetg(a + 1, t_VECSMALL);
! 1050: x = cgetg(a + 1, t_VECSMALL);
! 1051: y = cgetg(a + 1, t_VECSMALL);
! 1052: G = (GEN *) cgetg(a + 1, t_VECSMALL); /* Don't worry */
! 1053: W = matheadlong((GEN) td->PV[td->ordre[n]],td->ladic);
! 1054: for (cx = 1, i = 1, j = 1; cx <= a; cx++, i++)
! 1055: {
! 1056: x[cx] = (i != d) ? 0 : t;
! 1057: y[cx] = 1;
! 1058: G[cx] = (GEN) F[((long **) B)[j][i]]; /* Be happy */
! 1059: if (i == d)
! 1060: {
! 1061: i = 0;
! 1062: j++;
! 1063: }
! 1064: } /* Be happy now! */
! 1065: NN = divis(gpowgs(stoi(b), c * (d - 1)),cut);
! 1066: if (DEBUGLEVEL >= 4)
! 1067: fprintferr("GaloisConj:I will try %Z permutations\n", NN);
! 1068: N1=100000;
! 1069: NQ=dvmdis(NN,N1,&NR);
! 1070: if (cmpis(NQ,1000000000)>0)
! 1071: {
! 1072: avma=avm;
! 1073: err(warner,"Combinatorics too hard : would need %Z tests!\n I'll skip it but you will get a partial result...",NN);
! 1074: return permidentity(n);
! 1075: }
! 1076: N2=itos(NQ); R1=itos(NR);
! 1077: for (l2 = 0; l2 <= N2; l2++)
! 1078: {
! 1079: long nbiter = (l2<N2) ? N1: R1;
! 1080: if (DEBUGLEVEL >= 2 && N2)
! 1081: fprintferr("%d%% ", muldiv(l2,100,N2));
! 1082: for (l1 = 0; l1 < nbiter; l1++)
! 1083: {
! 1084: if (start)
! 1085: {
! 1086: for (i = 1, j = d; i < a;)
! 1087: {
! 1088: y[i] = 1;
! 1089: if ((++(x[i])) != b)
! 1090: break;
! 1091: x[i++] = 0;
! 1092: if (i == j)
! 1093: {
! 1094: y[i++] = 1;
! 1095: j += d;
! 1096: }
! 1097: }
! 1098: y[i + 1] = 1;
! 1099: }
! 1100: else start=1;
! 1101: for (p1 = 1, p5 = d; p1 <= a; p1++, p5++)
! 1102: if (y[p1])
! 1103: {
! 1104: if (p5 == d)
! 1105: {
! 1106: p5 = 0;
! 1107: p4 = 0;
! 1108: }
! 1109: else
! 1110: p4 = x[p1 - 1];
! 1111: if (p5 == d - 1)
! 1112: p6 = p1 + 1 - d;
! 1113: else
! 1114: p6 = p1 + 1;
! 1115: y[p1] = 0;
! 1116: V = 0;
! 1117: for (p2 = 1 + p4, p3 = 1 + x[p1]; p2 <= b; p2++)
! 1118: {
! 1119: V += mael(W,G[p6][p3],G[p1][p2]);
! 1120: p3 += s;
! 1121: if (p3 > b)
! 1122: p3 -= b;
! 1123: }
! 1124: p3 = 1 + x[p1] - s;
! 1125: if (p3 <= 0)
! 1126: p3 += b;
! 1127: for (p2 = p4; p2 >= 1; p2--)
! 1128: {
! 1129: V += mael(W,G[p6][p3],G[p1][p2]);
! 1130: p3 -= s;
! 1131: if (p3 <= 0)
! 1132: p3 += b;
! 1133: }
! 1134: ar[p1]=V;
! 1135: }
! 1136: V = 0;
! 1137: for (p1 = 1; p1 <= a; p1++)
! 1138: V += ar[p1];
! 1139: if (labs(V)<=n)
! 1140: {
! 1141: for (p1 = 1, p5 = d; p1 <= a; p1++, p5++)
! 1142: {
! 1143: if (p5 == d)
! 1144: {
! 1145: p5 = 0;
! 1146: p4 = 0;
! 1147: }
! 1148: else
! 1149: p4 = x[p1 - 1];
! 1150: if (p5 == d - 1)
! 1151: p6 = p1 + 1 - d;
! 1152: else
! 1153: p6 = p1 + 1;
! 1154: for (p2 = 1 + p4, p3 = 1 + x[p1]; p2 <= b; p2++)
! 1155: {
! 1156: pf[G[p1][p2]] = G[p6][p3];
! 1157: p3 += s;
! 1158: if (p3 > b)
! 1159: p3 -= b;
! 1160: }
! 1161: p3 = 1 + x[p1] - s;
! 1162: if (p3 <= 0)
! 1163: p3 += b;
! 1164: for (p2 = p4; p2 >= 1; p2--)
! 1165: {
! 1166: pf[G[p1][p2]] = G[p6][p3];
! 1167: p3 -= s;
! 1168: if (p3 <= 0)
! 1169: p3 += b;
! 1170: }
! 1171: }
! 1172: if (verifietest(pf, td))
! 1173: {
! 1174: if (DEBUGLEVEL >= 1)
! 1175: {
! 1176: GEN nb=addis(mulss(l2,N1),l1);
! 1177: msgtimer("testpermutation(%Z)", nb);
! 1178: if (DEBUGLEVEL >= 2 && hop)
! 1179: fprintferr("GaloisConj:%d hop sur %Z iterations\n", hop, nb);
! 1180: }
! 1181: avma = av;
! 1182: return pf;
! 1183: }
! 1184: else
! 1185: hop++;
! 1186: }
! 1187: }
! 1188: }
! 1189: if (DEBUGLEVEL >= 1)
! 1190: {
! 1191: msgtimer("testpermutation(%Z)", NN);
! 1192: if (DEBUGLEVEL >= 2 && hop)
! 1193: fprintferr("GaloisConj:%d hop sur %Z iterations\n", hop, NN);
! 1194: }
! 1195: avma = avm;
! 1196: return gzero;
! 1197: }
! 1198: /* Compute generators for the subgroup of (Z/nZ)* given in HNF.
! 1199: * I apologize for the following spec:
! 1200: * If zns=znstar(2) then
! 1201: * zn2=gtovecsmall((GEN)zns[2]);
! 1202: * zn3=lift((GEN)zns[3]);
! 1203: * gen and ord : VECSMALL of length lg(zn3).
! 1204: * the result is in gen.
! 1205: * ord contains the relatives orders of the generators.
! 1206: */
! 1207: GEN hnftogeneratorslist(long n, GEN zn2, GEN zn3, GEN lss, GEN gen, GEN ord)
! 1208: {
! 1209: ulong ltop=avma;
! 1210: int j,h;
! 1211: GEN m=stoi(n);
! 1212: for (j = 1; j < lg(gen); j++)
! 1213: {
! 1214: gen[j] = 1;
! 1215: for (h = 1; h < lg(lss); h++)
! 1216: gen[j] = mulssmod(gen[j], itos(powmodulo((GEN)zn3[h],gmael(lss,j,h),m)),n);
! 1217: ord[j] = zn2[j] / itos(gmael(lss,j,j));
! 1218: }
! 1219: avma=ltop;
! 1220: return gen;
! 1221: }
! 1222: /*in place sort. Return V for convenience only*/
! 1223: GEN sortvecsmall(GEN V)
! 1224: {
! 1225: long i,j,k,l,m;
! 1226: for(l=1;l<lg(V);l<<=1)
! 1227: for(j=1;(j>>1)<lg(V);j+=l<<1)
! 1228: for(k=j+l,i=j; i<k && k<j+(l<<1) && k<lg(V);)
! 1229: if (V[i]>V[k])
! 1230: {
! 1231: long z=V[k];
! 1232: for (m=k;m>i;m--)
! 1233: V[m]=V[m-1];
! 1234: V[i]=z;
! 1235: k++;
! 1236: }
! 1237: else
! 1238: i++;
! 1239: return V ;
! 1240: }
! 1241: GEN uniqvecsmall(GEN V)
! 1242: {
! 1243: ulong ltop=avma;
! 1244: GEN W;
! 1245: long i,j;
! 1246: if ( lg(V) == 1 ) return gcopy(V);
! 1247: W=cgetg(lg(V),t_VECSMALL);
! 1248: W[1]=V[1];
! 1249: for(i=2,j=1;i<lg(V);i++)
! 1250: if (V[i]!=W[j])
! 1251: W[++j]=V[i];
! 1252: setlg(W,j+1);
! 1253: return gerepileupto(ltop,W);
! 1254: }
! 1255: int pari_compare_lg(GEN x, GEN y)
! 1256: {
! 1257: return lg(x)-lg(y);
! 1258: }
! 1259: /* Compute the elements of the subgroup of (Z/nZ)* given in HNF.
! 1260: * I apologize for the following spec:
! 1261: * If zns=znstar(2) then
! 1262: * zn2=gtovecsmall((GEN)zns[2]);
! 1263: * zn3=lift((GEN)zns[3]);
! 1264: * card=cardinal of the subgroup(i.e phi(n)/det(lss))
! 1265: */
! 1266: GEN
! 1267: hnftoelementslist(long n, GEN zn2, GEN zn3, GEN lss, long card)
! 1268: {
! 1269: ulong ltop;
! 1270: GEN sg, gen, ord;
! 1271: int k, j, l;
! 1272: sg = cgetg(1 + card, t_VECSMALL);
! 1273: ltop=avma;
! 1274: gen = cgetg(lg(zn3), t_VECSMALL);
! 1275: ord = cgetg(lg(zn3), t_VECSMALL);
! 1276: sg[1] = 1;
! 1277: hnftogeneratorslist(n,zn2,zn3,lss,gen,ord);
! 1278: for (j = 1, l = 1; j < lg(gen); j++)
! 1279: {
! 1280: int c = l * (ord[j] - 1);
! 1281: for (k = 1; k <= c; k++) /* I like it */
! 1282: sg[++l] = mulssmod(sg[k], gen[j], n);
! 1283: }
! 1284: sortvecsmall(sg);
! 1285: avma=ltop;
! 1286: return sg;
! 1287: }
! 1288:
! 1289: /*
! 1290: * Calcule la liste des sous groupes de \ZZ/m\ZZ d'ordre divisant p et
! 1291: * retourne la liste de leurs elements
! 1292: */
! 1293: GEN
! 1294: listsousgroupes(long m, long p)
! 1295: {
! 1296: ulong ltop = avma;
! 1297: GEN zns, lss, sg, res, zn2, zn3;
! 1298: int k, card, i, phi;
! 1299: if (m == 2)
! 1300: {
! 1301: res = cgetg(2, t_VEC);
! 1302: sg = cgetg(2, t_VECSMALL);
! 1303: res[1] = (long) sg;
! 1304: sg[1] = 1;
! 1305: return res;
! 1306: }
! 1307: zns = znstar(stoi(m));
! 1308: phi = itos((GEN) zns[1]);
! 1309: zn2 = gtovecsmall((GEN)zns[2]);
! 1310: zn3 = lift((GEN)zns[3]);
! 1311: lss = subgrouplist((GEN) zns[2], 0);
! 1312: res = cgetg(lg(lss), t_VEC);
! 1313: for (k = 1, i = lg(lss) - 1; i >= 1; i--)
! 1314: {
! 1315: long av;
! 1316: av = avma;
! 1317: card = phi / itos(det((GEN) lss[i]));
! 1318: avma = av;
! 1319: if (p % card == 0)
! 1320: res[k++] = (long) hnftoelementslist(m,zn2,zn3,(GEN)lss[i],card);
! 1321: }
! 1322: setlg(res,k);
! 1323: res=gen_sort(res,0,&pari_compare_lg);
! 1324: return gerepileupto(ltop,res);
! 1325: }
! 1326:
! 1327: /* Compute the orbits decomposition of a permutation
! 1328: * or of a set of permutations given by a vector.
! 1329: * v must not be an empty vector, becaude there is no way then
! 1330: * to tell the length of the permutation.
! 1331: */
! 1332:
! 1333: GEN
! 1334: permorbite(GEN v)
! 1335: {
! 1336: ulong ltop = avma, lbot;
! 1337: int i, j, k, l, m, n, o, p, flag;
! 1338: GEN bit, cycle, cy, u;
! 1339: if (typ(v) == t_VECSMALL)
! 1340: {
! 1341: u = cgetg(2, t_VEC);
! 1342: u[1] = (long) v;
! 1343: v = u;
! 1344: }
! 1345: n = lg(v[1]);
! 1346: cycle = cgetg(n, t_VEC);
! 1347: bit = cgetg(n, t_VECSMALL);
! 1348: for (i = 1; i < n; i++)
! 1349: bit[i] = 0;
! 1350: for (k = 1, l = 1; k < n;)
! 1351: {
! 1352: for (j = 1; bit[j]; j++);
! 1353: cy = cgetg(n, t_VECSMALL);
! 1354: m = 1;
! 1355: cy[m++] = j;
! 1356: bit[j] = 1;
! 1357: k++;
! 1358: do
! 1359: {
! 1360: flag = 0;
! 1361: for (o = 1; o < lg(v); o++)
! 1362: {
! 1363: for (p = 1; p < m; p++) /* m varie! */
! 1364: {
! 1365: j = ((long **) v)[o][cy[p]];
! 1366: if (!bit[j])
! 1367: {
! 1368: flag = 1;
! 1369: bit[j] = 1;
! 1370: k++;
! 1371: cy[m++] = j;
! 1372: }
! 1373: }
! 1374: }
! 1375: }
! 1376: while (flag);
! 1377: setlg(cy, m);
! 1378: cycle[l++] = (long) cy;
! 1379: }
! 1380: setlg(cycle, l);
! 1381: lbot = avma;
! 1382: cycle = gcopy(cycle);
! 1383: return gerepile(ltop, lbot, cycle);
! 1384: }
! 1385:
! 1386: GEN
! 1387: fixedfieldpolsigma(GEN sigma, GEN p, GEN Tp, GEN sym, GEN deg, long g)
! 1388: {
! 1389: ulong ltop=avma;
! 1390: long i, j, npows;
! 1391: GEN s, f, pows;
! 1392: sigma=lift(gmul(sigma,gmodulsg(1,p)));
! 1393: f=polx[varn(sigma)];
! 1394: s=zeropol(varn(sigma));
! 1395: for(j=1;j<lg(sym);j++)
! 1396: if(sym[j])
! 1397: {
! 1398: s=FpX_add(s,FpX_Fp_mul(FpXQ_pow(f,stoi(deg[j]),Tp,p),stoi(sym[j]),p),p);
! 1399: }
! 1400: npows = brent_kung_optpow(lgef(Tp)-4,g-1);
! 1401: pows = FpXQ_powers(sigma,npows,Tp,p);
! 1402: for(i=2; i<=g;i++)
! 1403: {
! 1404: f=FpX_FpXQV_compo(f,pows,Tp,p);
! 1405: for(j=1;j<lg(sym);j++)
! 1406: if(sym[j])
! 1407: {
! 1408: s=FpX_add(s,FpX_Fp_mul(FpXQ_pow(f,stoi(deg[j]),Tp,p),stoi(sym[j]),p),p);
! 1409: }
! 1410: }
! 1411: return gerepileupto(ltop, s);
! 1412: }
! 1413:
! 1414: GEN
! 1415: fixedfieldfactmod(GEN Sp, GEN p, GEN Tmod)
! 1416: {
! 1417: long i;
! 1418: long l=lg(Tmod);
! 1419: GEN F=cgetg(l,t_VEC);
! 1420: for(i=1;i<l;i++)
! 1421: F[i]=(long)FpXQ_minpoly(Sp, (GEN) Tmod[i],p);
! 1422: return F;
! 1423: }
! 1424:
! 1425: GEN
! 1426: fixedfieldnewtonsumaut(GEN sigma, GEN p, GEN Tp, GEN e, long g)
! 1427: {
! 1428: ulong ltop=avma;
! 1429: long i;
! 1430: GEN s,f,V;
! 1431: long rt;
! 1432: sigma=lift(gmul(sigma,gmodulsg(1,p)));
! 1433: f=polx[varn(sigma)];
! 1434: rt=brent_kung_optpow(lgef(Tp)-4,g-1);
! 1435: V=FpXQ_powers(sigma,rt,Tp,p);
! 1436: s=FpXQ_pow(f,e,Tp,p);
! 1437: for(i=2; i<=g;i++)
! 1438: {
! 1439: f=FpX_FpXQV_compo(f,V,Tp,p);
! 1440: s=FpX_add(s,FpXQ_pow(f,e,Tp,p),p);
! 1441: }
! 1442: return gerepileupto(ltop, s);
! 1443: }
! 1444:
! 1445: GEN
! 1446: fixedfieldnewtonsum(GEN O, GEN L, GEN mod, GEN e)
! 1447: {
! 1448: long f,g;
! 1449: long i,j;
! 1450: GEN PL;
! 1451: f=lg(O)-1;
! 1452: g=lg(O[1])-1;
! 1453: PL=cgetg(lg(O), t_COL);
! 1454: for(i=1; i<=f; i++)
! 1455: {
! 1456: ulong ltop=avma;
! 1457: GEN s=gzero;
! 1458: for(j=1; j<=g; j++)
! 1459: s=addii(s,powmodulo((GEN)L[mael(O,i,j)],e,mod));
! 1460: PL[i]=lpileupto(ltop,modii(s,mod));
! 1461: }
! 1462: return PL;
! 1463: }
! 1464:
! 1465: GEN
! 1466: fixedfieldpol(GEN O, GEN L, GEN mod, GEN sym, GEN deg)
! 1467: {
! 1468: ulong ltop=avma;
! 1469: long i;
! 1470: GEN S=gzero;
! 1471: for(i=1;i<lg(sym);i++)
! 1472: if (sym[i])
! 1473: S=gadd(S,gmulsg(sym[i],fixedfieldnewtonsum(O, L, mod, stoi(deg[i]))));
! 1474: S=FpV_red(S,mod);
! 1475: return gerepileupto(ltop, S);
! 1476: }
! 1477:
! 1478: static long
! 1479: fixedfieldtests(GEN LN, long n)
! 1480: {
! 1481: long i,j,k;
! 1482: for (i=1;i<lg(LN[1]);i++)
! 1483: for(j=i+1;j<lg(LN[1]);j++)
! 1484: {
! 1485: for(k=1;k<=n;k++)
! 1486: if (cmpii(gmael(LN,k,j),gmael(LN,k,i)))
! 1487: break;
! 1488: if (k>n)
! 1489: return 0;
! 1490: }
! 1491: return 1;
! 1492: }
! 1493:
! 1494: static long
! 1495: fixedfieldtest(GEN V)
! 1496: {
! 1497: long i,j;
! 1498: for (i=1;i<lg(V);i++)
! 1499: for(j=i+1;j<lg(V);j++)
! 1500: if (!cmpii((GEN)V[i],(GEN)V[j]))
! 1501: return 0;
! 1502: return 1;
! 1503: }
! 1504:
! 1505: void
! 1506: debug_surmer(char *s,GEN S, long n)
! 1507: {
! 1508: long l=lg(S);
! 1509: setlg(S,n+1);
! 1510: fprintferr(s,S);
! 1511: setlg(S,l);
! 1512: }
! 1513:
! 1514: /*We try hard to find a polynomial R squarefree mod p.
! 1515: Unfortunately, it may be too much asked.
! 1516: A less bugprone way would be to only ask it to be squarefree mod p^e
! 1517: with e not too big. Most of the code is here, but I lack the theoretical
! 1518: knowledge to make it to work smoothly.B.A.
! 1519: */
! 1520: static GEN
! 1521: fixedfieldsurmer(GEN O, GEN L, GEN mod, GEN l, GEN p, GEN S, GEN deg, long v, GEN LN,long n)
! 1522: {
! 1523: long i;
! 1524: GEN V;
! 1525: V=(GEN)LN[n];
! 1526: for (i=1;i<lg(S);i++)
! 1527: S[i]=0;
! 1528: S[n]=1;
! 1529: for (i=0;i<2*n-1;i++)
! 1530: {
! 1531: long k;
! 1532: if (fixedfieldtest(V))
! 1533: {
! 1534: ulong av=avma;
! 1535: GEN s=fixedfieldpol(O,L,mod,S,deg);
! 1536: GEN P=FpV_roots_to_pol(s,mod,v);
! 1537: P=FpX_center(P,mod);
! 1538: if (p==gun || FpX_is_squarefree(P,p))
! 1539: {
! 1540: GEN V;
! 1541: if (DEBUGLEVEL>=4)
! 1542: debug_surmer("FixedField: Sym: %Z\n",S,n);
! 1543: V=cgetg(3,t_VEC);
! 1544: V[1]=lcopy(s);/*do not swap*/
! 1545: V[2]=lcopy(P);
! 1546: return V;
! 1547: }
! 1548: else
! 1549: {
! 1550: if (DEBUGLEVEL>=6)
! 1551: debug_surmer("FixedField: bad mod: %Z\n",S,n);
! 1552: avma=av;
! 1553: }
! 1554: }
! 1555: else
! 1556: {
! 1557: if (DEBUGLEVEL>=6)
! 1558: debug_surmer("FixedField: Tested: %Z\n",S,n);
! 1559: }
! 1560: k=1+(i%n);
! 1561: S[k]++;
! 1562: V=FpV_red(gadd(V,(GEN)LN[k]),l);
! 1563: }
! 1564: return NULL;
! 1565: }
! 1566:
! 1567: GEN
! 1568: fixedfieldsympol(GEN O, GEN L, GEN mod, GEN l, GEN p, GEN S, GEN deg, long v)
! 1569: {
! 1570: ulong ltop=avma;
! 1571: GEN V=NULL;
! 1572: GEN LN=cgetg(lg(L),t_VEC);
! 1573: GEN Ll=FpV_red(L,l);
! 1574: long n,i;
! 1575: for(i=1;i<lg(deg);i++)
! 1576: deg[i]=0;
! 1577: for(n=1,i=1;i<lg(L);i++)
! 1578: {
! 1579: long j;
! 1580: LN[n]=(long)fixedfieldnewtonsum(O,Ll,l,stoi(i));
! 1581: for(j=2;j<lg(LN[n]);j++)
! 1582: if(cmpii(gmael(LN,n,j),gmael(LN,n,1)))
! 1583: break;
! 1584: if(j==lg(LN[n]) && j>2)
! 1585: continue;
! 1586: if (DEBUGLEVEL>=6) fprintferr("FixedField: LN[%d]=%Z\n",n,LN[n]);
! 1587: deg[n]=i;
! 1588: if (fixedfieldtests(LN,n))
! 1589: {
! 1590: V=fixedfieldsurmer(O,L,mod,l,p,S,deg,v,LN,n);
! 1591: if (V)
! 1592: break;
! 1593: }
! 1594: n++;
! 1595: }
! 1596: if (DEBUGLEVEL>=4)
! 1597: {
! 1598: long l=lg(deg);
! 1599: setlg(deg,n);
! 1600: fprintferr("FixedField: Computed degrees: %Z\n",deg);
! 1601: setlg(deg,l);
! 1602: }
! 1603: if (!V)
! 1604: err(talker, "prime too small in fixedfield");
! 1605: return gerepileupto(ltop,V);
! 1606: }
! 1607: /*
! 1608: * Calcule l'inclusion de R dans T i.e. S telque T|R\compo S
! 1609: * Ne recopie pas PL.
! 1610: */
! 1611: GEN
! 1612: fixedfieldinclusion(GEN O, GEN PL)
! 1613: {
! 1614: GEN S;
! 1615: int i, j;
! 1616: S = cgetg((lg(O) - 1) * (lg(O[1]) - 1) + 1, t_COL);
! 1617: for (i = 1; i < lg(O); i++)
! 1618: for (j = 1; j < lg(O[i]); j++)
! 1619: S[((GEN *) O)[i][j]] = PL[i];
! 1620: return S;
! 1621: }
! 1622:
! 1623: /*Usually mod is bigger than than den so there is no need to reduce it.*/
! 1624: GEN
! 1625: vandermondeinversemod(GEN L, GEN T, GEN den, GEN mod)
! 1626: {
! 1627: ulong av;
! 1628: int i, j, n = lg(L);
! 1629: long x = varn(T);
! 1630: GEN M, P, Tp;
! 1631: M = cgetg(n, t_MAT);
! 1632: av=avma;
! 1633: Tp = gclone(FpX_red(deriv(T, x),mod)); /*clone*/
! 1634: avma=av;
! 1635: for (i = 1; i < n; i++)
! 1636: {
! 1637: GEN z;
! 1638: av = avma;
! 1639: z = mpinvmod(FpX_eval(Tp, (GEN) L[i],mod),mod);
! 1640: z = modii(mulii(den,z),mod);
! 1641: P = FpX_Fp_mul(FpX_div(T, deg1pol(gun,negi((GEN) L[i]),x),mod), z, mod);
! 1642: M[i] = lgetg(n, t_COL);
! 1643: for (j = 1; j < n; j++)
! 1644: mael(M,i,j) = lcopy((GEN) P[1 + j]);
! 1645: M[i] = lpileupto(av,(GEN)M[i]);
! 1646: }
! 1647: gunclone(Tp); /*unclone*/
! 1648: return M;
! 1649: }
! 1650: /* Calcule le polynome associe a un vecteur conjugue v*/
! 1651: static GEN
! 1652: vectopol(GEN v, GEN M, GEN den , GEN mod, long x)
! 1653: {
! 1654: long n = lg(v),i,k,av;
! 1655: GEN z = cgetg(n+1,t_POL),p1,p2,mod2;
! 1656: av=avma;
! 1657: mod2=gclone(shifti(mod,-1));/*clone*/
! 1658: avma=av;
! 1659: z[1]=evalsigne(1)|evalvarn(x)|evallgef(n+1);
! 1660: for (i=2; i<=n; i++)
! 1661: {
! 1662: p1=gzero; av=avma;
! 1663: for (k=1; k<n; k++)
! 1664: {
! 1665: p2=mulii(gcoeff(M,i-1,k),(GEN)v[k]);
! 1666: p1=addii(p1,p2);
! 1667: }
! 1668: p1=modii(p1,mod);
! 1669: if (cmpii(p1,mod2)>0) p1=subii(p1,mod);
! 1670: p1=gdiv(p1,den);
! 1671: z[i]=lpileupto(av,p1);
! 1672: }
! 1673: gunclone(mod2);/*unclone*/
! 1674: return normalizepol_i(z,n+1);
! 1675: }
! 1676: /* Calcule le polynome associe a une permutation des racines*/
! 1677: static GEN
! 1678: permtopol(GEN p, GEN L, GEN M, GEN den, GEN mod, long x)
! 1679: {
! 1680: long n = lg(L),i,k,av;
! 1681: GEN z = cgetg(n+1,t_POL),p1,p2,mod2;
! 1682: if (lg(p) != n) err(talker,"incorrect permutation in permtopol");
! 1683: av=avma;
! 1684: mod2=gclone(shifti(mod,-1)); /*clone*/
! 1685: avma=av;
! 1686: z[1]=evalsigne(1)|evalvarn(x)|evallgef(n+1);
! 1687: for (i=2; i<=n; i++)
! 1688: {
! 1689: p1=gzero; av=avma;
! 1690: for (k=1; k<n; k++)
! 1691: {
! 1692: p2=mulii(gcoeff(M,i-1,k),(GEN)L[p[k]]);
! 1693: p1=addii(p1,p2);
! 1694: }
! 1695: p1=modii(p1,mod);
! 1696: if (cmpii(p1,mod2)>0) p1=subii(p1,mod);
! 1697: p1=gdiv(p1,den);
! 1698: z[i]=lpileupto(av,p1);
! 1699: }
! 1700: gunclone(mod2); /*unclone*/
! 1701: return normalizepol_i(z,n+1);
! 1702: }
! 1703:
! 1704: GEN
! 1705: galoisgrouptopol( GEN res, GEN L, GEN M, GEN den, GEN mod, long v)
! 1706: {
! 1707: GEN aut = cgetg(lg(res), t_COL);
! 1708: long i;
! 1709: for (i = 1; i < lg(res); i++)
! 1710: {
! 1711: if (DEBUGLEVEL>=6)
! 1712: fprintferr("%d ",i);
! 1713: aut[i] = (long) permtopol((GEN) res[i], L, M, den, mod, v);
! 1714: }
! 1715: return aut;
! 1716: }
! 1717: /* No more used and ugly.
! 1718: * However adding corepartial to GP seem a good idea.
! 1719: *
! 1720: * factorise partiellement n sous la forme n=d*u*f^2 ou d est un
! 1721: * discriminant fondamental et u n'est pas un carre parfait et
! 1722: * retourne u*f
! 1723: * Note: Parfois d n'est pas un discriminant (congru a 3 mod 4).
! 1724: * Cela se produit si u est congrue a 3 mod 4.
! 1725: */
! 1726: GEN
! 1727: corediscpartial(GEN n)
! 1728: {
! 1729: ulong av = avma;
! 1730: long i,r,s;
! 1731: GEN fa, p1, p2, p3, res = gun, res2 = gun, res3=gun;
! 1732: /*d=res,f=res2,u=res3*/
! 1733: if (gcmp1(n))
! 1734: return gun;
! 1735: fa = auxdecomp(n, 0);
! 1736: p1 = (GEN) fa[1];
! 1737: p2 = (GEN) fa[2];
! 1738: for (i = 1; i < lg(p1) - 1; i++)
! 1739: {
! 1740: p3 = (GEN) p2[i];
! 1741: if (mod2(p3))
! 1742: res = mulii(res, (GEN) p1[i]);
! 1743: if (!gcmp1(p3))
! 1744: res2 = mulii(res2, gpow((GEN) p1[i], shifti(p3, -1), 0));
! 1745: }
! 1746: p3 = (GEN) p2[i];
! 1747: if (mod2(p3)) /* impaire: verifions */
! 1748: {
! 1749: if (!gcmp1(p3))
! 1750: res2 = mulii(res2, gpow((GEN) p1[i], shifti(p3, -1), 0));
! 1751: if (isprime((GEN) p1[i]))
! 1752: res = mulii(res, (GEN) p1[i]);
! 1753: else
! 1754: res3 = (GEN)p1[i];
! 1755: }
! 1756: else /* paire:OK */
! 1757: res2 = mulii(res2, gpow((GEN) p1[i], shifti(p3, -1), 0));
! 1758: r = mod4(res);
! 1759: if (signe(res) < 0)
! 1760: r = 4 - r;
! 1761: s = mod4(res3);/*res2 est >0*/
! 1762: if (r == 3 && s!=3)
! 1763: /*d est n'est pas un discriminant mais res2 l'est: corrige*/
! 1764: res2 = gmul2n((GEN) res2, -1);
! 1765: return gerepileupto(av,gmul(res2,res3));
! 1766: }
! 1767: GEN respm(GEN x,GEN y,GEN pm);
! 1768: GEN ZX_disc(GEN x);
! 1769:
! 1770: GEN
! 1771: indexpartial(GEN P, GEN DP)
! 1772: {
! 1773: ulong av = avma;
! 1774: long i, nb;
! 1775: GEN fa, p1, res = gun, dP;
! 1776: dP = derivpol(P);
! 1777: if(DEBUGLEVEL>=5) gentimer(3);
! 1778: if(!DP) DP = ZX_disc(P);
! 1779: DP = mpabs(DP);
! 1780: if(DEBUGLEVEL>=5) genmsgtimer(3,"IndexPartial: discriminant");
! 1781: fa = auxdecomp(DP, 0);
! 1782: if(DEBUGLEVEL>=5) genmsgtimer(3,"IndexPartial: factorization");
! 1783: nb = lg(fa[1]);
! 1784: for (i = 1; i < nb; i++)
! 1785: {
! 1786: GEN p=gmael(fa,1,i);
! 1787: GEN e=gmael(fa,2,i);
! 1788: p1 = powgi(p,shifti(e,-1));
! 1789: if ( i==nb-1 )
! 1790: {
! 1791: if ( mod2(e) && !isprime(p) )
! 1792: p1 = mulii(p1,p);
! 1793: }
! 1794: else if ( cmpis(e,4)>=0 )
! 1795: {
! 1796: if(DEBUGLEVEL>=5) fprintferr("IndexPartial: factor %Z ",p1);
! 1797: p1 = mppgcd(p1, respm(P,dP,p1));
! 1798: if(DEBUGLEVEL>=5)
! 1799: {
! 1800: fprintferr("--> %Z : ",p1);
! 1801: genmsgtimer(3,"");
! 1802: }
! 1803: }
! 1804: res=mulii(res,p1);
! 1805: }
! 1806: return gerepileupto(av,res);
! 1807: }
! 1808:
! 1809: /* Calcule la puissance exp d'une permutation decompose en cycle */
! 1810: GEN
! 1811: permcyclepow(GEN O, long exp)
! 1812: {
! 1813: int j, k, n;
! 1814: GEN p;
! 1815: for (n = 0, j = 1; j < lg(O); j++)
! 1816: n += lg(O[j]) - 1;
! 1817: p = cgetg(n + 1, t_VECSMALL);
! 1818: for (j = 1; j < lg(O); j++)
! 1819: {
! 1820: n = lg(O[j]) - 1;
! 1821: for (k = 1; k <= n; k++)
! 1822: p[mael(O,j,k)] = mael(O,j,1 + (k + exp - 1) % n);
! 1823: }
! 1824: return p;
! 1825: }
! 1826:
! 1827: /*
! 1828: * Casse l'orbite O en sous-orbite d'ordre premier correspondant a des
! 1829: * puissance de l'element
! 1830: */
! 1831: GEN
! 1832: splitorbite(GEN O)
! 1833: {
! 1834: ulong ltop = avma, lbot;
! 1835: int i, n;
! 1836: GEN F, fc, res;
! 1837: n = lg(O[1]) - 1;
! 1838: F = factor(stoi(n));
! 1839: fc = cgetg(lg(F[1]), t_VECSMALL);
! 1840: for (i = 1; i < lg(fc); i++)
! 1841: fc[i] = itos(powgi(gmael(F,1,i), gmael(F,2,i)));
! 1842: lbot = avma;
! 1843: res = cgetg(3, t_VEC);
! 1844: res[1] = lgetg(lg(fc), t_VEC);
! 1845: res[2] = lgetg(lg(fc), t_VECSMALL);
! 1846: for (i = 1; i < lg(fc); i++)
! 1847: {
! 1848: ((GEN **) res)[1][lg(fc) - i] = permcyclepow(O, n / fc[i]);
! 1849: ((GEN *) res)[2][lg(fc) - i] = fc[i];
! 1850: }
! 1851: if ( DEBUGLEVEL>=4 )
! 1852: fprintferr("GaloisConj:splitorbite: %Z\n",res);
! 1853: return gerepile(ltop, lbot, res);
! 1854: }
! 1855:
! 1856: /*
! 1857: * structure qui contient l'analyse du groupe de Galois par etude des degres
! 1858: * de Frobenius:
! 1859: *
! 1860: * p: nombre premier a relever deg: degre du lift à effectuer pgp: plus grand
! 1861: * diviseur premier du degre de T.
! 1862: * l: un nombre premier tel que T est totalement decompose
! 1863: * modulo l ppp: plus petit diviseur premier du degre de T. primepointer:
! 1864: * permet de calculer les premiers suivant p.
! 1865: */
! 1866: struct galois_analysis
! 1867: {
! 1868: long p;
! 1869: long deg;
! 1870: long ord;
! 1871: long l;
! 1872: long ppp;
! 1873: long p4;
! 1874: enum {ga_all_normal=1,ga_ext_2=2,ga_non_wss=4} group;
! 1875: byteptr primepointer;
! 1876: };
! 1877: void
! 1878: galoisanalysis(GEN T, struct galois_analysis *ga, long calcul_l, long karma_type)
! 1879: {
! 1880: ulong ltop=avma;
! 1881: long n,p;
! 1882: long i;
! 1883: enum {k_amoeba=0,k_snake=1,k_fish=2,k_bird=4,k_rodent=6,k_dog=8,k_human=9,k_cat=12} karma;
! 1884: long group,linf;
! 1885: /*TODO: complete the table to at least 200*/
! 1886: const int prim_nonss_orders[]={36,48,56,60,72,75,80,96,108,120,132,0};
! 1887: GEN F,Fp,Fe,Fpe,O;
! 1888: long np;
! 1889: long order,phi_order;
! 1890: long plift,nbmax,nbtest,deg;
! 1891: byteptr primepointer,pp;
! 1892: if (DEBUGLEVEL >= 1)
! 1893: timer2();
! 1894: n = degpol(T);
! 1895: if (!karma_type) karma_type=n;
! 1896: else err(warner,"entering black magic computation");
! 1897: O = cgetg(n+1,t_VECSMALL);
! 1898: for(i=1;i<=n;i++) O[i]=0;
! 1899: F = factor(stoi(n));
! 1900: Fp=gtovecsmall((GEN)F[1]);
! 1901: Fe=gtovecsmall((GEN)F[2]);
! 1902: np=lg(Fp)-1;
! 1903: Fpe=cgetg(np+1, t_VECSMALL);
! 1904: for (i = 1; i < lg(Fpe); i++)
! 1905: Fpe[i] = itos(powgi(gmael(F,1,i), gmael(F,2,i)));
! 1906: /*In this part, we study the cardinal of the group to have an information
! 1907: about the orders, so if we are unlucky we can continue.*/
! 1908:
! 1909: /*Are there non WSS groups of this order ?*/
! 1910: group=0;
! 1911: for(i=0;prim_nonss_orders[i];i++)
! 1912: if (n%prim_nonss_orders[i] == 0)
! 1913: {
! 1914: group |= ga_non_wss;
! 1915: break;
! 1916: }
! 1917: if ( n>12 && n%12 == 0 )
! 1918: {
! 1919: /*We need to know the greatest prime dividing n/12*/
! 1920: if ( Fp[np] == 3 && Fe[np] == 1 )
! 1921: group |= ga_ext_2;
! 1922: }
! 1923: phi_order = 1;
! 1924: order = 1;
! 1925: for (i = np; i > 0; i--)
! 1926: {
! 1927: p = Fp[i];
! 1928: if (phi_order % p != 0)
! 1929: {
! 1930: order *= p;
! 1931: phi_order *= p - 1;
! 1932: }
! 1933: else
! 1934: {
! 1935: group |= ga_all_normal;
! 1936: break;
! 1937: }
! 1938: if (Fe[i]>1)
! 1939: break;
! 1940: }
! 1941: /*Now, we study the orders of the Frobenius elements*/
! 1942: plift = 0;
! 1943: nbmax = 8+(n>>1);
! 1944: nbtest = 0; karma = k_amoeba;
! 1945: deg = Fp[np];
! 1946: for (p = 0, pp = primepointer = diffptr;
! 1947: (plift == 0
! 1948: || (nbtest < nbmax && (nbtest <=8 || order < (n>>1)))
! 1949: || (n == 24 && O[6] == 0 && O[4] == 0)
! 1950: || ((group&ga_non_wss) && order == Fp[np]))
! 1951: && (nbtest < 3 * nbmax || !(group&ga_non_wss)) ;)
! 1952: {
! 1953: ulong av;
! 1954: long prime_incr;
! 1955: GEN ip,FS,p1;
! 1956: long o,norm_o=1;
! 1957: prime_incr = *primepointer++;
! 1958: if (!prime_incr)
! 1959: err(primer1);
! 1960: p += prime_incr;
! 1961: /*discard small primes*/
! 1962: if (p <= (n << 1))
! 1963: continue;
! 1964: ip=stoi(p);
! 1965: if (!FpX_is_squarefree(T,ip))
! 1966: continue;
! 1967: nbtest++;
! 1968: av=avma;
! 1969: FS=(GEN)simplefactmod(T,ip)[1];
! 1970: p1=(GEN)FS[1];
! 1971: for(i=2;i<lg(FS);i++)
! 1972: if (cmpii(p1,(GEN)FS[i]))
! 1973: break;
! 1974: if (i<lg(FS))
! 1975: {
! 1976: avma = ltop;
! 1977: if (DEBUGLEVEL >= 2)
! 1978: fprintferr("GaloisAnalysis:non Galois for p=%ld\n", p);
! 1979: ga->p = p;
! 1980: ga->deg = 0;
! 1981: return; /* Not a Galois polynomial */
! 1982: }
! 1983: o=n/(lg(FS)-1);
! 1984: avma=av;
! 1985: if (!O[o])
! 1986: O[o]=p;
! 1987: if (o % deg == 0)
! 1988: {
! 1989: /*We try to find a power of the Frobenius which generate
! 1990: a normal subgroup just by looking at the order.*/
! 1991: if (o * Fp[1] >= n)
! 1992: /*Subgroup of smallest index are normal*/
! 1993: norm_o = o;
! 1994: else
! 1995: {
! 1996: norm_o = 1;
! 1997: for (i = np; i > 0; i--)
! 1998: {
! 1999: if (o % Fpe[i] == 0)
! 2000: norm_o *= Fpe[i];
! 2001: else
! 2002: break;
! 2003: }
! 2004: }
! 2005: if (norm_o != 1)
! 2006: {
! 2007: if (!(group&ga_all_normal) || o > order ||
! 2008: (o == order && (plift == 0 || norm_o > deg
! 2009: || (norm_o == deg && cgcd(p-1,karma_type)>karma ))))
! 2010: {
! 2011: deg = norm_o;
! 2012: order = o;
! 2013: plift = p;
! 2014: karma=cgcd(p-1,karma_type);
! 2015: pp = primepointer;
! 2016: group |= ga_all_normal;
! 2017: }
! 2018: }
! 2019: else if (!(group&ga_all_normal) && (plift == 0 || o > order
! 2020: || ( o == order && cgcd(p-1,karma_type)>karma )))
! 2021: {
! 2022: order = o;
! 2023: plift = p;
! 2024: karma=cgcd(p-1,karma_type);
! 2025: pp = primepointer;
! 2026: }
! 2027: }
! 2028: if (DEBUGLEVEL >= 5)
! 2029: fprintferr("GaloisAnalysis:Nbtest=%ld,p=%ld,o=%ld,n_o=%d,best p=%ld,ord=%ld,k=%ld\n",
! 2030: nbtest, p, o, norm_o, plift, order,karma);
! 2031: }
! 2032: /* This is to avoid looping on non-wss group.
! 2033: To be checked for large groups. */
! 2034: /*Would it be better to disable this check if we are in a good case ?
! 2035: * (ga_all_normal and !(ga_ext_2) (e.g. 60)) ?*/
! 2036: if (plift == 0 || ((group&ga_non_wss) && order == Fp[np]))
! 2037: {
! 2038: deg = 0;
! 2039: err(warner, "Galois group almost certainly not weakly super solvable");
! 2040: }
! 2041: /*linf=(n*(n-1))>>2;*/
! 2042: linf=n;
! 2043: if (calcul_l && O[1]<=linf)
! 2044: {
! 2045: ulong av;
! 2046: long prime_incr;
! 2047: long l=0;
! 2048: /*we need a totally splited prime l*/
! 2049: av = avma;
! 2050: while (l == 0)
! 2051: {
! 2052: long nb;
! 2053: prime_incr = *primepointer++;
! 2054: if (!prime_incr)
! 2055: err(primer1);
! 2056: p += prime_incr;
! 2057: if (p<=linf) continue;
! 2058: nb=FpX_nbroots(T,stoi(p));
! 2059: if (nb == n)
! 2060: l = p;
! 2061: else if (nb && FpX_is_squarefree(T,stoi(p)))
! 2062: {
! 2063: avma = ltop;
! 2064: if (DEBUGLEVEL >= 2)
! 2065: fprintferr("GaloisAnalysis:non Galois for p=%ld\n", p);
! 2066: ga->p = p;
! 2067: ga->deg = 0;
! 2068: return; /* Not a Galois polynomial */
! 2069: }
! 2070: avma = av;
! 2071: }
! 2072: O[1]=l;
! 2073: }
! 2074: ga->p = plift;
! 2075: ga->group = group;
! 2076: ga->deg = deg;
! 2077: ga->ord = order;
! 2078: ga->l = O[1];
! 2079: ga->primepointer = pp;
! 2080: ga->ppp = Fp[1];
! 2081: ga->p4 = O[4];
! 2082: if (DEBUGLEVEL >= 4)
! 2083: fprintferr("GaloisAnalysis:p=%ld l=%ld group=%ld deg=%ld ord=%ld\n",
! 2084: plift, O[1], group, deg, order);
! 2085: if (DEBUGLEVEL >= 1)
! 2086: msgtimer("galoisanalysis()");
! 2087: avma = ltop;
! 2088: }
! 2089:
! 2090: /* Groupe A4 */
! 2091: GEN
! 2092: a4galoisgen(GEN T, struct galois_test *td)
! 2093: {
! 2094: ulong ltop = avma, av, av2;
! 2095: long i, j, k;
! 2096: long n;
! 2097: long N, hop = 0;
! 2098: long **O;
! 2099: GEN *ar, **mt; /* tired of casting */
! 2100: GEN t, u;
! 2101: GEN res, orb, ry;
! 2102: GEN pft, pfu, pfv;
! 2103: n = degpol(T);
! 2104: res = cgetg(3, t_VEC);
! 2105: ry = cgetg(4, t_VEC);
! 2106: res[1] = (long) ry;
! 2107: pft = cgetg(n + 1, t_VECSMALL);
! 2108: pfu = cgetg(n + 1, t_VECSMALL);
! 2109: pfv = cgetg(n + 1, t_VECSMALL);
! 2110: ry[1] = (long) pft;
! 2111: ry[2] = (long) pfu;
! 2112: ry[3] = (long) pfv;
! 2113: ry = cgetg(4, t_VECSMALL);
! 2114: ry[1] = 2;
! 2115: ry[2] = 2;
! 2116: ry[3] = 3;
! 2117: res[2] = (long) ry;
! 2118: av = avma;
! 2119: ar = (GEN *) alloue_ardoise(n, 1 + lg(td->ladic));
! 2120: mt = (GEN **) td->PV[td->ordre[n]];
! 2121: t = cgetg(n + 1, t_VECSMALL) + 1; /* Sorry for this hack */
! 2122: u = cgetg(n + 1, t_VECSMALL) + 1; /* too lazy to correct */
! 2123: av2 = avma;
! 2124: N = itos(gdiv(mpfact(n), mpfact(n >> 1))) >> (n >> 1);
! 2125: if (DEBUGLEVEL >= 4)
! 2126: fprintferr("A4GaloisConj:I will test %ld permutations\n", N);
! 2127: avma = av2;
! 2128: for (i = 0; i < n; i++)
! 2129: t[i] = i + 1;
! 2130: for (i = 0; i < N; i++)
! 2131: {
! 2132: GEN g;
! 2133: int a, x, y;
! 2134: if (i == 0)
! 2135: {
! 2136: gaffect(gzero, ar[(n - 2) >> 1]);
! 2137: for (k = n - 2; k > 2; k -= 2)
! 2138: addiiz(ar[k >> 1], addii(mt[k + 1][k + 2], mt[k + 2][k + 1]),
! 2139: ar[(k >> 1) - 1]);
! 2140: }
! 2141: else
! 2142: {
! 2143: x = i;
! 2144: y = 1;
! 2145: do
! 2146: {
! 2147: y += 2;
! 2148: a = x%y;
! 2149: x = x/y;
! 2150: }
! 2151: while (!a);
! 2152: switch (y)
! 2153: {
! 2154: case 3:
! 2155: x = t[2];
! 2156: if (a == 1)
! 2157: {
! 2158: t[2] = t[1];
! 2159: t[1] = x;
! 2160: }
! 2161: else
! 2162: {
! 2163: t[2] = t[0];
! 2164: t[0] = x;
! 2165: }
! 2166: break;
! 2167: case 5:
! 2168: x = t[0];
! 2169: t[0] = t[2];
! 2170: t[2] = t[1];
! 2171: t[1] = x;
! 2172: x = t[4];
! 2173: t[4] = t[4 - a];
! 2174: t[4 - a] = x;
! 2175: addiiz(ar[2], addii(mt[t[4]][t[5]], mt[t[5]][t[4]]), ar[1]);
! 2176: break;
! 2177: case 7:
! 2178: x = t[0];
! 2179: t[0] = t[4];
! 2180: t[4] = t[3];
! 2181: t[3] = t[1];
! 2182: t[1] = t[2];
! 2183: t[2] = x;
! 2184: x = t[6];
! 2185: t[6] = t[6 - a];
! 2186: t[6 - a] = x;
! 2187: addiiz(ar[3], addii(mt[t[6]][t[7]], mt[t[7]][t[6]]), ar[2]);
! 2188: addiiz(ar[2], addii(mt[t[4]][t[5]], mt[t[5]][t[4]]), ar[1]);
! 2189: break;
! 2190: case 9:
! 2191: x = t[0];
! 2192: t[0] = t[6];
! 2193: t[6] = t[5];
! 2194: t[5] = t[3];
! 2195: t[3] = x;
! 2196: x = t[4];
! 2197: t[4] = t[1];
! 2198: t[1] = x;
! 2199: x = t[8];
! 2200: t[8] = t[8 - a];
! 2201: t[8 - a] = x;
! 2202: addiiz(ar[4], addii(mt[t[8]][t[9]], mt[t[9]][t[8]]), ar[3]);
! 2203: addiiz(ar[3], addii(mt[t[6]][t[7]], mt[t[7]][t[6]]), ar[2]);
! 2204: addiiz(ar[2], addii(mt[t[4]][t[5]], mt[t[5]][t[4]]), ar[1]);
! 2205: break;
! 2206: default:
! 2207: y--;
! 2208: x = t[0];
! 2209: t[0] = t[2];
! 2210: t[2] = t[1];
! 2211: t[1] = x;
! 2212: for (k = 4; k < y; k += 2)
! 2213: {
! 2214: int j;
! 2215: x = t[k];
! 2216: for (j = k; j > 0; j--)
! 2217: t[j] = t[j - 1];
! 2218: t[0] = x;
! 2219: }
! 2220: x = t[y];
! 2221: t[y] = t[y - a];
! 2222: t[y - a] = x;
! 2223: for (k = y; k > 2; k -= 2)
! 2224: addiiz(ar[k >> 1],
! 2225: addii(mt[t[k]][t[k + 1]], mt[t[k + 1]][t[k]]),
! 2226: ar[(k >> 1) - 1]);
! 2227: }
! 2228: }
! 2229: g = addii(ar[1], addii(addii(mt[t[0]][t[1]], mt[t[1]][t[0]]),
! 2230: addii(mt[t[2]][t[3]], mt[t[3]][t[2]])));
! 2231: if (padicisint(g, td))
! 2232: {
! 2233: for (k = 0; k < n; k += 2)
! 2234: {
! 2235: pft[t[k]] = t[k + 1];
! 2236: pft[t[k + 1]] = t[k];
! 2237: }
! 2238: if (verifietest(pft, td))
! 2239: break;
! 2240: else
! 2241: hop++;
! 2242: }
! 2243: avma = av2;
! 2244: }
! 2245: if (i == N)
! 2246: {
! 2247: avma = ltop;
! 2248: if (DEBUGLEVEL >= 1 && hop)
! 2249: fprintferr("A4GaloisConj: %ld hop sur %ld iterations\n", hop, N);
! 2250: return gzero;
! 2251: }
! 2252: if (DEBUGLEVEL >= 1 && hop)
! 2253: fprintferr("A4GaloisConj: %ld hop sur %ld iterations\n", hop, N);
! 2254: N = itos(gdiv(mpfact(n >> 1), mpfact(n >> 2))) >> 1;
! 2255: avma = av2;
! 2256: if (DEBUGLEVEL >= 4)
! 2257: fprintferr("A4GaloisConj:sigma=%Z \n", pft);
! 2258: for (i = 0; i < N; i++)
! 2259: {
! 2260: GEN g;
! 2261: int a, x, y;
! 2262: if (i == 0)
! 2263: {
! 2264: for (k = 0; k < n; k += 4)
! 2265: {
! 2266: u[k + 3] = t[k + 3];
! 2267: u[k + 2] = t[k + 1];
! 2268: u[k + 1] = t[k + 2];
! 2269: u[k] = t[k];
! 2270: }
! 2271: }
! 2272: else
! 2273: {
! 2274: x = i;
! 2275: y = -2;
! 2276: do
! 2277: {
! 2278: y += 4;
! 2279: a = x%y;
! 2280: x = x/y;
! 2281: }
! 2282: while (!a);
! 2283: x = u[2];
! 2284: u[2] = u[0];
! 2285: u[0] = x;
! 2286: switch (y)
! 2287: {
! 2288: case 2:
! 2289: break;
! 2290: case 6:
! 2291: x = u[4];
! 2292: u[4] = u[6];
! 2293: u[6] = x;
! 2294: if (!(a & 1))
! 2295: {
! 2296: a = 4 - (a >> 1);
! 2297: x = u[6];
! 2298: u[6] = u[a];
! 2299: u[a] = x;
! 2300: x = u[4];
! 2301: u[4] = u[a - 2];
! 2302: u[a - 2] = x;
! 2303: }
! 2304: break;
! 2305: case 10:
! 2306: x = u[6];
! 2307: u[6] = u[3];
! 2308: u[3] = u[2];
! 2309: u[2] = u[4];
! 2310: u[4] = u[1];
! 2311: u[1] = u[0];
! 2312: u[0] = x;
! 2313: if (a >= 3)
! 2314: a += 2;
! 2315: a = 8 - a;
! 2316: x = u[10];
! 2317: u[10] = u[a];
! 2318: u[a] = x;
! 2319: x = u[8];
! 2320: u[8] = u[a - 2];
! 2321: u[a - 2] = x;
! 2322: break;
! 2323: }
! 2324: }
! 2325: g = gzero;
! 2326: for (k = 0; k < n; k += 2)
! 2327: g = addii(g, addii(mt[u[k]][u[k + 1]], mt[u[k + 1]][u[k]]));
! 2328: if (padicisint(g, td))
! 2329: {
! 2330: for (k = 0; k < n; k += 2)
! 2331: {
! 2332: pfu[u[k]] = u[k + 1];
! 2333: pfu[u[k + 1]] = u[k];
! 2334: }
! 2335: if (verifietest(pfu, td))
! 2336: break;
! 2337: else
! 2338: hop++;
! 2339: }
! 2340: avma = av2;
! 2341: }
! 2342: if (i == N)
! 2343: {
! 2344: avma = ltop;
! 2345: return gzero;
! 2346: }
! 2347: if (DEBUGLEVEL >= 1 && hop)
! 2348: fprintferr("A4GaloisConj: %ld hop sur %ld iterations\n", hop, N);
! 2349: if (DEBUGLEVEL >= 4)
! 2350: fprintferr("A4GaloisConj:tau=%Z \n", pfu);
! 2351: avma = av2;
! 2352: orb = cgetg(3, t_VEC);
! 2353: orb[1] = (long) pft;
! 2354: orb[2] = (long) pfu;
! 2355: if (DEBUGLEVEL >= 4)
! 2356: fprintferr("A4GaloisConj:orb=%Z \n", orb);
! 2357: O = (long**)permorbite(orb);
! 2358: if (DEBUGLEVEL >= 4)
! 2359: fprintferr("A4GaloisConj:O=%Z \n", O);
! 2360: av2 = avma;
! 2361: for (j = 0; j < 2; j++)
! 2362: {
! 2363: pfv[O[1][1]] = O[2][1];
! 2364: pfv[O[1][2]] = O[2][3 + j];
! 2365: pfv[O[1][3]] = O[2][4 - (j << 1)];
! 2366: pfv[O[1][4]] = O[2][2 + j];
! 2367: for (i = 0; i < 4; i++)
! 2368: {
! 2369: long x;
! 2370: GEN g;
! 2371: switch (i)
! 2372: {
! 2373: case 0:
! 2374: break;
! 2375: case 1:
! 2376: x = O[3][1];
! 2377: O[3][1] = O[3][2];
! 2378: O[3][2] = x;
! 2379: x = O[3][3];
! 2380: O[3][3] = O[3][4];
! 2381: O[3][4] = x;
! 2382: break;
! 2383: case 2:
! 2384: x = O[3][1];
! 2385: O[3][1] = O[3][4];
! 2386: O[3][4] = x;
! 2387: x = O[3][2];
! 2388: O[3][2] = O[3][3];
! 2389: O[3][3] = x;
! 2390: break;
! 2391: case 3:
! 2392: x = O[3][1];
! 2393: O[3][1] = O[3][2];
! 2394: O[3][2] = x;
! 2395: x = O[3][3];
! 2396: O[3][3] = O[3][4];
! 2397: O[3][4] = x;
! 2398: }
! 2399: pfv[O[2][1]] = O[3][1];
! 2400: pfv[O[2][3 + j]] = O[3][4 - j];
! 2401: pfv[O[2][4 - (j << 1)]] = O[3][2 + (j << 1)];
! 2402: pfv[O[2][2 + j]] = O[3][3 - j];
! 2403: pfv[O[3][1]] = O[1][1];
! 2404: pfv[O[3][4 - j]] = O[1][2];
! 2405: pfv[O[3][2 + (j << 1)]] = O[1][3];
! 2406: pfv[O[3][3 - j]] = O[1][4];
! 2407: g = gzero;
! 2408: for (k = 1; k <= n; k++)
! 2409: g = addii(g, mt[k][pfv[k]]);
! 2410: if (padicisint(g, td) && verifietest(pfv, td))
! 2411: {
! 2412: avma = av;
! 2413: if (DEBUGLEVEL >= 1)
! 2414: fprintferr("A4GaloisConj:%ld hop sur %d iterations max\n",
! 2415: hop, 10395 + 68);
! 2416: return res;
! 2417: }
! 2418: else
! 2419: hop++;
! 2420: avma = av2;
! 2421: }
! 2422: }
! 2423: /* Echec? */
! 2424: avma = ltop;
! 2425: return gzero;
! 2426: }
! 2427:
! 2428: /* Groupe S4 */
! 2429: static void
! 2430: s4makelift(GEN u, struct galois_lift *gl, GEN liftpow)
! 2431: {
! 2432: int i;
! 2433: liftpow[1] = (long) automorphismlift(u, gl, NULL);
! 2434: for (i = 2; i < lg(liftpow); i++)
! 2435: liftpow[i] = (long) FpXQ_mul((GEN) liftpow[i - 1], (GEN) liftpow[1],gl->TQ,gl->Q);
! 2436: }
! 2437: static long
! 2438: s4test(GEN u, GEN liftpow, struct galois_lift *gl, GEN phi)
! 2439: {
! 2440: ulong ltop=avma;
! 2441: GEN res;
! 2442: int bl,i,d = lgef(u)-2;
! 2443: if (DEBUGLEVEL >= 6)
! 2444: timer2();
! 2445: if ( !d ) return 0;
! 2446: res=(GEN)u[2];
! 2447: for (i = 1; i < d; i++)
! 2448: {
! 2449: if (lgef(liftpow[i])>2)
! 2450: res=addii(res,mulii(gmael(liftpow,i,2), (GEN) u[i + 2]));
! 2451: }
! 2452: res=modii(mulii(res,gl->den),gl->Q);
! 2453: if (cmpii(res,gl->gb->bornesol)>0
! 2454: && cmpii(res,subii(gl->Q,gl->gb->bornesol))<0)
! 2455: {
! 2456: avma=ltop;
! 2457: return 0;
! 2458: }
! 2459: res = (GEN) scalarpol((GEN)u[2],varn(u));
! 2460: for (i = 1; i < d ; i++)
! 2461: {
! 2462: GEN z;
! 2463: z = FpX_Fp_mul((GEN) liftpow[i], (GEN) u[i + 2],NULL);
! 2464: res = FpX_add(res,z ,gl->Q);
! 2465: }
! 2466: res = FpX_center(FpX_Fp_mul(res,gl->den,gl->Q), gl->Q);
! 2467: if (DEBUGLEVEL >= 6)
! 2468: msgtimer("s4test()");
! 2469: bl = poltopermtest(res, gl, phi);
! 2470: avma=ltop;
! 2471: return bl;
! 2472: }
! 2473: static GEN
! 2474: s4releveauto(GEN misom,GEN Tmod,GEN Tp,GEN p,long a1,long a2,long a3,long a4,long a5,long a6)
! 2475: {
! 2476: ulong ltop=avma;
! 2477: GEN u1,u2,u3,u4,u5;
! 2478: GEN pu1,pu2,pu3,pu4;
! 2479: pu1=FpX_mul( (GEN) Tmod[a2], (GEN) Tmod[a1],p);
! 2480: u1 = FpX_chinese_coprime(gmael(misom,a1,a2),gmael(misom,a2,a1),
! 2481: (GEN) Tmod[a2], (GEN) Tmod[a1],pu1,p);
! 2482: pu2=FpX_mul( (GEN) Tmod[a4], (GEN) Tmod[a3],p);
! 2483: u2 = FpX_chinese_coprime(gmael(misom,a3,a4),gmael(misom,a4,a3),
! 2484: (GEN) Tmod[a4], (GEN) Tmod[a3],pu2,p);
! 2485: pu3=FpX_mul( (GEN) Tmod[a6], (GEN) Tmod[a5],p);
! 2486: u3 = FpX_chinese_coprime(gmael(misom,a5,a6),gmael(misom,a6,a5),
! 2487: (GEN) Tmod[a6], (GEN) Tmod[a5],pu3,p);
! 2488: pu4=FpX_mul(pu1,pu2,p);
! 2489: u4 = FpX_chinese_coprime(u1,u2,pu1,pu2,pu4,p);
! 2490: u5 = FpX_chinese_coprime(u4,u3,pu4,pu3,Tp,p);
! 2491: return gerepileupto(ltop,u5);
! 2492: }
! 2493: GEN
! 2494: s4galoisgen(struct galois_lift *gl)
! 2495: {
! 2496: struct galois_testlift gt;
! 2497: ulong ltop = avma, av, ltop2;
! 2498: GEN Tmod, isom, isominv, misom;
! 2499: int i, j;
! 2500: GEN sg;
! 2501: GEN sigma, tau, phi;
! 2502: GEN res, ry;
! 2503: GEN pj;
! 2504: GEN p,Q,TQ,Tp;
! 2505: GEN bezoutcoeff, pauto, liftpow, aut;
! 2506:
! 2507: p = gl->p;
! 2508: Q = gl->Q;
! 2509: res = cgetg(3, t_VEC);
! 2510: ry = cgetg(5, t_VEC);
! 2511: res[1] = (long) ry;
! 2512: for (i = 1; i < lg(ry); i++)
! 2513: ry[i] = lgetg(lg(gl->L), t_VECSMALL);
! 2514: ry = cgetg(5, t_VECSMALL);
! 2515: res[2] = (long) ry;
! 2516: ry[1] = 2;
! 2517: ry[2] = 2;
! 2518: ry[3] = 3;
! 2519: ry[4] = 2;
! 2520: ltop2 = avma;
! 2521: sg = cgetg(7, t_VECSMALL);
! 2522: pj = cgetg(7, t_VECSMALL);
! 2523: sigma = cgetg(lg(gl->L), t_VECSMALL);
! 2524: tau = cgetg(lg(gl->L), t_VECSMALL);
! 2525: phi = cgetg(lg(gl->L), t_VECSMALL);
! 2526: for (i = 1; i < lg(sg); i++)
! 2527: sg[i] = i;
! 2528: Tp = FpX_red(gl->T,p);
! 2529: TQ = gl->TQ;
! 2530: Tmod = lift((GEN) factmod(gl->T, p)[1]);
! 2531: isom = cgetg(lg(Tmod), t_VEC);
! 2532: isominv = cgetg(lg(Tmod), t_VEC);
! 2533: misom = cgetg(lg(Tmod), t_MAT);
! 2534: aut=galoisdolift(gl, NULL);
! 2535: inittestlift(aut,Tmod, gl, >);
! 2536: bezoutcoeff = gt.bezoutcoeff;
! 2537: pauto = gt.pauto;
! 2538: for (i = 1; i < lg(pj); i++)
! 2539: pj[i] = 0;
! 2540: for (i = 1; i < lg(isom); i++)
! 2541: {
! 2542: misom[i] = lgetg(lg(Tmod), t_COL);
! 2543: isom[i] = (long) Fp_isom((GEN) Tmod[1], (GEN) Tmod[i], p);
! 2544: if (DEBUGLEVEL >= 6)
! 2545: fprintferr("S4GaloisConj:Computing isomorphisms %d:%Z\n", i,
! 2546: (GEN) isom[i]);
! 2547: isominv[i] = (long) Fp_inv_isom((GEN) isom[i], (GEN) Tmod[i],p);
! 2548: }
! 2549: for (i = 1; i < lg(isom); i++)
! 2550: for (j = 1; j < lg(isom); j++)
! 2551: mael(misom,i,j) = (long) FpX_FpXQ_compo((GEN) isominv[i],(GEN) isom[j],
! 2552: (GEN) Tmod[j],p);
! 2553: liftpow = cgetg(24, t_VEC);
! 2554: av = avma;
! 2555: for (i = 0; i < 3; i++)
! 2556: {
! 2557: ulong av2;
! 2558: GEN u;
! 2559: int j1, j2, j3;
! 2560: ulong avm1, avm2;
! 2561: GEN u1, u2, u3;
! 2562: if (i)
! 2563: {
! 2564: int x;
! 2565: x = sg[3];
! 2566: if (i == 1)
! 2567: {
! 2568: sg[3] = sg[2];
! 2569: sg[2] = x;
! 2570: }
! 2571: else
! 2572: {
! 2573: sg[3] = sg[1];
! 2574: sg[1] = x;
! 2575: }
! 2576: }
! 2577: u=s4releveauto(misom,Tmod,Tp,p,sg[1],sg[2],sg[3],sg[4],sg[5],sg[6]);
! 2578: s4makelift(u, gl, liftpow);
! 2579: av2 = avma;
! 2580: for (j1 = 0; j1 < 4; j1++)
! 2581: {
! 2582: u1 = FpX_add(FpXQ_mul((GEN) bezoutcoeff[sg[5]],
! 2583: (GEN) pauto[1 + j1],TQ,Q),
! 2584: FpXQ_mul((GEN) bezoutcoeff[sg[6]],
! 2585: (GEN) pauto[((-j1) & 3) + 1],TQ,Q),Q);
! 2586: avm1 = avma;
! 2587: for (j2 = 0; j2 < 4; j2++)
! 2588: {
! 2589: u2 = FpX_add(u1, FpXQ_mul((GEN) bezoutcoeff[sg[3]],
! 2590: (GEN) pauto[1 + j2],TQ,Q),NULL);
! 2591: u2 = FpX_add(u2, FpXQ_mul((GEN) bezoutcoeff[sg[4]], (GEN)
! 2592: pauto[((-j2) & 3) + 1],TQ,Q),Q);
! 2593: avm2 = avma;
! 2594: for (j3 = 0; j3 < 4; j3++)
! 2595: {
! 2596: u3 = FpX_add(u2, FpXQ_mul((GEN) bezoutcoeff[sg[1]],
! 2597: (GEN) pauto[1 + j3],TQ,Q),NULL);
! 2598: u3 = FpX_add(u3, FpXQ_mul((GEN) bezoutcoeff[sg[2]], (GEN)
! 2599: pauto[((-j3) & 3) + 1],TQ,Q),Q);
! 2600: if (DEBUGLEVEL >= 4)
! 2601: fprintferr("S4GaloisConj:Testing %d/3:%d/4:%d/4:%d/4:%Z\n",
! 2602: i, j1,j2, j3, sg);
! 2603: if (s4test(u3, liftpow, gl, sigma))
! 2604: {
! 2605: pj[1] = j3;
! 2606: pj[2] = j2;
! 2607: pj[3] = j1;
! 2608: goto suites4;
! 2609: }
! 2610: avma = avm2;
! 2611: }
! 2612: avma = avm1;
! 2613: }
! 2614: avma = av2;
! 2615: }
! 2616: avma = av;
! 2617: }
! 2618: avma = ltop;
! 2619: return gzero;
! 2620: suites4:
! 2621: if (DEBUGLEVEL >= 4)
! 2622: fprintferr("S4GaloisConj:sigma=%Z\n", sigma);
! 2623: if (DEBUGLEVEL >= 4)
! 2624: fprintferr("S4GaloisConj:pj=%Z\n", pj);
! 2625: avma = av;
! 2626: for (j = 1; j <= 3; j++)
! 2627: {
! 2628: ulong av2;
! 2629: GEN u;
! 2630: int w, l;
! 2631: int z;
! 2632: z = sg[1]; sg[1] = sg[3]; sg[3] = sg[5]; sg[5] = z;
! 2633: z = sg[2]; sg[2] = sg[4]; sg[4] = sg[6]; sg[6] = z;
! 2634: z = pj[1]; pj[1] = pj[2]; pj[2] = pj[3]; pj[3] = z;
! 2635: for (l = 0; l < 2; l++)
! 2636: {
! 2637: u=s4releveauto(misom,Tmod,Tp,p,sg[1],sg[3],sg[2],sg[4],sg[5],sg[6]);
! 2638: s4makelift(u, gl, liftpow);
! 2639: av2 = avma;
! 2640: for (w = 0; w < 4; w += 2)
! 2641: {
! 2642: GEN uu;
! 2643: long av3;
! 2644: pj[6] = (w + pj[3]) & 3;
! 2645: uu =FpX_add(FpXQ_mul((GEN) bezoutcoeff[sg[5]],
! 2646: (GEN) pauto[(pj[6] & 3) + 1],TQ,Q),
! 2647: FpXQ_mul((GEN) bezoutcoeff[sg[6]],
! 2648: (GEN) pauto[((-pj[6]) & 3) + 1],TQ,Q),Q);
! 2649: av3 = avma;
! 2650: for (i = 0; i < 4; i++)
! 2651: {
! 2652: GEN u;
! 2653: pj[4] = i;
! 2654: pj[5] = (i + pj[2] - pj[1]) & 3;
! 2655: if (DEBUGLEVEL >= 4)
! 2656: fprintferr("S4GaloisConj:Testing %d/3:%d/2:%d/2:%d/4:%Z:%Z\n",
! 2657: j - 1, w >> 1, l, i, sg, pj);
! 2658: u = FpX_add(uu, FpXQ_mul((GEN) pauto[(pj[4] & 3) + 1],
! 2659: (GEN) bezoutcoeff[sg[1]],TQ,Q),Q);
! 2660: u = FpX_add(u, FpXQ_mul((GEN) pauto[((-pj[4]) & 3) + 1],
! 2661: (GEN) bezoutcoeff[sg[3]],TQ,Q),Q);
! 2662: u = FpX_add(u, FpXQ_mul((GEN) pauto[(pj[5] & 3) + 1],
! 2663: (GEN) bezoutcoeff[sg[2]],TQ,Q),Q);
! 2664: u = FpX_add(u, FpXQ_mul((GEN) pauto[((-pj[5]) & 3) + 1],
! 2665: (GEN) bezoutcoeff[sg[4]],TQ,Q),Q);
! 2666: if (s4test(u, liftpow, gl, tau))
! 2667: goto suites4_2;
! 2668: avma = av3;
! 2669: }
! 2670: avma = av2;
! 2671: }
! 2672: z = sg[4];
! 2673: sg[4] = sg[3];
! 2674: sg[3] = z;
! 2675: pj[2] = (-pj[2]) & 3;
! 2676: avma = av;
! 2677: }
! 2678: }
! 2679: avma = ltop;
! 2680: return gzero;
! 2681: suites4_2:
! 2682: avma = av;
! 2683: {
! 2684: int abc, abcdef;
! 2685: GEN u;
! 2686: ulong av2;
! 2687: abc = (pj[1] + pj[2] + pj[3]) & 3;
! 2688: abcdef = (((abc + pj[4] + pj[5] - pj[6]) & 3) >> 1);
! 2689: u = s4releveauto(misom,Tmod,Tp,p,sg[1],sg[4],sg[2],sg[5],sg[3],sg[6]);
! 2690: s4makelift(u, gl, liftpow);
! 2691: av2 = avma;
! 2692: for (j = 0; j < 8; j++)
! 2693: {
! 2694: int h, g, i;
! 2695: h = j & 3;
! 2696: g = abcdef + ((j & 4) >> 1);
! 2697: i = h + abc - g;
! 2698: u = FpXQ_mul((GEN) pauto[(g & 3) + 1],
! 2699: (GEN) bezoutcoeff[sg[1]],TQ,Q);
! 2700: u = FpX_add(u, FpXQ_mul((GEN) pauto[((-g) & 3) + 1],
! 2701: (GEN) bezoutcoeff[sg[4]],TQ,Q),NULL);
! 2702: u = FpX_add(u, FpXQ_mul((GEN) pauto[(h & 3) + 1],
! 2703: (GEN) bezoutcoeff[sg[2]],TQ,Q),NULL);
! 2704: u = FpX_add(u, FpXQ_mul((GEN) pauto[((-h) & 3) + 1],
! 2705: (GEN) bezoutcoeff[sg[5]],TQ,Q),NULL);
! 2706: u = FpX_add(u, FpXQ_mul((GEN) pauto[(i & 3) + 1],
! 2707: (GEN) bezoutcoeff[sg[3]],TQ,Q),NULL);
! 2708: u = FpX_add(u, FpXQ_mul((GEN) pauto[((-i) & 3) + 1],
! 2709: (GEN) bezoutcoeff[sg[6]],TQ,Q),Q);
! 2710: if (DEBUGLEVEL >= 4)
! 2711: fprintferr("S4GaloisConj:Testing %d/8 %d:%d:%d\n",
! 2712: j, g & 3, h & 3, i & 3);
! 2713: if (s4test(u, liftpow, gl, phi))
! 2714: break;
! 2715: avma = av2;
! 2716: }
! 2717: }
! 2718: if (j == 8)
! 2719: {
! 2720: avma = ltop;
! 2721: return gzero;
! 2722: }
! 2723: for (i = 1; i < lg(gl->L); i++)
! 2724: {
! 2725: ((GEN **) res)[1][1][i] = sigma[tau[i]];
! 2726: ((GEN **) res)[1][2][i] = phi[sigma[tau[phi[i]]]];
! 2727: ((GEN **) res)[1][3][i] = phi[sigma[i]];
! 2728: ((GEN **) res)[1][4][i] = sigma[i];
! 2729: }
! 2730: avma = ltop2;
! 2731: return res;
! 2732: }
! 2733: struct galois_frobenius
! 2734: {
! 2735: long p;
! 2736: long fp;
! 2737: long deg;
! 2738: GEN Tmod;
! 2739: GEN psi;
! 2740: };
! 2741:
! 2742: /*Warning : the output of this function is not gerepileupto
! 2743: * compatible...*/
! 2744: static GEN
! 2745: galoisfindgroups(GEN lo, GEN sg, long f)
! 2746: {
! 2747: ulong ltop=avma;
! 2748: GEN V,W;
! 2749: long i,j,k;
! 2750: V=cgetg(lg(lo),t_VEC);
! 2751: for(j=1,i=1;i<lg(lo);i++)
! 2752: {
! 2753: ulong av=avma;
! 2754: GEN U;
! 2755: W=cgetg(lg(lo[i]),t_VECSMALL);
! 2756: for(k=1;k<lg(lo[i]);k++)
! 2757: W[k]=mael(lo,i,k)%f;
! 2758: sortvecsmall(W);
! 2759: U=uniqvecsmall(W);
! 2760: if (gegal(U, sg))
! 2761: {
! 2762: cgiv(U);
! 2763: V[j++]=lo[i];
! 2764: }
! 2765: else
! 2766: avma=av;
! 2767: }
! 2768: setlg(V,j);
! 2769: /*warning components of V point to W*/
! 2770: return gerepileupto(ltop,V);
! 2771: }
! 2772:
! 2773: static long
! 2774: galoisfrobeniustest(GEN aut, struct galois_lift *gl, GEN frob)
! 2775: {
! 2776: ulong ltop=avma;
! 2777: GEN tlift = FpX_center(FpX_Fp_mul(aut,gl->den,gl->Q), gl->Q);
! 2778: long res=poltopermtest(tlift, gl, frob);
! 2779: avma=ltop;
! 2780: return res;
! 2781: }
! 2782:
! 2783: static GEN
! 2784: galoismakepsiab(long g)
! 2785: {
! 2786: GEN psi=cgetg(g+1,t_VECSMALL);
! 2787: long i;
! 2788: for (i = 1; i <= g; i++)
! 2789: psi[i] = 1;
! 2790: return psi;
! 2791: }
! 2792:
! 2793: static GEN
! 2794: galoismakepsi(long g, GEN sg, GEN pf)
! 2795: {
! 2796: GEN psi=cgetg(g+1,t_VECSMALL);
! 2797: long i;
! 2798: for (i = 1; i < g; i++)
! 2799: psi[i] = sg[pf[i]];
! 2800: psi[g]=sg[1];
! 2801: return psi;
! 2802: }
! 2803:
! 2804: static GEN
! 2805: galoisfrobeniuslift(GEN T, GEN den, GEN L, GEN Lden, long gmask,
! 2806: struct galois_frobenius *gf, struct galois_borne *gb)
! 2807: {
! 2808: ulong ltop=avma,av2;
! 2809: struct galois_testlift gt;
! 2810: struct galois_lift gl;
! 2811: GEN res;
! 2812: long i,j,k;
! 2813: long n=lg(L)-1, deg=1, g=lg(gf->Tmod)-1;
! 2814: GEN F,Fp,Fe;
! 2815: GEN ip=stoi(gf->p), aut, frob;
! 2816: if (DEBUGLEVEL >= 4)
! 2817: fprintferr("GaloisConj:p=%ld deg=%ld fp=%ld\n", gf->p, deg, gf->fp);
! 2818: res = cgetg(lg(L), t_VECSMALL);
! 2819: gf->psi = galoismakepsiab(g);
! 2820: av2=avma;
! 2821: initlift(T, den, ip, L, Lden, gb, &gl);
! 2822: aut = galoisdolift(&gl, res);
! 2823: if (!aut || galoisfrobeniustest(aut,&gl,res))
! 2824: {
! 2825: avma=av2;
! 2826: gf->deg = gf->fp;
! 2827: return res;
! 2828: }
! 2829: inittestlift(aut,gf->Tmod, &gl, >);
! 2830: gt.C=cgetg(gf->fp+1,t_VEC);
! 2831: for (i = 1; i <= gf->fp; i++)
! 2832: {
! 2833: gt.C[i]=lgetg(gt.g+1,t_VECSMALL);
! 2834: for(j = 1; j <= gt.g; j++)
! 2835: mael(gt.C,i,j)=0;
! 2836: }
! 2837: gt.Cd=gcopy(gt.C);
! 2838:
! 2839: F=factor(stoi(gf->fp));
! 2840: Fp=gtovecsmall((GEN)F[1]);
! 2841: Fe=gtovecsmall((GEN)F[2]);
! 2842: frob = cgetg(lg(L), t_VECSMALL);
! 2843: for(k=lg(Fp)-1;k>=1;k--)
! 2844: {
! 2845: long btop=avma;
! 2846: GEN psi=NULL,fres=NULL,sg;
! 2847: long el=gf->fp, dg=1, dgf=1;
! 2848: long e,pr;
! 2849: sg=permidentity(1);
! 2850: for(e=1;e<=Fe[k];e++)
! 2851: {
! 2852: long l;
! 2853: GEN lo;
! 2854: GEN pf;
! 2855: dg *= Fp[k]; el /= Fp[k];
! 2856: if ( DEBUGLEVEL>=4 )
! 2857: fprintferr("Trying degre %d.\n",dg);
! 2858: if (galoisfrobeniustest((GEN)gt.pauto[el+1],&gl,frob))
! 2859: {
! 2860: dgf = dg;
! 2861: psi = galoismakepsiab(g);
! 2862: fres= gcopy(frob);
! 2863: continue;
! 2864: }
! 2865: disable_dbg(0);
! 2866: lo = listsousgroupes(dg, n / gf->fp);
! 2867: disable_dbg(-1);
! 2868: if (e!=1)
! 2869: lo = galoisfindgroups(lo, sg, dgf);
! 2870: if (DEBUGLEVEL >= 4)
! 2871: fprintferr("Galoisconj:Subgroups list:%Z\n", lo);
! 2872: for (l = 1; l < lg(lo); l++)
! 2873: if ( lg(lo[l])>2 &&
! 2874: frobeniusliftall((GEN)lo[l], el, &pf, &gl, >, frob))
! 2875: {
! 2876: sg = gcopy((GEN)lo[l]);
! 2877: psi = galoismakepsi(g,sg,pf);
! 2878: dgf = dg;
! 2879: fres=gcopy(frob);
! 2880: break;
! 2881: }
! 2882: if ( l == lg(lo) )
! 2883: break;
! 2884: }
! 2885: if (dgf==1) { avma=btop; continue; }
! 2886: pr=deg*dgf;
! 2887: if (deg==1)
! 2888: {
! 2889: for(i=1;i<lg(res);i++) res[i]=fres[i];
! 2890: for(i=1;i<lg(psi);i++) gf->psi[i]=psi[i];
! 2891: }
! 2892: else
! 2893: {
! 2894: GEN cp=permapply(res,fres);
! 2895: for(i=1;i<lg(res);i++) res[i]=cp[i];
! 2896: for(i=1;i<lg(psi);i++) gf->psi[i]=(dgf*gf->psi[i]+deg*psi[i])%pr;
! 2897: }
! 2898: deg=pr;
! 2899: avma=btop;
! 2900: }
! 2901: for (i = 1; i <= gf->fp; i++)
! 2902: for (j = 1; j <= gt.g; j++)
! 2903: if (mael(gt.C,i,j))
! 2904: gunclone(gmael(gt.C,i,j));
! 2905: if (DEBUGLEVEL>=4 && res)
! 2906: fprintferr("Best lift: %d\n",deg);
! 2907: if (deg==1)
! 2908: {
! 2909: avma=ltop;
! 2910: return NULL;
! 2911: }
! 2912: else
! 2913: {
! 2914: /*We need to normalise result so that psi[g]=1*/
! 2915: long im=itos(mpinvmod(stoi(gf->psi[g]),stoi(gf->fp)));
! 2916: GEN cp=permcyclepow(permorbite(res), im);
! 2917: for(i=1;i<lg(res);i++) res[i]=cp[i];
! 2918: for(i=1;i<lg(gf->psi);i++) gf->psi[i]=mulssmod(im,gf->psi[i],gf->fp);
! 2919: avma=av2;
! 2920: gf->deg=deg;
! 2921: return res;
! 2922: }
! 2923: }
! 2924:
! 2925: static GEN
! 2926: galoisfindfrobenius(GEN T, GEN L, GEN M, GEN den, struct galois_frobenius *gf,
! 2927: struct galois_borne *gb, const struct galois_analysis *ga)
! 2928: {
! 2929: ulong ltop=avma,lbot;
! 2930: long try=0;
! 2931: long n = degpol(T), deg, gmask;
! 2932: byteptr primepointer = ga->primepointer;
! 2933: GEN Lden,frob;
! 2934: Lden=makeLden(L,den,gb);
! 2935: gf->deg=ga->deg;gf->p=ga->p; deg=ga->deg;
! 2936: gmask=(ga->group&ga_ext_2)?3:1;
! 2937: for (;;)
! 2938: {
! 2939: ulong av = avma;
! 2940: long isram;
! 2941: long i,c;
! 2942: GEN ip,Tmod;
! 2943: ip = stoi(gf->p);
! 2944: Tmod = lift((GEN) factmod(T, ip));
! 2945: isram = 0;
! 2946: for (i = 1; i < lg(Tmod[2]) && !isram; i++)
! 2947: if (!gcmp1(gmael(Tmod,2,i)))
! 2948: isram = 1;
! 2949: if (isram == 0)
! 2950: {
! 2951: gf->fp = degpol(gmael(Tmod,1,1));
! 2952: for (i = 2; i < lg(Tmod[1]); i++)
! 2953: if (degpol(gmael(Tmod,1,i)) != gf->fp)
! 2954: {
! 2955: avma = ltop;
! 2956: return NULL; /* Not Galois polynomial */
! 2957: }
! 2958: lbot=avma;
! 2959: gf->Tmod=gcopy((GEN)Tmod[1]);
! 2960: if ( ((gmask&1) && gf->fp % deg == 0) || ((gmask&2) && gf->fp % 2== 0) )
! 2961: {
! 2962: frob=galoisfrobeniuslift(T, den, L, Lden, gmask, gf, gb);
! 2963: if (frob)
! 2964: {
! 2965: GEN *gptr[3];
! 2966: gptr[0]=&gf->Tmod;
! 2967: gptr[1]=&gf->psi;
! 2968: gptr[2]=&frob;
! 2969: gerepilemanysp(ltop,lbot,gptr,3);
! 2970: return frob;
! 2971: }
! 2972: if ((ga->group&ga_all_normal) && gf->fp % deg == 0)
! 2973: gmask&=~1;
! 2974: /*The first prime degree is always divisible by deg, so we don't
! 2975: * have to worry about ext_2 being used before regular supersolvable*/
! 2976: if (!gmask)
! 2977: {
! 2978: avma = ltop;
! 2979: return NULL;
! 2980: }
! 2981: try++;
! 2982: if ( (ga->group&ga_non_wss) && try > n )
! 2983: err(warner, "galoisconj _may_ hang up for this polynomial");
! 2984: }
! 2985: }
! 2986: c = *primepointer++;
! 2987: if (!c)
! 2988: err(primer1);
! 2989: gf->p += c;
! 2990: if (DEBUGLEVEL >= 4)
! 2991: fprintferr("GaloisConj:next p=%ld\n", gf->p);
! 2992: avma = av;
! 2993: }
! 2994: }
! 2995:
! 2996: GEN
! 2997: galoisgen(GEN T, GEN L, GEN M, GEN den, struct galois_borne *gb,
! 2998: const struct galois_analysis *ga)
! 2999: {
! 3000: struct galois_test td;
! 3001: struct galois_frobenius gf;
! 3002: ulong ltop = avma, lbot, ltop2;
! 3003: long n, p, deg;
! 3004: long fp,gp;
! 3005: long x;
! 3006: int i, j;
! 3007: GEN Lden, sigma;
! 3008: GEN Tmod, res, pf = gzero, split, ip;
! 3009: GEN frob;
! 3010: GEN O;
! 3011: GEN P, PG, PL, Pden, PM, Pmod, Pp, Pladicabs;
! 3012: n = degpol(T);
! 3013: if (!ga->deg)
! 3014: return gzero;
! 3015: x = varn(T);
! 3016: if (DEBUGLEVEL >= 9)
! 3017: fprintferr("GaloisConj:denominator:%Z\n", den);
! 3018: if (n == 12 && ga->ord==3) /* A4 is very probable,so test it first */
! 3019: {
! 3020: long av = avma;
! 3021: if (DEBUGLEVEL >= 4)
! 3022: fprintferr("GaloisConj:Testing A4 first\n");
! 3023: inittest(L, M, gb->bornesol, gb->ladicsol, &td);
! 3024: lbot = avma;
! 3025: PG = a4galoisgen(T, &td);
! 3026: freetest(&td);
! 3027: if (PG != gzero)
! 3028: return gerepile(ltop, lbot, PG);
! 3029: avma = av;
! 3030: }
! 3031: if (n == 24 && ga->ord==3) /* S4 is very probable,so test it first */
! 3032: {
! 3033: long av = avma;
! 3034: struct galois_lift gl;
! 3035: if (DEBUGLEVEL >= 4)
! 3036: fprintferr("GaloisConj:Testing S4 first\n");
! 3037: lbot = avma;
! 3038: Lden=makeLden(L,den,gb);
! 3039: initlift(T, den, stoi(ga->p4), L, Lden, gb, &gl);
! 3040: PG = s4galoisgen(&gl);
! 3041: if (PG != gzero)
! 3042: return gerepile(ltop, lbot, PG);
! 3043: avma = av;
! 3044: }
! 3045: frob=galoisfindfrobenius(T, L, M, den, &gf, gb, ga);
! 3046: if (!frob)
! 3047: {
! 3048: ltop=avma;
! 3049: return gzero;
! 3050: }
! 3051: p=gf.p;deg=gf.deg;
! 3052: ip=stoi(p);
! 3053: Tmod=gf.Tmod;
! 3054: fp=gf.fp; gp=n/fp;
! 3055: O = permorbite(frob);
! 3056: split = splitorbite(O);
! 3057: deg=lg(O[1])-1;
! 3058: sigma = permtopol(frob, L, M, den, gb->ladicabs, x);
! 3059: if (DEBUGLEVEL >= 9)
! 3060: fprintferr("GaloisConj:Orbite:%Z\n", O);
! 3061: if (deg == n) /* Cyclique */
! 3062: {
! 3063: lbot = avma;
! 3064: res = cgetg(3, t_VEC);
! 3065: res[1] = lgetg(lg(split[1]), t_VEC);
! 3066: res[2] = lgetg(lg(split[2]), t_VECSMALL);
! 3067: for (i = 1; i < lg(split[1]); i++)
! 3068: {
! 3069: mael(res,1,i) = lcopy(gmael(split,1,i));
! 3070: mael(res,2,i) = mael(split,2,i);
! 3071: }
! 3072: return gerepile(ltop, lbot, res);
! 3073: }
! 3074: if (DEBUGLEVEL >= 9)
! 3075: fprintferr("GaloisConj:Frobenius:%Z\n", sigma);
! 3076: {
! 3077: GEN V, Tp, Sp, sym, dg;
! 3078: sym=cgetg(lg(L),t_VECSMALL);
! 3079: dg=cgetg(lg(L),t_VECSMALL);
! 3080: V = fixedfieldsympol(O, L, gb->ladicabs, gb->l, ip, sym, dg, x);
! 3081: P=(GEN)V[2];
! 3082: PL=(GEN)V[1];
! 3083: Tp = FpX_red(T,ip);
! 3084: Sp = fixedfieldpolsigma(sigma,ip,Tp,sym,dg,deg);
! 3085: Pmod = fixedfieldfactmod(Sp,ip,Tmod);
! 3086: Pp = FpX_red(P,ip);
! 3087: if (DEBUGLEVEL >= 4)
! 3088: {
! 3089: fprintferr("GaloisConj:P=%Z \n", P);
! 3090: fprintferr("GaloisConj:psi=%Z\n", gf.psi);
! 3091: fprintferr("GaloisConj:Sp=%Z\n", Sp);
! 3092: fprintferr("GaloisConj:Pmod=%Z\n", Pmod);
! 3093: fprintferr("GaloisConj:Tmod=%Z\n", Tmod);
! 3094: }
! 3095: if ( n == (deg<<1))
! 3096: {
! 3097: PG=cgetg(3,t_VEC);
! 3098: PG[1]=lgetg(2,t_VEC);
! 3099: PG[2]=lgetg(2,t_VECSMALL);
! 3100: mael(PG,1,1)=lgetg(3,t_VECSMALL);
! 3101: mael(PG,2,1)=2;
! 3102: mael3(PG,1,1,1)=2;
! 3103: mael3(PG,1,1,2)=1;
! 3104: Pden=NULL; PM=NULL; Pladicabs=NULL;/*make KB happy*/
! 3105: }
! 3106: else
! 3107: {
! 3108: struct galois_analysis Pga;
! 3109: struct galois_borne Pgb;
! 3110: galoisanalysis(P, &Pga, 0, 0);
! 3111: if (Pga.deg == 0)
! 3112: {
! 3113: avma = ltop;
! 3114: return gzero; /* Avoid computing the discriminant */
! 3115: }
! 3116: Pgb.l = gb->l;
! 3117: Pden = galoisborne(P, NULL, &Pgb, Pga.ppp);
! 3118: Pladicabs=Pgb.ladicabs;
! 3119: if (Pgb.valabs > gb->valabs)
! 3120: {
! 3121: if (DEBUGLEVEL>=4)
! 3122: fprintferr("GaloisConj:increase prec of p-adic roots of %ld.\n"
! 3123: ,Pgb.valabs-gb->valabs);
! 3124: PL = rootpadicliftroots(P,PL,gb->l,Pgb.valabs);
! 3125: }
! 3126: PM = vandermondeinversemod(PL, P, Pden, Pgb.ladicabs);
! 3127: PG = galoisgen(P, PL, PM, Pden, &Pgb, &Pga);
! 3128: }
! 3129: }
! 3130: if (DEBUGLEVEL >= 4)
! 3131: fprintferr("GaloisConj:Retour sur Terre:%Z\n", PG);
! 3132: if (PG == gzero)
! 3133: {
! 3134: avma = ltop;
! 3135: return gzero;
! 3136: }
! 3137: inittest(L, M, gb->bornesol, gb->ladicsol, &td);
! 3138: lbot = avma;
! 3139: res = cgetg(3, t_VEC);
! 3140: res[1] = lgetg(lg(PG[1]) + lg(split[1]) - 1, t_VEC);
! 3141: res[2] = lgetg(lg(PG[1]) + lg(split[1]) - 1, t_VECSMALL);
! 3142: for (i = 1; i < lg(split[1]); i++)
! 3143: {
! 3144: ((GEN **) res)[1][i] = gcopy(((GEN **) split)[1][i]);
! 3145: ((GEN *) res)[2][i] = ((GEN *) split)[2][i];
! 3146: }
! 3147: for (i = lg(split[1]); i < lg(res[1]); i++)
! 3148: ((GEN **) res)[1][i] = cgetg(n + 1, t_VECSMALL);
! 3149: ltop2 = avma;
! 3150: for (j = 1; j < lg(PG[1]); j++)
! 3151: {
! 3152: GEN B, tau;
! 3153: long t, g;
! 3154: long w, sr, dss;
! 3155: if (DEBUGLEVEL >= 6)
! 3156: fprintferr("GaloisConj:G[%d]=%Z d'ordre relatif %d\n", j,
! 3157: ((GEN **) PG)[1][j], ((long **) PG)[2][j]);
! 3158: B = permorbite(((GEN **) PG)[1][j]);
! 3159: if (DEBUGLEVEL >= 6)
! 3160: fprintferr("GaloisConj:B=%Z\n", B);
! 3161: if (n == (deg<<1))
! 3162: tau=deg1pol(stoi(-1),negi((GEN)P[3]),x);
! 3163: else
! 3164: tau = permtopol(gmael(PG,1,j), PL, PM, Pden, Pladicabs, x);
! 3165: if (DEBUGLEVEL >= 9)
! 3166: fprintferr("GaloisConj:Paut=%Z\n", tau);
! 3167: /*tau not in Z[X]*/
! 3168: tau = lift(gmul(tau,gmodulcp(gun,ip)));
! 3169: tau = FpX_FpXQ_compo((GEN) Pmod[gp], tau,Pp,ip);
! 3170: tau = FpX_gcd(Pp, tau,ip);
! 3171: /*normalize gcd*/
! 3172: tau = FpX_Fp_mul(tau,mpinvmod((GEN) tau[lgef(tau) - 1],ip),ip);
! 3173: if (DEBUGLEVEL >= 6)
! 3174: fprintferr("GaloisConj:tau=%Z\n", tau);
! 3175: for (g = 1; g < lg(Pmod); g++)
! 3176: if (gegal(tau, (GEN) Pmod[g]))
! 3177: break;
! 3178: if (DEBUGLEVEL >= 6)
! 3179: fprintferr("GaloisConj:g=%ld\n", g);
! 3180: if (g == lg(Pmod))
! 3181: {
! 3182: freetest(&td);
! 3183: avma = ltop;
! 3184: return gzero;
! 3185: }
! 3186: w = gf.psi[g];
! 3187: dss = deg / cgcd(deg, w - 1);
! 3188: sr = 1;
! 3189: for (i = 1; i < lg(B[1]) - 1; i++)
! 3190: sr = (1 + sr * w) % deg;
! 3191: sr = cgcd(sr, deg);
! 3192: if (DEBUGLEVEL >= 6)
! 3193: fprintferr("GaloisConj:w=%ld [%ld] sr=%ld dss=%ld\n", w, deg, sr, dss);
! 3194: for (t = 0; t < sr; t += dss)
! 3195: {
! 3196: pf = testpermutation(O, B, w, t, sr, &td);
! 3197: if (pf != gzero)
! 3198: break;
! 3199: }
! 3200: if (pf == gzero)
! 3201: {
! 3202: freetest(&td);
! 3203: avma = ltop;
! 3204: return gzero;
! 3205: }
! 3206: else
! 3207: {
! 3208: int i;
! 3209: for (i = 1; i <= n; i++)
! 3210: ((GEN **) res)[1][lg(split[1]) + j - 1][i] = pf[i];
! 3211: ((GEN *) res)[2][lg(split[1]) + j - 1] = ((GEN *) PG)[2][j];
! 3212: }
! 3213: avma = ltop2;
! 3214: }
! 3215: if (DEBUGLEVEL >= 4)
! 3216: fprintferr("GaloisConj:Fini!\n");
! 3217: freetest(&td);
! 3218: return gerepile(ltop, lbot, res);
! 3219: }
! 3220:
! 3221: /*
! 3222: * T: polynome ou nf den optionnel: si présent doit etre un multiple du
! 3223: * dénominateur commun des solutions Si T est un nf et den n'est pas présent,
! 3224: * le denominateur de la base d'entiers est pris pour den
! 3225: */
! 3226: GEN
! 3227: galoisconj4(GEN T, GEN den, long flag, long karma)
! 3228: {
! 3229: ulong ltop = avma;
! 3230: GEN G, L, M, res, aut, grp=NULL;/*keep gcc happy on the wall*/
! 3231: struct galois_analysis ga;
! 3232: struct galois_borne gb;
! 3233: int n, i, j, k;
! 3234: if (typ(T) != t_POL)
! 3235: {
! 3236: T = checknf(T);
! 3237: if (!den)
! 3238: den = gcoeff((GEN) T[8], degpol(T[1]), degpol(T[1]));
! 3239: T = (GEN) T[1];
! 3240: }
! 3241: n = degpol(T);
! 3242: if (n <= 0)
! 3243: err(constpoler, "galoisconj4");
! 3244: for (k = 2; k <= n + 2; k++)
! 3245: if (typ(T[k]) != t_INT)
! 3246: err(talker, "polynomial not in Z[X] in galoisconj4");
! 3247: if (!gcmp1((GEN) T[n + 2]))
! 3248: err(talker, "non-monic polynomial in galoisconj4");
! 3249: n = degpol(T);
! 3250: if (n == 1) /* Too easy! */
! 3251: {
! 3252: if (!flag)
! 3253: {
! 3254: res = cgetg(2, t_COL);
! 3255: res[1] = (long) polx[varn(T)];
! 3256: return res;
! 3257: }
! 3258: ga.l = 3;
! 3259: ga.deg = 1;
! 3260: ga.ppp = 1;
! 3261: den = gun;
! 3262: }
! 3263: else
! 3264: galoisanalysis(T, &ga, 1, karma);
! 3265: if (ga.deg == 0)
! 3266: {
! 3267: avma = ltop;
! 3268: return stoi(ga.p); /* Avoid computing the discriminant */
! 3269: }
! 3270: if (den)
! 3271: {
! 3272: if (typ(den) != t_INT)
! 3273: err(talker, "Second arg. must be integer in galoisconj4");
! 3274: den = absi(den);
! 3275: }
! 3276: gb.l = stoi(ga.l);
! 3277: if (DEBUGLEVEL >= 1)
! 3278: timer2();
! 3279: den = galoisborne(T, den, &gb, ga.ppp);
! 3280: if (DEBUGLEVEL >= 1)
! 3281: msgtimer("galoisborne()");
! 3282: L = rootpadicfast(T, gb.l, gb.valabs);
! 3283: if (DEBUGLEVEL >= 1)
! 3284: msgtimer("rootpadicfast()");
! 3285: M = vandermondeinversemod(L, T, den, gb.ladicabs);
! 3286: if (DEBUGLEVEL >= 1)
! 3287: msgtimer("vandermondeinversemod()");
! 3288: if (n == 1)
! 3289: {
! 3290: G = cgetg(3, t_VEC);
! 3291: G[1] = lgetg(1, t_VECSMALL);
! 3292: G[2] = lgetg(1, t_VECSMALL);
! 3293: }
! 3294: else
! 3295: G = galoisgen(T, L, M, den, &gb, &ga);
! 3296: if (DEBUGLEVEL >= 6)
! 3297: fprintferr("GaloisConj:%Z\n", G);
! 3298: if (G == gzero)
! 3299: {
! 3300: avma = ltop;
! 3301: return gzero;
! 3302: }
! 3303: if (DEBUGLEVEL >= 1)
! 3304: timer2();
! 3305: if (flag)
! 3306: {
! 3307: grp = cgetg(9, t_VEC);
! 3308: grp[1] = (long) gcopy(T);
! 3309: grp[2] = lgetg(4,t_VEC); /*Make K.B. happy(8 components)*/
! 3310: mael(grp,2,1) = lstoi(ga.l);
! 3311: mael(grp,2,2) = lstoi(gb.valabs);
! 3312: mael(grp,2,3) = licopy(gb.ladicabs);
! 3313: grp[3] = (long) gcopy(L);
! 3314: grp[4] = (long) gcopy(M);
! 3315: grp[5] = (long) gcopy(den);
! 3316: grp[7] = (long) gcopy((GEN) G[1]);
! 3317: grp[8] = (long) gcopy((GEN) G[2]);
! 3318: }
! 3319: res = cgetg(n + 1, t_VEC);
! 3320: res[1] = (long) permidentity(n);
! 3321: k = 1;
! 3322: for (i = 1; i < lg(G[1]); i++)
! 3323: {
! 3324: int c = k * (((long **) G)[2][i] - 1);
! 3325: for (j = 1; j <= c; j++) /* I like it */
! 3326: res[++k] = (long) permapply((GEN) res[j], ((GEN **) G)[1][i]);
! 3327: }
! 3328: if (flag)
! 3329: {
! 3330: grp[6] = (long) res;
! 3331: return gerepileupto(ltop, grp);
! 3332: }
! 3333: aut = galoisgrouptopol(res,L,M,den,gb.ladicsol, varn(T));
! 3334: if (DEBUGLEVEL >= 1)
! 3335: msgtimer("Calcul polynomes");
! 3336: return gerepileupto(ltop, aut);
! 3337: }
! 3338:
! 3339: /* Calcule le nombre d'automorphisme de T avec forte probabilité */
! 3340: /* pdepart premier nombre premier a tester */
! 3341: long
! 3342: numberofconjugates(GEN T, long pdepart)
! 3343: {
! 3344: ulong ltop = avma, ltop2;
! 3345: long n, p, nbmax, nbtest;
! 3346: long card;
! 3347: byteptr primepointer;
! 3348: int i;
! 3349: GEN L;
! 3350: n = degpol(T);
! 3351: card = sturm(T);
! 3352: card = cgcd(card, n - card);
! 3353: nbmax = (n<<1) + 1;
! 3354: if (nbmax < 20) nbmax=20;
! 3355: nbtest = 0;
! 3356: L = cgetg(n + 1, t_VECSMALL);
! 3357: ltop2 = avma;
! 3358: for (p = 0, primepointer = diffptr; nbtest < nbmax && card > 1;)
! 3359: {
! 3360: long s, c;
! 3361: long isram;
! 3362: GEN S;
! 3363: c = *primepointer++;
! 3364: if (!c)
! 3365: err(primer1);
! 3366: p += c;
! 3367: if (p < pdepart)
! 3368: continue;
! 3369: S = (GEN) simplefactmod(T, stoi(p));
! 3370: isram = 0;
! 3371: for (i = 1; i < lg(S[2]) && !isram; i++)
! 3372: if (!gcmp1(((GEN **) S)[2][i]))
! 3373: isram = 1;
! 3374: if (isram == 0)
! 3375: {
! 3376: for (i = 1; i <= n; i++)
! 3377: L[i] = 0;
! 3378: for (i = 1; i < lg(S[1]) && !isram; i++)
! 3379: L[itos(((GEN **) S)[1][i])]++;
! 3380: s = L[1];
! 3381: for (i = 2; i <= n; i++)
! 3382: s = cgcd(s, L[i] * i);
! 3383: card = cgcd(s, card);
! 3384: }
! 3385: if (DEBUGLEVEL >= 6)
! 3386: fprintferr("NumberOfConjugates:Nbtest=%ld,card=%ld,p=%ld\n", nbtest,
! 3387: card, p);
! 3388: nbtest++;
! 3389: avma = ltop2;
! 3390: }
! 3391: if (DEBUGLEVEL >= 2)
! 3392: fprintferr("NumberOfConjugates:card=%ld,p=%ld\n", card, p);
! 3393: avma = ltop;
! 3394: return card;
! 3395: }
! 3396:
! 3397: GEN
! 3398: galoisconj0(GEN nf, long flag, GEN d, long prec)
! 3399: {
! 3400: ulong ltop;
! 3401: GEN G, T;
! 3402: long card;
! 3403: if (typ(nf) != t_POL)
! 3404: {
! 3405: nf = checknf(nf);
! 3406: T = (GEN) nf[1];
! 3407: }
! 3408: else
! 3409: T = nf;
! 3410: switch (flag)
! 3411: {
! 3412: case 0:
! 3413: ltop = avma;
! 3414: G = galoisconj4(nf, d, 0, 0);
! 3415: if (typ(G) != t_INT) /* Success */
! 3416: return G;
! 3417: else
! 3418: {
! 3419: card = numberofconjugates(T, G == gzero ? 2 : itos(G));
! 3420: avma = ltop;
! 3421: if (card != 1)
! 3422: {
! 3423: if (typ(nf) == t_POL)
! 3424: {
! 3425: G = galoisconj2pol(nf, card, prec);
! 3426: if (lg(G) <= card)
! 3427: err(warner, "conjugates list may be incomplete in nfgaloisconj");
! 3428: return G;
! 3429: }
! 3430: else
! 3431: return galoisconj(nf);
! 3432: }
! 3433: }
! 3434: break; /* Failure */
! 3435: case 1:
! 3436: return galoisconj(nf);
! 3437: case 2:
! 3438: return galoisconj2(nf, degpol(T), prec);
! 3439: case 4:
! 3440: G = galoisconj4(nf, d, 0, 0);
! 3441: if (typ(G) != t_INT)
! 3442: return G;
! 3443: break; /* Failure */
! 3444: default:
! 3445: err(flagerr, "nfgaloisconj");
! 3446: }
! 3447: G = cgetg(2, t_COL);
! 3448: G[1] = (long) polx[varn(T)];
! 3449: return G; /* Failure */
! 3450: }
! 3451:
! 3452:
! 3453:
! 3454: /******************************************************************************/
! 3455: /* Isomorphism between number fields */
! 3456: /******************************************************************************/
! 3457: long
! 3458: isomborne(GEN P, GEN den, GEN p)
! 3459: {
! 3460: ulong ltop=avma;
! 3461: struct galois_borne gb;
! 3462: gb.l=p;
! 3463: galoisborne(P,den,&gb,degpol(P));
! 3464: avma=ltop;
! 3465: return gb.valsol;
! 3466: }
! 3467:
! 3468:
! 3469:
! 3470: /******************************************************************************/
! 3471: /* Galois theory related algorithms */
! 3472: /******************************************************************************/
! 3473: GEN
! 3474: checkgal(GEN gal)
! 3475: {
! 3476: if (typ(gal) == t_POL)
! 3477: err(talker, "please apply galoisinit first");
! 3478: if (typ(gal) != t_VEC || lg(gal) != 9)
! 3479: err(talker, "Not a Galois field in a Galois related function");
! 3480: return gal;
! 3481: }
! 3482:
! 3483: GEN
! 3484: galoisinit(GEN nf, GEN den, long karma)
! 3485: {
! 3486: GEN G;
! 3487: G = galoisconj4(nf, den, 1, karma);
! 3488: if (typ(G) == t_INT)
! 3489: err(talker,
! 3490: "galoisinit: field not Galois or Galois group not weakly super solvable");
! 3491: return G;
! 3492: }
! 3493:
! 3494: GEN
! 3495: galoispermtopol(GEN gal, GEN perm)
! 3496: {
! 3497: GEN v;
! 3498: long t = typ(perm);
! 3499: int i;
! 3500: gal = checkgal(gal);
! 3501: switch (t)
! 3502: {
! 3503: case t_VECSMALL:
! 3504: return permtopol(perm, (GEN) gal[3], (GEN) gal[4], (GEN) gal[5],
! 3505: gmael(gal,2,3), varn((GEN) gal[1]));
! 3506: case t_VEC:
! 3507: case t_COL:
! 3508: case t_MAT:
! 3509: v = cgetg(lg(perm), t);
! 3510: for (i = 1; i < lg(v); i++)
! 3511: v[i] = (long) galoispermtopol(gal, (GEN) perm[i]);
! 3512: return v;
! 3513: }
! 3514: err(typeer, "galoispermtopol");
! 3515: return NULL; /* not reached */
! 3516: }
! 3517:
! 3518:
! 3519: GEN
! 3520: galoiscosets(GEN O, GEN perm)
! 3521: {
! 3522: long i,j,k,u;
! 3523: long o = lg(O)-1, f = lg(O[1])-1;
! 3524: GEN C = cgetg(lg(O),t_VECSMALL);
! 3525: ulong av=avma;
! 3526: GEN RC=cgetg(lg(perm),t_VECSMALL);
! 3527: for(i=1;i<lg(RC);i++)
! 3528: RC[i]=0;
! 3529: u=mael(O,1,1);
! 3530: for(i=1,j=1;j<=o;i++)
! 3531: {
! 3532: if (RC[mael(perm,i,u)])
! 3533: continue;
! 3534: for(k=1;k<=f;k++)
! 3535: RC[mael(perm,i,mael(O,1,k))]=1;
! 3536: C[j++]=i;
! 3537: }
! 3538: avma=av;
! 3539: return C;
! 3540: }
! 3541:
! 3542: GEN
! 3543: fixedfieldfactor(GEN L, GEN O, GEN perm, GEN M, GEN den, GEN mod,
! 3544: long x,long y)
! 3545: {
! 3546: ulong ltop=avma;
! 3547: GEN F,G,V,res,cosets;
! 3548: int i, j, k;
! 3549: F=cgetg(lg(O[1])+1,t_COL);
! 3550: F[lg(O[1])]=un;
! 3551: G=cgetg(lg(O),t_VEC);
! 3552: for (i = 1; i < lg(O); i++)
! 3553: {
! 3554: GEN Li = cgetg(lg(O[i]),t_VEC);
! 3555: for (j = 1; j < lg(O[i]); j++)
! 3556: Li[j] = L[mael(O,i,j)];
! 3557: G[i]=(long)FpV_roots_to_pol(Li,mod,x);
! 3558: }
! 3559:
! 3560: cosets=galoiscosets(O,perm);
! 3561: if (DEBUGLEVEL>=4) fprintferr("GaloisFixedField:cosets=%Z \n",cosets);
! 3562: V=cgetg(lg(O),t_COL);
! 3563: if (DEBUGLEVEL>=6) fprintferr("GaloisFixedField:den=%Z mod=%Z \n",den,mod);
! 3564: res=cgetg(lg(O),t_VEC);
! 3565: for (i = 1; i < lg(O); i++)
! 3566: {
! 3567: ulong av=avma;
! 3568: long ii,jj;
! 3569: GEN G=cgetg(lg(O),t_VEC);
! 3570: for (ii = 1; ii < lg(O); ii++)
! 3571: {
! 3572: GEN Li = cgetg(lg(O[ii]),t_VEC);
! 3573: for (jj = 1; jj < lg(O[ii]); jj++)
! 3574: Li[jj] = L[mael(perm,cosets[i],mael(O,ii,jj))];
! 3575: G[ii]=(long)FpV_roots_to_pol(Li,mod,x);
! 3576: }
! 3577: for (j = 1; j < lg(O[1]); j++)
! 3578: {
! 3579: for(k = 1; k < lg(O); k++)
! 3580: V[k]=mael(G,k,j+1);
! 3581: F[j]=(long) vectopol(V, M, den, mod, y);
! 3582: }
! 3583: res[i]=(long)gerepileupto(av,gtopolyrev(F,x));
! 3584: }
! 3585: return gerepileupto(ltop,res);
! 3586: }
! 3587:
! 3588: GEN
! 3589: galoisfixedfield(GEN gal, GEN perm, long flag, long y)
! 3590: {
! 3591: ulong ltop = avma, lbot;
! 3592: GEN L, P, S, PL, O, res, mod;
! 3593: long x, n;
! 3594: int i;
! 3595: gal = checkgal(gal);
! 3596: x = varn((GEN) gal[1]);
! 3597: L = (GEN) gal[3]; n=lg(L)-1;
! 3598: mod = gmael(gal,2,3);
! 3599: if (flag<0 || flag>2)
! 3600: err(flagerr, "galoisfixedfield");
! 3601: if (typ(perm) == t_VEC)
! 3602: {
! 3603: if (lg(perm)==1)
! 3604: perm=permidentity(n);
! 3605: else
! 3606: for (i = 1; i < lg(perm); i++)
! 3607: if (typ(perm[i]) != t_VECSMALL || lg(perm[i])!=n+1)
! 3608: err(typeer, "galoisfixedfield");
! 3609: }
! 3610: else if (typ(perm) != t_VECSMALL || lg(perm)!=n+1 )
! 3611: err(typeer, "galoisfixedfield");
! 3612: O = permorbite(perm);
! 3613: {
! 3614: GEN sym=cgetg(lg(L),t_VECSMALL);
! 3615: GEN dg=cgetg(lg(L),t_VECSMALL);
! 3616: GEN V;
! 3617: V = fixedfieldsympol(O, L, mod, gmael(gal,2,1), gun, sym, dg, x);
! 3618: P=(GEN)V[2];
! 3619: PL=(GEN)V[1];
! 3620: }
! 3621: if (flag==1)
! 3622: return gerepileupto(ltop,P);
! 3623: S = fixedfieldinclusion(O, PL);
! 3624: S = vectopol(S, (GEN) gal[4], (GEN) gal[5], mod, x);
! 3625: if (flag==0)
! 3626: {
! 3627: lbot = avma;
! 3628: res = cgetg(3, t_VEC);
! 3629: res[1] = (long) gcopy(P);
! 3630: res[2] = (long) gmodulcp(S, (GEN) gal[1]);
! 3631: return gerepile(ltop, lbot, res);
! 3632: }
! 3633: else
! 3634: {
! 3635: GEN PM,Pden;
! 3636: {
! 3637: struct galois_borne Pgb;
! 3638: long val=itos(gmael(gal,2,2));
! 3639: Pgb.l = gmael(gal,2,1);
! 3640: Pden = galoisborne(P, NULL, &Pgb, 1);
! 3641: if (Pgb.valabs > val)
! 3642: {
! 3643: if (DEBUGLEVEL>=4)
! 3644: fprintferr("GaloisConj:increase prec of p-adic roots of %ld.\n"
! 3645: ,Pgb.valabs-val);
! 3646: PL = rootpadicliftroots(P,PL,Pgb.l,Pgb.valabs);
! 3647: L = rootpadicliftroots((GEN) gal[1],L,Pgb.l,Pgb.valabs);
! 3648: mod = Pgb.ladicabs;
! 3649: }
! 3650: }
! 3651: PM = vandermondeinversemod(PL, P, Pden, mod);
! 3652: lbot = avma;
! 3653: if (y==-1)
! 3654: y = fetch_user_var("y");
! 3655: if (y<=x)
! 3656: err(talker,"priority of optional variable too high in galoisfixedfield");
! 3657: res = cgetg(4, t_VEC);
! 3658: res[1] = (long) gcopy(P);
! 3659: res[2] = (long) gmodulcp(S, (GEN) gal[1]);
! 3660: res[3] = (long) fixedfieldfactor(L,O,(GEN)gal[6],
! 3661: PM,Pden,mod,x,y);
! 3662: return gerepile(ltop, lbot, res);
! 3663: }
! 3664: }
! 3665:
! 3666: /*return 1 if gal is abelian, else 0*/
! 3667: long
! 3668: galoistestabelian(GEN gal)
! 3669: {
! 3670: ulong ltop=avma;
! 3671: long i,j;
! 3672: gal = checkgal(gal);
! 3673: for(i=2;i<lg(gal[7]);i++)
! 3674: for(j=1;j<i;j++)
! 3675: {
! 3676: long test=gegal(permapply(gmael(gal,7,i),gmael(gal,7,j)),
! 3677: permapply(gmael(gal,7,j),gmael(gal,7,i)));
! 3678: avma=ltop;
! 3679: if (!test) return 0;
! 3680: }
! 3681: return 1;
! 3682: }
! 3683:
! 3684: GEN galoisisabelian(GEN gal, long flag)
! 3685: {
! 3686: long i, j;
! 3687: long test, n=lg(gal[7]);
! 3688: GEN M;
! 3689: gal = checkgal(gal);
! 3690: test=galoistestabelian(gal);
! 3691: if (!test) return gzero;
! 3692: if (flag==1) return gun;
! 3693: if (flag) err(flagerr,"galoisisabelian");
! 3694: M=cgetg(n,t_MAT);
! 3695: for(i=1;i<n;i++)
! 3696: {
! 3697: ulong btop;
! 3698: GEN P;
! 3699: long k;
! 3700: M[i]=lgetg(n,t_COL);
! 3701: btop=avma;
! 3702: P=permcyclepow(permorbite(gmael(gal,7,i)),mael(gal,8,i));
! 3703: for(j=1;j<lg(gal[6]);j++)
! 3704: if (gegal(P,gmael(gal,6,j)))
! 3705: break;
! 3706: avma=btop;
! 3707: if (j==lg(gal[6])) err(talker,"wrong argument in galoisisabelian");
! 3708: j--;
! 3709: for(k=1;k<i;k++)
! 3710: {
! 3711: mael(M,i,k)=lstoi(j%mael(gal,8,k));
! 3712: j/=mael(gal,8,k);
! 3713: }
! 3714: mael(M,i,k++)=lstoi(mael(gal,8,i));
! 3715: for( ;k<n;k++)
! 3716: mael(M,i,k)=zero;
! 3717: }
! 3718: return M;
! 3719: }
! 3720: /* Calcule les orbites d'un sous-groupe de Z/nZ donne par un
! 3721: * generateur ou d'un ensemble de generateur donne par un vecteur.
! 3722: */
! 3723: GEN
! 3724: subgroupcoset(long n, GEN v)
! 3725: {
! 3726: ulong ltop = avma, lbot;
! 3727: int i, j, k = 1, l, m, o, p, flag;
! 3728: GEN bit, cycle, cy;
! 3729: cycle = cgetg(n, t_VEC);
! 3730: bit = cgetg(n, t_VECSMALL);
! 3731: for (i = 1; i < n; i++)
! 3732: if (cgcd(i,n)==1)
! 3733: bit[i] = 0;
! 3734: else
! 3735: {
! 3736: bit[i] = -1;
! 3737: k++;
! 3738: }
! 3739: for (l = 1; k < n;)
! 3740: {
! 3741: for (j = 1; bit[j]; j++);
! 3742: cy = cgetg(n, t_VECSMALL);
! 3743: m = 1;
! 3744: cy[m++] = j;
! 3745: bit[j] = 1;
! 3746: k++;
! 3747: do
! 3748: {
! 3749: flag = 0;
! 3750: for (o = 1; o < lg(v); o++)
! 3751: {
! 3752: for (p = 1; p < m; p++) /* m varie! */
! 3753: {
! 3754: j = mulssmod(v[o],cy[p],n);
! 3755: if (!bit[j])
! 3756: {
! 3757: flag = 1;
! 3758: bit[j] = 1;
! 3759: k++;
! 3760: cy[m++] = j;
! 3761: }
! 3762: }
! 3763: }
! 3764: }
! 3765: while (flag);
! 3766: setlg(cy, m);
! 3767: cycle[l++] = (long) cy;
! 3768: }
! 3769: setlg(cycle, l);
! 3770: lbot = avma;
! 3771: cycle = gcopy(cycle);
! 3772: return gerepile(ltop, lbot, cycle);
! 3773: }
! 3774: /* Calcule les elements d'un sous-groupe H de Z/nZ donne par un
! 3775: * generateur ou d'un ensemble de generateur donne par un vecteur (v).
! 3776: *
! 3777: * cy liste des elements VECSMALL de longueur au moins card H.
! 3778: * bit bitmap des elements VECSMALL de longueur au moins n.
! 3779: * retourne le nombre d'elements+1
! 3780: */
! 3781: long
! 3782: sousgroupeelem(long n, GEN v, GEN cy, GEN bit)
! 3783: {
! 3784: int j, m, o, p, flag;
! 3785: for(j=1;j<n;j++)
! 3786: bit[j]=0;
! 3787: m = 1;
! 3788: bit[m] = 1;
! 3789: cy[m++] = 1;
! 3790: do
! 3791: {
! 3792: flag = 0;
! 3793: for (o = 1; o < lg(v); o++)
! 3794: {
! 3795: for (p = 1; p < m; p++) /* m varie! */
! 3796: {
! 3797: j = mulssmod(v[o],cy[p],n);
! 3798: if (!bit[j])
! 3799: {
! 3800: flag = 1;
! 3801: bit[j] = 1;
! 3802: cy[m++] = j;
! 3803: }
! 3804: }
! 3805: }
! 3806: }
! 3807: while (flag);
! 3808: return m;
! 3809: }
! 3810: /* n,v comme precedemment.
! 3811: * Calcule le conducteur et retourne le nouveau groupe de congruence dans V
! 3812: * V doit etre un t_VECSMALL de taille n+1 au moins.
! 3813: */
! 3814: long znconductor(long n, GEN v, GEN V)
! 3815: {
! 3816: ulong ltop;
! 3817: int i,j;
! 3818: long m;
! 3819: GEN F,W;
! 3820: W = cgetg(n, t_VECSMALL);
! 3821: ltop=avma;
! 3822: m = sousgroupeelem(n,v,V,W);
! 3823: setlg(V,m);
! 3824: if (DEBUGLEVEL>=6)
! 3825: fprintferr("SubCyclo:elements:%Z\n",V);
! 3826: F = factor(stoi(n));
! 3827: for(i=lg((GEN)F[1])-1;i>0;i--)
! 3828: {
! 3829: long p,e,q;
! 3830: p=itos(gcoeff(F,i,1));
! 3831: e=itos(gcoeff(F,i,2));
! 3832: if (DEBUGLEVEL>=4)
! 3833: fprintferr("SubCyclo:testing %ld^%ld\n",p,e);
! 3834: while (e>=1)
! 3835: {
! 3836: int z = 1;
! 3837: q=n/p;
! 3838: for(j=1;j<p;j++)
! 3839: {
! 3840: z += q;
! 3841: if (!W[z] && z%p) break;
! 3842: }
! 3843: if (j<p)
! 3844: {
! 3845: if (DEBUGLEVEL>=4)
! 3846: fprintferr("SubCyclo:%ld not found\n",z);
! 3847: break;
! 3848: }
! 3849: e--;
! 3850: n=q;
! 3851: if (DEBUGLEVEL>=4)
! 3852: fprintferr("SubCyclo:new conductor:%ld\n",n);
! 3853: m=sousgroupeelem(n,v,V,W);
! 3854: setlg(V,m);
! 3855: if (DEBUGLEVEL>=6)
! 3856: fprintferr("SubCyclo:elements:%Z\n",V);
! 3857: }
! 3858: }
! 3859: if (DEBUGLEVEL>=6)
! 3860: fprintferr("SubCyclo:conductor:%ld\n",n);
! 3861: avma=ltop;
! 3862: return n;
! 3863: }
! 3864:
! 3865: static GEN gscycloconductor(GEN g, long n, long flag)
! 3866: {
! 3867: if (flag==2)
! 3868: {
! 3869: GEN V=cgetg(3,t_VEC);
! 3870: V[1]=lcopy(g);
! 3871: V[2]=lstoi(n);
! 3872: return V;
! 3873: }
! 3874: return g;
! 3875: }
! 3876: static GEN
! 3877: lift_check_modulus(GEN H, GEN n)
! 3878: {
! 3879: long t=typ(H), l=lg(H);
! 3880: long i;
! 3881: GEN V;
! 3882: switch(t)
! 3883: {
! 3884: case t_INTMOD:
! 3885: if (cmpii(n,(GEN)H[1]))
! 3886: err(talker,"wrong modulus in galoissubcyclo");
! 3887: H = (GEN)H[2];
! 3888: case t_INT:
! 3889: if (!is_pm1(mppgcd(H,n)))
! 3890: err(talker,"generators must be prime to conductor in galoissubcyclo");
! 3891: return modii(H,n);
! 3892: case t_VEC: case t_COL:
! 3893: V=cgetg(l,t);
! 3894: for(i=1;i<l;i++)
! 3895: V[i] = (long)lift_check_modulus((GEN)H[i],n);
! 3896: return V;
! 3897: case t_VECSMALL:
! 3898: return H;
! 3899: }
! 3900: err(talker,"wrong type in galoissubcyclo");
! 3901: return NULL;/*not reached*/
! 3902: }
! 3903:
! 3904: GEN
! 3905: galoissubcyclo(long n, GEN H, GEN Z, long v, long flag)
! 3906: {
! 3907: ulong ltop=avma,av;
! 3908: GEN l,borne,le,powz,z,V;
! 3909: long i;
! 3910: long e,val;
! 3911: long u,o,j;
! 3912: GEN O,g;
! 3913: if (flag<0 || flag>2) err(flagerr,"galoisubcyclo");
! 3914: if ( v==-1 ) v=0;
! 3915: if ( n<1 ) err(arither2);
! 3916: if ( n>=VERYBIGINT)
! 3917: err(impl,"galoissubcyclo for huge conductor");
! 3918: if ( typ(H)==t_MAT )
! 3919: {
! 3920: GEN zn2, zn3, gen, ord;
! 3921: if (lg(H) == 1 || lg(H) != lg(H[1]))
! 3922: err(talker,"not a HNF matrix in galoissubcyclo");
! 3923: if (!Z)
! 3924: Z=znstar(stoi(n));
! 3925: else if (typ(Z)!=t_VEC || lg(Z)!=4)
! 3926: err(talker,"Optionnal parameter must be as output by znstar in galoissubcyclo");
! 3927: zn2 = gtovecsmall((GEN)Z[2]);
! 3928: zn3 = lift((GEN)Z[3]);
! 3929: if ( lg(zn2) != lg(H) || lg(zn3) != lg(H))
! 3930: err(talker,"Matrix of wrong dimensions in galoissubcyclo");
! 3931: gen = cgetg(lg(zn3), t_VECSMALL);
! 3932: ord = cgetg(lg(zn3), t_VECSMALL);
! 3933: hnftogeneratorslist(n,zn2,zn3,H,gen,ord);
! 3934: H=gen;
! 3935: }
! 3936: else
! 3937: {
! 3938: H=lift_check_modulus(H,stoi(n));
! 3939: H=gtovecsmall(H);
! 3940: for (i=1;i<lg(H);i++)
! 3941: if (H[i]<0)
! 3942: H[i]=mulssmod(-H[i],n-1,n);
! 3943: /*Should check components are prime to n, but it is costly*/
! 3944: }
! 3945: V = cgetg(n, t_VECSMALL);
! 3946: if (DEBUGLEVEL >= 1)
! 3947: timer2();
! 3948: n = znconductor(n,H,V);
! 3949: if (flag==1) {avma=ltop;return stoi(n);}
! 3950: if (DEBUGLEVEL >= 1)
! 3951: msgtimer("znconductor.");
! 3952: H = V;
! 3953: O = subgroupcoset(n,H);
! 3954: if (DEBUGLEVEL >= 1)
! 3955: msgtimer("subgroupcoset.");
! 3956: if (DEBUGLEVEL >= 6)
! 3957: fprintferr("Subcyclo: orbit=%Z\n",O);
! 3958: if (lg(O)==1 || lg(O[1])==2)
! 3959: {
! 3960: avma=ltop;
! 3961: return gscycloconductor(cyclo(n,v),n,flag);
! 3962: }
! 3963: u=lg(O)-1;o=lg(O[1])-1;
! 3964: if (DEBUGLEVEL >= 4)
! 3965: fprintferr("Subcyclo: %ld orbits with %ld elements each\n",u,o);
! 3966: l=stoi(n+1);e=1;
! 3967: while(!isprime(l))
! 3968: {
! 3969: l=addis(l,n);
! 3970: e++;
! 3971: }
! 3972: if (DEBUGLEVEL >= 4)
! 3973: fprintferr("Subcyclo: prime l=%Z\n",l);
! 3974: av=avma;
! 3975: /*Borne utilise':
! 3976: Vecmax(Vec((x+o)^u)=max{binome(u,i)*o^i ;1<=i<=u}
! 3977: */
! 3978: i=u-(1+u)/(1+o);
! 3979: borne=mulii(binome(stoi(u),i),gpowgs(stoi(o),i));
! 3980: if (DEBUGLEVEL >= 4)
! 3981: fprintferr("Subcyclo: borne=%Z\n",borne);
! 3982: val=mylogint(shifti(borne,1),l);
! 3983: avma=av;
! 3984: if (DEBUGLEVEL >= 4)
! 3985: fprintferr("Subcyclo: val=%ld\n",val);
! 3986: le=gpowgs(l,val);
! 3987: z=lift(gpowgs(gener(l),e));
! 3988: z=padicsqrtnlift(gun,stoi(n),z,l,val);
! 3989: if (DEBUGLEVEL >= 1)
! 3990: msgtimer("padicsqrtnlift.");
! 3991: powz = cgetg(n,t_VEC); powz[1] = (long)z;
! 3992: for (i=2; i<n; i++) powz[i] = lmodii(mulii(z,(GEN)powz[i-1]),le);
! 3993: if (DEBUGLEVEL >= 1)
! 3994: msgtimer("computing roots.");
! 3995: g=cgetg(u+1,t_VEC);
! 3996: for(i=1;i<=u;i++)
! 3997: {
! 3998: GEN s;
! 3999: s=gzero;
! 4000: for(j=1;j<=o;j++)
! 4001: s=addii(s,(GEN)powz[mael(O,i,j)]);
! 4002: g[i]=(long)modii(negi(s),le);
! 4003: }
! 4004: if (DEBUGLEVEL >= 1)
! 4005: msgtimer("computing new roots.");
! 4006: g=FpV_roots_to_pol(g,le,v);
! 4007: if (DEBUGLEVEL >= 1)
! 4008: msgtimer("computing products.");
! 4009: g=FpX_center(g,le);
! 4010: return gerepileupto(ltop,gscycloconductor(g,n,flag));
! 4011: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to galconj.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 / basemath