Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - modules - krasner.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.2 lcov report (development 29115-f22e516b23) Lines: 440 449 98.0 %
Date: 2024-03-28 08:06:56 Functions: 34 34 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2009  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : #define DEBUGLEVEL DEBUGLEVEL_padicfields
      19             : 
      20             : /************************************************************/
      21             : /*     Computation of all the extensions of a given         */
      22             : /*               degree of a p-adic field                   */
      23             : /* Xavier Roblot                                            */
      24             : /************************************************************/
      25             : /* cf. Math. Comp, vol. 70, No. 236, pp. 1641-1659 for more details.
      26             :    Note that n is now e (since the e from the paper is always = 1) and l
      27             :    is now f */
      28             : /* Structure for a given type of extension */
      29             : typedef struct {
      30             :   GEN p;
      31             :   long e, f, j;
      32             :   long a, b, c;
      33             :   long v;     /* auxiliary variable */
      34             :   long r;     /* pr = p^r */
      35             :   GEN pr;     /* p-adic precision for poly. reduction */
      36             :   GEN q, qm1; /* p^f, q - 1 */
      37             :   GEN T;      /* polynomial defining K^ur */
      38             :   GEN frob;   /* Frobenius acting of the root of T (mod pr) */
      39             :   GEN u;      /* suitable root of unity (expressed in terms of the root of T) */
      40             :   GEN nbext;  /* number of extensions */
      41             :   GEN roottable; /* table of roots of polynomials over the residue field */
      42             :   GEN pk;     /* powers of p: [p^1, p^2, ..., p^c] */
      43             : } KRASNER_t;
      44             : 
      45             : /* Structure containing the field data (constructed with FieldData) */
      46             : typedef struct {
      47             :   GEN p;
      48             :   GEN top;  /* absolute polynomial of a = z + pi (mod pr) */
      49             :   GEN topr; /* top mod p */
      50             :   GEN z;    /* root of T in terms of a (mod pr) */
      51             :   GEN eis;  /* relative polynomial of pi in terms of z (mod pr)  */
      52             :   GEN pi;   /* prime element in terms of a */
      53             :   GEN ipi;  /* p/pi in terms of a (mod pr) (used to divide by pi) */
      54             :   GEN pik;  /* [1, pi, ..., pi^(e-1), pi^e/p] in terms of a (mod pr).
      55             :                Note the last one is different! */
      56             :   long cj;  /* number of conjugate fields */
      57             : } FAD_t;
      58             : 
      59             : #undef CHECK
      60             : 
      61             : /* Eval P(x) assuming P has positive coefficients and the result is a long */
      62             : static ulong
      63      169393 : ZX_z_eval(GEN P, ulong x)
      64             : {
      65      169393 :   long i, l = lg(P);
      66             :   ulong z;
      67             : 
      68      169393 :   if (typ(P) == t_INT) return itou(P);
      69       64463 :   if (l == 2) return 0;
      70       36995 :   z = itou(gel(P, l-1));
      71       95340 :   for (i = l-2; i > 1; i--) z = z*x + itou(gel(P,i));
      72       36995 :   return z;
      73             : }
      74             : 
      75             : /* Eval P(x, y) assuming P has positive coefficients and the result is a long */
      76             : static ulong
      77       40495 : ZXY_z_eval(GEN P, ulong x, ulong y)
      78             : {
      79       40495 :   long i, l = lg(P);
      80             :   ulong z;
      81             : 
      82       40495 :   if (l == 2) return 0;
      83       40495 :   z = ZX_z_eval(gel(P, l-1), y);
      84      169393 :   for (i = l-2; i > 1; i--) z = z*x + ZX_z_eval(gel(P, i), y);
      85       40495 :   return z;
      86             : }
      87             : 
      88             : /* P an Fq[X], where Fq = Fp[Y]/(T(Y)), a an FpX representing the automorphism
      89             :  * y -> a(y). Return a(P), applying a() coefficientwise. */
      90             : static GEN
      91       36806 : FqX_FpXQ_eval(GEN P, GEN a, GEN T, GEN p)
      92             : {
      93             :   long i, l;
      94       36806 :   GEN Q = cgetg_copy(P, &l);
      95             : 
      96       36806 :   Q[1] = P[1];
      97      166992 :   for (i = 2; i < l; i++)
      98             :   {
      99      130186 :     GEN cf = gel(P, i);
     100      130186 :     if (typ(cf) == t_POL) {
     101       87094 :       switch(degpol(cf))
     102             :       {
     103         945 :         case -1: cf = gen_0; break;
     104         511 :         case 0:  cf = gel(cf,2); break;
     105       85638 :         default:
     106       85638 :           cf = FpX_FpXQ_eval(cf, a, T, p);
     107       85638 :           cf = simplify_shallow(cf);
     108       85638 :           break;
     109             :       }
     110       43092 :     }
     111      130186 :     gel(Q, i) = cf;
     112             :   }
     113             : 
     114       36806 :   return Q;
     115             : }
     116             : 
     117             : /* Sieving routines */
     118             : static GEN
     119         497 : InitSieve(long l) { return zero_F2v(l); }
     120             : static long
     121         777 : NextZeroValue(GEN sv, long i)
     122             : {
     123        1890 :   for(; i <= sv[1]; i++)
     124        1393 :     if (!F2v_coeff(sv, i)) return i;
     125         497 :   return 0; /* sieve is full or out of the sieve! */
     126             : }
     127             : static void
     128        1113 : SetSieveValue(GEN sv, long i)
     129        1113 : { if (i >= 1 && i <= sv[1]) F2v_set(sv, i); }
     130             : 
     131             : /* return 1 if the data verify Ore condition and 0 otherwise */
     132             : static long
     133          21 : VerifyOre(GEN p, long e, long j)
     134             : {
     135             :   long nv, b, vb, nb;
     136             : 
     137          21 :   if (j < 0) return 0;
     138          21 :   nv = e * u_pval(e, p);
     139          21 :   b  = j%e;
     140          21 :   if (b == 0) return (j == nv);
     141          14 :   if (j > nv) return 0;
     142             :   /* j < nv */
     143          14 :   vb = u_pval(b, p);
     144          14 :   nb = vb*e;
     145          14 :   return (nb <= j);
     146             : }
     147             : 
     148             : /* Given [K:Q_p] = m and disc(K/Q_p) = p^d, return all decompositions K/K^ur/Q_p
     149             :  * as [e, f, j] with
     150             :  *   K^ur/Q_p unramified of degree f,
     151             :  *   K/K^ur totally ramified of degree e and discriminant p^(e+j-1);
     152             :  * thus d = f*(e+j-1) and j > 0 iff ramification is wild */
     153             : static GEN
     154           7 : possible_efj_by_d(GEN p, long m, long d)
     155             : {
     156             :   GEN rep, div;
     157             :   long i, ctr, l;
     158             : 
     159           7 :   if (d == 0) return mkvec(mkvecsmall3(1, m, 0)); /* unramified */
     160           7 :   div = divisorsu(ugcd(m, d));
     161           7 :   l = lg(div); ctr = 1;
     162           7 :   rep = cgetg(l, t_VEC);
     163          28 :   for (i = 1; i < l; i++)
     164             :   {
     165          21 :     long f = div[i], e = m/f, j = d/f - e + 1;
     166          21 :     if (VerifyOre(p, e, j)) gel(rep, ctr++) = mkvecsmall3(e, f, j);
     167             :   }
     168           7 :   setlg(rep, ctr); return rep;
     169             : }
     170             : 
     171             : /* return the number of extensions corresponding to (e,f,j) */
     172             : static GEN
     173        1743 : NumberExtensions(KRASNER_t *data)
     174             : {
     175             :   ulong pp, pa;
     176             :   long e, a, b;
     177             :   GEN p, q, s0, p1;
     178             : 
     179        1743 :   e = data->e;
     180        1743 :   p = data->p;
     181        1743 :   q = data->q;
     182        1743 :   a = data->a; /* floor(j/e) <= v_p(e), hence p^a | e */
     183        1743 :   b = data->b; /* j % e */
     184        1743 :   if (is_bigint(p)) /* implies a = 0 */
     185          70 :     return b == 0? utoipos(e): mulsi(e, data->qm1);
     186             : 
     187        1673 :   pp = p[2];
     188        1673 :   switch(a)
     189             :   {
     190        1617 :     case 0:  pa = 1;  s0 = utoipos(e); break;
     191          49 :     case 1:  pa = pp; s0 = mului(e, powiu(q, e / pa)); break;
     192           7 :     default:
     193           7 :       pa = upowuu(pp, a); /* p^a */
     194           7 :       s0 = mulsi(e, powiu(q, (e / pa) * ((pa - 1) / (pp - 1))));
     195             :   }
     196             :   /* s0 = e * q^(e*sum(p^-i, i=1...a)) */
     197        1673 :   if (b == 0) return s0;
     198             : 
     199             :   /* q^floor((b-1)/p^(a+1)) */
     200          98 :   p1 = powiu(q, sdivsi(b-1, muluu(pa, pp)));
     201          98 :   return mulii(mulii(data->qm1, s0), p1);
     202             : }
     203             : 
     204             : static GEN
     205        3213 : get_topx(KRASNER_t *data, GEN eis)
     206             : {
     207             :   GEN p1, p2, rpl, c;
     208             :   pari_sp av;
     209             :   long j;
     210             :   /* top poly. is the minimal polynomial of root(T) + root(eis) */
     211        3213 :   if (data->f == 1) return eis;
     212        1953 :   c = FpX_neg(pol_x(data->v),data->pr);
     213        1953 :   rpl = FqX_translate(eis, c, data->T, data->pr);
     214        1953 :   p1 = p2 = rpl; av = avma;
     215       21525 :   for (j = 1; j < data->f; j++)
     216             :   {
     217             :     /* compute conjugate polynomials using the Frobenius */
     218       19572 :     p1 = FqX_FpXQ_eval(p1, data->frob, data->T, data->pr);
     219       19572 :     p2 = FqX_mul(p2, p1, data->T, data->pr);
     220       19572 :     if (gc_needed(av,2)) gerepileall(av, 2, &p1,&p2);
     221             :   }
     222        1953 :   return simplify_shallow(p2); /* ZX */
     223             : }
     224             : 
     225             : /* eis (ZXY): Eisenstein polynomial over the field defined by T.
     226             :  * topx (ZX): absolute equation of root(T) + root(eis).
     227             :  * Return the struct FAD corresponding to the field it defines (GENs created
     228             :  * as clones). Assume e > 1. */
     229             : static void
     230         854 : FieldData(KRASNER_t *data, FAD_t *fdata, GEN eis, GEN topx)
     231             : {
     232             :   GEN p1, p2, p3, z, ipi, cipi, dipi, t, Q;
     233             : 
     234         854 :   fdata->p = data->p;
     235         854 :   t = leafcopy(topx); setvarn(t, data->v);
     236         854 :   fdata->top  = t;
     237         854 :   fdata->topr = FpX_red(t, data->pr);
     238             : 
     239         854 :   if (data->f == 1) z = gen_0;
     240             :   else
     241             :   { /* Compute a root of T in K(top) using Hensel's lift */
     242         238 :     z = pol_x(data->v);
     243         238 :     p1 = FpX_deriv(data->T, data->p);
     244             :     /* First lift to a root mod p */
     245             :     for (;;) {
     246         560 :       p2 = FpX_FpXQ_eval(data->T, z, fdata->top, data->p);
     247         560 :       if (gequal0(p2)) break;
     248         322 :       p3 = FpX_FpXQ_eval(p1, z, fdata->top, data->p);
     249         322 :       z  = FpX_sub(z, FpXQ_div(p2, p3, fdata->top, data->p), data->p);
     250             :     }
     251             :     /* Then a root mod p^r */
     252         238 :     z = ZpX_ZpXQ_liftroot(data->T, z, fdata->top, data->p, data->r);
     253             :   }
     254             : 
     255         854 :   fdata->z  = z;
     256         854 :   fdata->eis = eis;
     257         854 :   fdata->pi  = Fq_sub(pol_x(data->v), fdata->z,
     258             :                       FpX_red(fdata->top, data->p), data->p);
     259         854 :   ipi = RgXQ_inv(fdata->pi, fdata->top);
     260         854 :   ipi = Q_remove_denom(ipi, &dipi);
     261         854 :   Q = mulii(data->pr, data->p);
     262         854 :   cipi = Fp_inv(diviiexact(dipi, data->p), Q);
     263         854 :   fdata->ipi = FpX_Fp_mul(ipi, cipi, Q); /* p/pi mod p^(pr+1) */
     264             : 
     265             :   /* Last one is set to pi^e/p (so we compute pi^e with one extra precision) */
     266         854 :   p2 = mulii(data->pr, data->p);
     267         854 :   p1 = FpXQ_powers(fdata->pi, data->e, fdata->topr, p2);
     268         854 :   gel(p1, data->e+1) = ZX_Z_divexact(gel(p1, data->e+1), data->p);
     269         854 :   fdata->pik  = p1;
     270         854 : }
     271             : 
     272             : /* return pol*ipi/p (mod top, pp) if it has integral coefficient, NULL
     273             :    otherwise. The result is reduced mod top, pp */
     274             : static GEN
     275      307468 : DivideByPi(FAD_t *fdata, GEN pp, GEN ppp, GEN pol)
     276             : {
     277             :   GEN P;
     278             :   long i, l;
     279      307468 :   pari_sp av = avma;
     280             : 
     281             :   /* "pol" is a t_INT or an FpX: signe() works for both */
     282      307468 :   if (!signe(pol)) return pol;
     283             : 
     284      300384 :   P = Fq_mul(pol, fdata->ipi, fdata->top, ppp); /* FpX */
     285      300384 :   l  = lg(P);
     286     1562631 :   for (i = 2; i < l; i++)
     287             :   {
     288             :     GEN r;
     289     1343377 :     gel(P,i) = dvmdii(gel(P,i), fdata->p, &r);
     290     1343377 :     if (r != gen_0) return gc_NULL(av);
     291             :   }
     292      219254 :   return FpX_red(P, pp);
     293             : }
     294             : 
     295             : /* return pol# = pol~/pi^vpi(pol~) mod pp where pol~(x) = pol(pi.x + beta)
     296             :  * has coefficients in the field defined by top and pi is a prime element */
     297             : static GEN
     298       40565 : GetSharp(FAD_t *fdata, GEN pp, GEN ppp, GEN pol, GEN beta, long *pl)
     299             : {
     300             :   GEN p1, p2;
     301       40565 :   long i, v, l, d = degpol(pol);
     302       40565 :   pari_sp av = avma;
     303             : 
     304       40565 :   if (!gequal0(beta))
     305       31927 :     p1 = FqX_translate(pol, beta, fdata->topr, pp);
     306             :   else
     307        8638 :     p1 = shallowcopy(pol);
     308       40565 :   p2 = p1;
     309      132930 :   for (v = 0;; v++)
     310             :   {
     311      334103 :     for (i = 0; i <= v; i++)
     312             :     {
     313      241738 :       GEN c = gel(p2, i+2);
     314      241738 :       c = DivideByPi(fdata, pp, ppp, c);
     315      241738 :       if (!c) break;
     316      201173 :       gel(p2, i+2) = c;
     317             :     }
     318      132930 :     if (i <= v) break;
     319       92365 :     p1 = shallowcopy(p2);
     320             :   }
     321       40565 :   if (!v) pari_err_BUG("GetSharp [no division]");
     322             : 
     323       65730 :   for (l = v; l >= 0; l--)
     324             :   {
     325       65730 :     GEN c = gel(p1, l+2);
     326       65730 :     c = DivideByPi(fdata, pp, ppp, c);
     327       65730 :     if (!c) { break; }
     328             :   }
     329             : 
     330       40565 :   *pl = l; if (l < 2) return NULL;
     331             : 
     332             :   /* adjust powers */
     333       31955 :   for (i = v+1; i <= d; i++)
     334       10458 :     gel(p1, i+2) = Fq_mul(gel(p1, i+2),
     335       10458 :                           gel(fdata->pik, i-v+1), fdata->topr, pp);
     336             : 
     337       21497 :   return gerepilecopy(av, normalizepol(p1));
     338             : }
     339             : 
     340             : /* Compute roots of pol in the residue field. Use table look-up if possible */
     341             : static GEN
     342       40495 : Quick_FqX_roots(KRASNER_t *data, GEN pol)
     343             : {
     344             :   GEN rts;
     345       40495 :   ulong ind = 0;
     346             : 
     347       40495 :   if (data->f == 1)
     348       22414 :     pol = FpXY_evalx(pol, gen_0, data->p);
     349             :   else
     350       18081 :     pol = FqX_red(pol, data->T, data->p);
     351       40495 :   if (data->roottable)
     352             :   {
     353       40495 :     ind = ZXY_z_eval(pol, data->q[2], data->p[2]);
     354       40495 :     if (gel(data->roottable, ind)) return gel(data->roottable, ind);
     355             :   }
     356        2856 :   rts = FqX_roots(pol, data->T, data->p);
     357        2856 :   if (ind) gel(data->roottable, ind) = gclone(rts);
     358        2856 :   return rts;
     359             : }
     360             : 
     361             : static void
     362         133 : FreeRootTable(GEN T)
     363             : {
     364         133 :   if (T)
     365             :   {
     366         133 :     long j, l = lg(T);
     367      590373 :     for (j = 1; j < l; j++) guncloneNULL(gel(T,j));
     368             :   }
     369         133 : }
     370             : 
     371             : /* Return the number of roots of pol congruent to alpha modulo pi working
     372             :    modulo pp (all roots if alpha is NULL); if flag is nonzero, return 1
     373             :    as soon as a root is found. (Note: ppp = pp*p for DivideByPi) */
     374             : static long
     375       59563 : RootCongruents(KRASNER_t *data, FAD_t *fdata, GEN pol, GEN alpha, GEN pp, GEN ppp, long sl, long flag)
     376             : {
     377             :   GEN R;
     378             :   long s, i;
     379             : 
     380       59563 :   if (alpha)
     381             :   {
     382             :     long l;
     383       40565 :     pol = GetSharp(fdata, pp, ppp, pol, alpha, &l);
     384       40565 :     if (l <= 1) return l;
     385             :     /* decrease precision if sl gets bigger than a multiple of e */
     386       21497 :     sl += l;
     387       21497 :     if (sl >= data->e)
     388             :     {
     389       18249 :       sl -= data->e;
     390       18249 :       ppp = pp;
     391       18249 :       pp = diviiexact(pp, data->p);
     392             :     }
     393             :   }
     394             : 
     395       40495 :   R  = Quick_FqX_roots(data, pol);
     396       75390 :   for (s = 0, i = 1; i < lg(R); i++)
     397             :   {
     398       40565 :     s += RootCongruents(data, fdata, pol, gel(R, i), pp, ppp, sl, flag);
     399       40565 :     if (flag && s) return 1;
     400             :   }
     401       34825 :   return s;
     402             : }
     403             : 
     404             : /* pol is a ZXY defining a polynomial over the field defined by fdata
     405             :    If flag != 0, return 1 as soon as a root is found. Computations are done with
     406             :    a precision of pr. */
     407             : static long
     408       18998 : RootCountingAlgorithm(KRASNER_t *data, FAD_t *fdata, GEN pol, long flag)
     409             : {
     410             :   long j, l, d;
     411       18998 :   GEN P = cgetg_copy(pol, &l);
     412             : 
     413       18998 :   P[1] = pol[1];
     414       18998 :   d = l-3;
     415       84693 :   for (j = 0; j < d; j++)
     416             :   {
     417       65695 :     GEN cf = gel(pol, j+2);
     418       65695 :     if (typ(cf) == t_INT)
     419       38430 :       cf = diviiexact(cf, data->p);
     420             :     else
     421       27265 :       cf = ZX_Z_divexact(cf, data->p);
     422       65695 :     gel(P, j+2) = Fq_mul(cf, gel(fdata->pik, j+1), fdata->topr, data->pr);
     423             :   }
     424       18998 :   gel(P, d+2) = gel(fdata->pik, d+1); /* pik[d] = pi^d/p */
     425             : 
     426       18998 :   return RootCongruents(data, fdata, P, NULL, diviiexact(data->pr, data->p), data->pr, 0, flag);
     427             : }
     428             : 
     429             : /* Return nonzero if the field given by fdata defines a field isomorphic to
     430             :  * the one defined by pol */
     431             : static long
     432       11662 : IsIsomorphic(KRASNER_t *data, FAD_t *fdata, GEN pol)
     433             : {
     434             :   long j, nb;
     435       11662 :   pari_sp av = avma;
     436             : 
     437       11662 :   if (RgX_is_ZX(pol)) return RootCountingAlgorithm(data, fdata, pol, 1);
     438             : 
     439       13545 :   for (j = 1; j <= data->f; j++)
     440             :   {
     441       10108 :     GEN p1 = FqX_FpXQ_eval(pol, fdata->z, fdata->top, data->pr);
     442       10108 :     nb = RootCountingAlgorithm(data, fdata, p1, 1);
     443       10108 :     if (nb) return gc_long(av, nb);
     444        9513 :     if (j < data->f)
     445        6076 :       pol = FqX_FpXQ_eval(pol, data->frob, data->T, data->pr);
     446             :   }
     447        3437 :   return gc_long(av,0);
     448             : }
     449             : 
     450             : /* Compute the number of conjugates fields of the field given by fdata */
     451             : static void
     452         854 : NbConjugateFields(KRASNER_t *data, FAD_t *fdata)
     453             : {
     454         854 :   GEN pol = fdata->eis;
     455             :   long j, nb;
     456         854 :   pari_sp av = avma;
     457             : 
     458         854 :   if (RgX_is_ZX(pol)) { /* split for efficiency; contains the case f = 1 */
     459         616 :     fdata->cj = data->e / RootCountingAlgorithm(data, fdata, pol, 0);
     460         616 :     set_avma(av); return;
     461             :   }
     462             : 
     463         238 :   nb = 0;
     464         882 :   for (j = 1; j <= data->f; j++)
     465             :   { /* Transform to pol. in z to pol. in a */
     466         644 :     GEN p1 = FqX_FpXQ_eval(pol, fdata->z, fdata->top, data->pr);
     467         644 :     nb += RootCountingAlgorithm(data, fdata, p1, 0);
     468             :     /* Look at the roots of conjugates polynomials */
     469         644 :     if (j < data->f)
     470         406 :       pol = FqX_FpXQ_eval(pol, data->frob, data->T, data->pr);
     471             :   }
     472         238 :   set_avma(av);
     473         238 :   fdata->cj = data->e * data->f / nb;
     474         238 :   return;
     475             : }
     476             : 
     477             : /* return a minimal list of polynomials generating all the totally
     478             :    ramified extensions of K^ur of degree e and discriminant
     479             :    p^{e + j - 1} in the tamely ramified case */
     480             : static GEN
     481        1582 : TamelyRamifiedCase(KRASNER_t *data)
     482             : {
     483        1582 :   long av = avma;
     484        1582 :   long g = ugcdui(data->e, data->qm1); /* (e, q-1) */
     485        1582 :   GEN rep, eis, Xe = gpowgs(pol_x(0), data->e), m = stoi(data->e / g);
     486             : 
     487        1582 :   rep = zerovec(g);
     488        1582 :   eis = gadd(Xe, data->p);
     489        1582 :   gel(rep, 1) = mkvec2(get_topx(data,eis), m);
     490        1582 :   if (g > 1)
     491             :   {
     492         497 :     ulong pmodg = umodiu(data->p, g);
     493         497 :     long r = 1, ct = 1;
     494         497 :     GEN sv = InitSieve(g-1);
     495             :     /* let Frobenius act to get a minimal set of polynomials over Q_p */
     496        1274 :     while (r)
     497             :     {
     498             :       long gr;
     499        1554 :       GEN p1 = (typ(data->u) == t_INT)?
     500         777 :         mulii(Fp_powu(data->u, r, data->p), data->p):
     501         588 :         ZX_Z_mul(FpXQ_powu(data->u, r, data->T, data->p), data->p);
     502         777 :       eis = gadd(Xe, p1); /* add cst coef */
     503         777 :       ct++;
     504         777 :       gel(rep, ct) = mkvec2(get_topx(data,eis), m);
     505         777 :       gr = r;
     506        1113 :       do { SetSieveValue(sv, gr); gr = Fl_mul(gr, pmodg, g); } while (gr != r);
     507         777 :       r  = NextZeroValue(sv, r);
     508             :     }
     509         497 :     setlg(rep, ct+1);
     510             :   }
     511        1582 :   return gerepilecopy(av, rep);
     512             : }
     513             : 
     514             : static long
     515         399 : function_l(GEN p, long a, long b, long i)
     516             : {
     517         399 :   long l = 1 + a - z_pval(i, p);
     518         399 :   if (i < b) l++;
     519         399 :   return (l < 1)? 1: l;
     520             : }
     521             : 
     522             : /* Structure of the coefficients set Omega. Each coefficient is [start, zr]
     523             :  * meaning all the numbers of the form:
     524             :  *   zeta_0 * p^start + ... + zeta_s * p^c (s = c - start)
     525             :  * with zeta_i roots of unity (powers of zq + zero), zeta_0 = 0 is
     526             :  * possible iff zr = 1 */
     527             : static GEN
     528         133 : StructureOmega(KRASNER_t *data, GEN *pnbpol)
     529             : {
     530         133 :   GEN nbpol, O = cgetg(data->e + 1, t_VEC);
     531             :   long i;
     532             : 
     533             :   /* i = 0 */
     534         133 :   gel(O, 1) = mkvecsmall2(1, 0);
     535         133 :   nbpol = mulii(data->qm1, powiu(data->q, data->c - 1));
     536         532 :   for (i = 1; i < data->e; i++)
     537             :   {
     538         399 :     long v_start, zero = 0;
     539             :     GEN nbcf, p1;
     540         399 :     v_start = function_l(data->p, data->a, data->b, i);
     541         399 :     p1 = powiu(data->q, data->c - v_start);
     542         399 :     if (i == data->b)
     543          98 :       nbcf = mulii(p1, data->qm1);
     544             :     else
     545             :     {
     546         301 :       nbcf = mulii(p1, data->q);
     547         301 :       zero = 1;
     548             :     }
     549         399 :     gel(O, i+1) = mkvecsmall2(v_start, zero);
     550         399 :     nbpol = mulii(nbpol, nbcf);
     551             :   }
     552         133 :   *pnbpol = nbpol; return O;
     553             : }
     554             : 
     555             : /* a random element of the finite field */
     556             : static GEN
     557       37450 : RandomFF(KRASNER_t *data)
     558             : {
     559       37450 :   long i, l = data->f + 2, p = itou(data->p);
     560       37450 :   GEN c = cgetg(l, t_POL);
     561       37450 :   c[1] = evalvarn(data->v);
     562       86919 :   for (i = 2; i < l; i++) gel(c, i) = utoi(random_Fl(p));
     563       37450 :   return ZX_renormalize(c, l);
     564             : }
     565             : static GEN
     566        2709 : RandomPol(KRASNER_t *data, GEN Omega)
     567             : {
     568        2709 :   long i, j, l = data->e + 3, end = data->c;
     569        2709 :   GEN pol = cgetg(l, t_POL);
     570        2709 :   pol[1] = evalsigne(1) | evalvarn(0);
     571       13230 :   for (i = 1; i <= data->e; i++)
     572             :   {
     573       10521 :     GEN c, cf = gel(Omega, i);
     574       10521 :     long st = cf[1], zr = cf[2];
     575             :     /* c = sum_{st <= j <= end} x_j p^j, where x_j are random Fq mod (p,upl)
     576             :      * If (!zr), insist on x_st != 0 */
     577             :     for (;;) {
     578       13027 :       c = RandomFF(data);
     579       13027 :       if (zr || signe(c)) break;
     580             :     }
     581       34944 :     for (j = 1; j <= end-st; j++)
     582       24423 :       c = ZX_add(c, ZX_Z_mul(RandomFF(data), gel(data->pk, j)));
     583       10521 :     c = ZX_Z_mul(c, gel(data->pk, st));
     584       10521 :     c = FpX_red(c, data->pr);
     585       10521 :     switch(degpol(c))
     586             :     {
     587         672 :       case -1: c = gen_0; break;
     588        7770 :       case  0: c = gel(c,2); break;
     589             :     }
     590       10521 :     gel(pol, i+1) = c;
     591             :   }
     592        2709 :   gel(pol, i+1) = gen_1; /* monic */
     593        2709 :   return pol;
     594             : }
     595             : 
     596             : static void
     597         854 : CloneFieldData(FAD_t *fdata)
     598             : {
     599         854 :  fdata->top = gclone(fdata->top);
     600         854 :  fdata->topr= gclone(fdata->topr);
     601         854 :  fdata->z   = gclone(fdata->z);
     602         854 :  fdata->eis = gclone(fdata->eis);
     603         854 :  fdata->pi  = gclone(fdata->pi);
     604         854 :  fdata->ipi = gclone(fdata->ipi);
     605         854 :  fdata->pik = gclone(fdata->pik);
     606         854 : }
     607             : static void
     608         854 : FreeFieldData(FAD_t *fdata)
     609             : {
     610         854 :   gunclone(fdata->top);
     611         854 :   gunclone(fdata->topr);
     612         854 :   gunclone(fdata->z);
     613         854 :   gunclone(fdata->eis);
     614         854 :   gunclone(fdata->pi);
     615         854 :   gunclone(fdata->ipi);
     616         854 :   gunclone(fdata->pik);
     617         854 : }
     618             : 
     619             : static GEN
     620         133 : WildlyRamifiedCase(KRASNER_t *data)
     621             : {
     622         133 :   long nbext, ct, fd, nb = 0, j;
     623             :   GEN nbpol, rpl, rep, Omega;
     624             :   FAD_t **vfd;
     625             :   pari_timer T;
     626         133 :   pari_sp av = avma, av2;
     627             : 
     628             :   /* Omega = vector giving the structure of the set Omega */
     629             :   /* nbpol = number of polynomials to consider = |Omega| */
     630         133 :   Omega = StructureOmega(data, &nbpol);
     631         133 :   nbext = itos_or_0(data->nbext);
     632         133 :   if (!nbext || (nbext & ~LGBITS))
     633           0 :     pari_err_OVERFLOW("padicfields [too many extensions]");
     634             : 
     635         133 :   if (DEBUGLEVEL>0) {
     636           0 :     err_printf("There are %ld extensions to find and %Ps polynomials to consider\n", nbext, nbpol);
     637           0 :     timer_start(&T);
     638             :   }
     639             : 
     640         133 :   vfd = (FAD_t**)new_chunk(nbext);
     641        2541 :   for (j = 0; j < nbext; j++) vfd[j] = (FAD_t*)new_chunk(sizeof(FAD_t));
     642             : 
     643         133 :   ct = 0;
     644         133 :   fd = 0;
     645         133 :   av2 = avma;
     646             : 
     647        2842 :   while (fd < nbext)
     648             :   { /* Jump randomly among the polynomials : seems best... */
     649        2709 :     rpl = RandomPol(data, Omega);
     650        2709 :     if (DEBUGLEVEL>3) err_printf("considering polynomial %Ps\n", rpl);
     651       12516 :     for (j = 0; j < ct; j++)
     652             :     {
     653       11662 :       nb = IsIsomorphic(data, vfd[j], rpl);
     654       11662 :       if (nb) break;
     655             :     }
     656        2709 :     if (!nb)
     657             :     {
     658         854 :       GEN topx = get_topx(data, rpl);
     659         854 :       FAD_t *f = (FAD_t*)vfd[ct];
     660         854 :       FieldData(data, f, rpl, topx);
     661         854 :       CloneFieldData(f);
     662         854 :       NbConjugateFields(data, f);
     663         854 :       nb = f->cj;
     664         854 :       fd += nb;
     665         854 :       ct++;
     666         854 :       if (DEBUGLEVEL>1)
     667           0 :         err_printf("%ld more extension%s\t(%ld/%ld, %ldms)\n",
     668             :                    nb, (nb == 1)? "": "s", fd, nbext, timer_delay(&T));
     669             :     }
     670        2709 :     set_avma(av2);
     671             :   }
     672             : 
     673         133 :   rep = cgetg(ct+1, t_VEC);
     674         987 :   for (j = 0; j < ct; j++)
     675             :   {
     676         854 :     FAD_t *f = (FAD_t*)vfd[j];
     677         854 :     GEN topx = ZX_copy(f->top);
     678         854 :     setvarn(topx, 0);
     679         854 :     gel(rep, j+1) = mkvec2(topx, stoi(f->cj));
     680         854 :     FreeFieldData(f);
     681             :   }
     682         133 :   FreeRootTable(data->roottable);
     683         133 :   return gerepileupto(av, rep);
     684             : }
     685             : 
     686             : /* return the minimal polynomial T of a generator of K^ur and the expression (mod pr)
     687             :  * in terms of T of a root of unity u such that u is l-maximal for all primes l
     688             :  * dividing g = (e,q-1). */
     689             : static void
     690        1715 : setUnramData(KRASNER_t *d)
     691             : {
     692        1715 :   if (d->f == 1)
     693             :   {
     694         560 :     d->T = pol_x(d->v);
     695         560 :     d->u = pgener_Fp(d->p);
     696         560 :     d->frob = pol_x(d->v);
     697             :   }
     698             :   else
     699             :   {
     700        1155 :     GEN L, z, T, p = d->p;
     701        1155 :     d->T = T = init_Fq(p, d->f, d->v);
     702        1155 :     L = gel(factoru(ugcdui(d->e, d->qm1)), 1);
     703        1155 :     z = gener_FpXQ_local(T, p, zv_to_ZV(L));
     704        1155 :     d->u = ZpXQ_sqrtnlift(pol_1(d->v), d->qm1, z, T, p, d->r);
     705        1155 :     d->frob = ZpX_ZpXQ_liftroot(T, FpXQ_pow(pol_x(d->v), p, T, p), T, p, d->r);
     706             :   }
     707        1715 : }
     708             : 
     709             : /* return [ p^1, p^2, ..., p^c ] */
     710             : static GEN
     711         133 : get_pk(KRASNER_t *data)
     712             : {
     713         133 :   long i, l = data->c + 1;
     714         133 :   GEN pk = cgetg(l, t_VEC), p = data->p;
     715         133 :   gel(pk, 1) = p;
     716         455 :   for (i = 2; i <= data->c; i++) gel(pk, i) = mulii(gel(pk, i-1), p);
     717         133 :   return pk;
     718             : }
     719             : 
     720             : /* Compute an absolute polynomial for all the totally ramified
     721             :    extensions of Q_p(z) of degree e and discriminant p^{e + j - 1}
     722             :    where z is a root of upl defining an unramified extension of Q_p */
     723             : /* See padicfields for the meaning of flag */
     724             : static GEN
     725        1743 : GetRamifiedPol(GEN p, GEN efj, long flag)
     726             : {
     727        1743 :   long e = efj[1], f = efj[2], j = efj[3], i, l;
     728        1743 :   const long v = 1;
     729             :   GEN L, pols;
     730             :   KRASNER_t data;
     731        1743 :   pari_sp av = avma;
     732             : 
     733        1743 :   if (DEBUGLEVEL>1)
     734           0 :     err_printf("  Computing extensions with decomposition [e, f, j] = [%ld, %ld, %ld]\n", e,f,j);
     735        1743 :   data.p   = p;
     736        1743 :   data.e   = e;
     737        1743 :   data.f   = f;
     738        1743 :   data.j   = j;
     739        1743 :   data.a   = j/e;
     740        1743 :   data.b   = j%e;
     741        1743 :   data.c   = (e+2*j)/e+1;
     742        1743 :   data.q   = powiu(p, f);
     743        1743 :   data.qm1 = subiu(data.q, 1);
     744        1743 :   data.v   = v;
     745        1743 :   data.r   = 1 + (long)ceildivuu(2*j + 3, e); /* enough precision */
     746        1743 :   data.pr  = powiu(p, data.r);
     747        1743 :   data.nbext = NumberExtensions(&data);
     748             : 
     749        1743 :   if (flag == 2) return data.nbext;
     750             : 
     751        1715 :   setUnramData(&data);
     752        1715 :   if (DEBUGLEVEL>1)
     753           0 :     err_printf("  Unramified part %Ps\n", data.T);
     754        1715 :   data.roottable = NULL;
     755        1715 :   if (j)
     756             :   {
     757         133 :     if (lgefint(data.q) == 3)
     758             :     {
     759         133 :       ulong npol = upowuu(data.q[2], e+1);
     760         133 :       if (npol && npol < (1<<19)) data.roottable = const_vec(npol, NULL);
     761             :     }
     762         133 :     data.pk = get_pk(&data);
     763         133 :     L = WildlyRamifiedCase(&data);
     764             :   }
     765             :   else
     766        1582 :     L = TamelyRamifiedCase(&data);
     767             : 
     768        1715 :   pols = cgetg_copy(L, &l);
     769        1715 :   if (flag == 1)
     770             :   {
     771        1631 :     GEN E = utoipos(e), F = utoipos(f), D = utoi(f*(e+j-1));
     772        4578 :     for (i = 1; i < l; i++)
     773             :     {
     774        2947 :       GEN T = gel(L,i);
     775        2947 :       gel(pols, i) = mkvec5(ZX_copy(gel(T, 1)), E,F,D, icopy(gel(T, 2)));
     776             :     }
     777             :   }
     778             :   else
     779             :   {
     780         350 :     for (i = 1; i < l; i++)
     781             :     {
     782         266 :       GEN T = gel(L,i);
     783         266 :       gel(pols, i) = ZX_copy(gel(T,1));
     784             :     }
     785             :   }
     786        1715 :   return gerepileupto(av, pols);
     787             : }
     788             : 
     789             : static GEN
     790         483 : possible_efj(GEN p, long m)
     791             : { /* maximal possible discriminant valuation d <= m * (1+v_p(m)) - 1 */
     792             :   /* 1) [j = 0, tame] d = m - f with f | m and v_p(f) = v_p(m), or
     793             :    * 2) [j > 0, wild] d >= m, j <= v_p(e)*e   */
     794         483 :   ulong m1, pve, pp = p[2]; /* pp used only if v > 0 */
     795         483 :   long ve, v = u_pvalrem(m, p, &m1);
     796         483 :   GEN L, D = divisorsu(m1);
     797         483 :   long i, taum1 = lg(D)-1, nb = 0;
     798             : 
     799         483 :   if (v) {
     800          49 :     for (pve = 1,ve = 1; ve <= v; ve++) { pve *= pp; nb += pve * ve; }
     801          21 :     nb = itos_or_0(muluu(nb, zv_sum(D)));
     802          21 :     if (!nb || is_bigint( mului(pve, sqru(v)) ) )
     803           0 :       pari_err_OVERFLOW("padicfields [too many ramification possibilities]");
     804             :   }
     805         483 :   nb += taum1; /* upper bound for the number of possible triples [e,f,j] */
     806             : 
     807         483 :   L = cgetg(nb + 1, t_VEC);
     808             :   /* 1) tame */
     809        2093 :   for (nb=1, i=1; i < lg(D); i++)
     810             :   {
     811        1610 :     long e = D[i], f = m / e;
     812        1610 :     gel(L, nb++) = mkvecsmall3(e, f, 0);
     813             :   }
     814             :   /* 2) wild */
     815             :   /* Ore's condition: either
     816             :    * 1) j = v_p(e) * e, or
     817             :    * 2) j = a e + b, with 0 < b < e and v_p(b) <= a < v_p(e) */
     818         511 :   for (pve = 1, ve = 1; ve <= v; ve++)
     819             :   {
     820          28 :     pve *= pp; /* = p^ve */
     821          56 :     for (i = 1; i < lg(D); i++)
     822             :     {
     823          28 :       long a,b, e = D[i] * pve, f = m / e;
     824          98 :       for (b = 1; b < e; b++)
     825         154 :         for (a = u_lval(b, pp); a < ve; a++)
     826          84 :           gel(L, nb++) = mkvecsmall3(e, f,  a*e + b);
     827          28 :       gel(L, nb++) = mkvecsmall3(e, f, ve*e);
     828             :     }
     829             :   }
     830         483 :   setlg(L, nb); return L;
     831             : }
     832             : 
     833             : #ifdef CHECK
     834             : static void
     835             : checkpols(GEN p, GEN EFJ, GEN pols)
     836             : {
     837             :   GEN pol, p1, p2;
     838             :   long c1, c2, e, f, j, i, l = lg(pols);
     839             : 
     840             :   if (typ(pols) == t_INT) return;
     841             : 
     842             :   e = EFJ[1];
     843             :   f = EFJ[2];
     844             :   j = EFJ[3];
     845             : 
     846             :   for (i = 1; i < l; i++)
     847             :   {
     848             :     pol = gel(pols, i);
     849             :     if (typ(pol) == t_VEC) pol = gel(pol, 1);
     850             :     if (!polisirreducible(pol)) pari_err_BUG("Polynomial is reducible");
     851             :     p1 = poldisc0(pol, -1);
     852             :     if (gvaluation(p1, p) != f*(e+j-1)) pari_err_BUG("Discriminant is incorrect");
     853             :     /* only compute a p-maximal order */
     854             :     p1 = nfinit0(mkvec2(pol, mkvec(p)), 0, DEFAULTPREC);
     855             :     p2 = idealprimedec(p1, p);
     856             :     if(lg(p2) > 2) pari_err_BUG("Prime p is split");
     857             :     p2 = gel(p2, 1);
     858             :     if (cmpis(gel(p2, 3), e)) pari_err_BUG("Ramification index is incorrect");
     859             :     if (cmpis(gel(p2, 4), f)) pari_err_BUG("inertia degree is incorrect");
     860             :   }
     861             : 
     862             :   if (l == 2) return;
     863             :   if (e*f > 20) return;
     864             : 
     865             :   /* TODO: check that (random) distinct polynomials give nonisomorphic extensions */
     866             :   for (i = 1; i < 3*l; i++)
     867             :   {
     868             :     c1 = random_Fl(l-1)+1;
     869             :     c2 = random_Fl(l-1)+1;
     870             :     if (c1 == c2) continue;
     871             :     p1 = gel(pols, c1);
     872             :     if (typ(p1) == t_VEC) p1 = gel(p1, 1);
     873             :     p2 = gel(pols, c2);
     874             :     if (typ(p2) == t_VEC) p2 = gel(p2, 1);
     875             :     pol = polcompositum0(p1, p2, 0);
     876             :     pol = gel(pol, 1);
     877             :     if (poldegree(pol, -1) > 100) continue;
     878             :     p1 = factorpadic(pol, p, 2);
     879             :     p1 = gmael(p1, 1, 1);
     880             :     if (poldegree(p1, -1) == e*f) pari_err_BUG("fields are isomorphic");
     881             :     /*
     882             :       p1 = nfinit0(mkvec2(pol, mkvec(p)), 0, DEFAULTPREC);
     883             :       p2 = idealprimedec_galois(p1, p);
     884             :       if (!cmpis(mulii(gel(p2, 3), gel(p2, 4)), e*f)) pari_err_BUG("fields are isomorphic");
     885             :     */
     886             :   }
     887             : }
     888             : #endif
     889             : 
     890             : static GEN
     891         490 : pols_from_efj(pari_sp av, GEN EFJ, GEN p, long flag)
     892             : {
     893             :   long i, l;
     894         490 :   GEN L = cgetg_copy(EFJ, &l);
     895         490 :   if (l == 1) { set_avma(av); return flag == 2? gen_0: cgetg(1, t_VEC); }
     896        2233 :   for (i = 1; i < l; i++)
     897             :   {
     898        1743 :     gel(L,i) = GetRamifiedPol(p, gel(EFJ,i), flag);
     899             : #ifdef CHECK
     900             :     checkpols(p, gel(EFJ, i), gel(L, i));
     901             : #endif
     902             :   }
     903         490 :   if (flag == 2) return gerepileuptoint(av, ZV_sum(L));
     904         483 :   return gerepilecopy(av, shallowconcat1(L));
     905             : }
     906             : 
     907             : /* return a minimal list of polynomials generating all the extensions of
     908             :    Q_p of given degree N; if N is a t_VEC [n,d], return extensions of degree n
     909             :    and discriminant p^d. */
     910             : /* Return only the polynomials if flag = 0 (default), also the ramification
     911             :    degree, the residual degree and the discriminant if flag = 1 and only the
     912             :    number of extensions if flag = 2 */
     913             : GEN
     914         490 : padicfields0(GEN p, GEN N, long flag)
     915             : {
     916         490 :   pari_sp av = avma;
     917         490 :   long m = 0, d = -1;
     918             :   GEN L;
     919             : 
     920         490 :   if (typ(p) != t_INT) pari_err_TYPE("padicfields",p);
     921             :   /* be nice to silly users */
     922         490 :   if (!BPSW_psp(p)) pari_err_PRIME("padicfields",p);
     923         490 :   switch(typ(N))
     924             :   {
     925           7 :     case t_VEC:
     926           7 :       if (lg(N) != 3) pari_err_TYPE("padicfields",N);
     927           7 :       d = gtos(gel(N,2));
     928           7 :       N = gel(N,1); /* fall through */
     929         490 :     case t_INT:
     930         490 :       m = itos(N);
     931         490 :       if (m <= 0) pari_err_DOMAIN("padicfields", "degree", "<=", gen_0,N);
     932         490 :       break;
     933           0 :     default:
     934           0 :       pari_err_TYPE("padicfields",N);
     935             :   }
     936         490 :   if (d >= 0) return padicfields(p, m, d, flag);
     937         483 :   L = possible_efj(p, m);
     938         483 :   return pols_from_efj(av, L, p, flag);
     939             : }
     940             : 
     941             : GEN
     942           7 : padicfields(GEN p, long m, long d, long flag)
     943             : {
     944           7 :   long av = avma;
     945           7 :   GEN L = possible_efj_by_d(p, m, d);
     946           7 :   return pols_from_efj(av, L, p, flag);
     947             : }

Generated by: LCOV version 1.14