[BACK]Return to galconj.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / pari / src / basemath

Annotation of OpenXM_contrib/pari/src/basemath/galconj.c, Revision 1.1.1.1

1.1       maekawa     1: /*************************************************************************/
                      2: /**                                                                    **/
                      3: /**                           GALOIS CONJUGATES                                **/
                      4: /**                                                                    **/
                      5: /*************************************************************************/
                      6: /* $Id: galconj.c,v 1.8 1999/09/23 17:50:56 karim Exp $ */
                      7: #include "pari.h"
                      8: GEN
                      9: galoisconj(GEN nf)
                     10: {
                     11:   GEN x, y, z;
                     12:   long i, lz, lx, v, av = avma;
                     13:   nf = checknf(nf);
                     14:   x = (GEN) nf[1];
                     15:   v = varn(x);
                     16:   lx = lgef(x) - 2;
                     17:   if (v == 0)
                     18:     nf = gsubst(nf, 0, polx[MAXVARN]);
                     19:   else
                     20:   {
                     21:     x = dummycopy(x);
                     22:     setvarn(x, 0);
                     23:   }
                     24:   z = nfroots(nf, x);
                     25:   lz = lg(z);
                     26:   y = cgetg(lz, t_COL);
                     27:   for (i = 1; i < lz; i++)
                     28:   {
                     29:     GEN p1 = lift((GEN) z[i]);
                     30:     setvarn(p1, v);
                     31:     y[i] = (long) p1;
                     32:   }
                     33:   return gerepileupto(av, y);
                     34: }
                     35: /* nbmax: maximum number of possible conjugates */
                     36: GEN
                     37: galoisconj2pol(GEN x, long nbmax, long prec)
                     38: {
                     39:   long av = avma, i, n, v, nbauto;
                     40:   GEN y, w, polr, p1, p2;
                     41:   n = lgef(x) - 3;
                     42:   if (n <= 0)
                     43:     return cgetg(1, t_VEC);
                     44:   if (gisirreducible(x) == gzero)
                     45:     err(redpoler, "galoisconj2pol");
                     46:   polr = roots(x, prec);
                     47:   p1 = (GEN) polr[1];
                     48:   nbauto = 1;
                     49:   prec = (long) (bit_accuracy(prec) * L2SL10 * 0.75);
                     50:   w = cgetg(n + 2, t_VEC);
                     51:   w[1] = un;
                     52:   for (i = 2; i <= n; i++)
                     53:     w[i] = lmul(p1, (GEN) w[i - 1]);
                     54:   v = varn(x);
                     55:   y = cgetg(nbmax + 1, t_COL);
                     56:   y[1] = (long) polx[v];
                     57:   for (i = 2; i <= n && nbauto < nbmax; i++)
                     58:   {
                     59:     w[n + 1] = polr[i];
                     60:     p1 = lindep2(w, prec);
                     61:     if (signe(p1[n + 1]))
                     62:     {
                     63:       setlg(p1, n + 1);
                     64:       p2 = gdiv(gtopolyrev(p1, v), negi((GEN) p1[n + 1]));
                     65:       if (gdivise(poleval(x, p2), x))
                     66:       {
                     67:        y[++nbauto] = (long) p2;
                     68:        if (DEBUGLEVEL > 1)
                     69:          fprintferr("conjugate %ld: %Z\n", i, y[nbauto]);
                     70:       }
                     71:     }
                     72:   }
                     73:   setlg(y, 1 + nbauto);
                     74:   return gerepileupto(av, gen_sort(y, 0, cmp_pol));
                     75: }
                     76: GEN
                     77: galoisconj2(GEN nf, long nbmax, long prec)
                     78: {
                     79:   long av = avma, i, j, n, r1, ru, nbauto;
                     80:   GEN x, y, w, polr, p1, p2;
                     81:   if (typ(nf) == t_POL)
                     82:     return galoisconj2pol(nf, nbmax, prec);
                     83:   nf = checknf(nf);
                     84:   x = (GEN) nf[1];
                     85:   n = lgef(x) - 3;
                     86:   if (n <= 0)
                     87:     return cgetg(1, t_VEC);
                     88:   r1 = itos(gmael(nf, 2, 1));
                     89:   p1 = (GEN) nf[6];
                     90:   prec = precision((GEN) p1[1]);
                     91:   /* accuracy in decimal digits */
                     92:   prec = (long) (bit_accuracy(prec) * L2SL10 * 0.75);
                     93:   ru = (n + r1) >> 1;
                     94:   nbauto = 1;
                     95:   polr = cgetg(n + 1, t_VEC);
                     96:   for (i = 1; i <= r1; i++)
                     97:     polr[i] = p1[i];
                     98:   for (j = i; i <= ru; i++)
                     99:   {
                    100:     polr[j++] = p1[i];
                    101:     polr[j++] = lconj((GEN) p1[i]);
                    102:   }
                    103:   p2 = gmael(nf, 5, 1);
                    104:   w = cgetg(n + 2, t_VEC);
                    105:   for (i = 1; i <= n; i++)
                    106:     w[i] = coeff(p2, 1, i);
                    107:   y = cgetg(nbmax + 1, t_COL);
                    108:   y[1] = (long) polx[varn(x)];
                    109:   for (i = 2; i <= n && nbauto < nbmax; i++)
                    110:   {
                    111:     w[n + 1] = polr[i];
                    112:     p1 = lindep2(w, prec);
                    113:     if (signe(p1[n + 1]))
                    114:     {
                    115:       setlg(p1, n + 1);
                    116:       settyp(p1, t_COL);
                    117:       p2 = gdiv(gmul((GEN) nf[7], p1), negi((GEN) p1[n + 1]));
                    118:       if (gdivise(poleval(x, p2), x))
                    119:       {
                    120:        y[++nbauto] = (long) p2;
                    121:        if (DEBUGLEVEL > 1)
                    122:          fprintferr("conjugate %ld: %Z\n", i, y[nbauto]);
                    123:       }
                    124:     }
                    125:   }
                    126:   setlg(y, 1 + nbauto);
                    127:   return gerepileupto(av, gen_sort(y, 0, cmp_pol));
                    128: }
                    129: /*************************************************************************/
                    130: /**                                                                    **/
                    131: /**                           GALOIS4                                  **/
                    132: /**                                                                    **/
                    133: /**                                                                     **/
                    134: /**                                                                     **/
                    135: /*************************************************************************/
                    136: struct galois_lift
                    137: {
                    138:   GEN T;
                    139:   GEN den;
                    140:   GEN p;
                    141:   GEN borne;
                    142:   GEN L;
                    143:   long e;
                    144:   GEN Q;
                    145: };
                    146: /* Initialise la structure galois_lift */
                    147: void
                    148: initlift(GEN T, GEN den, GEN p, GEN borne, GEN L, struct galois_lift * gl)
                    149: {
                    150:   ulong ltop, lbot;
                    151:   gl->T = T;
                    152:   gl->den = den;
                    153:   gl->p = p;
                    154:   gl->borne = borne;
                    155:   gl->L = L;
                    156:   ltop = avma;
                    157:   gl->e = itos(gmax(gdeux, gceil(gdiv(glog(gmul2n(borne, 1), DEFAULTPREC),
                    158:                                      glog(p, DEFAULTPREC)))));
                    159:   lbot = avma;
                    160:   gl->Q = gpowgs(p, gl->e);
                    161:   gl->Q = gerepile(ltop, lbot, gl->Q);
                    162: }
                    163: /*
                    164:  * T est le polynome \sum t_i*X^i S est un Polmod T Retourne un vecteur Spow
                    165:  * verifiant la condition: i>=1 et t_i!=0 ==> Spow[i]=S^(i-1)*t_i Essaie
                    166:  * d'etre efficace sur les polynomes lacunaires
                    167:  */
                    168: GEN
                    169: compoTS(GEN T, GEN S)
                    170: {
                    171:   GEN Spow;
                    172:   int i;
                    173:   Spow = cgetg(lgef(T) - 2, t_VEC);
                    174:   for (i = 3; i < lg(Spow); i++)
                    175:     Spow[i] = 0;
                    176:   Spow[1] = un;
                    177:   Spow[2] = (long) S;
                    178:   for (i = 2; i < lg(Spow) - 1; i++)
                    179:   {
                    180:     if (!gcmp0((GEN) T[i + 3]))
                    181:       for (;;)
                    182:       {
                    183:        int k, l;
                    184:        for (k = 1; k <= (i >> 1); k++)
                    185:          if (Spow[k + 1] && Spow[i - k + 1])
                    186:            break;
                    187:        if ((k << 1) < i)
                    188:        {
                    189:          Spow[i + 1] = lmul((GEN) Spow[k + 1], (GEN) Spow[i - k + 1]);
                    190:          break;
                    191:        } else if ((k << 1) == i)       /* En esperant que sqr est plus
                    192:                                         * rapide... */
                    193:        {
                    194:          Spow[i + 1] = lsqr((GEN) Spow[k + 1]);
                    195:          break;
                    196:        }
                    197:        for (k = i - 1; k > 0; k--)
                    198:          if (Spow[k + 1])
                    199:            break;
                    200:        if ((k << 1) < i)
                    201:        {
                    202:          Spow[(k << 1) + 1] = lsqr((GEN) Spow[k + 1]);
                    203:          continue;
                    204:        }
                    205:        for (l = i - k; l > 0; l--)
                    206:          if (Spow[l + 1])
                    207:            break;
                    208:        if (Spow[i - l - k + 1])
                    209:          Spow[i - k + 1] = lmul((GEN) Spow[i - l - k + 1], (GEN) Spow[l + 1]);
                    210:        else
                    211:          Spow[l + k + 1] = lmul((GEN) Spow[k + 1], (GEN) Spow[l + 1]);
                    212:       }
                    213:   }
                    214:   for (i = 1; i < lg(Spow); i++)
                    215:     if (!gcmp0((GEN) T[i + 2]))
                    216:       Spow[i] = lmul((GEN) Spow[i], (GEN) T[i + 2]);
                    217:   return Spow;
                    218: }
                    219: /*
                    220:  * Calcule T(S) a l'aide du vecteur Spow
                    221:  */
                    222: static GEN
                    223: calcTS(GEN Spow, GEN T, GEN S)
                    224: {
                    225:   GEN res = gzero;
                    226:   int i;
                    227:   for (i = 1; i < lg(Spow); i++)
                    228:     if (!gcmp0((GEN) T[i + 2]))
                    229:       res = gadd(res, (GEN) Spow[i]);
                    230:   res = gmul(res, S);
                    231:   return gadd(res, (GEN) T[2]);
                    232: }
                    233: /*
                    234:  * Calcule T'(S) a l'aide du vecteur Spow
                    235:  */
                    236: static GEN
                    237: calcderivTS(GEN Spow, GEN T, GEN mod)
                    238: {
                    239:   GEN res = gzero;
                    240:   int i;
                    241:   for (i = 1; i < lg(Spow); i++)
                    242:     if (!gcmp0((GEN) T[i + 2]))
                    243:       res = gadd(res, gmul((GEN) Spow[i], stoi(i)));
                    244:   res = gmul(lift(lift(res)), mod);
                    245:   return gmodulcp(res, T);
                    246: }
                    247: /*
                    248:  * Verifie que S est une solution presque surement et calcule sa permutation
                    249:  */
                    250: static int
                    251: poltopermtest(GEN f, GEN L, GEN pf)
                    252: {
                    253:   ulong ltop;
                    254:   GEN fx, fp;
                    255:   int i, j;
                    256:   fp = cgetg(lg(L), t_VECSMALL);
                    257:   ltop = avma;
                    258:   for (i = 1; i < lg(fp); i++)
                    259:     fp[i] = 1;
                    260:   for (i = 1; i < lg(L); i++)
                    261:   {
                    262:     fx = gsubst(f, varn(f), (GEN) L[i]);
                    263:     for (j = 1; j < lg(L); j++)
                    264:       if (fp[j] && gegal(fx, (GEN) L[j]))
                    265:       {
                    266:        pf[i] = j;
                    267:        fp[j] = 0;
                    268:        break;
                    269:       }
                    270:     if (j == lg(L))
                    271:       return 0;
                    272:     avma = ltop;
                    273:   }
                    274:   return 1;
                    275: }
                    276: /*
                    277:  * Soit T un polynome de \ZZ[X] , p un nombre premier , S\in\FF_p[X]/(T) tel
                    278:  * que T(S)=0 [p,T] Relever S en S_0 tel que T(S_0)=0 [T,p^e]
                    279:  */
                    280: GEN
                    281: automorphismlift(GEN S, struct galois_lift * gl)
                    282: {
                    283:   ulong ltop, lbot;
                    284:   long x;
                    285:   GEN q, mod, modold;
                    286:   GEN W, Tr, Sr, Wr = gzero, Trold, Spow;
                    287:   int flag, init, ex;
                    288:   GEN *gptr[2];
                    289:   if (DEBUGLEVEL >= 1)
                    290:     timer2();
                    291:   x = varn(gl->T);
                    292:   Tr = gmul(gl->T, gmodulcp(gun, gl->p));
                    293:   W = ginv(gsubst(deriv(Tr, x), x, S));
                    294:   q = gl->p;
                    295:   modold = gl->p;
                    296:   Trold = Tr;
                    297:   ex = 1;
                    298:   flag = 1;
                    299:   init = 0;
                    300:   gptr[0] = &S;
                    301:   gptr[1] = &Wr;
                    302:   while (flag)
                    303:   {
                    304:     if (DEBUGLEVEL >= 1)
                    305:       timer2();
                    306:     q = gsqr(q);
                    307:     ex <<= 1;
                    308:     if (ex >= gl->e)
                    309:     {
                    310:       flag = 0;
                    311:       q = gl->Q;
                    312:       ex = gl->e;
                    313:     }
                    314:     mod = gmodulcp(gun, q);
                    315:     Tr = gmul(gl->T, mod);
                    316:     ltop = avma;
                    317:     Sr = gmodulcp(gmul(lift(lift(S)), mod), Tr);
                    318:     Spow = compoTS(Tr, Sr);
                    319:     if (init)
                    320:       W = gmul(Wr, gsub(gdeux, gmul(Wr, calcderivTS(Spow, Trold, modold))));
                    321:     else
                    322:       init = 1;
                    323:     Wr = gmodulcp(gmul(lift(lift(W)), mod), Tr);
                    324:     S = gmul(Wr, calcTS(Spow, Tr, Sr));
                    325:     lbot = avma;
                    326:     Wr = gcopy(Wr);
                    327:     S = gsub(Sr, S);
                    328:     gerepilemanysp(ltop, lbot, gptr, 2);
                    329:     modold = mod;
                    330:     Trold = Tr;
                    331:   }
                    332:   if (DEBUGLEVEL >= 1)
                    333:     msgtimer("automorphismlift()");
                    334:   return S;
                    335: }
                    336: struct galois_testlift
                    337: {
                    338:   long n;
                    339:   long f;
                    340:   long g;
                    341:   GEN bezoutcoeff;
                    342:   GEN pauto;
                    343: };
                    344: /*
                    345:  * Si polb est un polynomes de Z[X] et pola un facteur modulo p, retourne b*v
                    346:  * telqu' il existe u et a tel que: a=pola [mod p] a*b=polb [mod p^e]
                    347:  * b*v+a*u=1 [mod p^e]
                    348:  */
                    349: GEN
                    350: bezout_lift_fact(GEN pola, GEN polb, GEN p, GEN pev, long e)
                    351: {
                    352:   long ev;
                    353:   GEN ae, be, u, v, ae2, be2, s, t, pe, pe2, z, g;
                    354:   long ltop = avma, lbot;
                    355:   if (DEBUGLEVEL >= 1)
                    356:     timer2();
                    357:   ae = lift(pola);
                    358:   be = Fp_poldivres(polb, ae, p, NULL);
                    359:   g = (GEN) Fp_pol_extgcd(ae, be, p, &u, &v)[2];       /* deg g = 0 */
                    360:   if (!gcmp1(g))
                    361:   {
                    362:     g = mpinvmod(g, p);
                    363:     u = gmul(u, g);
                    364:     v = gmul(v, g);
                    365:   }
                    366:   for (pe = p, ev = 1; ev < e;)
                    367:   {
                    368:     ev <<= 1;
                    369:     pe2 = (ev >= e) ? pev : sqri(pe);
                    370:     g = gadd(polb, gneg_i(gmul(ae, be)));
                    371:     g = Fp_pol_red(g, pe2);
                    372:     g = gdivexact(g, pe);
                    373:     z = Fp_pol_red(gmul(v, g), pe);
                    374:     t = Fp_poldivres(z, ae, pe, &s);
                    375:     t = gadd(gmul(u, g), gmul(t, be));
                    376:     t = Fp_pol_red(t, pe);
                    377:     be2 = gadd(be, gmul(t, pe));
                    378:     ae2 = gadd(ae, gmul(s, pe));/* already reduced mod pe2 */
                    379:     g = gadd(gun, gneg_i(gadd(gmul(u, ae2), gmul(v, be2))));
                    380:     g = Fp_pol_red(g, pe2);
                    381:     g = gdivexact(g, pe);
                    382:     z = Fp_pol_red(gmul(v, g), pe);
                    383:     t = Fp_poldivres(z, ae, pe, &s);
                    384:     t = gadd(gmul(u, g), gmul(t, be));
                    385:     t = Fp_pol_red(t, pe);
                    386:     u = gadd(u, gmul(t, pe));
                    387:     v = gadd(v, gmul(s, pe));
                    388:     pe = pe2;
                    389:     ae = ae2;
                    390:     be = be2;
                    391:   }
                    392:   lbot = avma;
                    393:   g = gmul(v, be);
                    394:   g = gerepile(ltop, lbot, g);
                    395:   if (DEBUGLEVEL >= 1)
                    396:     msgtimer("bezout_lift_fact()");
                    397:   return g;
                    398: }
                    399: static long
                    400: inittestlift(GEN Tmod, long elift, struct galois_lift * gl, struct galois_testlift * gt, GEN frob)
                    401: {
                    402:   ulong ltop = avma;
                    403:   int i, j;
                    404:   long v;
                    405:   GEN pe, autpow, plift;
                    406:   GEN Tmodp, xmodp, modQ, TmodQ, xmodQ;
                    407:   GEN *gptr[2];
                    408:   v = varn(gl->T);
                    409:   gt->n = lg(gl->L) - 1;
                    410:   gt->g = lg(Tmod) - 1;
                    411:   gt->f = gt->n / gt->g;
                    412:   Tmodp = gmul(gl->T, gmodulcp(gun, gl->p));
                    413:   xmodp = gmodulcp(gmul(polx[v], gmodulcp(gun, gl->p)), Tmodp);
                    414:   pe = gpowgs(gl->p, elift);
                    415:   plift = automorphismlift(powgi(xmodp, pe), gl);
                    416:   if (frob)
                    417:   {
                    418:     GEN tlift;
                    419:     tlift = gdiv(centerlift(gmul((GEN) plift[2], gl->den)), gl->den);
                    420:     if (poltopermtest(tlift, gl->L, frob))
                    421:     {
                    422:       avma = ltop;
                    423:       return 1;
                    424:     }
                    425:   }
                    426:   modQ = gmodulcp(gun, gl->Q);
                    427:   TmodQ = gmul(gl->T, modQ);
                    428:   xmodQ = gmodulcp(gmul(polx[v], modQ), TmodQ);
                    429:   gt->bezoutcoeff = cgetg(gt->g + 1, t_VEC);
                    430:   for (i = 1; i <= gt->g; i++)
                    431:   {
                    432:     GEN blift;
                    433:     blift = bezout_lift_fact((GEN) Tmod[i], gl->T, gl->p, gl->Q, gl->e);
                    434:     gt->bezoutcoeff[i] = (long) gmodulcp(gmul(blift, modQ), TmodQ);
                    435:   }
                    436:   gt->pauto = cgetg(gt->f + 1, t_VEC);
                    437:   gt->pauto[1] = (long) xmodQ;
                    438:   gt->pauto[2] = (long) plift;
                    439:   if (gt->f > 2)
                    440:   {
                    441:     autpow = cgetg(gt->n, t_VEC);
                    442:     autpow[1] = (long) plift;
                    443:     for (i = 2; i < gt->n; i++)
                    444:       autpow[i] = lmul((GEN) autpow[i - 1], plift);
                    445:     for (i = 3; i <= gt->f; i++)
                    446:     {
                    447:       GEN s, P;
                    448:       P = ((GEN **) gt->pauto)[i - 1][2];
                    449:       s = (GEN) P[2];
                    450:       for (j = 1; j < lgef(P) - 2; j++)
                    451:        s = gadd(s, gmul((GEN) autpow[j], (GEN) P[j + 2]));
                    452:       gt->pauto[i] = (long) s;
                    453:     }
                    454:   }
                    455:   gptr[0] = &gt->bezoutcoeff;
                    456:   gptr[1] = &gt->pauto;
                    457:   gerepilemany(ltop, gptr, 2);
                    458:   return 0;
                    459: }
                    460: /*
                    461:  *
                    462:  */
                    463: long
                    464: frobeniusliftall(GEN sg, GEN * psi, struct galois_lift * gl, struct galois_testlift * gt, GEN frob)
                    465: {
                    466:   ulong ltop = avma, av, ltop2;
                    467:   long d, N, z, m, c;
                    468:   int i, j, k;
                    469:   GEN pf, u, v;
                    470:   m = gt->g;
                    471:   c = lg(sg) - 1;
                    472:   d = m / c;
                    473:   pf = cgetg(m, t_VECSMALL);
                    474:   *psi = pf;
                    475:   ltop2 = avma;
                    476:   N = itos(gdiv(mpfact(m), gmul(stoi(c), gpowgs(mpfact(d), c))));
                    477:   avma = ltop2;
                    478:   v = gmul((GEN) gt->pauto[2], (GEN) gt->bezoutcoeff[m]);
                    479:   for (i = 1; i < m; i++)
                    480:     pf[i] = 1 + i / d;
                    481:   av = avma;
                    482:   for (i = 0;; i++)
                    483:   {
                    484:     u = v;
                    485:     if (DEBUGLEVEL >= 4)
                    486:     {
                    487:       fprintferr("GaloisConj:Testing %Z:%d:%Z:", sg, i, pf);
                    488:       timer2();
                    489:     }
                    490:     for (j = 1; j < m; j++)
                    491:       u = gadd(u, gmul((GEN) gt->pauto[sg[pf[j]] + 1], (GEN) gt->bezoutcoeff[j]));
                    492:     u = gdiv(centerlift(gmul((GEN) u[2], gl->den)), gl->den);
                    493:     if (poltopermtest(u, gl->L, frob))
                    494:     {
                    495:       if (DEBUGLEVEL >= 4)
                    496:        msgtimer("");
                    497:       avma = ltop2;
                    498:       return 1;
                    499:     }
                    500:     if (DEBUGLEVEL >= 4)
                    501:       msgtimer("");
                    502:     avma = av;
                    503:     if (i == N - 1)
                    504:       break;
                    505:     for (j = 2; j < m && pf[j - 1] >= pf[j]; j++);
                    506:     for (k = 1; k < j - k && pf[k] != pf[j - k]; k++)
                    507:     {
                    508:       z = pf[k];
                    509:       pf[k] = pf[j - k];
                    510:       pf[j - k] = z;
                    511:     }
                    512:     for (k = j - 1; pf[k] >= pf[j]; k--);
                    513:     z = pf[j];
                    514:     pf[j] = pf[k];
                    515:     pf[k] = z;
                    516:   }
                    517:   avma = ltop;
                    518:   *psi = NULL;
                    519:   return 0;
                    520: }
                    521: /*
                    522:  * applique une permutation t a un vecteur s sans duplication
                    523:  */
                    524: static GEN
                    525: applyperm(GEN s, GEN t)
                    526: {
                    527:   GEN u;
                    528:   int i;
                    529:   if (lg(s) < lg(t))
                    530:     err(talker, "First permutation shorter than second in applyperm");
                    531:   u = cgetg(lg(s), typ(s));
                    532:   for (i = 1; i < lg(s); i++)
                    533:     u[i] = s[t[i]];
                    534:   return u;
                    535: }
                    536: /*
                    537:  * alloue une ardoise pour n entiers de longueurs pour les test de
                    538:  * permutation
                    539:  */
                    540: GEN
                    541: alloue_ardoise(long n, long s)
                    542: {
                    543:   GEN ar;
                    544:   long i;
                    545:   ar = cgetg(n + 1, t_VEC);
                    546:   for (i = 1; i <= n; i++)
                    547:     ar[i] = lgetg(s, t_INT);
                    548:   return ar;
                    549: }
                    550: /*
                    551:  * structure contenant toutes les données pour le tests des permutations:
                    552:  *
                    553:  * ordre :ordre des tests pour verifie_test ordre[lg(ordre)]: numero du test
                    554:  * principal borne : borne sur les coefficients a trouver ladic: modulo
                    555:  * l-adique des racines lborne:ladic-borne TM:vecteur des ligne de M
                    556:  * PV:vecteur des clones des matrices de test (Vmatrix) (ou NULL si non
                    557:  * calcule) L,M comme d'habitude (voir plus bas)
                    558:  */
                    559: struct galois_test
                    560: {
                    561:   GEN ordre;
                    562:   GEN borne, lborne, ladic;
                    563:   GEN PV, TM;
                    564:   GEN L, M;
                    565: };
                    566: /* Calcule la matrice de tests correspondant a la n-ieme ligne de V */
                    567: static GEN
                    568: Vmatrix(long n, struct galois_test * td)
                    569: {
                    570:   ulong ltop = avma, lbot;
                    571:   GEN V;
                    572:   long i;
                    573:   V = cgetg(lg(td->L), t_VEC);
                    574:   for (i = 1; i < lg(V); i++)
                    575:     V[i] = ((GEN **) td->M)[i][n][2];
                    576:   V = centerlift(gmul(td->L, V));
                    577:   lbot = avma;
                    578:   V = gmod(V, td->ladic);
                    579:   return gerepile(ltop, lbot, V);
                    580: }
                    581: /*
                    582:  * Initialise la structure galois_test
                    583:  */
                    584: static void
                    585: inittest(GEN L, GEN M, GEN borne, GEN ladic, struct galois_test * td)
                    586: {
                    587:   ulong ltop;
                    588:   long i;
                    589:   int n = lg(L) - 1;
                    590:   if (DEBUGLEVEL >= 8)
                    591:     fprintferr("GaloisConj:Entree Init Test\n");
                    592:   td->ordre = cgetg(n + 1, t_VECSMALL);
                    593:   for (i = 1; i <= n - 2; i++)
                    594:     td->ordre[i] = i + 2;
                    595:   for (; i <= n; i++)
                    596:     td->ordre[i] = i - n + 2;
                    597:   td->borne = borne;
                    598:   td->lborne = gsub(ladic, borne);
                    599:   td->ladic = ladic;
                    600:   td->L = L;
                    601:   td->M = M;
                    602:   td->PV = cgetg(n + 1, t_VECSMALL);   /* peut-etre t_VEC est correct ici */
                    603:   for (i = 1; i <= n; i++)
                    604:     td->PV[i] = 0;
                    605:   ltop = avma;
                    606:   td->PV[td->ordre[n]] = (long) gclone(Vmatrix(td->ordre[n], td));
                    607:   avma = ltop;
                    608:   td->TM = gtrans(M);
                    609:   settyp(td->TM, t_VEC);
                    610:   for (i = 1; i < lg(td->TM); i++)
                    611:     settyp(td->TM[i], t_VEC);
                    612:   if (DEBUGLEVEL >= 8)
                    613:     fprintferr("GaloisConj:Sortie Init Test\n");
                    614: }
                    615: /*
                    616:  * liberer les clones de la structure galois_test
                    617:  *
                    618:  * Reservé a l'accreditation ultra-violet:Liberez les clones! Paranoia(tm)
                    619:  */
                    620: static void
                    621: freetest(struct galois_test * td)
                    622: {
                    623:   int i;
                    624:   for (i = 1; i < lg(td->PV); i++)
                    625:     if (td->PV[i])
                    626:     {
                    627:       gunclone((GEN) td->PV[i]);
                    628:       td->PV[i] = 0;
                    629:     }
                    630: }
                    631: /*
                    632:  * Test si le nombre padique P est proche d'un entier inferieur a td->borne
                    633:  * en valeur absolue.
                    634:  */
                    635: long
                    636: padicisint(GEN P, struct galois_test * td)
                    637: {
                    638:   long r;
                    639:   ulong ltop = avma;
                    640:   GEN U;
                    641:   if (typ(P) != t_INT)
                    642:     err(typeer, "padicisint");
                    643:   U = modii(P, td->ladic);
                    644:   r = gcmp(U, (GEN) td->borne) <= 0 || gcmp(U, (GEN) td->lborne) >= 0;
                    645:   avma = ltop;
                    646:   return r;
                    647: }
                    648: /*
                    649:  * Verifie si pf est une vrai solution et non pas un "hop"
                    650:  */
                    651: static long
                    652: verifietest(GEN pf, struct galois_test * td)
                    653: {
                    654:   ulong av = avma;
                    655:   GEN P, V;
                    656:   int i, j;
                    657:   int n = lg(td->L) - 1;
                    658:   if (DEBUGLEVEL >= 8)
                    659:     fprintferr("GaloisConj:Entree Verifie Test\n");
                    660:   P = applyperm(td->L, pf);
                    661:   for (i = 1; i < n; i++)
                    662:   {
                    663:     long ord;
                    664:     GEN PW;
                    665:     ord = td->ordre[i];
                    666:     PW = (GEN) td->PV[ord];
                    667:     if (PW)
                    668:     {
                    669:       V = ((GEN **) PW)[1][pf[1]];
                    670:       for (j = 2; j <= n; j++)
                    671:        V = gadd(V, ((GEN **) PW)[j][pf[j]]);
                    672:     } else
                    673:       V = centerlift(gmul((GEN) td->TM[ord], P));
                    674:     if (!padicisint(V, td))
                    675:       break;
                    676:   }
                    677:   if (i == n)
                    678:   {
                    679:     if (DEBUGLEVEL >= 8)
                    680:       fprintferr("GaloisConj:Sortie Verifie Test:1\n");
                    681:     avma = av;
                    682:     return 1;
                    683:   }
                    684:   if (!td->PV[td->ordre[i]])
                    685:   {
                    686:     td->PV[td->ordre[i]] = (long) gclone(Vmatrix(td->ordre[i], td));
                    687:     if (DEBUGLEVEL >= 2)
                    688:       fprintferr("M");
                    689:   }
                    690:   if (DEBUGLEVEL >= 2)
                    691:     fprintferr("%d.", i);
                    692:   if (i > 1)
                    693:   {
                    694:     long z;
                    695:     z = td->ordre[i];
                    696:     for (j = i; j > 1; j--)
                    697:       td->ordre[j] = td->ordre[j - 1];
                    698:     td->ordre[1] = z;
                    699:     if (DEBUGLEVEL >= 6)
                    700:       fprintferr("%Z", td->ordre);
                    701:   }
                    702:   if (DEBUGLEVEL >= 8)
                    703:     fprintferr("GaloisConj:Sortie Verifie Test:0\n");
                    704:   avma = av;
                    705:   return 0;
                    706: }
                    707: /*
                    708:  * F est la decomposition en cycle de sigma, B est la decomposition en cycle
                    709:  * de cl(tau) Teste toutes les permutations pf pouvant correspondre a tau tel
                    710:  * que : tau*sigma*tau^-1=sigma^s et tau^d=sigma^t  ou d est l'ordre de
                    711:  * cl(tau) --------- x: vecteur des choix y: vecteur des mises a jour G:
                    712:  * vecteur d'acces a F linéaire
                    713:  */
                    714: GEN
                    715: testpermutation(GEN F, GEN B, long s, long t, long cut, struct galois_test * td)
                    716: {
                    717:   ulong av, avm = avma, av2;
                    718:   long a, b, c, d, n;
                    719:   GEN pf, x, ar, y, *G;
                    720:   int N, cx, p1, p2, p3, p4, p5, p6;
                    721:   int i, j, hop = 0;
                    722:   GEN V, W;
                    723:   if (DEBUGLEVEL >= 1)
                    724:     timer2();
                    725:   a = lg(F) - 1;
                    726:   b = lg(F[1]) - 1;
                    727:   c = lg(B) - 1;
                    728:   d = lg(B[1]) - 1;
                    729:   n = a * b;
                    730:   s = (b + s) % b;
                    731:   pf = cgetg(n + 1, t_VECSMALL);
                    732:   av = avma;
                    733:   ar = alloue_ardoise(a, 1 + lg(td->ladic));
                    734:   x = cgetg(a + 1, t_VECSMALL);
                    735:   y = cgetg(a + 1, t_VECSMALL);
                    736:   G = (GEN *) cgetg(a + 1, t_VECSMALL);        /* Don't worry */
                    737:   av2 = avma;
                    738:   W = (GEN) td->PV[td->ordre[n]];
                    739:   for (cx = 1, i = 1, j = 1; cx <= a; cx++, i++)
                    740:   {
                    741:     x[cx] = (i != d) ? 0 : t;
                    742:     y[cx] = 1;
                    743:     G[cx] = (GEN) F[((long **) B)[j][i]];      /* Be happy */
                    744:     if (i == d)
                    745:     {
                    746:       i = 0;
                    747:       j++;
                    748:     }
                    749:   }                            /* Be happy now! */
                    750:   N = itos(gpowgs(stoi(b), c * (d - 1))) / cut;
                    751:   avma = av2;
                    752:   if (DEBUGLEVEL >= 4)
                    753:     fprintferr("GaloisConj:I will try %d permutations\n", N);
                    754:   for (cx = 0; cx < N; cx++)
                    755:   {
                    756:     if (DEBUGLEVEL >= 2 && cx % 1000 == 999)
                    757:       fprintferr("%d%% ", (100 * cx) / N);
                    758:     if (cx)
                    759:     {
                    760:       for (i = 1, j = d; i < a;)
                    761:       {
                    762:        y[i] = 1;
                    763:        if ((++(x[i])) != b)
                    764:          break;
                    765:        x[i++] = 0;
                    766:        if (i == j)
                    767:        {
                    768:          y[i++] = 1;
                    769:          j += d;
                    770:        }
                    771:       }
                    772:       y[i + 1] = 1;
                    773:     }
                    774:     for (p1 = 1, p5 = d; p1 <= a; p1++, p5++)
                    775:       if (y[p1])
                    776:       {
                    777:        if (p5 == d)
                    778:        {
                    779:          p5 = 0;
                    780:          p4 = 0;
                    781:        } else
                    782:          p4 = x[p1 - 1];
                    783:        if (p5 == d - 1)
                    784:          p6 = p1 + 1 - d;
                    785:        else
                    786:          p6 = p1 + 1;
                    787:        y[p1] = 0;
                    788:        V = gzero;
                    789:        for (p2 = 1 + p4, p3 = 1 + x[p1]; p2 <= b; p2++)
                    790:        {
                    791:          V = gadd(V, ((GEN **) W)[G[p6][p3]][G[p1][p2]]);
                    792:          p3 += s;
                    793:          if (p3 > b)
                    794:            p3 -= b;
                    795:        }
                    796:        p3 = 1 + x[p1] - s;
                    797:        if (p3 <= 0)
                    798:          p3 += b;
                    799:        for (p2 = p4; p2 >= 1; p2--)
                    800:        {
                    801:          V = gadd(V, ((GEN **) W)[G[p6][p3]][G[p1][p2]]);
                    802:          p3 -= s;
                    803:          if (p3 <= 0)
                    804:            p3 += b;
                    805:        }
                    806:        gaffect((GEN) V, (GEN) ar[p1]);
                    807:       }
                    808:     V = gzero;
                    809:     for (p1 = 1; p1 <= a; p1++)
                    810:       V = gadd(V, (GEN) ar[p1]);
                    811:     if (padicisint(V, td))
                    812:     {
                    813:       for (p1 = 1, p5 = d; p1 <= a; p1++, p5++)
                    814:       {
                    815:        if (p5 == d)
                    816:        {
                    817:          p5 = 0;
                    818:          p4 = 0;
                    819:        } else
                    820:          p4 = x[p1 - 1];
                    821:        if (p5 == d - 1)
                    822:          p6 = p1 + 1 - d;
                    823:        else
                    824:          p6 = p1 + 1;
                    825:        for (p2 = 1 + p4, p3 = 1 + x[p1]; p2 <= b; p2++)
                    826:        {
                    827:          pf[G[p1][p2]] = G[p6][p3];
                    828:          p3 += s;
                    829:          if (p3 > b)
                    830:            p3 -= b;
                    831:        }
                    832:        p3 = 1 + x[p1] - s;
                    833:        if (p3 <= 0)
                    834:          p3 += b;
                    835:        for (p2 = p4; p2 >= 1; p2--)
                    836:        {
                    837:          pf[G[p1][p2]] = G[p6][p3];
                    838:          p3 -= s;
                    839:          if (p3 <= 0)
                    840:            p3 += b;
                    841:        }
                    842:       }
                    843:       if (verifietest(pf, td))
                    844:       {
                    845:        avma = av;
                    846:        if (DEBUGLEVEL >= 1)
                    847:          msgtimer("testpermutation(%d)", cx + 1);
                    848:        if (DEBUGLEVEL >= 2 && hop)
                    849:          fprintferr("GaloisConj:%d hop sur %d iterations\n", hop, cx + 1);
                    850:        return pf;
                    851:       } else
                    852:        hop++;
                    853:     }
                    854:     avma = av2;
                    855:   }
                    856:   avma = avm;
                    857:   if (DEBUGLEVEL >= 1)
                    858:     msgtimer("testpermutation(%d)", N);
                    859:   if (DEBUGLEVEL >= 2 && hop)
                    860:     fprintferr("GaloisConj:%d hop sur %d iterations\n", hop, N);
                    861:   return gzero;
                    862: }
                    863: /*
                    864:  * Calcule la liste des sous groupes de \ZZ/m\ZZ d'ordre divisant p et
                    865:  * retourne la liste de leurs elements
                    866:  */
                    867: GEN
                    868: listsousgroupes(long m, long p)
                    869: {
                    870:   ulong ltop = avma, lbot;
                    871:   GEN zns, lss, res, sg, gen, ord, res1;
                    872:   int k, card, i, j, l, n, phi, h;
                    873:   if (m == 2)
                    874:   {
                    875:     res = cgetg(2, t_VEC);
                    876:     sg = cgetg(2, t_VECSMALL);
                    877:     res[1] = (long) sg;
                    878:     sg[1] = 1;
                    879:     return res;
                    880:   }
                    881:   zns = znstar(stoi(m));
                    882:   phi = itos((GEN) zns[1]);
                    883:   lss = subgrouplist((GEN) zns[2], 0);
                    884:   gen = cgetg(lg(zns[3]), t_VECSMALL);
                    885:   ord = cgetg(lg(zns[3]), t_VECSMALL);
                    886:   res1 = cgetg(lg(lss), t_VECSMALL);
                    887:   lbot = avma;
                    888:   for (k = 1, i = lg(lss) - 1; i >= 1; i--)
                    889:   {
                    890:     long av;
                    891:     av = avma;
                    892:     card = phi / itos(det((GEN) lss[i]));
                    893:     avma = av;
                    894:     if (p % card == 0)
                    895:     {
                    896:       sg = cgetg(1 + card, t_VECSMALL);
                    897:       sg[1] = 1;
                    898:       av = avma;
                    899:       for (j = 1; j < lg(gen); j++)
                    900:       {
                    901:        gen[j] = 1;
                    902:        for (h = 1; h < lg(lss[i]); h++)
                    903:          gen[j] = (gen[j] * itos(lift(powgi(((GEN **) zns)[3][h],
                    904:                                           ((GEN ***) lss)[i][j][h])))) % m;
                    905:        ord[j] = itos(((GEN **) zns)[2][j]) / itos(((GEN ***) lss)[i][j][j]);
                    906:       }
                    907:       avma = av;
                    908:       for (j = 1, l = 1; j < lg(gen); j++)
                    909:       {
                    910:        int c = l * (ord[j] - 1);
                    911:        for (n = 1; n <= c; n++)/* I like it */
                    912:          sg[++l] = (sg[n] * gen[j]) % m;
                    913:       }
                    914:       res1[k++] = (long) sg;
                    915:     }
                    916:   }
                    917:   res = cgetg(k, t_VEC);
                    918:   for (i = 1; i < k; i++)
                    919:     res[i] = res1[i];
                    920:   return gerepile(ltop, lbot, res);
                    921: }
                    922: /* retourne la permutation identite */
                    923: GEN
                    924: permidentity(long l)
                    925: {
                    926:   GEN perm;
                    927:   int i;
                    928:   perm = cgetg(l + 1, t_VECSMALL);
                    929:   for (i = 1; i <= l; i++)
                    930:     perm[i] = i;
                    931:   return perm;
                    932: }
                    933: /* retourne la decomposition en cycle */
                    934: GEN
                    935: permtocycle(GEN p)
                    936: {
                    937:   long ltop = avma, lbot;
                    938:   int i, j, k, l, m;
                    939:   GEN bit, cycle, cy;
                    940:   cycle = cgetg(lg(p), t_VEC);
                    941:   bit = cgetg(lg(p), t_VECSMALL);
                    942:   for (i = 1; i < lg(p); i++)
                    943:     bit[i] = 0;
                    944:   for (k = 1, l = 1; k < lg(p);)
                    945:   {
                    946:     for (j = 1; bit[j]; j++);
                    947:     cy = cgetg(lg(p), t_VECSMALL);
                    948:     m = 1;
                    949:     do
                    950:     {
                    951:       bit[j] = 1;
                    952:       k++;
                    953:       cy[m++] = j;
                    954:       j = p[j];
                    955:     } while (!bit[j]);
                    956:     setlg(cy, m);
                    957:     cycle[l++] = (long) cy;
                    958:   }
                    959:   setlg(cycle, l);
                    960:   lbot = avma;
                    961:   cycle = gcopy(cycle);
                    962:   return gerepile(ltop, lbot, cycle);
                    963: }
                    964: /* Calcule les orbites d'un ensemble de permutations */
                    965: GEN
                    966: vecpermorbite(GEN v)
                    967: {
                    968:   ulong ltop = avma, lbot;
                    969:   int i, j, k, l, m, n, o, p, flag;
                    970:   GEN bit, cycle, cy;
                    971:   n = lg(v[1]);
                    972:   cycle = cgetg(n, t_VEC);
                    973:   bit = cgetg(n, t_VECSMALL);
                    974:   for (i = 1; i < n; i++)
                    975:     bit[i] = 0;
                    976:   for (k = 1, l = 1; k < n;)
                    977:   {
                    978:     for (j = 1; bit[j]; j++);
                    979:     cy = cgetg(n, t_VECSMALL);
                    980:     m = 1;
                    981:     cy[m++] = j;
                    982:     bit[j] = 1;
                    983:     k++;
                    984:     do
                    985:     {
                    986:       flag = 0;
                    987:       for (o = 1; o < lg(v); o++)
                    988:       {
                    989:        for (p = 1; p < m; p++) /* m varie! */
                    990:        {
                    991:          j = ((long **) v)[o][cy[p]];
                    992:          if (!bit[j])
                    993:          {
                    994:            flag = 1;
                    995:            bit[j] = 1;
                    996:            k++;
                    997:            cy[m++] = j;
                    998:          }
                    999:        }
                   1000:       }
                   1001:     } while (flag);
                   1002:     setlg(cy, m);
                   1003:     cycle[l++] = (long) cy;
                   1004:   }
                   1005:   setlg(cycle, l);
                   1006:   lbot = avma;
                   1007:   cycle = gcopy(cycle);
                   1008:   return gerepile(ltop, lbot, cycle);
                   1009: }
                   1010: /*
                   1011:  * Calcule la permutation associe a un polynome f des racines L
                   1012:  */
                   1013: GEN
                   1014: poltoperm(GEN f, GEN L)
                   1015: {
                   1016:   ulong ltop, ltop2;
                   1017:   GEN pf, fx, fp;
                   1018:   int i, j;
                   1019:   pf = cgetg(lg(L), t_VECSMALL);
                   1020:   ltop = avma;
                   1021:   fp = cgetg(lg(L), t_VECSMALL);
                   1022:   ltop2 = avma;
                   1023:   for (i = 1; i < lg(fp); i++)
                   1024:     fp[i] = 1;
                   1025:   for (i = 1; i < lg(L); i++)
                   1026:   {
                   1027:     fx = gsubst(f, varn(f), (GEN) L[i]);
                   1028:     for (j = 1; j < lg(L); j++)
                   1029:       if (fp[j] && gegal(fx, (GEN) L[j]))
                   1030:       {
                   1031:        pf[i] = j;
                   1032:        fp[j] = 0;
                   1033:        break;
                   1034:       }
                   1035:     avma = ltop2;
                   1036:   }
                   1037:   avma = ltop;
                   1038:   return pf;
                   1039: }
                   1040: /*
                   1041:  * Calcule un polynome R definissant  le corps fixe de T pour les orbites O
                   1042:  * tel que R est sans-facteur carre modulo mod et l retourne dans U les
                   1043:  * racines de R
                   1044:  */
                   1045: GEN
                   1046: corpsfixeorbitemod(GEN L, GEN O, long x, GEN mod, GEN l, GEN * U)
                   1047: {
                   1048:   ulong ltop = avma, lbot, av, av2;
                   1049:   GEN g, p, PL, *gptr[2], gmod, gl, modl;
                   1050:   GEN id;
                   1051:   int i, j, d, dmax;
                   1052:   PL = cgetg(lg(O), t_COL);
                   1053:   modl = gmodulcp(gun, l);
                   1054:   av2 = avma;
                   1055:   dmax = lg(L) + 1;
                   1056:   d = 0;
                   1057:   do
                   1058:   {
                   1059:     avma = av2;
                   1060:     id = stoi(d);
                   1061:     g = polun[x];
                   1062:     for (i = 1; i < lg(O); i++)
                   1063:     {
                   1064:       p = gadd(id, (GEN) L[((GEN *) O)[i][1]]);
                   1065:       for (j = 2; j < lg(O[i]); j++)
                   1066:        p = gmul(p, gadd(id, (GEN) L[((GEN *) O)[i][j]]));
                   1067:       PL[i] = (long) p;
                   1068:       g = gmul(g, gsub(polx[x], p));
                   1069:     }
                   1070:     lbot = avma;
                   1071:     g = centerlift(g);
                   1072:     av = avma;
                   1073:     gmod = gmul(g, mod);
                   1074:     gl = gmul(g, modl);
                   1075:     if (DEBUGLEVEL >= 6)
                   1076:       fprintferr("GaloisConj:corps fixe:%d:%Z\n", d, g);
                   1077:     d = (d > 0 ? -d : 1 - d);
                   1078:     if (d > dmax)
                   1079:       err(talker, "prime too small in corpsfixeorbitemod");
                   1080:   } while (degree(ggcd(deriv(gl, x), gl)) || degree(ggcd(deriv(gmod, x), gmod)));
                   1081:   avma = av;
                   1082:   *U = gcopy(PL);
                   1083:   gptr[0] = &g;
                   1084:   gptr[1] = U;
                   1085:   gerepilemanysp(ltop, lbot, gptr, 2);
                   1086:   return g;
                   1087: }
                   1088: /*
                   1089:  * Calcule l'inclusion de R dans T i.e. S telque T|R\compo S
                   1090:  */
                   1091: GEN
                   1092: corpsfixeinclusion(GEN O, GEN PL)
                   1093: {
                   1094:   GEN S;
                   1095:   int i, j;
                   1096:   S = cgetg((lg(O) - 1) * (lg(O[1]) - 1) + 1, t_COL);
                   1097:   for (i = 1; i < lg(O); i++)
                   1098:     for (j = 1; j < lg(O[i]); j++)
                   1099:       S[((GEN *) O)[i][j]] = lcopy((GEN) PL[i]);
                   1100:   return S;
                   1101: }
                   1102: /* Calcule l'inverse de la matrice de van der Monde de T multiplie par den */
                   1103: GEN
                   1104: vandermondeinverse(GEN L, GEN T, GEN den)
                   1105: {
                   1106:   ulong ltop = avma, lbot;
                   1107:   int i, j, n = lg(L);
                   1108:   long x = varn(T);
                   1109:   GEN M, P, Tp;
                   1110:   if (DEBUGLEVEL >= 1)
                   1111:     timer2();
                   1112:   M = cgetg(n, t_MAT);
                   1113:   Tp = deriv(T, x);
                   1114:   for (i = 1; i < n; i++)
                   1115:   {
                   1116:     M[i] = lgetg(n, t_COL);
                   1117:     P = gdiv(gdeuc(T, gsub(polx[x], (GEN) L[i])), gsubst(Tp, x, (GEN) L[i]));
                   1118:     for (j = 1; j < n; j++)
                   1119:       ((GEN *) M)[i][j] = P[1 + j];
                   1120:   }
                   1121:   if (DEBUGLEVEL >= 1)
                   1122:     msgtimer("vandermondeinverse");
                   1123:   lbot = avma;
                   1124:   M = gmul(den, M);
                   1125:   return gerepile(ltop, lbot, M);
                   1126: }
                   1127: /* Calcule le polynome associe a un vecteur conjugue v */
                   1128: static GEN
                   1129: vectopol(GEN v, GEN L, GEN M, GEN den, long x)
                   1130: {
                   1131:   ulong ltop = avma, lbot;
                   1132:   GEN res;
                   1133:   res = gmul(M, v);
                   1134:   res = gtopolyrev(centerlift(res), x);
                   1135:   lbot = avma;
                   1136:   res = gdiv(res, den);
                   1137:   return gerepile(ltop, lbot, res);
                   1138: }
                   1139: /* Calcule le polynome associe a une permuation p */
                   1140: static GEN
                   1141: permtopol(GEN p, GEN L, GEN M, GEN den, long x)
                   1142: {
                   1143:   ulong ltop = avma, lbot;
                   1144:   GEN res;
                   1145:   res = gmul(M, applyperm(L, p));
                   1146:   res = gtopolyrev(centerlift(res), x);
                   1147:   lbot = avma;
                   1148:   res = gdiv(res, den);
                   1149:   return gerepile(ltop, lbot, res);
                   1150: }
                   1151: /*
                   1152:  * factorise partiellement n sous la forme n=d*u*f^2 ou d est sans facteur
                   1153:  * carre et u n'est pas un carre parfait et retourne u*f
                   1154:  */
                   1155: GEN
                   1156: mycoredisc(GEN n)
                   1157: {
                   1158:   ulong av = avma, tetpil, r;
                   1159:   GEN y, p1, p2, ret;
                   1160:   {
                   1161:     long av = avma, tetpil, i;
                   1162:     GEN fa, p1, p2, p3, res = gun, res2 = gun, y;
                   1163:     fa = auxdecomp(n, 0);
                   1164:     p1 = (GEN) fa[1];
                   1165:     p2 = (GEN) fa[2];
                   1166:     for (i = 1; i < lg(p1) - 1; i++)
                   1167:     {
                   1168:       p3 = (GEN) p2[i];
                   1169:       if (mod2(p3))
                   1170:        res = mulii(res, (GEN) p1[i]);
                   1171:       if (!gcmp1(p3))
                   1172:        res2 = mulii(res2, gpui((GEN) p1[i], shifti(p3, -1), 0));
                   1173:     }
                   1174:     p3 = (GEN) p2[i];
                   1175:     if (mod2(p3))              /* impaire: verifions */
                   1176:     {
                   1177:       if (!gcmp1(p3))
                   1178:        res2 = mulii(res2, gpui((GEN) p1[i], shifti(p3, -1), 0));
                   1179:       if (isprime((GEN) p1[i]))
                   1180:        res = mulii(res, (GEN) p1[i]);
                   1181:       else
                   1182:        res2 = mulii(res2, (GEN) p1[i]);
                   1183:     } else                     /* paire:OK */
                   1184:       res2 = mulii(res2, gpui((GEN) p1[i], shifti(p3, -1), 0));
                   1185:     tetpil = avma;
                   1186:     y = cgetg(3, t_VEC);
                   1187:     y[1] = (long) icopy(res);
                   1188:     y[2] = (long) icopy(res2);
                   1189:     ret = gerepile(av, tetpil, y);
                   1190:   }
                   1191:   p2 = ret;
                   1192:   p1 = (GEN) p2[1];
                   1193:   r = mod4(p1);
                   1194:   if (signe(p1) < 0)
                   1195:     r = 4 - r;
                   1196:   if (r == 1 || r == 4)
                   1197:     return (GEN) p2[2];
                   1198:   tetpil = avma;
                   1199:   y = gmul2n((GEN) p2[2], -1);
                   1200:   return gerepile(av, tetpil, y);
                   1201: }
                   1202: /* Calcule la puissance exp d'une permutation decompose en cycle */
                   1203: GEN
                   1204: permcyclepow(GEN O, long exp)
                   1205: {
                   1206:   int j, k, n;
                   1207:   GEN p;
                   1208:   for (n = 0, j = 1; j < lg(O); j++)
                   1209:     n += lg(O[j]) - 1;
                   1210:   p = cgetg(n + 1, t_VECSMALL);
                   1211:   for (j = 1; j < lg(O); j++)
                   1212:   {
                   1213:     n = lg(O[j]) - 1;
                   1214:     for (k = 1; k <= n; k++)
                   1215:       p[((long **) O)[j][k]] = ((long **) O)[j][1 + (k + exp - 1) % n];
                   1216:   }
                   1217:   return p;
                   1218: }
                   1219: /*
                   1220:  * Casse l'orbite O en sous-orbite d'ordre premier correspondant a des
                   1221:  * puissance de l'element
                   1222:  */
                   1223: GEN
                   1224: splitorbite(GEN O)
                   1225: {
                   1226:   ulong ltop = avma, lbot;
                   1227:   int i, n;
                   1228:   GEN F, fc, res;
                   1229:   n = lg(O[1]) - 1;
                   1230:   F = factor(stoi(n));
                   1231:   fc = cgetg(lg(F[1]), t_VECSMALL);
                   1232:   for (i = 1; i < lg(fc); i++)
                   1233:     fc[i] = itos(gpow(((GEN **) F)[1][i], ((GEN **) F)[2][i], DEFAULTPREC));
                   1234:   lbot = avma;
                   1235:   res = cgetg(lg(fc), t_VEC);
                   1236:   for (i = 1; i < lg(fc); i++)
                   1237:   {
                   1238:     GEN v;
                   1239:     v = cgetg(3, t_VEC);
                   1240:     res[i] = (long) v;
                   1241:     v[1] = (long) permcyclepow(O, n / fc[i]);
                   1242:     v[2] = (long) stoi(fc[i]);
                   1243:   }
                   1244:   return gerepile(ltop, lbot, res);
                   1245: }
                   1246: /*
                   1247:  * structure qui contient l'analyse du groupe de Galois par etude des degres
                   1248:  * de Frobenius:
                   1249:  *
                   1250:  * p: nombre premier a relever deg: degre du lift à effectuer pgp: plus grand
                   1251:  * diviseur premier du degre de T. exception: 1 si le pgp-Sylow est
                   1252:  * non-cyclique. l: un nombre premier tel que T est totalement decompose
                   1253:  * modulo l ppp:  plus petit diviseur premier du degre de T. primepointer:
                   1254:  * permet de calculer les premiers suivant p.
                   1255:  */
                   1256: struct galois_analysis
                   1257: {
                   1258:   long p;
                   1259:   long deg;
                   1260:   long exception;
                   1261:   long ord;
                   1262:   long l;
                   1263:   long ppp;
                   1264:   long p4;
                   1265:   byteptr primepointer;
                   1266: };
                   1267: /* peut-etre on peut accelerer(distinct degree factorisation */
                   1268: void
                   1269: galoisanalysis(GEN T, struct galois_analysis * ga, long calcul_l)
                   1270: {
                   1271:   ulong ltop = avma;
                   1272:   long n, p, ex, plift, nbmax, nbtest, exist6 = 0, p4 = 0;
                   1273:   GEN F, fc;
                   1274:   byteptr primepointer, pp = 0;
                   1275:   long pha, ord, deg, ppp, pgp, ordmax = 0, l = 0;
                   1276:   int i;
                   1277:   /* Etude du cardinal: */
                   1278:   /* Calcul de l'ordre garanti d'un sous-groupe cyclique */
                   1279:   /* Tout les groupes d'ordre n sont cyclique ssi gcd(n,phi(n))==1 */
                   1280:   if (DEBUGLEVEL >= 1)
                   1281:     timer2();
                   1282:   n = degree(T);
                   1283:   F = factor(stoi(n));
                   1284:   fc = cgetg(lg(F[1]), t_VECSMALL);
                   1285:   for (i = 1; i < lg(fc); i++)
                   1286:     fc[i] = itos(gpow(((GEN **) F)[1][i], ((GEN **) F)[2][i], DEFAULTPREC));
                   1287:   ppp = itos(((GEN **) F)[1][1]);      /* Plus Petit diviseur Premier */
                   1288:   pgp = itos(((GEN **) F)[1][lg(F[1]) - 1]);   /* Plus Grand diviseur
                   1289:                                                 * Premier */
                   1290:   pha = 1;
                   1291:   ord = 1;
                   1292:   ex = 0;
                   1293:   for (i = lg(F[1]) - 1; i > 0; i--)
                   1294:   {
                   1295:     p = itos(((GEN **) F)[1][i]);
                   1296:     if (pha % p != 0)
                   1297:     {
                   1298:       ord *= p;
                   1299:       pha *= p - 1;
                   1300:     } else
                   1301:     {
                   1302:       ex = 1;
                   1303:       break;
                   1304:     }
                   1305:     if (!gcmp1(((GEN **) F)[2][i]))
                   1306:       break;
                   1307:   }
                   1308:   plift = 0;
                   1309:   /* Etude des ordres des Frobenius */
                   1310:   nbmax = max(4, n / 2);       /* Nombre maxi d'éléments à tester */
                   1311:   nbtest = 0;
                   1312:   deg = 0;
                   1313:   for (p = 0, primepointer = diffptr; plift == 0 || (nbtest < nbmax && ord != n) || (n == 24 && exist6 == 0 && p4 == 0);)
                   1314:   {
                   1315:     long u, s, c;
                   1316:     long isram;
                   1317:     GEN S;
                   1318:     c = *primepointer++;
                   1319:     if (!c)
                   1320:       err(primer1);
                   1321:     p += c;
                   1322:     if (p <= (n << 1))
                   1323:       continue;
                   1324:     S = (GEN) simplefactmod(T, stoi(p));
                   1325:     isram = 0;
                   1326:     for (i = 1; i < lg(S[2]) && !isram; i++)
                   1327:       if (!gcmp1(((GEN **) S)[2][i]))
                   1328:        isram = 1;
                   1329:     if (isram == 0)
                   1330:     {
                   1331:       s = itos(((GEN **) S)[1][1]);
                   1332:       for (i = 2; i < lg(S[1]) && !isram; i++)
                   1333:        if (itos(((GEN **) S)[1][i]) != s)
                   1334:        {
                   1335:          avma = ltop;
                   1336:          if (DEBUGLEVEL >= 2)
                   1337:            fprintferr("GaloisAnalysis:non Galois for p=%d\n", p);
                   1338:          ga->p = p;
                   1339:          ga->deg = 0;
                   1340:          return;               /* Not a Galois polynomial */
                   1341:        }
                   1342:       if (l == 0 && s == 1)
                   1343:        l = p;
                   1344:       nbtest++;
                   1345:       if (nbtest > (3 * nbmax) && (n == 60 || n == 75))
                   1346:        break;
                   1347:       if (s % 6 == 0)
                   1348:        exist6 = 1;
                   1349:       if (p4 == 0 && s == 4)
                   1350:        p4 = p;
                   1351:       if (DEBUGLEVEL >= 6)
                   1352:        fprintferr("GaloisAnalysis:Nbtest=%d,plift=%d,p=%d,s=%d,ord=%d\n", nbtest, plift, p, s, ord);
                   1353:       if (s > ordmax)
                   1354:        ordmax = s;
                   1355:       if (s >= ord)            /* Calcul de l'exposant distinguant minimal
                   1356:                                 * garanti */
                   1357:       {
                   1358:        if (s * ppp == n)       /* Tout sous groupe d'indice ppp est distingué */
                   1359:          u = s;
                   1360:        else                    /* Théorème de structure des groupes
                   1361:                                 * hyper-résoluble */
                   1362:        {
                   1363:          u = 1;
                   1364:          for (i = lg(fc) - 1; i > 0; i--)
                   1365:          {
                   1366:            if (s % fc[i] == 0)
                   1367:              u *= fc[i];
                   1368:            else
                   1369:              break;
                   1370:          }
                   1371:        }
                   1372:        if (u != 1)
                   1373:        {
                   1374:          if (!ex || s > ord || (s == ord && (plift == 0 || u > deg)))
                   1375:          {
                   1376:            deg = u;
                   1377:            ord = s;
                   1378:            plift = p;
                   1379:            pp = primepointer;
                   1380:            ex = 1;
                   1381:          }
                   1382:        } else if (!ex && (plift == 0 || s > ord))
                   1383:          /*
                   1384:           * J'ai besoin de plus de connaissance sur les p-groupes, surtout
                   1385:           * pour p=2;
                   1386:           */
                   1387:        {
                   1388:          deg = pgp;
                   1389:          ord = s;
                   1390:          plift = p;
                   1391:          pp = primepointer;
                   1392:          ex = 0;
                   1393:        }
                   1394:       }
                   1395:     }
                   1396:   }
                   1397:   if (plift == 0 || ((n == 36 || n == 48) && !exist6) || (n == 56 && ordmax == 7) || (n == 60 && ordmax == 5) || (n == 72 && !exist6) || (n == 80 && ordmax == 5))
                   1398:   {
                   1399:     deg = 0;
                   1400:     err(warner, "Galois group almost certainly not weakly super solvable");
                   1401:   }
                   1402:   avma = ltop;
                   1403:   if (calcul_l)
                   1404:   {
                   1405:     ulong av;
                   1406:     GEN x;
                   1407:     long c;
                   1408:     x = cgetg(3, t_POLMOD);
                   1409:     x[2] = lpolx[varn(T)];
                   1410:     av = avma;
                   1411:     while (l == 0)
                   1412:     {
                   1413:       c = *primepointer++;
                   1414:       if (!c)
                   1415:        err(primer1);
                   1416:       p += c;
                   1417:       x[1] = lmul(T, gmodulss(1, p));
                   1418:       if (gegal(gpuigs(x, p), x))
                   1419:        l = p;
                   1420:       avma = av;
                   1421:     }
                   1422:     avma = ltop;
                   1423:   }
                   1424:   ga->p = plift;
                   1425:   ga->exception = ex;
                   1426:   ga->deg = deg;
                   1427:   ga->ord = ord;
                   1428:   ga->l = l;
                   1429:   ga->primepointer = pp;
                   1430:   ga->ppp = ppp;
                   1431:   ga->p4 = p4;
                   1432:   if (DEBUGLEVEL >= 4)
                   1433:     fprintferr("GaloisAnalysis:p=%d l=%d exc=%d deg=%d ord=%d ppp=%d\n", p, l, ex, deg, ord, ppp);
                   1434:   if (DEBUGLEVEL >= 1)
                   1435:     msgtimer("galoisanalysis()");
                   1436: }
                   1437: /* Calcule les bornes sur les coefficients a chercher */
                   1438: struct galois_borne
                   1439: {
                   1440:   GEN l;
                   1441:   long valsol;
                   1442:   long valabs;
                   1443:   GEN bornesol;
                   1444:   GEN ladicsol;
                   1445:   GEN ladicabs;
                   1446: };
                   1447: /*
                   1448:  * recalcule L avec une precision superieur
                   1449:  */
                   1450: GEN
                   1451: recalculeL(GEN T, GEN Tp, GEN L, struct galois_borne * gb, struct galois_borne * ngb)
                   1452: {
                   1453:   ulong ltop = avma, lbot = avma;
                   1454:   GEN L2, y, z, lad;
                   1455:   long i, val;
                   1456:   if (DEBUGLEVEL >= 1)
                   1457:     timer2();
                   1458:   L2 = cgetg(lg(L), typ(L));
                   1459:   for (i = 1; i < lg(L); i++)
                   1460:   {
                   1461:     ltop = avma;
                   1462:     val = gb->valabs;
                   1463:     lad = gb->ladicabs;
                   1464:     z = (GEN) L[i];
                   1465:     while (val < ngb->valabs)
                   1466:     {
                   1467:       val <<= 1;
                   1468:       if (val >= ngb->valabs)
                   1469:       {
                   1470:        val = ngb->valabs;
                   1471:        lad = ngb->ladicabs;
                   1472:       } else
                   1473:        lad = gsqr(lad);
                   1474:       y = gmodulcp((GEN) z[2], lad);
                   1475:       z = gdiv(poleval(T, y), poleval(Tp, y));
                   1476:       lbot = avma;
                   1477:       z = gsub(y, z);
                   1478:     }
                   1479:     L2[i] = (long) gerepile(ltop, lbot, z);
                   1480:   }
                   1481:   return L2;
                   1482: }
                   1483: void
                   1484: initborne(GEN T, GEN disc, struct galois_borne * gb, long ppp)
                   1485: {
                   1486:   ulong ltop = avma, lbot, av2;
                   1487:   GEN borne, borneroots, borneabs;
                   1488:   int i, j;
                   1489:   int n;
                   1490:   GEN L, M, z;
                   1491:   L = roots(T, DEFAULTPREC);
                   1492:   n = lg(L) - 1;
                   1493:   for (i = 1; i <= n; i++)
                   1494:   {
                   1495:     z = (GEN) L[i];
                   1496:     if (signe(z[2]))
                   1497:       break;
                   1498:     L[i] = z[1];
                   1499:   }
                   1500:   M = vandermondeinverse(L, gmul(T, realun(DEFAULTPREC)), disc);
                   1501:   borne = gzero;
                   1502:   for (i = 1; i <= n; i++)
                   1503:   {
                   1504:     z = gzero;
                   1505:     for (j = 1; j <= n; j++)
                   1506:       z = gadd(z, gabs(((GEN **) M)[j][i], DEFAULTPREC));
                   1507:     if (gcmp(z, borne) > 0)
                   1508:       borne = z;
                   1509:   }
                   1510:   borneroots = gzero;
                   1511:   for (i = 1; i <= n; i++)
                   1512:   {
                   1513:     z = gabs((GEN) L[i], DEFAULTPREC);
                   1514:     if (gcmp(z, borneroots) > 0)
                   1515:       borneroots = z;
                   1516:   }
                   1517:   borneabs = addsr(1, gpowgs(addsr(n, borneroots), n / ppp));
                   1518:   lbot = avma;
                   1519:   borneroots = addsr(1, gmul(borne, borneroots));
                   1520:   av2 = avma;
                   1521:   borneabs = gmul2n(gmul(borne, borneabs), 4);
                   1522:   gb->valsol = itos(gceil(gdiv(glog(gmul2n(borneroots, 4 + (n >> 1)), DEFAULTPREC), glog(gb->l, DEFAULTPREC))));
                   1523:   if (DEBUGLEVEL >= 4)
                   1524:     fprintferr("GaloisConj:val1=%d\n", gb->valsol);
                   1525:   gb->valabs = max(gb->valsol, itos(gceil(gdiv(glog(borneabs, DEFAULTPREC), glog(gb->l, DEFAULTPREC)))));
                   1526:   if (DEBUGLEVEL >= 4)
                   1527:     fprintferr("GaloisConj:val2=%d\n", gb->valabs);
                   1528:   avma = av2;
                   1529:   gb->bornesol = gerepile(ltop, lbot, borneroots);
                   1530:   gb->ladicsol = gpowgs(gb->l, gb->valsol);
                   1531:   gb->ladicabs = gpowgs(gb->l, gb->valabs);
                   1532: }
                   1533: /* Groupe A4 */
                   1534: GEN
                   1535: a4galoisgen(GEN T, struct galois_test * td)
                   1536: {
                   1537:   int ltop = avma, av, av2;
                   1538:   int i, j, k;
                   1539:   int n;
                   1540:   int N, hop = 0;
                   1541:   GEN *ar, **mt;               /* tired of casting */
                   1542:   GEN t, u, O;
                   1543:   GEN res, orb, ry;
                   1544:   GEN pft, pfu, pfv;
                   1545:   n = degree(T);
                   1546:   res = cgetg(4, t_VEC);
                   1547:   ry = cgetg(3, t_VEC);
                   1548:   res[1] = (long) ry;
                   1549:   pft = cgetg(n + 1, t_VECSMALL);
                   1550:   ry[1] = (long) pft;
                   1551:   ry[2] = (long) stoi(2);
                   1552:   ry = cgetg(3, t_VEC);
                   1553:   pfu = cgetg(n + 1, t_VECSMALL);
                   1554:   ry[1] = (long) pfu;
                   1555:   ry[2] = (long) stoi(2);
                   1556:   res[2] = (long) ry;
                   1557:   ry = cgetg(3, t_VEC);
                   1558:   pfv = cgetg(n + 1, t_VECSMALL);
                   1559:   ry[1] = (long) pfv;
                   1560:   ry[2] = (long) stoi(3);
                   1561:   res[3] = (long) ry;
                   1562:   av = avma;
                   1563:   ar = (GEN *) alloue_ardoise(n, 1 + lg(td->ladic));
                   1564:   mt = (GEN **) td->PV[td->ordre[n]];
                   1565:   t = cgetg(n + 1, t_VECSMALL) + 1;    /* Sorry for this hack */
                   1566:   u = cgetg(n + 1, t_VECSMALL) + 1;    /* too lazy to correct */
                   1567:   av2 = avma;
                   1568:   N = itos(gdiv(mpfact(n), mpfact(n >> 1))) >> (n >> 1);
                   1569:   if (DEBUGLEVEL >= 4)
                   1570:     fprintferr("A4GaloisConj:I will test %d permutations\n", N);
                   1571:   avma = av2;
                   1572:   for (i = 0; i < n; i++)
                   1573:     t[i] = i + 1;
                   1574:   for (i = 0; i < N; i++)
                   1575:   {
                   1576:     GEN g;
                   1577:     int a, x, y;
                   1578:     if (i == 0)
                   1579:     {
                   1580:       gaffect(gzero, ar[(n - 2) >> 1]);
                   1581:       for (k = n - 2; k > 2; k -= 2)
                   1582:        gaddz(ar[k >> 1], gadd(mt[k + 1][k + 2], mt[k + 2][k + 1]), ar[(k >> 1) - 1]);
                   1583:     } else
                   1584:     {
                   1585:       x = i;
                   1586:       y = 1;
                   1587:       do
                   1588:       {
                   1589:        hiremainder = 0;
                   1590:        y += 2;
                   1591:        x = divll(x, y);
                   1592:        a = hiremainder;
                   1593:       }
                   1594:       while (!a);
                   1595:       switch (y)
                   1596:       {
                   1597:        case 3:
                   1598:        x = t[2];
                   1599:        if (a == 1)
                   1600:        {
                   1601:          t[2] = t[1];
                   1602:          t[1] = x;
                   1603:        } else
                   1604:        {
                   1605:          t[2] = t[0];
                   1606:          t[0] = x;
                   1607:        }
                   1608:        break;
                   1609:        case 5:
                   1610:        x = t[0];
                   1611:        t[0] = t[2];
                   1612:        t[2] = t[1];
                   1613:        t[1] = x;
                   1614:        x = t[4];
                   1615:        t[4] = t[4 - a];
                   1616:        t[4 - a] = x;
                   1617:        gaddz(ar[2], gadd(mt[t[4]][t[5]], mt[t[5]][t[4]]), ar[1]);
                   1618:        break;
                   1619:        case 7:
                   1620:        x = t[0];
                   1621:        t[0] = t[4];
                   1622:        t[4] = t[3];
                   1623:        t[3] = t[1];
                   1624:        t[1] = t[2];
                   1625:        t[2] = x;
                   1626:        x = t[6];
                   1627:        t[6] = t[6 - a];
                   1628:        t[6 - a] = x;
                   1629:        gaddz(ar[3], gadd(mt[t[6]][t[7]], mt[t[7]][t[6]]), ar[2]);
                   1630:        gaddz(ar[2], gadd(mt[t[4]][t[5]], mt[t[5]][t[4]]), ar[1]);
                   1631:        break;
                   1632:        case 9:
                   1633:        x = t[0];
                   1634:        t[0] = t[6];
                   1635:        t[6] = t[5];
                   1636:        t[5] = t[3];
                   1637:        t[3] = x;
                   1638:        x = t[4];
                   1639:        t[4] = t[1];
                   1640:        t[1] = x;
                   1641:        x = t[8];
                   1642:        t[8] = t[8 - a];
                   1643:        t[8 - a] = x;
                   1644:        gaddz(ar[4], gadd(mt[t[8]][t[9]], mt[t[9]][t[8]]), ar[3]);
                   1645:        gaddz(ar[3], gadd(mt[t[6]][t[7]], mt[t[7]][t[6]]), ar[2]);
                   1646:        gaddz(ar[2], gadd(mt[t[4]][t[5]], mt[t[5]][t[4]]), ar[1]);
                   1647:        break;
                   1648:        default:
                   1649:        y--;
                   1650:        x = t[0];
                   1651:        t[0] = t[2];
                   1652:        t[2] = t[1];
                   1653:        t[1] = x;
                   1654:        for (k = 4; k < y; k += 2)
                   1655:        {
                   1656:          int j;
                   1657:          x = t[k];
                   1658:          for (j = k; j > 0; j--)
                   1659:            t[j] = t[j - 1];
                   1660:          t[0] = x;
                   1661:        }
                   1662:        x = t[y];
                   1663:        t[y] = t[y - a];
                   1664:        t[y - a] = x;
                   1665:        for (k = y; k > 2; k -= 2)
                   1666:          gaddz(ar[k >> 1], gadd(mt[t[k]][t[k + 1]], mt[t[k + 1]][t[k]]), ar[(k >> 1) - 1]);
                   1667:       }
                   1668:     }
                   1669:     g = gadd(ar[1], gadd(gadd(mt[t[0]][t[1]], mt[t[1]][t[0]]),
                   1670:                         gadd(mt[t[2]][t[3]], mt[t[3]][t[2]])));
                   1671:     if (padicisint(g, td))
                   1672:     {
                   1673:       for (k = 0; k < n; k += 2)
                   1674:       {
                   1675:        pft[t[k]] = t[k + 1];
                   1676:        pft[t[k + 1]] = t[k];
                   1677:       }
                   1678:       if (verifietest(pft, td))
                   1679:        break;
                   1680:       else
                   1681:        hop++;
                   1682:     }
                   1683:     avma = av2;
                   1684:   }
                   1685:   if (i == N)
                   1686:   {
                   1687:     avma = ltop;
                   1688:     if (DEBUGLEVEL >= 1 && hop)
                   1689:       fprintferr("A4GaloisConj: %d hop sur %d iterations\n", hop, N);
                   1690:     return gzero;
                   1691:   }
                   1692:   if (DEBUGLEVEL >= 1 && hop)
                   1693:     fprintferr("A4GaloisConj: %d hop sur %d iterations\n", hop, N);
                   1694:   N = itos(gdiv(mpfact(n >> 1), mpfact(n >> 2))) >> 1;
                   1695:   avma = av2;
                   1696:   if (DEBUGLEVEL >= 4)
                   1697:     fprintferr("A4GaloisConj:sigma=%Z \n", pft);
                   1698:   for (i = 0; i < N; i++)
                   1699:   {
                   1700:     GEN g;
                   1701:     int a, x, y;
                   1702:     if (i == 0)
                   1703:     {
                   1704:       for (k = 0; k < n; k += 4)
                   1705:       {
                   1706:        u[k + 3] = t[k + 3];
                   1707:        u[k + 2] = t[k + 1];
                   1708:        u[k + 1] = t[k + 2];
                   1709:        u[k] = t[k];
                   1710:       }
                   1711:     } else
                   1712:     {
                   1713:       x = i;
                   1714:       y = -2;
                   1715:       do
                   1716:       {
                   1717:        hiremainder = 0;
                   1718:        y += 4;
                   1719:        x = divll(x, y);
                   1720:        a = hiremainder;
                   1721:       } while (!a);
                   1722:       x = u[2];
                   1723:       u[2] = u[0];
                   1724:       u[0] = x;
                   1725:       switch (y)
                   1726:       {
                   1727:        case 2:
                   1728:        break;
                   1729:        case 6:
                   1730:        x = u[4];
                   1731:        u[4] = u[6];
                   1732:        u[6] = x;
                   1733:        if (a % 2 == 0)
                   1734:        {
                   1735:          a = 4 - a / 2;
                   1736:          x = u[6];
                   1737:          u[6] = u[a];
                   1738:          u[a] = x;
                   1739:          x = u[4];
                   1740:          u[4] = u[a - 2];
                   1741:          u[a - 2] = x;
                   1742:        }
                   1743:        break;
                   1744:        case 10:
                   1745:        x = u[6];
                   1746:        u[6] = u[3];
                   1747:        u[3] = u[2];
                   1748:        u[2] = u[4];
                   1749:        u[4] = u[1];
                   1750:        u[1] = u[0];
                   1751:        u[0] = x;
                   1752:        if (a >= 3)
                   1753:          a += 2;
                   1754:        a = 8 - a;
                   1755:        x = u[10];
                   1756:        u[10] = u[a];
                   1757:        u[a] = x;
                   1758:        x = u[8];
                   1759:        u[8] = u[a - 2];
                   1760:        u[a - 2] = x;
                   1761:        break;
                   1762:       }
                   1763:     }
                   1764:     g = gzero;
                   1765:     for (k = 0; k < n; k += 2)
                   1766:       g = gadd(g, gadd(mt[u[k]][u[k + 1]], mt[u[k + 1]][u[k]]));
                   1767:     if (padicisint(g, td))
                   1768:     {
                   1769:       for (k = 0; k < n; k += 2)
                   1770:       {
                   1771:        pfu[u[k]] = u[k + 1];
                   1772:        pfu[u[k + 1]] = u[k];
                   1773:       }
                   1774:       if (verifietest(pfu, td))
                   1775:        break;
                   1776:       else
                   1777:        hop++;
                   1778:     }
                   1779:     avma = av2;
                   1780:   }
                   1781:   if (i == N)
                   1782:   {
                   1783:     avma = ltop;
                   1784:     return gzero;
                   1785:   }
                   1786:   if (DEBUGLEVEL >= 1 && hop)
                   1787:     fprintferr("A4GaloisConj: %d hop sur %d iterations\n", hop, N);
                   1788:   if (DEBUGLEVEL >= 4)
                   1789:     fprintferr("A4GaloisConj:tau=%Z \n", u);
                   1790:   avma = av2;
                   1791:   orb = cgetg(3, t_VEC);
                   1792:   orb[1] = (long) pft;
                   1793:   orb[2] = (long) pfu;
                   1794:   if (DEBUGLEVEL >= 4)
                   1795:     fprintferr("A4GaloisConj:orb=%Z \n", orb);
                   1796:   O = vecpermorbite(orb);
                   1797:   if (DEBUGLEVEL >= 4)
                   1798:     fprintferr("A4GaloisConj:O=%Z \n", O);
                   1799:   av2 = avma;
                   1800:   for (j = 0; j < 2; j++)
                   1801:   {
                   1802:     pfv[((long **) O)[1][1]] = ((long **) O)[2][1];
                   1803:     pfv[((long **) O)[1][2]] = ((long **) O)[2][3 + j];
                   1804:     pfv[((long **) O)[1][3]] = ((long **) O)[2][4 - (j << 1)];
                   1805:     pfv[((long **) O)[1][4]] = ((long **) O)[2][2 + j];
                   1806:     for (i = 0; i < 4; i++)
                   1807:     {
                   1808:       long x;
                   1809:       GEN g;
                   1810:       switch (i)
                   1811:       {
                   1812:        case 0:
                   1813:        break;
                   1814:        case 1:
                   1815:        x = ((long **) O)[3][1];
                   1816:        ((long **) O)[3][1] = ((long **) O)[3][2];
                   1817:        ((long **) O)[3][2] = x;
                   1818:        x = ((long **) O)[3][3];
                   1819:        ((long **) O)[3][3] = ((long **) O)[3][4];
                   1820:        ((long **) O)[3][4] = x;
                   1821:        break;
                   1822:        case 2:
                   1823:        x = ((long **) O)[3][1];
                   1824:        ((long **) O)[3][1] = ((long **) O)[3][4];
                   1825:        ((long **) O)[3][4] = x;
                   1826:        x = ((long **) O)[3][2];
                   1827:        ((long **) O)[3][2] = ((long **) O)[3][3];
                   1828:        ((long **) O)[3][3] = x;
                   1829:        break;
                   1830:        case 3:
                   1831:        x = ((long **) O)[3][1];
                   1832:        ((long **) O)[3][1] = ((long **) O)[3][2];
                   1833:        ((long **) O)[3][2] = x;
                   1834:        x = ((long **) O)[3][3];
                   1835:        ((long **) O)[3][3] = ((long **) O)[3][4];
                   1836:        ((long **) O)[3][4] = x;
                   1837:       }
                   1838:       pfv[((long **) O)[2][1]] = ((long **) O)[3][1];
                   1839:       pfv[((long **) O)[2][3 + j]] = ((long **) O)[3][4 - j];
                   1840:       pfv[((long **) O)[2][4 - (j << 1)]] = ((long **) O)[3][2 + (j << 1)];
                   1841:       pfv[((long **) O)[2][2 + j]] = ((long **) O)[3][3 - j];
                   1842:       pfv[((long **) O)[3][1]] = ((long **) O)[1][1];
                   1843:       pfv[((long **) O)[3][4 - j]] = ((long **) O)[1][2];
                   1844:       pfv[((long **) O)[3][2 + (j << 1)]] = ((long **) O)[1][3];
                   1845:       pfv[((long **) O)[3][3 - j]] = ((long **) O)[1][4];
                   1846:       g = gzero;
                   1847:       for (k = 1; k <= n; k++)
                   1848:        g = gadd(g, mt[k][pfv[k]]);
                   1849:       if (padicisint(g, td) && verifietest(pfv, td))
                   1850:       {
                   1851:        avma = av;
                   1852:        if (DEBUGLEVEL >= 1)
                   1853:          fprintferr("A4GaloisConj:%d hop sur %d iterations max\n", hop, 10395 + 68);
                   1854:        return res;
                   1855:       } else
                   1856:        hop++;
                   1857:       avma = av2;
                   1858:     }
                   1859:   }
                   1860:   /* Echec? */
                   1861:   avma = ltop;
                   1862:   return gzero;
                   1863: }
                   1864: /* Groupe S4 */
                   1865: GEN
                   1866: Fpisom(GEN p, GEN T1, GEN T2)
                   1867: {
                   1868:   ulong ltop = avma, lbot;
                   1869:   GEN T, res;
                   1870:   long v;
                   1871:   if (T1 == T2)
                   1872:     return gmodulcp(polx[varn(T1)], T1);
                   1873:   v = varn(T2);
                   1874:   setvarn(T2, MAXVARN);
                   1875:   T = (GEN) factmod9(T1, p, T2)[1];
                   1876:   setvarn(T2, v);
                   1877:   lbot = avma;
                   1878:   res = gneg(((GEN **) T)[1][2]);
                   1879:   setvarn(res[1], v);
                   1880:   setvarn(res[2], v);
                   1881:   return gerepile(ltop, lbot, res);
                   1882: }
                   1883: GEN
                   1884: Fpinvisom(GEN S)
                   1885: {
                   1886:   ulong ltop = avma, lbot;
                   1887:   GEN M, U, V;
                   1888:   int n, i, j;
                   1889:   long x;
                   1890:   x = varn(S[1]);
                   1891:   U = gmodulcp(polun[x], (GEN) S[1]);
                   1892:   n = degree((GEN) S[1]);
                   1893:   M = cgetg(n + 1, t_MAT);
                   1894:   for (i = 1; i <= n; i++)
                   1895:   {
                   1896:     int d;
                   1897:     if (i > 1)
                   1898:       U = gmul(U, S);
                   1899:     M[i] = lgetg(n + 1, t_COL);
                   1900:     d = degree((GEN) U[2]);
                   1901:     for (j = 1; j <= d + 1; j++)
                   1902:       ((GEN *) M)[i][j] = ((GEN *) U)[2][1 + j];
                   1903:     for (j = d + 2; j <= n; j++)
                   1904:       ((GEN *) M)[i][j] = zero;
                   1905:   }
                   1906:   V = cgetg(n + 1, t_COL);
                   1907:   for (i = 1; i <= n; i++)
                   1908:     V[i] = zero;
                   1909:   V[2] = un;
                   1910:   V = gauss(M, V);
                   1911:   lbot = avma;
                   1912:   V = gtopolyrev(V, x);
                   1913:   return gerepile(ltop, lbot, V);
                   1914: }
                   1915: static void
                   1916: makelift(GEN u, struct galois_lift * gl, GEN liftpow)
                   1917: {
                   1918:   int i;
                   1919:   liftpow[1] = (long) automorphismlift(u, gl);
                   1920:   for (i = 2; i < lg(liftpow); i++)
                   1921:     liftpow[i] = lmul((GEN) liftpow[i - 1], (GEN) liftpow[1]);
                   1922: }
                   1923: static long
                   1924: s4test(GEN u, GEN liftpow, struct galois_lift * gl, GEN phi)
                   1925: {
                   1926:   GEN res;
                   1927:   int i;
                   1928:   if (DEBUGLEVEL >= 6)
                   1929:     timer2();
                   1930:   u = (GEN) u[2];
                   1931:   res = (GEN) u[2];
                   1932:   for (i = 1; i < lgef(u) - 2; i++)
                   1933:     res = gadd(res, gmul((GEN) liftpow[i], (GEN) u[i + 2]));
                   1934:   res = gdiv(centerlift(gmul((GEN) res[2], gl->den)), gl->den);
                   1935:   if (DEBUGLEVEL >= 6)
                   1936:     msgtimer("s4test()");
                   1937:   return poltopermtest(res, gl->L, phi);
                   1938: }
                   1939: GEN
                   1940: s4galoisgen(struct galois_lift * gl)
                   1941: {
                   1942:   struct galois_testlift gt;
                   1943:   ulong ltop = avma, av, ltop2;
                   1944:   GEN Tmod, isom, isominv, misom;
                   1945:   int i, j;
                   1946:   GEN sg;
                   1947:   GEN sigma, tau, phi;
                   1948:   GEN res, ry;
                   1949:   GEN pj;
                   1950:   GEN p;
                   1951:   GEN bezoutcoeff, pauto, liftpow;
                   1952:   long v;
                   1953:   p = gl->p;
                   1954:   v = varn(gl->T);
                   1955:   res = cgetg(5, t_VEC);
                   1956:   ry = cgetg(3, t_VEC);
                   1957:   res[1] = (long) ry;
                   1958:   ry[1] = lgetg(lg(gl->L), t_VECSMALL);
                   1959:   ry[2] = (long) stoi(2);
                   1960:   ry = cgetg(3, t_VEC);
                   1961:   res[2] = (long) ry;
                   1962:   ry[1] = lgetg(lg(gl->L), t_VECSMALL);
                   1963:   ry[2] = (long) stoi(2);
                   1964:   ry = cgetg(3, t_VEC);
                   1965:   res[3] = (long) ry;
                   1966:   ry[1] = lgetg(lg(gl->L), t_VECSMALL);
                   1967:   ry[2] = (long) stoi(3);
                   1968:   ry = cgetg(3, t_VEC);
                   1969:   res[4] = (long) ry;
                   1970:   ry[1] = lgetg(lg(gl->L), t_VECSMALL);
                   1971:   ry[2] = (long) stoi(2);
                   1972:   ltop2 = avma;
                   1973:   sg = cgetg(7, t_VECSMALL);
                   1974:   pj = cgetg(7, t_VECSMALL);
                   1975:   sigma = cgetg(lg(gl->L), t_VECSMALL);
                   1976:   tau = cgetg(lg(gl->L), t_VECSMALL);
                   1977:   phi = cgetg(lg(gl->L), t_VECSMALL);
                   1978:   for (i = 1; i < lg(sg); i++)
                   1979:     sg[i] = i;
                   1980:   Tmod = (GEN) factmod(gl->T, p)[1];
                   1981:   isom = cgetg(lg(Tmod), t_VEC);
                   1982:   isominv = cgetg(lg(Tmod), t_VEC);
                   1983:   misom = cgetg(lg(Tmod), t_MAT);
                   1984:   inittestlift(Tmod, 1, gl, &gt, NULL);
                   1985:   bezoutcoeff = gt.bezoutcoeff;
                   1986:   pauto = gt.pauto;
                   1987:   for (i = 1; i < lg(pj); i++)
                   1988:     pj[i] = 0;
                   1989:   for (i = 1; i < lg(isom); i++)
                   1990:   {
                   1991:     misom[i] = lgetg(lg(Tmod), t_COL);
                   1992:     isom[i] = (long) Fpisom(p, (GEN) Tmod[1], (GEN) Tmod[i]);
                   1993:     if (DEBUGLEVEL >= 4)
                   1994:       fprintferr("S4GaloisConj:Computing isomorphisms %d:%Z\n", i, (GEN) isom[i]);
                   1995:     isominv[i] = (long) lift(Fpinvisom((GEN) isom[i]));
                   1996:   }
                   1997:   for (i = 1; i < lg(isom); i++)
                   1998:     for (j = 1; j < lg(isom); j++)
                   1999:       ((GEN **) misom)[i][j] = gsubst((GEN) isominv[i], v, (GEN) isom[j]);
                   2000:   liftpow = cgetg(24, t_VEC);
                   2001:   av = avma;
                   2002:   for (i = 0; i < 3; i++)
                   2003:   {
                   2004:     ulong av2;
                   2005:     GEN u;
                   2006:     int j1, j2, j3;
                   2007:     ulong avm1, avm2;
                   2008:     GEN u1, u2, u3;
                   2009:     if (i)
                   2010:     {
                   2011:       int x;
                   2012:       x = sg[3];
                   2013:       if (i == 1)
                   2014:       {
                   2015:        sg[3] = sg[2];
                   2016:        sg[2] = x;
                   2017:       } else
                   2018:       {
                   2019:        sg[3] = sg[1];
                   2020:        sg[1] = x;
                   2021:       }
                   2022:     }
                   2023:     u = chinois(((GEN **) misom)[sg[1]][sg[2]],
                   2024:                ((GEN **) misom)[sg[2]][sg[1]]);
                   2025:     u = chinois(u, ((GEN **) misom)[sg[3]][sg[4]]);
                   2026:     u = chinois(u, ((GEN **) misom)[sg[4]][sg[3]]);
                   2027:     u = chinois(u, ((GEN **) misom)[sg[5]][sg[6]]);
                   2028:     u = chinois(u, ((GEN **) misom)[sg[6]][sg[5]]);
                   2029:     makelift(u, gl, liftpow);
                   2030:     av2 = avma;
                   2031:     for (j1 = 0; j1 < 4; j1++)
                   2032:     {
                   2033:       u1 = gadd(gmul((GEN) bezoutcoeff[sg[5]], (GEN) pauto[1 + j1]),
                   2034:              gmul((GEN) bezoutcoeff[sg[6]], (GEN) pauto[((-j1) & 3) + 1]));
                   2035:       avm1 = avma;
                   2036:       for (j2 = 0; j2 < 4; j2++)
                   2037:       {
                   2038:        u2 = gadd(u1, gmul((GEN) bezoutcoeff[sg[3]], (GEN) pauto[1 + j2]));
                   2039:        u2 = gadd(u2, gmul((GEN) bezoutcoeff[sg[4]], (GEN) pauto[((-j2) & 3) + 1]));
                   2040:        avm2 = avma;
                   2041:        for (j3 = 0; j3 < 4; j3++)
                   2042:        {
                   2043:          u3 = gadd(u2, gmul((GEN) bezoutcoeff[sg[1]], (GEN) pauto[1 + j3]));
                   2044:          u3 = gadd(u3, gmul((GEN) bezoutcoeff[sg[2]], (GEN) pauto[((-j3) & 3) + 1]));
                   2045:          if (DEBUGLEVEL >= 4)
                   2046:            fprintferr("S4GaloisConj:Testing %d/3:%d/4:%d/4:%d/4:%Z\n", i, j1, j2, j3, sg);
                   2047:          if (s4test(u3, liftpow, gl, sigma))
                   2048:          {
                   2049:            pj[1] = j3;
                   2050:            pj[2] = j2, pj[3] = j1;
                   2051:            goto suites4;
                   2052:          }
                   2053:          avma = avm2;
                   2054:        }
                   2055:        avma = avm1;
                   2056:       }
                   2057:       avma = av2;
                   2058:     }
                   2059:     avma = av;
                   2060:   }
                   2061:   avma = ltop;
                   2062:   return gzero;
                   2063: suites4:
                   2064:   if (DEBUGLEVEL >= 4)
                   2065:     fprintferr("S4GaloisConj:sigma=%Z\n", sigma);
                   2066:   if (DEBUGLEVEL >= 4)
                   2067:     fprintferr("S4GaloisConj:pj=%Z\n", pj);
                   2068:   avma = av;
                   2069:   for (j = 1; j <= 3; j++)
                   2070:   {
                   2071:     ulong av2;
                   2072:     GEN u;
                   2073:     int w, l;
                   2074:     int z;
                   2075:     z = sg[1];
                   2076:     sg[1] = sg[3];
                   2077:     sg[3] = sg[5];
                   2078:     sg[5] = z;
                   2079:     z = sg[2];
                   2080:     sg[2] = sg[4];
                   2081:     sg[4] = sg[6];
                   2082:     sg[6] = z;
                   2083:     z = pj[1];
                   2084:     pj[1] = pj[2];
                   2085:     pj[2] = pj[3];
                   2086:     pj[3] = z;
                   2087:     for (l = 0; l < 2; l++)
                   2088:     {
                   2089:       u = chinois(((GEN **) misom)[sg[1]][sg[3]],
                   2090:                  ((GEN **) misom)[sg[3]][sg[1]]);
                   2091:       u = chinois(u, ((GEN **) misom)[sg[2]][sg[4]]);
                   2092:       u = chinois(u, ((GEN **) misom)[sg[4]][sg[2]]);
                   2093:       u = chinois(u, ((GEN **) misom)[sg[5]][sg[6]]);
                   2094:       u = chinois(u, ((GEN **) misom)[sg[6]][sg[5]]);
                   2095:       makelift(u, gl, liftpow);
                   2096:       av2 = avma;
                   2097:       for (w = 0; w < 4; w += 2)
                   2098:       {
                   2099:        GEN uu;
                   2100:        long av3;
                   2101:        pj[6] = (w + pj[3]) & 3;
                   2102:        uu = gadd(gmul((GEN) bezoutcoeff[sg[5]], (GEN) pauto[(pj[6] & 3) + 1]),
                   2103:           gmul((GEN) bezoutcoeff[sg[6]], (GEN) pauto[((-pj[6]) & 3) + 1]));
                   2104:        av3 = avma;
                   2105:        for (i = 0; i < 4; i++)
                   2106:        {
                   2107:          GEN u;
                   2108:          pj[4] = i;
                   2109:          pj[5] = (i + pj[2] - pj[1]) & 3;
                   2110:          if (DEBUGLEVEL >= 4)
                   2111:            fprintferr("S4GaloisConj:Testing %d/3:%d/2:%d/2:%d/4:%Z:%Z\n", j - 1, w >> 1, l, i, sg, pj);
                   2112:          u = gadd(uu, gmul((GEN) pauto[(pj[4] & 3) + 1], (GEN) bezoutcoeff[sg[1]]));
                   2113:          u = gadd(u, gmul((GEN) pauto[((-pj[4]) & 3) + 1], (GEN) bezoutcoeff[sg[3]]));
                   2114:          u = gadd(u, gmul((GEN) pauto[(pj[5] & 3) + 1], (GEN) bezoutcoeff[sg[2]]));
                   2115:          u = gadd(u, gmul((GEN) pauto[((-pj[5]) & 3) + 1], (GEN) bezoutcoeff[sg[4]]));
                   2116:          if (s4test(u, liftpow, gl, tau))
                   2117:            goto suites4_2;
                   2118:          avma = av3;
                   2119:        }
                   2120:        avma = av2;
                   2121:       }
                   2122:       z = sg[4];
                   2123:       sg[4] = sg[3];
                   2124:       sg[3] = z;
                   2125:       pj[2] = (-pj[2]) & 3;
                   2126:       avma = av;
                   2127:     }
                   2128:   }
                   2129:   avma = ltop;
                   2130:   return gzero;
                   2131: suites4_2:
                   2132:   avma = av;
                   2133:   {
                   2134:     int abc, abcdef;
                   2135:     GEN u;
                   2136:     ulong av2;
                   2137:     abc = (pj[1] + pj[2] + pj[3]) & 3;
                   2138:     abcdef = (((abc + pj[4] + pj[5] - pj[6]) & 3) >> 1);
                   2139:     u = chinois(((GEN **) misom)[sg[1]][sg[4]],
                   2140:                ((GEN **) misom)[sg[4]][sg[1]]);
                   2141:     u = chinois(u, ((GEN **) misom)[sg[2]][sg[5]]);
                   2142:     u = chinois(u, ((GEN **) misom)[sg[5]][sg[2]]);
                   2143:     u = chinois(u, ((GEN **) misom)[sg[3]][sg[6]]);
                   2144:     u = chinois(u, ((GEN **) misom)[sg[6]][sg[3]]);
                   2145:     makelift(u, gl, liftpow);
                   2146:     av2 = avma;
                   2147:     for (j = 0; j < 8; j++)
                   2148:     {
                   2149:       int h, g, i;
                   2150:       h = j & 3;
                   2151:       g = abcdef + ((j & 4) >> 1);
                   2152:       i = h + abc - g;
                   2153:       u = gmul((GEN) pauto[(g & 3) + 1], (GEN) bezoutcoeff[sg[1]]);
                   2154:       u = gadd(u, gmul((GEN) pauto[((-g) & 3) + 1], (GEN) bezoutcoeff[sg[4]]));
                   2155:       u = gadd(u, gmul((GEN) pauto[(h & 3) + 1], (GEN) bezoutcoeff[sg[2]]));
                   2156:       u = gadd(u, gmul((GEN) pauto[((-h) & 3) + 1], (GEN) bezoutcoeff[sg[5]]));
                   2157:       u = gadd(u, gmul((GEN) pauto[(i & 3) + 1], (GEN) bezoutcoeff[sg[3]]));
                   2158:       u = gadd(u, gmul((GEN) pauto[((-i) & 3) + 1], (GEN) bezoutcoeff[sg[6]]));
                   2159:       if (DEBUGLEVEL >= 4)
                   2160:        fprintferr("S4GaloisConj:Testing %d/8 %d:%d:%d\n", j, g & 3, h & 3, i & 3);
                   2161:       if (s4test(u, liftpow, gl, phi))
                   2162:        break;
                   2163:       avma = av2;
                   2164:     }
                   2165:   }
                   2166:   if (j == 8)
                   2167:   {
                   2168:     avma = ltop;
                   2169:     return gzero;
                   2170:   }
                   2171:   for (i = 1; i < lg(gl->L); i++)
                   2172:   {
                   2173:     ((GEN **) res)[1][1][i] = sigma[tau[i]];
                   2174:     ((GEN **) res)[2][1][i] = phi[sigma[tau[phi[i]]]];
                   2175:     ((GEN **) res)[3][1][i] = phi[sigma[i]];
                   2176:     ((GEN **) res)[4][1][i] = sigma[i];
                   2177:   }
                   2178:   avma = ltop2;
                   2179:   return res;
                   2180: }
                   2181: GEN
                   2182: galoisgen(GEN T, GEN L, GEN M, GEN den, struct galois_borne * gb, const struct galois_analysis * ga)
                   2183: {
                   2184:   struct galois_analysis Pga;
                   2185:   struct galois_borne Pgb;
                   2186:   struct galois_test td;
                   2187:   struct galois_testlift gt;
                   2188:   struct galois_lift gl;
                   2189:   ulong ltop = avma, lbot;
                   2190:   long n, p, deg, ex;
                   2191:   byteptr primepointer;
                   2192:   long sg, Pm = 0, fp;
                   2193:   long x = varn(T);
                   2194:   int i, j;
                   2195:   GEN Tmod, res, pf = gzero, split, psi, ip, mod, ppsi;
                   2196:   GEN frob;
                   2197:   GEN O;
                   2198:   GEN P, PG, PL, Pden, PM, Pmod, Pp;
                   2199:   GEN *lo;                     /* tired of casting */
                   2200:   n = degree(T);
                   2201:   if (!ga->deg)
                   2202:     return gzero;
                   2203:   p = ga->p;
                   2204:   ex = ga->exception;
                   2205:   deg = ga->deg;
                   2206:   primepointer = ga->primepointer;
                   2207:   if (DEBUGLEVEL >= 9)
                   2208:     fprintferr("GaloisConj:denominator:%Z\n", den);
                   2209:   if (n == 12 && ga->ord == 3) /* A4 is very probable,so test it first */
                   2210:   {
                   2211:     long av = avma;
                   2212:     if (DEBUGLEVEL >= 4)
                   2213:       fprintferr("GaloisConj:Testing A4 first\n");
                   2214:     inittest(L, M, gb->bornesol, gb->ladicsol, &td);
                   2215:     lbot = avma;
                   2216:     PG = a4galoisgen(T, &td);
                   2217:     freetest(&td);
                   2218:     if (PG != gzero)
                   2219:     {
                   2220:       return gerepile(ltop, lbot, PG);
                   2221:     }
                   2222:     avma = av;
                   2223:   }
                   2224:   if (n == 24 && ga->ord == 3) /* S4 is very probable,so test it first */
                   2225:   {
                   2226:     long av = avma;
                   2227:     if (DEBUGLEVEL >= 4)
                   2228:       fprintferr("GaloisConj:Testing S4 first\n");
                   2229:     lbot = avma;
                   2230:     initlift(T, den, stoi(ga->p4), gb->bornesol, L, &gl);
                   2231:     PG = s4galoisgen(&gl);
                   2232:     if (PG != gzero)
                   2233:     {
                   2234:       return gerepile(ltop, lbot, PG);
                   2235:     }
                   2236:     avma = av;
                   2237:   }
                   2238:   frob = cgetg(lg(L), t_VECSMALL);
                   2239:   for (;;)
                   2240:   {
                   2241:     ulong av = avma;
                   2242:     long isram;
                   2243:     long c;
                   2244:     ip = stoi(p);
                   2245:     Tmod = (GEN) factmod(T, ip);
                   2246:     isram = 0;
                   2247:     for (i = 1; i < lg(Tmod[2]) && !isram; i++)
                   2248:       if (!gcmp1(((GEN **) Tmod)[2][i]))
                   2249:        isram = 1;
                   2250:     if (isram == 0)
                   2251:     {
                   2252:       fp = degree(((GEN **) Tmod)[1][1]);
                   2253:       for (i = 2; i < lg(Tmod[1]); i++)
                   2254:        if (degree(((GEN **) Tmod)[1][i]) != fp)
                   2255:        {
                   2256:          avma = ltop;
                   2257:          return gzero;         /* Not Galois polynomial */
                   2258:        }
                   2259:       if (fp % deg == 0)
                   2260:       {
                   2261:        if (DEBUGLEVEL >= 4)
                   2262:          fprintferr("Galoisconj:p=%d deg=%d fp=%d\n", p, deg, fp);
                   2263:        lo = (GEN *) listsousgroupes(deg, n / fp);
                   2264:        initlift(T, den, ip, gb->bornesol, L, &gl);
                   2265:        if (inittestlift((GEN) Tmod[1], fp / deg, &gl, &gt, frob))
                   2266:        {
                   2267:          int k;
                   2268:          sg = lg(lo) - 1;
                   2269:          psi = cgetg(gt.g + 1, t_VECSMALL);
                   2270:          for (k = 1; k <= gt.g; k++)
                   2271:            psi[k] = 1;
                   2272:          goto suite;
                   2273:        }
                   2274:        if (DEBUGLEVEL >= 4)
                   2275:          fprintferr("Galoisconj:Subgroups list:%Z\n", (GEN) lo);
                   2276:        if (deg == fp && cgcd(deg, n / deg) == 1)       /* Pretty sure it's the
                   2277:                                                         * biggest ? */
                   2278:        {
                   2279:          for (sg = 1; sg < lg(lo) - 1; sg++)
                   2280:            if (frobeniusliftall(lo[sg], &psi, &gl, &gt, frob))
                   2281:              goto suite;
                   2282:        } else
                   2283:        {
                   2284:          for (sg = lg(lo) - 2; sg >= 1; sg--)  /* else start with the
                   2285:                                                 * fastest! */
                   2286:            if (frobeniusliftall(lo[sg], &psi, &gl, &gt, frob))
                   2287:              goto suite;
                   2288:        }
                   2289:        if (ex == 1 && (n == 12 || n % 12 != 0))
                   2290:        {
                   2291:          avma = ltop;
                   2292:          return gzero;
                   2293:        } else
                   2294:        {
                   2295:          ex--;
                   2296:          if (n % 12 == 0)
                   2297:          {
                   2298:            if (n % 36 == 0)
                   2299:              deg = 5 - deg;
                   2300:            else
                   2301:              deg = 2;          /* For group like Z/2ZxA4 */
                   2302:          }
                   2303:          if (n % 12 == 0 && ex == -5)
                   2304:            err(warner, "galoisconj _may_ hang up for this polynomial");
                   2305:        }
                   2306:       }
                   2307:     }
                   2308:     c = *primepointer++;
                   2309:     if (!c)
                   2310:       err(primer1);
                   2311:     p += c;
                   2312:     if (DEBUGLEVEL >= 4)
                   2313:       fprintferr("GaloisConj:next p=%d\n", p);
                   2314:     avma = av;
                   2315:   }
                   2316: suite:                         /* Djikstra probably hates me. (Linus
                   2317:                                 * Torvalds linux/kernel/sched.c) */
                   2318:   if (deg == n)                        /* Cyclique */
                   2319:   {
                   2320:     lbot = avma;
                   2321:     res = cgetg(2, t_VEC);
                   2322:     res[1] = lgetg(3, t_VEC);
                   2323:     ((GEN **) res)[1][1] = gcopy(frob);
                   2324:     ((GEN **) res)[1][2] = stoi(deg);
                   2325:     return gerepile(ltop, lbot, res);
                   2326:   }
                   2327:   if (DEBUGLEVEL >= 9)
                   2328:     fprintferr("GaloisConj:Frobenius:%Z\n", frob);
                   2329:   O = permtocycle(frob);
                   2330:   if (DEBUGLEVEL >= 9)
                   2331:     fprintferr("GaloisConj:Orbite:%Z\n", O);
                   2332:   {
                   2333:     GEN S, Tp, Fi, Sp;
                   2334:     long gp = n / fp;
                   2335:     ppsi = cgetg(gp + 1, t_VECSMALL);
                   2336:     mod = gmodulcp(gun, ip);
                   2337:     P = corpsfixeorbitemod(L, O, x, mod, gb->l, &PL);
                   2338:     S = corpsfixeinclusion(O, PL);
                   2339:     S = vectopol(S, L, M, den, x);
                   2340:     if (DEBUGLEVEL >= 6)
                   2341:       fprintferr("GaloisConj:Inclusion:%Z\n", S);
                   2342:     Pmod = (GEN) factmod(P, ip)[1];
                   2343:     Tp = gmul(T, mod);
                   2344:     Sp = gmul(S, mod);
                   2345:     Pp = gmul(P, mod);
                   2346:     if (DEBUGLEVEL >= 4)
                   2347:       fprintferr("GaloisConj:psi=%Z\n", psi);
                   2348:     if (DEBUGLEVEL >= 4)
                   2349:       fprintferr("GaloisConj:Pmod=%Z\n", Pmod);
                   2350:     if (DEBUGLEVEL >= 4)
                   2351:       fprintferr("GaloisConj:Tmod=%Z\n", Tmod);
                   2352:     for (i = 1; i <= gp; i++)
                   2353:     {
                   2354:       Fi = ggcd(Tp, gsubst((GEN) Pmod[i], x, Sp));
                   2355:       Fi = gdiv(Fi, (GEN) Fi[lgef(Fi) - 1]);
                   2356:       if (DEBUGLEVEL >= 6)
                   2357:        fprintferr("GaloisConj:Fi=%Z  %d", Fi, i);
                   2358:       for (j = 1; j <= gp; j++)
                   2359:        if (gegal(Fi, ((GEN **) Tmod)[1][j]))
                   2360:          break;
                   2361:       if (DEBUGLEVEL >= 6)
                   2362:        fprintferr("-->%d\n", j);
                   2363:       if (j == gp + 1)
                   2364:       {
                   2365:        avma = ltop;
                   2366:        return gzero;
                   2367:       }
                   2368:       if (j == gp)
                   2369:       {
                   2370:        Pm = i;
                   2371:        ppsi[i] = 1;
                   2372:       } else
                   2373:        ppsi[i] = psi[j];
                   2374:     }
                   2375:     if (DEBUGLEVEL >= 4)
                   2376:       fprintferr("GaloisConj:Pm=%d   ppsi=%Z\n", Pm, ppsi);
                   2377:     galoisanalysis(P, &Pga, 0);
                   2378:     if (Pga.deg == 0)
                   2379:     {
                   2380:       avma = ltop;
                   2381:       return gzero;            /* Avoid computing the discriminant */
                   2382:     }
                   2383:     Pden = gabs(mycoredisc(discsr(P)), DEFAULTPREC);
                   2384:     Pgb.l = gb->l;
                   2385:     initborne(P, Pden, &Pgb, Pga.ppp);
                   2386:     if (Pgb.valabs > gb->valabs)
                   2387:       PL = recalculeL(P, deriv(P, x), PL, gb, &Pgb);
                   2388:     PM = vandermondeinverse(PL, P, Pden);
                   2389:     PG = galoisgen(P, PL, PM, Pden, &Pgb, &Pga);
                   2390:   }
                   2391:   if (DEBUGLEVEL >= 4)
                   2392:     fprintferr("GaloisConj:Retour sur Terre:%Z\n", PG);
                   2393:   if (PG == gzero)
                   2394:   {
                   2395:     avma = ltop;
                   2396:     return gzero;
                   2397:   }
                   2398:   inittest(L, M, gb->bornesol, gb->ladicsol, &td);
                   2399:   split = splitorbite(O);
                   2400:   lbot = avma;
                   2401:   res = cgetg(lg(PG) + lg(split) - 1, t_VEC);
                   2402:   for (i = 1; i < lg(split); i++)
                   2403:     res[i] = lcopy((GEN) split[i]);
                   2404:   for (j = 1; j < lg(PG); j++)
                   2405:   {
                   2406:     ulong lbot = 0, av = avma;
                   2407:     GEN B, tau;
                   2408:     long t, g;
                   2409:     long w, sr, dss;
                   2410:     if (DEBUGLEVEL >= 6)
                   2411:       fprintferr("GaloisConj:G[%d][1]=%Z\n", j, ((GEN **) PG)[j][1]);
                   2412:     B = permtocycle(((GEN **) PG)[j][1]);
                   2413:     if (DEBUGLEVEL >= 6)
                   2414:       fprintferr("GaloisConj:B=%Z\n", B);
                   2415:     tau = gmul(mod, permtopol(((GEN **) PG)[j][1], PL, PM, Pden, x));
                   2416:     tau = gsubst((GEN) Pmod[Pm], x, tau);
                   2417:     tau = ggcd(Pp, tau);
                   2418:     tau = gdiv(tau, (GEN) tau[lgef(tau) - 1]);
                   2419:     if (DEBUGLEVEL >= 6)
                   2420:       fprintferr("GaloisConj:tau=%Z\n", tau);
                   2421:     for (g = 1; g < lg(Pmod); g++)
                   2422:       if (gegal(tau, (GEN) Pmod[g]))
                   2423:        break;
                   2424:     if (g == lg(Pmod))
                   2425:     {
                   2426:       freetest(&td);
                   2427:       avma = ltop;
                   2428:       return gzero;
                   2429:     }
                   2430:     w = ((long **) lo)[sg][ppsi[g]];
                   2431:     dss = deg / cgcd(deg, w - 1);
                   2432:     sr = 1;
                   2433:     for (i = 1; i < lg(B[1]) - 1; i++)
                   2434:       sr = (1 + sr * w) % deg;
                   2435:     sr = cgcd(sr, deg);
                   2436:     if (DEBUGLEVEL >= 6)
                   2437:       fprintferr("GaloisConj:w=%d [%d] sr=%d dss=%d\n", w, deg, sr, dss);
                   2438:     for (t = 0; t < sr; t += dss)
                   2439:     {
                   2440:       lbot = avma;
                   2441:       pf = testpermutation(O, B, w, t, sr, &td);
                   2442:       if (pf != gzero)
                   2443:        break;
                   2444:     }
                   2445:     if (pf == gzero)
                   2446:     {
                   2447:       freetest(&td);
                   2448:       avma = ltop;
                   2449:       return gzero;
                   2450:     } else
                   2451:     {
                   2452:       GEN f;
                   2453:       f = cgetg(3, t_VEC);
                   2454:       f[1] = (long) pf;
                   2455:       f[2] = (long) gcopy(((GEN **) PG)[j][2]);
                   2456:       res[lg(split) + j - 1] = (long) gerepile(av, lbot, f);
                   2457:     }
                   2458:   }
                   2459:   if (DEBUGLEVEL >= 4)
                   2460:     fprintferr("GaloisConj:Fini!\n");
                   2461:   freetest(&td);
                   2462:   return gerepile(ltop, lbot, res);
                   2463: }
                   2464: /*
                   2465:  * T: polynome ou nf den optionnel: si présent doit etre un multiple du
                   2466:  * dénominateur commun des solutions Si T est un nf et den n'est pas présent,
                   2467:  * le denominateur de la base d'entiers est pris pour den
                   2468:  */
                   2469: GEN
                   2470: galoisconj4(GEN T, GEN den, long flag)
                   2471: {
                   2472:   ulong ltop = avma, lbot;
                   2473:   GEN G, L, M, res, aut, L2, M2, grp;
                   2474:   struct galois_analysis ga;
                   2475:   struct galois_borne gb;
                   2476:   int n, i, j, k;
                   2477:   if (typ(T) != t_POL)
                   2478:   {
                   2479:     T = checknf(T);
                   2480:     if (!den)
                   2481:       den = gcoeff((GEN) T[8], degree((GEN) T[1]), degree((GEN) T[1]));
                   2482:     T = (GEN) T[1];
                   2483:   }
                   2484:   n = lgef(T) - 3;
                   2485:   if (n <= 0)
                   2486:     err(constpoler, "galoisconj4");
                   2487:   for (k = 2; k <= n + 2; k++)
                   2488:     if (typ(T[k]) != t_INT)
                   2489:       err(talker, "polynomial not in Z[X] in galoisconj4");
                   2490:   if (!gcmp1((GEN) T[n + 2]))
                   2491:     err(talker, "non-monic polynomial in galoisconj4");
                   2492:   n = degree(T);
                   2493:   if (n == 1)                  /* Too easy! */
                   2494:   {
                   2495:     res = cgetg(2, t_COL);
                   2496:     res[1] = (long) polx[varn(T)];
                   2497:     return res;
                   2498:   }
                   2499:   galoisanalysis(T, &ga, 1);
                   2500:   if (ga.deg == 0)
                   2501:   {
                   2502:     avma = ltop;
                   2503:     return stoi(ga.p);         /* Avoid computing the discriminant */
                   2504:   }
                   2505:   if (!den)
                   2506:     den = absi(mycoredisc(discsr(T)));
                   2507:   else
                   2508:   {
                   2509:     if (typ(den) != t_INT)
                   2510:       err(talker, "Second arg. must be integer in galoisconj4");
                   2511:     den = absi(den);
                   2512:   }
                   2513:   gb.l = stoi(ga.l);
                   2514:   initborne(T, den, &gb, ga.ppp);
                   2515:   if (DEBUGLEVEL >= 1)
                   2516:     timer2();
                   2517:   {
                   2518:     GEN f = rootpadic(T, gb.l, gb.valabs);
                   2519:     GEN _p = gmael(f, 1, 2);
                   2520:     L = gmul(gtrunc(f), gmodulcp(gun, gb.ladicabs));
                   2521:     gunclone(_p);
                   2522:   }
                   2523:   if (DEBUGLEVEL >= 1)
                   2524:     msgtimer("rootpadic()");
                   2525:   M = vandermondeinverse(L, T, den);
                   2526:   G = galoisgen(T, L, M, den, &gb, &ga);
                   2527:   if (DEBUGLEVEL >= 6)
                   2528:     fprintferr("GaloisConj:%Z\n", G);
                   2529:   if (G == gzero)
                   2530:   {
                   2531:     avma = ltop;
                   2532:     return gzero;
                   2533:   }
                   2534:   if (DEBUGLEVEL >= 1)
                   2535:     timer2();
                   2536:   L2 = gmul(L, gmodulcp(gun, gb.ladicsol));
                   2537:   M2 = gmul(M, gmodulcp(gun, gb.ladicsol));
                   2538:   if (flag)
                   2539:   {
                   2540:     lbot = avma;
                   2541:     grp = cgetg(7, t_VEC);
                   2542:     grp[1] = (long) gcopy(T);
                   2543:     grp[2] = (long) stoi(ga.l);
                   2544:     grp[3] = (long) gcopy(L);
                   2545:     grp[4] = (long) gcopy(M);
                   2546:     grp[5] = (long) gcopy(den);
                   2547:   }
                   2548:   res = cgetg(n + 1, t_VEC);
                   2549:   res[1] = (long) permidentity(n);
                   2550:   k = 1;
                   2551:   for (i = 1; i < lg(G); i++)
                   2552:   {
                   2553:     int c = k * (itos(((GEN **) G)[i][2]) - 1);
                   2554:     for (j = 1; j <= c; j++)   /* I like it */
                   2555:       res[++k] = (long) applyperm((GEN) res[j], ((GEN **) G)[i][1]);
                   2556:   }
                   2557:   if (flag)
                   2558:   {
                   2559:     grp[6] = (long) res;
                   2560:     return gerepile(ltop, lbot, grp);
                   2561:   }
                   2562:   aut = cgetg(n + 1, t_COL);
                   2563:   for (i = 1; i <= n; i++)
                   2564:     aut[i] = (long) permtopol((GEN) res[i], L2, M2, den, varn(T));
                   2565:   if (DEBUGLEVEL >= 1)
                   2566:     msgtimer("Calcul polynomes");
                   2567:   return gerepileupto(ltop, aut);
                   2568: }
                   2569: /* Calcule le nombre d'automorphisme de T avec forte probabilité */
                   2570: /* pdepart premier nombre premier a tester */
                   2571: long
                   2572: numberofconjugates(GEN T, long pdepart)
                   2573: {
                   2574:   ulong ltop = avma, ltop2;
                   2575:   long n, p, nbmax, nbtest;
                   2576:   long card;
                   2577:   byteptr primepointer;
                   2578:   int i;
                   2579:   GEN L;
                   2580:   n = degree(T);
                   2581:   card = sturm(T);
                   2582:   card = cgcd(card, n - card);
                   2583:   nbmax = n + 1;
                   2584:   nbtest = 0;
                   2585:   L = cgetg(n + 1, t_VECSMALL);
                   2586:   ltop2 = avma;
                   2587:   for (p = 0, primepointer = diffptr; nbtest < nbmax && card > 1;)
                   2588:   {
                   2589:     long s, c;
                   2590:     long isram;
                   2591:     GEN S;
                   2592:     c = *primepointer++;
                   2593:     if (!c)
                   2594:       err(primer1);
                   2595:     p += c;
                   2596:     if (p < pdepart)
                   2597:       continue;
                   2598:     S = (GEN) simplefactmod(T, stoi(p));
                   2599:     isram = 0;
                   2600:     for (i = 1; i < lg(S[2]) && !isram; i++)
                   2601:       if (!gcmp1(((GEN **) S)[2][i]))
                   2602:        isram = 1;
                   2603:     if (isram == 0)
                   2604:     {
                   2605:       for (i = 1; i <= n; i++)
                   2606:        L[i] = 0;
                   2607:       for (i = 1; i < lg(S[1]) && !isram; i++)
                   2608:        L[itos(((GEN **) S)[1][i])]++;
                   2609:       s = L[1];
                   2610:       for (i = 2; i <= n; i++)
                   2611:        s = cgcd(s, L[i] * i);
                   2612:       card = cgcd(s, card);
                   2613:     }
                   2614:     if (DEBUGLEVEL >= 6)
                   2615:       fprintferr("NumberOfConjugates:Nbtest=%d,card=%d,p=%d\n", nbtest, card, p);
                   2616:     nbtest++;
                   2617:     avma = ltop2;
                   2618:   }
                   2619:   if (DEBUGLEVEL >= 2)
                   2620:     fprintferr("NumberOfConjugates:card=%d,p=%d\n", card, p);
                   2621:   avma = ltop;
                   2622:   return card;
                   2623: }
                   2624: GEN
                   2625: galoisconj0(GEN nf, long flag, GEN d, long prec)
                   2626: {
                   2627:   GEN G, T;
                   2628:   long card;
                   2629:   if (typ(nf) != t_POL)
                   2630:   {
                   2631:     nf = checknf(nf);
                   2632:     T = (GEN) nf[1];
                   2633:   } else
                   2634:     T = nf;
                   2635:   switch (flag)
                   2636:   {
                   2637:    case 0:
                   2638:     G = galoisconj4(nf, d, 0);
                   2639:     if (typ(G) != t_INT)       /* Success */
                   2640:       return G;
                   2641:     else
                   2642:     {
                   2643:       card = G == gzero ? degree(T) : numberofconjugates(T, itos(G));
                   2644:       if (card != 1)
                   2645:       {
                   2646:        if (typ(nf) == t_POL)
                   2647:          return galoisconj2pol(nf, card, prec);
                   2648:        else
                   2649:          return galoisconj(nf);
                   2650:       }
                   2651:     }
                   2652:     break;                     /* Failure */
                   2653:    case 1:
                   2654:     return galoisconj(nf);
                   2655:    case 2:
                   2656:     return galoisconj2(nf, degree(T), prec);
                   2657:    case 4:
                   2658:     G = galoisconj4(nf, d, 0);
                   2659:     if (typ(G) != t_INT)
                   2660:       return G;
                   2661:     break;                     /* Failure */
                   2662:    default:
                   2663:     err(flagerr, "nfgaloisconj");
                   2664:   }
                   2665:   G = cgetg(2, t_COL);
                   2666:   G[1] = (long) polx[varn(T)];
                   2667:   return G;                    /* Failure */
                   2668: }
                   2669: /******************************************************************************/
                   2670: /* Galois theory related algorithms                         */
                   2671: /******************************************************************************/
                   2672: GEN
                   2673: checkgal(GEN gal)
                   2674: {
                   2675:   if (typ(gal) == t_POL)
                   2676:     err(talker, "please apply galoisinit first");
                   2677:   if (typ(gal) != t_VEC || lg(gal) != 7)
                   2678:     err(talker, "Not a Galois field in a Galois related function");
                   2679:   return gal;
                   2680: }
                   2681: GEN
                   2682: galoisinit(GEN nf, GEN den)
                   2683: {
                   2684:   GEN G;
                   2685:   G = galoisconj4(nf, den, 1);
                   2686:   if (typ(G) == t_INT)
                   2687:     err(talker, "galoisinit: field not Galois or Galois group not weakly super solvable");
                   2688:   return G;
                   2689: }
                   2690: GEN
                   2691: galoispermtopol(GEN gal, GEN perm)
                   2692: {
                   2693:   gal = checkgal(gal);
                   2694:   if (typ(perm) != t_VEC)
                   2695:     err(typeer, "galoispermtopol:");
                   2696:   return permtopol(perm, (GEN) gal[3], (GEN) gal[4], (GEN) gal[5], varn((GEN) gal[1]));
                   2697: }
                   2698: GEN
                   2699: galoisfixedfield(GEN gal, GEN perm, GEN p)
                   2700: {
                   2701:   ulong ltop = avma, lbot;
                   2702:   GEN P, S, PL, O, res, mod;
                   2703:   long x;
                   2704:   x = varn((GEN) gal[1]);
                   2705:   gal = checkgal(gal);
                   2706:   if (typ(perm) != t_VEC)
                   2707:     err(typeer, "galoisfixedfield:");
                   2708:   if (p)
                   2709:   {
                   2710:     if (typ(p) != t_INT)
                   2711:       err(talker, "galoisfixedfield: optionnal argument must be an prime integer");
                   2712:     mod = gmodulcp(gun, p);
                   2713:   } else
                   2714:     mod = gun;
                   2715:   O = vecpermorbite(perm);
                   2716:   P = corpsfixeorbitemod((GEN) gal[3], O, x, mod, (GEN) gal[2], &PL);
                   2717:   S = corpsfixeinclusion(O, PL);
                   2718:   S = vectopol(S, (GEN) gal[3], (GEN) gal[4], (GEN) gal[5], x);
                   2719:   lbot = avma;
                   2720:   res = cgetg(3, t_VEC);
                   2721:   res[1] = (long) gcopy(P);
                   2722:   res[2] = (long) gmodulcp(S, (GEN) gal[1]);
                   2723:   return gerepile(ltop, lbot, res);
                   2724: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>