"Fossies" - the Fresh Open Source Software Archive

Member "pari-2.13.1/src/basemath/base4.c" (1 Nov 2020, 82392 Bytes) of package /linux/misc/pari-2.13.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "base4.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.13.0_vs_2.13.1.

    1 /* Copyright (C) 2000  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 /*******************************************************************/
   16 /*                                                                 */
   17 /*                       BASIC NF OPERATIONS                       */
   18 /*                           (continued)                           */
   19 /*                                                                 */
   20 /*******************************************************************/
   21 #include "pari.h"
   22 #include "paripriv.h"
   23 
   24 /*******************************************************************/
   25 /*                                                                 */
   26 /*                     IDEAL OPERATIONS                            */
   27 /*                                                                 */
   28 /*******************************************************************/
   29 
   30 /* A valid ideal is either principal (valid nf_element), or prime, or a matrix
   31  * on the integer basis in HNF.
   32  * A prime ideal is of the form [p,a,e,f,b], where the ideal is p.Z_K+a.Z_K,
   33  * p is a rational prime, a belongs to Z_K, e=e(P/p), f=f(P/p), and b
   34  * is Lenstra's constant, such that p.P^(-1)= p Z_K + b Z_K.
   35  *
   36  * An extended ideal is a couple [I,F] where I is an ideal and F is either an
   37  * algebraic number, or a factorization matrix attached to an algebraic number.
   38  * All routines work with either extended ideals or ideals (an omitted F is
   39  * assumed to be factor(1)). All ideals are output in HNF form. */
   40 
   41 /* types and conversions */
   42 
   43 long
   44 idealtyp(GEN *ideal, GEN *arch)
   45 {
   46   GEN x = *ideal;
   47   long t,lx,tx = typ(x);
   48 
   49   if (tx!=t_VEC || lg(x)!=3) *arch = NULL;
   50   else
   51   {
   52     GEN a = gel(x,2);
   53     if (typ(a) == t_MAT && lg(a) != 3)
   54     { /* allow [;] */
   55       if (lg(a) != 1) pari_err_TYPE("idealtyp [extended ideal]",x);
   56       a = trivial_fact();
   57     }
   58     *arch = a;
   59     x = gel(x,1); tx = typ(x);
   60   }
   61   switch(tx)
   62   {
   63     case t_MAT: lx = lg(x);
   64       if (lx == 1) { t = id_PRINCIPAL; x = gen_0; break; }
   65       if (lx != lgcols(x)) pari_err_TYPE("idealtyp [nonsquare t_MAT]",x);
   66       t = id_MAT;
   67       break;
   68 
   69     case t_VEC: if (lg(x)!=6) pari_err_TYPE("idealtyp",x);
   70       t = id_PRIME; break;
   71 
   72     case t_POL: case t_POLMOD: case t_COL:
   73     case t_INT: case t_FRAC:
   74       t = id_PRINCIPAL; break;
   75     default:
   76       pari_err_TYPE("idealtyp",x);
   77       return 0; /*LCOV_EXCL_LINE*/
   78   }
   79   *ideal = x; return t;
   80 }
   81 
   82 /* true nf; v = [a,x,...], a in Z. Return (a,x) */
   83 GEN
   84 idealhnf_two(GEN nf, GEN v)
   85 {
   86   GEN p = gel(v,1), pi = gel(v,2), m = zk_scalar_or_multable(nf, pi);
   87   if (typ(m) == t_INT) return scalarmat(gcdii(m,p), nf_get_degree(nf));
   88   return ZM_hnfmodid(m, p);
   89 }
   90 /* true nf */
   91 GEN
   92 pr_hnf(GEN nf, GEN pr)
   93 {
   94   GEN p = pr_get_p(pr), m;
   95   if (pr_is_inert(pr)) return scalarmat(p, nf_get_degree(nf));
   96   m = zk_scalar_or_multable(nf, pr_get_gen(pr));
   97   return ZM_hnfmodprime(m, p);
   98 }
   99 
  100 GEN
  101 idealhnf_principal(GEN nf, GEN x)
  102 {
  103   GEN cx;
  104   x = nf_to_scalar_or_basis(nf, x);
  105   switch(typ(x))
  106   {
  107     case t_COL: break;
  108     case t_INT:  if (!signe(x)) return cgetg(1,t_MAT);
  109       return scalarmat(absi_shallow(x), nf_get_degree(nf));
  110     case t_FRAC:
  111       return scalarmat(Q_abs_shallow(x), nf_get_degree(nf));
  112     default: pari_err_TYPE("idealhnf",x);
  113   }
  114   x = Q_primitive_part(x, &cx);
  115   RgV_check_ZV(x, "idealhnf");
  116   x = zk_multable(nf, x);
  117   x = ZM_hnfmodid(x, zkmultable_capZ(x));
  118   return cx? ZM_Q_mul(x,cx): x;
  119 }
  120 
  121 /* x integral ideal in t_MAT form, nx columns */
  122 static GEN
  123 vec_mulid(GEN nf, GEN x, long nx, long N)
  124 {
  125   GEN m = cgetg(nx*N + 1, t_MAT);
  126   long i, j, k;
  127   for (i=k=1; i<=nx; i++)
  128     for (j=1; j<=N; j++) gel(m, k++) = zk_ei_mul(nf, gel(x,i),j);
  129   return m;
  130 }
  131 /* true nf */
  132 GEN
  133 idealhnf_shallow(GEN nf, GEN x)
  134 {
  135   long tx = typ(x), lx = lg(x), N;
  136 
  137   /* cannot use idealtyp because here we allow nonsquare matrices */
  138   if (tx == t_VEC && lx == 3) { x = gel(x,1); tx = typ(x); lx = lg(x); }
  139   if (tx == t_VEC && lx == 6) return pr_hnf(nf,x); /* PRIME */
  140   switch(tx)
  141   {
  142     case t_MAT:
  143     {
  144       GEN cx;
  145       long nx = lx-1;
  146       N = nf_get_degree(nf);
  147       if (nx == 0) return cgetg(1, t_MAT);
  148       if (nbrows(x) != N) pari_err_TYPE("idealhnf [wrong dimension]",x);
  149       if (nx == 1) return idealhnf_principal(nf, gel(x,1));
  150 
  151       if (nx == N && RgM_is_ZM(x) && ZM_ishnf(x)) return x;
  152       x = Q_primitive_part(x, &cx);
  153       if (nx < N) x = vec_mulid(nf, x, nx, N);
  154       x = ZM_hnfmod(x, ZM_detmult(x));
  155       return cx? ZM_Q_mul(x,cx): x;
  156     }
  157     case t_QFI:
  158     case t_QFR:
  159     {
  160       pari_sp av = avma;
  161       GEN u, D = nf_get_disc(nf), T = nf_get_pol(nf), f = nf_get_index(nf);
  162       GEN A = gel(x,1), B = gel(x,2);
  163       N = nf_get_degree(nf);
  164       if (N != 2)
  165         pari_err_TYPE("idealhnf [Qfb for nonquadratic fields]", x);
  166       if (!equalii(qfb_disc(x), D))
  167         pari_err_DOMAIN("idealhnf [Qfb]", "disc(q)", "!=", D, x);
  168       /* x -> A Z + (-B + sqrt(D)) / 2 Z
  169          K = Q[t]/T(t), t^2 + ut + v = 0,  u^2 - 4v = Df^2
  170          => t = (-u + sqrt(D) f)/2
  171          => sqrt(D)/2 = (t + u/2)/f */
  172       u = gel(T,3);
  173       B = deg1pol_shallow(ginv(f),
  174                           gsub(gdiv(u, shifti(f,1)), gdiv(B,gen_2)),
  175                           varn(T));
  176       return gerepileupto(av, idealhnf_two(nf, mkvec2(A,B)));
  177     }
  178     default: return idealhnf_principal(nf, x); /* PRINCIPAL */
  179   }
  180 }
  181 GEN
  182 idealhnf(GEN nf, GEN x)
  183 {
  184   pari_sp av = avma;
  185   GEN y = idealhnf_shallow(checknf(nf), x);
  186   return (avma == av)? gcopy(y): gerepileupto(av, y);
  187 }
  188 
  189 /* GP functions */
  190 
  191 GEN
  192 idealtwoelt0(GEN nf, GEN x, GEN a)
  193 {
  194   if (!a) return idealtwoelt(nf,x);
  195   return idealtwoelt2(nf,x,a);
  196 }
  197 
  198 GEN
  199 idealpow0(GEN nf, GEN x, GEN n, long flag)
  200 {
  201   if (flag) return idealpowred(nf,x,n);
  202   return idealpow(nf,x,n);
  203 }
  204 
  205 GEN
  206 idealmul0(GEN nf, GEN x, GEN y, long flag)
  207 {
  208   if (flag) return idealmulred(nf,x,y);
  209   return idealmul(nf,x,y);
  210 }
  211 
  212 GEN
  213 idealdiv0(GEN nf, GEN x, GEN y, long flag)
  214 {
  215   switch(flag)
  216   {
  217     case 0: return idealdiv(nf,x,y);
  218     case 1: return idealdivexact(nf,x,y);
  219     default: pari_err_FLAG("idealdiv");
  220   }
  221   return NULL; /* LCOV_EXCL_LINE */
  222 }
  223 
  224 GEN
  225 idealaddtoone0(GEN nf, GEN arg1, GEN arg2)
  226 {
  227   if (!arg2) return idealaddmultoone(nf,arg1);
  228   return idealaddtoone(nf,arg1,arg2);
  229 }
  230 
  231 /* b not a scalar */
  232 static GEN
  233 hnf_Z_ZC(GEN nf, GEN a, GEN b) { return hnfmodid(zk_multable(nf,b), a); }
  234 /* b not a scalar */
  235 static GEN
  236 hnf_Z_QC(GEN nf, GEN a, GEN b)
  237 {
  238   GEN db;
  239   b = Q_remove_denom(b, &db);
  240   if (db) a = mulii(a, db);
  241   b = hnf_Z_ZC(nf,a,b);
  242   return db? RgM_Rg_div(b, db): b;
  243 }
  244 /* b not a scalar (not point in trying to optimize for this case) */
  245 static GEN
  246 hnf_Q_QC(GEN nf, GEN a, GEN b)
  247 {
  248   GEN da, db;
  249   if (typ(a) == t_INT) return hnf_Z_QC(nf, a, b);
  250   da = gel(a,2);
  251   a = gel(a,1);
  252   b = Q_remove_denom(b, &db);
  253   /* write da = d*A, db = d*B, gcd(A,B) = 1
  254    * gcd(a/(d A), b/(d B)) = gcd(a B, A b) / A B d = gcd(a B, b) / A B d */
  255   if (db)
  256   {
  257     GEN d = gcdii(da,db);
  258     if (!is_pm1(d)) db = diviiexact(db,d); /* B */
  259     if (!is_pm1(db))
  260     {
  261       a = mulii(a, db); /* a B */
  262       da = mulii(da, db); /* A B d = lcm(denom(a),denom(b)) */
  263     }
  264   }
  265   return RgM_Rg_div(hnf_Z_ZC(nf,a,b), da);
  266 }
  267 static GEN
  268 hnf_QC_QC(GEN nf, GEN a, GEN b)
  269 {
  270   GEN da, db, d, x;
  271   a = Q_remove_denom(a, &da);
  272   b = Q_remove_denom(b, &db);
  273   if (da) b = ZC_Z_mul(b, da);
  274   if (db) a = ZC_Z_mul(a, db);
  275   d = mul_denom(da, db);
  276   a = zk_multable(nf,a); da = zkmultable_capZ(a);
  277   b = zk_multable(nf,b); db = zkmultable_capZ(b);
  278   x = ZM_hnfmodid(shallowconcat(a,b), gcdii(da,db));
  279   return d? RgM_Rg_div(x, d): x;
  280 }
  281 static GEN
  282 hnf_Q_Q(GEN nf, GEN a, GEN b) {return scalarmat(Q_gcd(a,b), nf_get_degree(nf));}
  283 GEN
  284 idealhnf0(GEN nf, GEN a, GEN b)
  285 {
  286   long ta, tb;
  287   pari_sp av;
  288   GEN x;
  289   if (!b) return idealhnf(nf,a);
  290 
  291   /* HNF of aZ_K+bZ_K */
  292   av = avma; nf = checknf(nf);
  293   a = nf_to_scalar_or_basis(nf,a); ta = typ(a);
  294   b = nf_to_scalar_or_basis(nf,b); tb = typ(b);
  295   if (ta == t_COL)
  296     x = (tb==t_COL)? hnf_QC_QC(nf, a,b): hnf_Q_QC(nf, b,a);
  297   else
  298     x = (tb==t_COL)? hnf_Q_QC(nf, a,b): hnf_Q_Q(nf, a,b);
  299   return gerepileupto(av, x);
  300 }
  301 
  302 /*******************************************************************/
  303 /*                                                                 */
  304 /*                       TWO-ELEMENT FORM                          */
  305 /*                                                                 */
  306 /*******************************************************************/
  307 static GEN idealapprfact_i(GEN nf, GEN x, int nored);
  308 
  309 static int
  310 ok_elt(GEN x, GEN xZ, GEN y)
  311 {
  312   pari_sp av = avma;
  313   return gc_bool(av, ZM_equal(x, ZM_hnfmodid(y, xZ)));
  314 }
  315 
  316 static GEN
  317 addmul_col(GEN a, long s, GEN b)
  318 {
  319   long i,l;
  320   if (!s) return a? leafcopy(a): a;
  321   if (!a) return gmulsg(s,b);
  322   l = lg(a);
  323   for (i=1; i<l; i++)
  324     if (signe(gel(b,i))) gel(a,i) = addii(gel(a,i), mulsi(s, gel(b,i)));
  325   return a;
  326 }
  327 
  328 /* a <-- a + s * b, all coeffs integers */
  329 static GEN
  330 addmul_mat(GEN a, long s, GEN b)
  331 {
  332   long j,l;
  333   /* copy otherwise next call corrupts a */
  334   if (!s) return a? RgM_shallowcopy(a): a;
  335   if (!a) return gmulsg(s,b);
  336   l = lg(a);
  337   for (j=1; j<l; j++)
  338     (void)addmul_col(gel(a,j), s, gel(b,j));
  339   return a;
  340 }
  341 
  342 static GEN
  343 get_random_a(GEN nf, GEN x, GEN xZ)
  344 {
  345   pari_sp av;
  346   long i, lm, l = lg(x);
  347   GEN a, z, beta, mul;
  348 
  349   beta= cgetg(l, t_VEC);
  350   mul = cgetg(l, t_VEC); lm = 1; /* = lg(mul) */
  351   /* look for a in x such that a O/xZ = x O/xZ */
  352   for (i = 2; i < l; i++)
  353   {
  354     GEN xi = gel(x,i);
  355     GEN t = FpM_red(zk_multable(nf,xi), xZ); /* ZM, cannot be a scalar */
  356     if (gequal0(t)) continue;
  357     if (ok_elt(x,xZ, t)) return xi;
  358     gel(beta,lm) = xi;
  359     /* mul[i] = { canonical generators for x[i] O/xZ as Z-module } */
  360     gel(mul,lm) = t; lm++;
  361   }
  362   setlg(mul, lm);
  363   setlg(beta,lm);
  364   z = cgetg(lm, t_VECSMALL);
  365   for(av = avma;; set_avma(av))
  366   {
  367     for (a=NULL,i=1; i<lm; i++)
  368     {
  369       long t = random_bits(4) - 7; /* in [-7,8] */
  370       z[i] = t;
  371       a = addmul_mat(a, t, gel(mul,i));
  372     }
  373     /* a = matrix (NOT HNF) of ideal generated by beta.z in O/xZ */
  374     if (a && ok_elt(x,xZ, a)) break;
  375   }
  376   for (a=NULL,i=1; i<lm; i++)
  377     a = addmul_col(a, z[i], gel(beta,i));
  378   return a;
  379 }
  380 
  381 /* x square matrix, assume it is HNF */
  382 static GEN
  383 mat_ideal_two_elt(GEN nf, GEN x)
  384 {
  385   GEN y, a, cx, xZ;
  386   long N = nf_get_degree(nf);
  387   pari_sp av, tetpil;
  388 
  389   if (lg(x)-1 != N) pari_err_DIM("idealtwoelt");
  390   if (N == 2) return mkvec2copy(gcoeff(x,1,1), gel(x,2));
  391 
  392   y = cgetg(3,t_VEC); av = avma;
  393   cx = Q_content(x);
  394   xZ = gcoeff(x,1,1);
  395   if (gequal(xZ, cx)) /* x = (cx) */
  396   {
  397     gel(y,1) = cx;
  398     gel(y,2) = gen_0; return y;
  399   }
  400   if (equali1(cx)) cx = NULL;
  401   else
  402   {
  403     x = Q_div_to_int(x, cx);
  404     xZ = gcoeff(x,1,1);
  405   }
  406   if (N < 6)
  407     a = get_random_a(nf, x, xZ);
  408   else
  409   {
  410     const long FB[] = { _evallg(15+1) | evaltyp(t_VECSMALL),
  411       2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47
  412     };
  413     GEN P, E, a1 = Z_smoothen(xZ, (GEN)FB, &P, &E);
  414     if (!a1) /* factors completely */
  415       a = idealapprfact_i(nf, idealfactor(nf,x), 1);
  416     else if (lg(P) == 1) /* no small factors */
  417       a = get_random_a(nf, x, xZ);
  418     else /* general case */
  419     {
  420       GEN A0, A1, a0, u0, u1, v0, v1, pi0, pi1, t, u;
  421       a0 = diviiexact(xZ, a1);
  422       A0 = ZM_hnfmodid(x, a0); /* smooth part of x */
  423       A1 = ZM_hnfmodid(x, a1); /* cofactor */
  424       pi0 = idealapprfact_i(nf, idealfactor(nf,A0), 1);
  425       pi1 = get_random_a(nf, A1, a1);
  426       (void)bezout(a0, a1, &v0,&v1);
  427       u0 = mulii(a0, v0);
  428       u1 = mulii(a1, v1);
  429       if (typ(pi0) != t_COL) t = addmulii(u0, pi0, u1);
  430       else
  431       { t = ZC_Z_mul(pi0, u1); gel(t,1) = addii(gel(t,1), u0); }
  432       u = ZC_Z_mul(pi1, u0); gel(u,1) = addii(gel(u,1), u1);
  433       a = nfmuli(nf, centermod(u, xZ), centermod(t, xZ));
  434     }
  435   }
  436   if (cx)
  437   {
  438     a = centermod(a, xZ);
  439     tetpil = avma;
  440     if (typ(cx) == t_INT)
  441     {
  442       gel(y,1) = mulii(xZ, cx);
  443       gel(y,2) = ZC_Z_mul(a, cx);
  444     }
  445     else
  446     {
  447       gel(y,1) = gmul(xZ, cx);
  448       gel(y,2) = RgC_Rg_mul(a, cx);
  449     }
  450   }
  451   else
  452   {
  453     tetpil = avma;
  454     gel(y,1) = icopy(xZ);
  455     gel(y,2) = centermod(a, xZ);
  456   }
  457   gerepilecoeffssp(av,tetpil,y+1,2); return y;
  458 }
  459 
  460 /* Given an ideal x, returns [a,alpha] such that a is in Q,
  461  * x = a Z_K + alpha Z_K, alpha in K^*
  462  * a = 0 or alpha = 0 are possible, but do not try to determine whether
  463  * x is principal. */
  464 GEN
  465 idealtwoelt(GEN nf, GEN x)
  466 {
  467   pari_sp av;
  468   GEN z;
  469   long tx = idealtyp(&x,&z);
  470   nf = checknf(nf);
  471   if (tx == id_MAT) return mat_ideal_two_elt(nf,x);
  472   if (tx == id_PRIME) return mkvec2copy(gel(x,1), gel(x,2));
  473   /* id_PRINCIPAL */
  474   av = avma; x = nf_to_scalar_or_basis(nf, x);
  475   return gerepilecopy(av, typ(x)==t_COL? mkvec2(gen_0,x):
  476                                          mkvec2(Q_abs_shallow(x),gen_0));
  477 }
  478 
  479 /*******************************************************************/
  480 /*                                                                 */
  481 /*                         FACTORIZATION                           */
  482 /*                                                                 */
  483 /*******************************************************************/
  484 /* x integral ideal in HNF, Zval = v_p(x \cap Z) > 0; return v_p(Nx) */
  485 static long
  486 idealHNF_norm_pval(GEN x, GEN p, long Zval)
  487 {
  488   long i, v = Zval, l = lg(x);
  489   for (i = 2; i < l; i++) v += Z_pval(gcoeff(x,i,i), p);
  490   return v;
  491 }
  492 
  493 /* x integral in HNF, f0 = partial factorization of a multiple of
  494  * x[1,1] = x\cap Z */
  495 GEN
  496 idealHNF_Z_factor_i(GEN x, GEN f0, GEN *pvN, GEN *pvZ)
  497 {
  498   GEN P, E, vN, vZ, xZ = gcoeff(x,1,1), f = f0? f0: Z_factor(xZ);
  499   long i, l;
  500   P = gel(f,1); l = lg(P);
  501   E = gel(f,2);
  502   *pvN = vN = cgetg(l, t_VECSMALL);
  503   *pvZ = vZ = cgetg(l, t_VECSMALL);
  504   for (i = 1; i < l; i++)
  505   {
  506     GEN p = gel(P,i);
  507     vZ[i] = f0? Z_pval(xZ, p): (long) itou(gel(E,i));
  508     vN[i] = idealHNF_norm_pval(x,p, vZ[i]);
  509   }
  510   return P;
  511 }
  512 /* return P, primes dividing Nx and xZ = x\cap Z, set v_p(Nx), v_p(xZ);
  513  * x integral in HNF */
  514 GEN
  515 idealHNF_Z_factor(GEN x, GEN *pvN, GEN *pvZ)
  516 { return idealHNF_Z_factor_i(x, NULL, pvN, pvZ); }
  517 
  518 /* v_P(A)*f(P) <= Nval [e.g. Nval = v_p(Norm A)], Zval = v_p(A \cap Z).
  519  * Return v_P(A) */
  520 static long
  521 idealHNF_val(GEN A, GEN P, long Nval, long Zval)
  522 {
  523   long f = pr_get_f(P), vmax, v, e, i, j, k, l;
  524   GEN mul, B, a, y, r, p, pk, cx, vals;
  525   pari_sp av;
  526 
  527   if (Nval < f) return 0;
  528   p = pr_get_p(P);
  529   e = pr_get_e(P);
  530   /* v_P(A) <= max [ e * v_p(A \cap Z), floor[v_p(Nix) / f ] */
  531   vmax = minss(Zval * e, Nval / f);
  532   mul = pr_get_tau(P);
  533   l = lg(mul);
  534   B = cgetg(l,t_MAT);
  535   /* B[1] not needed: v_pr(A[1]) = v_pr(A \cap Z) is known already */
  536   gel(B,1) = gen_0; /* dummy */
  537   for (j = 2; j < l; j++)
  538   {
  539     GEN x = gel(A,j);
  540     gel(B,j) = y = cgetg(l, t_COL);
  541     for (i = 1; i < l; i++)
  542     { /* compute a = (x.t0)_i, A in HNF ==> x[j+1..l-1] = 0 */
  543       a = mulii(gel(x,1), gcoeff(mul,i,1));
  544       for (k = 2; k <= j; k++) a = addii(a, mulii(gel(x,k), gcoeff(mul,i,k)));
  545       /* p | a ? */
  546       gel(y,i) = dvmdii(a,p,&r); if (signe(r)) return 0;
  547     }
  548   }
  549   vals = cgetg(l, t_VECSMALL);
  550   /* vals[1] not needed */
  551   for (j = 2; j < l; j++)
  552   {
  553     gel(B,j) = Q_primitive_part(gel(B,j), &cx);
  554     vals[j] = cx? 1 + e * Q_pval(cx, p): 1;
  555   }
  556   pk = powiu(p, ceildivuu(vmax, e));
  557   av = avma; y = cgetg(l,t_COL);
  558   /* can compute mod p^ceil((vmax-v)/e) */
  559   for (v = 1; v < vmax; v++)
  560   { /* we know v_pr(Bj) >= v for all j */
  561     if (e == 1 || (vmax - v) % e == 0) pk = diviiexact(pk, p);
  562     for (j = 2; j < l; j++)
  563     {
  564       GEN x = gel(B,j); if (v < vals[j]) continue;
  565       for (i = 1; i < l; i++)
  566       {
  567         pari_sp av2 = avma;
  568         a = mulii(gel(x,1), gcoeff(mul,i,1));
  569         for (k = 2; k < l; k++) a = addii(a, mulii(gel(x,k), gcoeff(mul,i,k)));
  570         /* a = (x.t_0)_i; p | a ? */
  571         a = dvmdii(a,p,&r); if (signe(r)) return v;
  572         if (lgefint(a) > lgefint(pk)) a = remii(a, pk);
  573         gel(y,i) = gerepileuptoint(av2, a);
  574       }
  575       gel(B,j) = y; y = x;
  576       if (gc_needed(av,3))
  577       {
  578         if(DEBUGMEM>1) pari_warn(warnmem,"idealval");
  579         gerepileall(av,3, &y,&B,&pk);
  580       }
  581     }
  582   }
  583   return v;
  584 }
  585 /* true nf, x != 0 integral ideal in HNF, cx t_INT or NULL,
  586  * FA integer factorization matrix or NULL. Return partial factorization of
  587  * cx * x above primes in FA (complete factorization if !FA)*/
  588 static GEN
  589 idealHNF_factor_i(GEN nf, GEN x, GEN cx, GEN FA)
  590 {
  591   const long N = lg(x)-1;
  592   long i, j, k, l, v;
  593   GEN vN, vZ, vP, vE, vp = idealHNF_Z_factor_i(x, FA, &vN,&vZ);
  594 
  595   l = lg(vp);
  596   i = cx? expi(cx)+1: 1;
  597   vP = cgetg((l+i-2)*N+1, t_COL);
  598   vE = cgetg((l+i-2)*N+1, t_COL);
  599   for (i = k = 1; i < l; i++)
  600   {
  601     GEN L, p = gel(vp,i);
  602     long Nval = vN[i], Zval = vZ[i], vc = cx? Z_pvalrem(cx,p,&cx): 0;
  603     if (vc)
  604     {
  605       L = idealprimedec(nf,p);
  606       if (is_pm1(cx)) cx = NULL;
  607     }
  608     else
  609       L = idealprimedec_limit_f(nf,p,Nval);
  610     for (j = 1; Nval && j < lg(L); j++) /* !Nval => only cx contributes */
  611     {
  612       GEN P = gel(L,j);
  613       pari_sp av = avma;
  614       v = idealHNF_val(x, P, Nval, Zval);
  615       set_avma(av);
  616       Nval -= v*pr_get_f(P);
  617       v += vc * pr_get_e(P); if (!v) continue;
  618       gel(vP,k) = P;
  619       gel(vE,k) = utoipos(v); k++;
  620     }
  621     if (vc) for (; j<lg(L); j++)
  622     {
  623       GEN P = gel(L,j);
  624       gel(vP,k) = P;
  625       gel(vE,k) = utoipos(vc * pr_get_e(P)); k++;
  626     }
  627   }
  628   if (cx && !FA)
  629   { /* complete factorization */
  630     GEN f = Z_factor(cx), cP = gel(f,1), cE = gel(f,2);
  631     long lc = lg(cP);
  632     for (i=1; i<lc; i++)
  633     {
  634       GEN p = gel(cP,i), L = idealprimedec(nf,p);
  635       long vc = itos(gel(cE,i));
  636       for (j=1; j<lg(L); j++)
  637       {
  638         GEN P = gel(L,j);
  639         gel(vP,k) = P;
  640         gel(vE,k) = utoipos(vc * pr_get_e(P)); k++;
  641       }
  642     }
  643   }
  644   setlg(vP, k);
  645   setlg(vE, k); return mkmat2(vP, vE);
  646 }
  647 /* true nf, x integral ideal */
  648 static GEN
  649 idealHNF_factor(GEN nf, GEN x, ulong lim)
  650 {
  651   GEN cx, F = NULL;
  652   if (lim)
  653   {
  654     GEN P, E;
  655     long i;
  656     /* strict useless because of prime table */
  657     F = absZ_factor_limit(gcoeff(x,1,1), lim);
  658     P = gel(F,1);
  659     E = gel(F,2);
  660     /* filter out entries > lim */
  661     for (i = lg(P)-1; i; i--)
  662       if (cmpiu(gel(P,i), lim) <= 0) break;
  663     setlg(P, i+1);
  664     setlg(E, i+1);
  665   }
  666   x = Q_primitive_part(x, &cx);
  667   return idealHNF_factor_i(nf, x, cx, F);
  668 }
  669 /* c * vector(#L,i,L[i].e), assume results fit in ulong */
  670 static GEN
  671 prV_e_muls(GEN L, long c)
  672 {
  673   long j, l = lg(L);
  674   GEN z = cgetg(l, t_COL);
  675   for (j = 1; j < l; j++) gel(z,j) = stoi(c * pr_get_e(gel(L,j)));
  676   return z;
  677 }
  678 /* true nf, y in Q */
  679 static GEN
  680 Q_nffactor(GEN nf, GEN y, ulong lim)
  681 {
  682   GEN f, P, E;
  683   long l, i;
  684   if (typ(y) == t_INT)
  685   {
  686     if (!signe(y)) pari_err_DOMAIN("idealfactor", "ideal", "=",gen_0,y);
  687     if (is_pm1(y)) return trivial_fact();
  688   }
  689   y = Q_abs_shallow(y);
  690   if (!lim) f = Q_factor(y);
  691   else
  692   {
  693     f = Q_factor_limit(y, lim);
  694     P = gel(f,1);
  695     E = gel(f,2);
  696     for (i = lg(P)-1; i > 0; i--)
  697       if (abscmpiu(gel(P,i), lim) < 0) break;
  698     setlg(P,i+1); setlg(E,i+1);
  699   }
  700   P = gel(f,1); l = lg(P); if (l == 1) return f;
  701   E = gel(f,2);
  702   for (i = 1; i < l; i++)
  703   {
  704     gel(P,i) = idealprimedec(nf, gel(P,i));
  705     gel(E,i) = prV_e_muls(gel(P,i), itos(gel(E,i)));
  706   }
  707   settyp(P,t_VEC); P = shallowconcat1(P);
  708   settyp(E,t_VEC); E = shallowconcat1(E);
  709   gel(f,1) = P; settyp(P, t_COL);
  710   gel(f,2) = E; return f;
  711 }
  712 
  713 GEN
  714 idealfactor_partial(GEN nf, GEN x, GEN L)
  715 {
  716   pari_sp av = avma;
  717   long i, j, l;
  718   GEN P, E;
  719   if (!L) return idealfactor(nf, x);
  720   if (typ(L) == t_INT) return idealfactor_limit(nf, x, itou(L));
  721   l = lg(L); if (l == 1) return trivial_fact();
  722   P = cgetg(l, t_VEC);
  723   for (i = 1; i < l; i++)
  724   {
  725     GEN p = gel(L,i);
  726     gel(P,i) = typ(p) == t_INT? idealprimedec(nf, p): mkvec(p);
  727   }
  728   P = shallowconcat1(P); settyp(P, t_COL);
  729   P = gen_sort_uniq(P, (void*)&cmp_prime_ideal, &cmp_nodata);
  730   E = cgetg_copy(P, &l);
  731   for (i = j = 1; i < l; i++)
  732   {
  733     long v = idealval(nf, x, gel(P,i));
  734     if (v) { gel(P,j) = gel(P,i); gel(E,j) = stoi(v); j++; }
  735   }
  736   setlg(P,j);
  737   setlg(E,j); return gerepilecopy(av, mkmat2(P, E));
  738 }
  739 GEN
  740 idealfactor_limit(GEN nf, GEN x, ulong lim)
  741 {
  742   pari_sp av = avma;
  743   GEN fa, y;
  744   long tx = idealtyp(&x,&y);
  745 
  746   if (tx == id_PRIME)
  747   {
  748     if (lim && abscmpiu(pr_get_p(x), lim) >= 0) return trivial_fact();
  749     retmkmat2(mkcolcopy(x), mkcol(gen_1));
  750   }
  751   nf = checknf(nf);
  752   if (tx == id_PRINCIPAL)
  753   {
  754     y = nf_to_scalar_or_basis(nf, x);
  755     if (typ(y) != t_COL) return gerepilecopy(av, Q_nffactor(nf, y, lim));
  756   }
  757   y = idealnumden(nf, x);
  758   fa = idealHNF_factor(nf, gel(y,1), lim);
  759   if (!isint1(gel(y,2)))
  760     fa = famat_div_shallow(fa, idealHNF_factor(nf, gel(y,2), lim));
  761   fa = gerepilecopy(av, fa);
  762   return sort_factor(fa, (void*)&cmp_prime_ideal, &cmp_nodata);
  763 }
  764 GEN
  765 idealfactor(GEN nf, GEN x) { return idealfactor_limit(nf, x, 0); }
  766 GEN
  767 gpidealfactor(GEN nf, GEN x, GEN lim)
  768 {
  769   ulong L = 0;
  770   if (lim)
  771   {
  772     if (typ(lim) != t_INT || signe(lim) < 0) pari_err_FLAG("idealfactor");
  773     L = itou(lim);
  774   }
  775   return idealfactor_limit(nf, x, L);
  776 }
  777 
  778 static GEN
  779 ramified_root(GEN nf, GEN R, GEN A, long n)
  780 {
  781   GEN v, P = gel(idealfactor(nf, R), 1);
  782   long i, l = lg(P);
  783   v = cgetg(l, t_VECSMALL);
  784   for (i = 1; i < l; i++)
  785   {
  786     long w = idealval(nf, A, gel(P,i));
  787     if (w % n) return NULL;
  788     v[i] = w / n;
  789   }
  790   return idealfactorback(nf, P, v, 0);
  791 }
  792 static int
  793 ramified_root_simple(GEN nf, long n, GEN P, GEN v)
  794 {
  795   long i, l = lg(v);
  796   for (i = 1; i < l; i++) if (v[i])
  797   {
  798     GEN vpr = idealprimedec(nf, gel(P,i));
  799     long lpr = lg(vpr), j;
  800     for (j = 1; j < lpr; j++)
  801     {
  802       long e = pr_get_e(gel(vpr,j));
  803       if ((e * v[i]) % n) return 0;
  804     }
  805   }
  806   return 1;
  807 }
  808 /* true nf; A is assumed to be the n-th power of an integral ideal,
  809  * return its n-th root; n > 1 */
  810 static long
  811 idealsqrtn_int(GEN nf, GEN A, long n, GEN *pB)
  812 {
  813   GEN C, root;
  814   long i, l;
  815 
  816   if (typ(A) == t_INT) /* > 0 */
  817   {
  818     GEN P = nf_get_ramified_primes(nf), v, q;
  819     l = lg(P); v = cgetg(l, t_VECSMALL);
  820     for (i = 1; i < l; i++) v[i] = Z_pvalrem(A, gel(P,i), &A);
  821     C = gen_1;
  822     if (!isint1(A) && !Z_ispowerall(A, n, pB? &C: NULL)) return 0;
  823     if (!pB) return ramified_root_simple(nf, n, P, v);
  824     q = factorback2(P, v);
  825     root = ramified_root(nf, q, q, n);
  826     if (!root) return 0;
  827     if (!equali1(C)) root = isint1(root)? C: ZM_Z_mul(root, C);
  828     *pB = root; return 1;
  829   }
  830   /* compute valuations at ramified primes */
  831   root = ramified_root(nf, idealadd(nf, nf_get_diff(nf), A), A, n);
  832   if (!root) return 0;
  833   /* remove ramified primes */
  834   if (isint1(root))
  835     root = matid(nf_get_degree(nf));
  836   else
  837     A = idealdivexact(nf, A, idealpows(nf,root,n));
  838   A = Q_primitive_part(A, &C);
  839   if (C)
  840   {
  841     if (!Z_ispowerall(C,n,&C)) return 0;
  842     if (pB) root = ZM_Z_mul(root, C);
  843   }
  844 
  845   /* compute final n-th root, at most degree(nf)-1 iterations */
  846   for (i = 0;; i++)
  847   {
  848     GEN J, b, a = gcoeff(A,1,1); /* A \cap Z */
  849     if (is_pm1(a)) break;
  850     if (!Z_ispowerall(a,n,&b)) return 0;
  851     J = idealadd(nf, b, A);
  852     A = idealdivexact(nf, idealpows(nf,J,n), A);
  853     /* div and not divexact here */
  854     if (pB) root = odd(i)? idealdiv(nf, root, J): idealmul(nf, root, J);
  855   }
  856   if (pB) *pB = root;
  857   return 1;
  858 }
  859 
  860 /* A is assumed to be the n-th power of an ideal in nf
  861  returns its n-th root. */
  862 long
  863 idealispower(GEN nf, GEN A, long n, GEN *pB)
  864 {
  865   pari_sp av = avma;
  866   GEN v, N, D;
  867   nf = checknf(nf);
  868   if (n <= 0) pari_err_DOMAIN("idealispower", "n", "<=", gen_0, stoi(n));
  869   if (n == 1) { if (pB) *pB = idealhnf(nf,A); return 1; }
  870   v = idealnumden(nf,A);
  871   if (gequal0(gel(v,1))) { set_avma(av); if (pB) *pB = cgetg(1,t_MAT); return 1; }
  872   if (!idealsqrtn_int(nf, gel(v,1), n, pB? &N: NULL)) return 0;
  873   if (!idealsqrtn_int(nf, gel(v,2), n, pB? &D: NULL)) return 0;
  874   if (pB) *pB = gerepileupto(av, idealdiv(nf,N,D)); else set_avma(av);
  875   return 1;
  876 }
  877 
  878 /* x t_INT or integral nonzero ideal in HNF */
  879 static GEN
  880 idealredmodpower_i(GEN nf, GEN x, ulong k, ulong B)
  881 {
  882   GEN cx, y, U, N, F, Q;
  883   if (typ(x) == t_INT)
  884   {
  885     if (!signe(x) || is_pm1(x)) return gen_1;
  886     F = Z_factor_limit(x, B);
  887     gel(F,2) = gdiventgs(gel(F,2), k);
  888     return ginv(factorback(F));
  889   }
  890   N = gcoeff(x,1,1); if (is_pm1(N)) return gen_1;
  891   F = absZ_factor_limit_strict(N, B, &U);
  892   if (U)
  893   {
  894     GEN M = powii(gel(U,1), gel(U,2));
  895     y = hnfmodid(x, M); /* coprime part to B! */
  896     if (!idealispower(nf, y, k, &U)) U = NULL;
  897     x = hnfmodid(x, diviiexact(N, M));
  898   }
  899   /* x = B-smooth part of initial x */
  900   x = Q_primitive_part(x, &cx);
  901   F = idealHNF_factor_i(nf, x, cx, F);
  902   gel(F,2) = gdiventgs(gel(F,2), k);
  903   Q = idealfactorback(nf, gel(F,1), gel(F,2), 0);
  904   if (U) Q = idealmul(nf,Q,U);
  905   if (typ(Q) == t_INT) return Q;
  906   y = idealred_elt(nf, idealHNF_inv_Z(nf, Q));
  907   return gdiv(y, gcoeff(Q,1,1));
  908 }
  909 GEN
  910 idealredmodpower(GEN nf, GEN x, ulong n, ulong B)
  911 {
  912   pari_sp av = avma;
  913   GEN a, b;
  914   nf = checknf(nf);
  915   if (!n) pari_err_DOMAIN("idealredmodpower","n", "=", gen_0, gen_0);
  916   x = idealnumden(nf, x);
  917   a = gel(x,1);
  918   if (isintzero(a)) { set_avma(av); return gen_1; }
  919   a = idealredmodpower_i(nf, gel(x,1), n, B);
  920   b = idealredmodpower_i(nf, gel(x,2), n, B);
  921   if (!isint1(b)) a = nf_to_scalar_or_basis(nf, nfdiv(nf, a, b));
  922   return gerepilecopy(av, a);
  923 }
  924 
  925 /* P prime ideal in idealprimedec format. Return valuation(A) at P */
  926 long
  927 idealval(GEN nf, GEN A, GEN P)
  928 {
  929   pari_sp av = avma;
  930   GEN a, p, cA;
  931   long vcA, v, Zval, tx = idealtyp(&A,&a);
  932 
  933   if (tx == id_PRINCIPAL) return nfval(nf,A,P);
  934   checkprid(P);
  935   if (tx == id_PRIME) return pr_equal(P, A)? 1: 0;
  936   /* id_MAT */
  937   nf = checknf(nf);
  938   A = Q_primitive_part(A, &cA);
  939   p = pr_get_p(P);
  940   vcA = cA? Q_pval(cA,p): 0;
  941   if (pr_is_inert(P)) return gc_long(av,vcA);
  942   Zval = Z_pval(gcoeff(A,1,1), p);
  943   if (!Zval) v = 0;
  944   else
  945   {
  946     long Nval = idealHNF_norm_pval(A, p, Zval);
  947     v = idealHNF_val(A, P, Nval, Zval);
  948   }
  949   return gc_long(av, vcA? v + vcA*pr_get_e(P): v);
  950 }
  951 GEN
  952 gpidealval(GEN nf, GEN ix, GEN P)
  953 {
  954   long v = idealval(nf,ix,P);
  955   return v == LONG_MAX? mkoo(): stoi(v);
  956 }
  957 
  958 /* gcd and generalized Bezout */
  959 
  960 GEN
  961 idealadd(GEN nf, GEN x, GEN y)
  962 {
  963   pari_sp av = avma;
  964   long tx, ty;
  965   GEN z, a, dx, dy, dz;
  966 
  967   tx = idealtyp(&x,&z);
  968   ty = idealtyp(&y,&z); nf = checknf(nf);
  969   if (tx != id_MAT) x = idealhnf_shallow(nf,x);
  970   if (ty != id_MAT) y = idealhnf_shallow(nf,y);
  971   if (lg(x) == 1) return gerepilecopy(av,y);
  972   if (lg(y) == 1) return gerepilecopy(av,x); /* check for 0 ideal */
  973   dx = Q_denom(x);
  974   dy = Q_denom(y); dz = lcmii(dx,dy);
  975   if (is_pm1(dz)) dz = NULL; else {
  976     x = Q_muli_to_int(x, dz);
  977     y = Q_muli_to_int(y, dz);
  978   }
  979   a = gcdii(gcoeff(x,1,1), gcoeff(y,1,1));
  980   if (is_pm1(a))
  981   {
  982     long N = lg(x)-1;
  983     if (!dz) { set_avma(av); return matid(N); }
  984     return gerepileupto(av, scalarmat(ginv(dz), N));
  985   }
  986   z = ZM_hnfmodid(shallowconcat(x,y), a);
  987   if (dz) z = RgM_Rg_div(z,dz);
  988   return gerepileupto(av,z);
  989 }
  990 
  991 static GEN
  992 trivial_merge(GEN x)
  993 { return (lg(x) == 1 || !is_pm1(gcoeff(x,1,1)))? NULL: gen_1; }
  994 /* true nf */
  995 static GEN
  996 _idealaddtoone(GEN nf, GEN x, GEN y, long red)
  997 {
  998   GEN a;
  999   long tx = idealtyp(&x, &a/*junk*/);
 1000   long ty = idealtyp(&y, &a/*junk*/);
 1001   long ea;
 1002   if (tx != id_MAT) x = idealhnf_shallow(nf, x);
 1003   if (ty != id_MAT) y = idealhnf_shallow(nf, y);
 1004   if (lg(x) == 1)
 1005     a = trivial_merge(y);
 1006   else if (lg(y) == 1)
 1007     a = trivial_merge(x);
 1008   else
 1009     a = hnfmerge_get_1(x, y);
 1010   if (!a) pari_err_COPRIME("idealaddtoone",x,y);
 1011   if (red && (ea = gexpo(a)) > 10)
 1012   {
 1013     GEN b = (typ(a) == t_COL)? a: scalarcol_shallow(a, nf_get_degree(nf));
 1014     b = ZC_reducemodlll(b, idealHNF_mul(nf,x,y));
 1015     if (gexpo(b) < ea) a = b;
 1016   }
 1017   return a;
 1018 }
 1019 /* true nf */
 1020 GEN
 1021 idealaddtoone_i(GEN nf, GEN x, GEN y)
 1022 { return _idealaddtoone(nf, x, y, 1); }
 1023 /* true nf */
 1024 GEN
 1025 idealaddtoone_raw(GEN nf, GEN x, GEN y)
 1026 { return _idealaddtoone(nf, x, y, 0); }
 1027 
 1028 GEN
 1029 idealaddtoone(GEN nf, GEN x, GEN y)
 1030 {
 1031   GEN z = cgetg(3,t_VEC), a;
 1032   pari_sp av = avma;
 1033   nf = checknf(nf);
 1034   a = gerepileupto(av, idealaddtoone_i(nf,x,y));
 1035   gel(z,1) = a;
 1036   gel(z,2) = typ(a) == t_COL? Z_ZC_sub(gen_1,a): subui(1,a);
 1037   return z;
 1038 }
 1039 
 1040 /* assume elements of list are integral ideals */
 1041 GEN
 1042 idealaddmultoone(GEN nf, GEN list)
 1043 {
 1044   pari_sp av = avma;
 1045   long N, i, l, nz, tx = typ(list);
 1046   GEN H, U, perm, L;
 1047 
 1048   nf = checknf(nf); N = nf_get_degree(nf);
 1049   if (!is_vec_t(tx)) pari_err_TYPE("idealaddmultoone",list);
 1050   l = lg(list);
 1051   L = cgetg(l, t_VEC);
 1052   if (l == 1)
 1053     pari_err_DOMAIN("idealaddmultoone", "sum(ideals)", "!=", gen_1, L);
 1054   nz = 0; /* number of nonzero ideals in L */
 1055   for (i=1; i<l; i++)
 1056   {
 1057     GEN I = gel(list,i);
 1058     if (typ(I) != t_MAT) I = idealhnf_shallow(nf,I);
 1059     if (lg(I) != 1)
 1060     {
 1061       nz++; RgM_check_ZM(I,"idealaddmultoone");
 1062       if (lgcols(I) != N+1) pari_err_TYPE("idealaddmultoone [not an ideal]", I);
 1063     }
 1064     gel(L,i) = I;
 1065   }
 1066   H = ZM_hnfperm(shallowconcat1(L), &U, &perm);
 1067   if (lg(H) == 1 || !equali1(gcoeff(H,1,1)))
 1068     pari_err_DOMAIN("idealaddmultoone", "sum(ideals)", "!=", gen_1, L);
 1069   for (i=1; i<=N; i++)
 1070     if (perm[i] == 1) break;
 1071   U = gel(U,(nz-1)*N + i); /* (L[1]|...|L[nz]) U = 1 */
 1072   nz = 0;
 1073   for (i=1; i<l; i++)
 1074   {
 1075     GEN c = gel(L,i);
 1076     if (lg(c) == 1)
 1077       c = gen_0;
 1078     else {
 1079       c = ZM_ZC_mul(c, vecslice(U, nz*N + 1, (nz+1)*N));
 1080       nz++;
 1081     }
 1082     gel(L,i) = c;
 1083   }
 1084   return gerepilecopy(av, L);
 1085 }
 1086 
 1087 /* multiplication */
 1088 
 1089 /* x integral ideal (without archimedean component) in HNF form
 1090  * y = [a,alpha] corresponds to the integral ideal aZ_K+alpha Z_K, a in Z,
 1091  * alpha a ZV or a ZM (multiplication table). Multiply them */
 1092 static GEN
 1093 idealHNF_mul_two(GEN nf, GEN x, GEN y)
 1094 {
 1095   GEN m, a = gel(y,1), alpha = gel(y,2);
 1096   long i, N;
 1097 
 1098   if (typ(alpha) != t_MAT)
 1099   {
 1100     alpha = zk_scalar_or_multable(nf, alpha);
 1101     if (typ(alpha) == t_INT) /* e.g. y inert ? 0 should not (but may) occur */
 1102       return signe(a)? ZM_Z_mul(x, gcdii(a, alpha)): cgetg(1,t_MAT);
 1103   }
 1104   N = lg(x)-1; m = cgetg((N<<1)+1,t_MAT);
 1105   for (i=1; i<=N; i++) gel(m,i)   = ZM_ZC_mul(alpha,gel(x,i));
 1106   for (i=1; i<=N; i++) gel(m,i+N) = ZC_Z_mul(gel(x,i), a);
 1107   return ZM_hnfmodid(m, mulii(a, gcoeff(x,1,1)));
 1108 }
 1109 
 1110 /* Assume ix and iy are integral in HNF form [NOT extended]. Not memory clean.
 1111  * HACK: ideal in iy can be of the form [a,b], a in Z, b in Z_K */
 1112 GEN
 1113 idealHNF_mul(GEN nf, GEN x, GEN y)
 1114 {
 1115   GEN z;
 1116   if (typ(y) == t_VEC)
 1117     z = idealHNF_mul_two(nf,x,y);
 1118   else
 1119   { /* reduce one ideal to two-elt form. The smallest */
 1120     GEN xZ = gcoeff(x,1,1), yZ = gcoeff(y,1,1);
 1121     if (cmpii(xZ, yZ) < 0)
 1122     {
 1123       if (is_pm1(xZ)) return gcopy(y);
 1124       z = idealHNF_mul_two(nf, y, mat_ideal_two_elt(nf,x));
 1125     }
 1126     else
 1127     {
 1128       if (is_pm1(yZ)) return gcopy(x);
 1129       z = idealHNF_mul_two(nf, x, mat_ideal_two_elt(nf,y));
 1130     }
 1131   }
 1132   return z;
 1133 }
 1134 
 1135 /* operations on elements in factored form */
 1136 
 1137 GEN
 1138 famat_mul_shallow(GEN f, GEN g)
 1139 {
 1140   if (typ(f) != t_MAT) f = to_famat_shallow(f,gen_1);
 1141   if (typ(g) != t_MAT) g = to_famat_shallow(g,gen_1);
 1142   if (lgcols(f) == 1) return g;
 1143   if (lgcols(g) == 1) return f;
 1144   return mkmat2(shallowconcat(gel(f,1), gel(g,1)),
 1145                 shallowconcat(gel(f,2), gel(g,2)));
 1146 }
 1147 GEN
 1148 famat_mulpow_shallow(GEN f, GEN g, GEN e)
 1149 {
 1150   if (!signe(e)) return f;
 1151   return famat_mul_shallow(f, famat_pow_shallow(g, e));
 1152 }
 1153 
 1154 GEN
 1155 famat_mulpows_shallow(GEN f, GEN g, long e)
 1156 {
 1157   if (e==0) return f;
 1158   return famat_mul_shallow(f, famat_pows_shallow(g, e));
 1159 }
 1160 
 1161 GEN
 1162 famat_div_shallow(GEN f, GEN g)
 1163 { return famat_mul_shallow(f, famat_inv_shallow(g)); }
 1164 
 1165 GEN
 1166 to_famat(GEN x, GEN y) { retmkmat2(mkcolcopy(x), mkcolcopy(y)); }
 1167 GEN
 1168 to_famat_shallow(GEN x, GEN y) { return mkmat2(mkcol(x), mkcol(y)); }
 1169 
 1170 /* concat the single elt x; not gconcat since x may be a t_COL */
 1171 static GEN
 1172 append(GEN v, GEN x)
 1173 {
 1174   long i, l = lg(v);
 1175   GEN w = cgetg(l+1, typ(v));
 1176   for (i=1; i<l; i++) gel(w,i) = gcopy(gel(v,i));
 1177   gel(w,i) = gcopy(x); return w;
 1178 }
 1179 /* add x^1 to famat f */
 1180 static GEN
 1181 famat_add(GEN f, GEN x)
 1182 {
 1183   GEN h = cgetg(3,t_MAT);
 1184   if (lgcols(f) == 1)
 1185   {
 1186     gel(h,1) = mkcolcopy(x);
 1187     gel(h,2) = mkcol(gen_1);
 1188   }
 1189   else
 1190   {
 1191     gel(h,1) = append(gel(f,1), x);
 1192     gel(h,2) = gconcat(gel(f,2), gen_1);
 1193   }
 1194   return h;
 1195 }
 1196 /* add x^-1 to famat f */
 1197 static GEN
 1198 famat_sub(GEN f, GEN x)
 1199 {
 1200   GEN h = cgetg(3,t_MAT);
 1201   if (lgcols(f) == 1)
 1202   {
 1203     gel(h,1) = mkcolcopy(x);
 1204     gel(h,2) = mkcol(gen_m1);
 1205   }
 1206   else
 1207   {
 1208     gel(h,1) = append(gel(f,1), x);
 1209     gel(h,2) = gconcat(gel(f,2), gen_m1);
 1210   }
 1211   return h;
 1212 }
 1213 
 1214 GEN
 1215 famat_mul(GEN f, GEN g)
 1216 {
 1217   GEN h;
 1218   if (typ(g) != t_MAT) {
 1219     if (typ(f) == t_MAT) return famat_add(f, g);
 1220     h = cgetg(3, t_MAT);
 1221     gel(h,1) = mkcol2(gcopy(f), gcopy(g));
 1222     gel(h,2) = mkcol2(gen_1, gen_1);
 1223   }
 1224   if (typ(f) != t_MAT) return famat_add(g, f);
 1225   if (lgcols(f) == 1) return gcopy(g);
 1226   if (lgcols(g) == 1) return gcopy(f);
 1227   h = cgetg(3,t_MAT);
 1228   gel(h,1) = gconcat(gel(f,1), gel(g,1));
 1229   gel(h,2) = gconcat(gel(f,2), gel(g,2));
 1230   return h;
 1231 }
 1232 
 1233 GEN
 1234 famat_div(GEN f, GEN g)
 1235 {
 1236   GEN h;
 1237   if (typ(g) != t_MAT) {
 1238     if (typ(f) == t_MAT) return famat_sub(f, g);
 1239     h = cgetg(3, t_MAT);
 1240     gel(h,1) = mkcol2(gcopy(f), gcopy(g));
 1241     gel(h,2) = mkcol2(gen_1, gen_m1);
 1242   }
 1243   if (typ(f) != t_MAT) return famat_sub(g, f);
 1244   if (lgcols(f) == 1) return famat_inv(g);
 1245   if (lgcols(g) == 1) return gcopy(f);
 1246   h = cgetg(3,t_MAT);
 1247   gel(h,1) = gconcat(gel(f,1), gel(g,1));
 1248   gel(h,2) = gconcat(gel(f,2), gneg(gel(g,2)));
 1249   return h;
 1250 }
 1251 
 1252 GEN
 1253 famat_sqr(GEN f)
 1254 {
 1255   GEN h;
 1256   if (typ(f) != t_MAT) return to_famat(f,gen_2);
 1257   if (lgcols(f) == 1) return gcopy(f);
 1258   h = cgetg(3,t_MAT);
 1259   gel(h,1) = gcopy(gel(f,1));
 1260   gel(h,2) = gmul2n(gel(f,2),1);
 1261   return h;
 1262 }
 1263 
 1264 GEN
 1265 famat_inv_shallow(GEN f)
 1266 {
 1267   if (typ(f) != t_MAT) return to_famat_shallow(f,gen_m1);
 1268   if (lgcols(f) == 1) return f;
 1269   return mkmat2(gel(f,1), ZC_neg(gel(f,2)));
 1270 }
 1271 GEN
 1272 famat_inv(GEN f)
 1273 {
 1274   if (typ(f) != t_MAT) return to_famat(f,gen_m1);
 1275   if (lgcols(f) == 1) return gcopy(f);
 1276   retmkmat2(gcopy(gel(f,1)), ZC_neg(gel(f,2)));
 1277 }
 1278 GEN
 1279 famat_pow(GEN f, GEN n)
 1280 {
 1281   if (typ(f) != t_MAT) return to_famat(f,n);
 1282   if (lgcols(f) == 1) return gcopy(f);
 1283   retmkmat2(gcopy(gel(f,1)), ZC_Z_mul(gel(f,2),n));
 1284 }
 1285 GEN
 1286 famat_pow_shallow(GEN f, GEN n)
 1287 {
 1288   if (is_pm1(n)) return signe(n) > 0? f: famat_inv_shallow(f);
 1289   if (typ(f) != t_MAT) return to_famat_shallow(f,n);
 1290   if (lgcols(f) == 1) return f;
 1291   return mkmat2(gel(f,1), ZC_Z_mul(gel(f,2),n));
 1292 }
 1293 
 1294 GEN
 1295 famat_pows_shallow(GEN f, long n)
 1296 {
 1297   if (n==1) return f;
 1298   if (n==-1) return famat_inv_shallow(f);
 1299   if (typ(f) != t_MAT) return to_famat_shallow(f, stoi(n));
 1300   if (lgcols(f) == 1) return f;
 1301   return mkmat2(gel(f,1), ZC_z_mul(gel(f,2),n));
 1302 }
 1303 
 1304 GEN
 1305 famat_Z_gcd(GEN M, GEN n)
 1306 {
 1307   pari_sp av=avma;
 1308   long i, j, l=lgcols(M);
 1309   GEN F=cgetg(3,t_MAT);
 1310   gel(F,1)=cgetg(l,t_COL);
 1311   gel(F,2)=cgetg(l,t_COL);
 1312   for (i=1, j=1; i<l; i++)
 1313   {
 1314     GEN p = gcoeff(M,i,1);
 1315     GEN e = gminsg(Z_pval(n,p),gcoeff(M,i,2));
 1316     if (signe(e))
 1317     {
 1318       gcoeff(F,j,1)=p;
 1319       gcoeff(F,j,2)=e;
 1320       j++;
 1321     }
 1322   }
 1323   setlg(gel(F,1),j); setlg(gel(F,2),j);
 1324   return gerepilecopy(av,F);
 1325 }
 1326 
 1327 /* x assumed to be a t_MATs (factorization matrix), or compatible with
 1328  * the element_* functions. */
 1329 static GEN
 1330 ext_sqr(GEN nf, GEN x)
 1331 { return (typ(x)==t_MAT)? famat_sqr(x): nfsqr(nf, x); }
 1332 static GEN
 1333 ext_mul(GEN nf, GEN x, GEN y)
 1334 { return (typ(x)==t_MAT)? famat_mul(x,y): nfmul(nf, x, y); }
 1335 static GEN
 1336 ext_inv(GEN nf, GEN x)
 1337 { return (typ(x)==t_MAT)? famat_inv(x): nfinv(nf, x); }
 1338 static GEN
 1339 ext_pow(GEN nf, GEN x, GEN n)
 1340 { return (typ(x)==t_MAT)? famat_pow(x,n): nfpow(nf, x, n); }
 1341 
 1342 GEN
 1343 famat_to_nf(GEN nf, GEN f)
 1344 {
 1345   GEN t, x, e;
 1346   long i;
 1347   if (lgcols(f) == 1) return gen_1;
 1348   x = gel(f,1);
 1349   e = gel(f,2);
 1350   t = nfpow(nf, gel(x,1), gel(e,1));
 1351   for (i=lg(x)-1; i>1; i--)
 1352     t = nfmul(nf, t, nfpow(nf, gel(x,i), gel(e,i)));
 1353   return t;
 1354 }
 1355 
 1356 GEN
 1357 famat_idealfactor(GEN nf, GEN x)
 1358 {
 1359   long i, l;
 1360   GEN g = gel(x,1), e = gel(x,2), h = cgetg_copy(g, &l);
 1361   for (i = 1; i < l; i++) gel(h,i) = idealfactor(nf, gel(g,i));
 1362   h = famat_reduce(famatV_factorback(h,e));
 1363   return sort_factor(h, (void*)&cmp_prime_ideal, &cmp_nodata);
 1364 }
 1365 
 1366 GEN
 1367 famat_reduce(GEN fa)
 1368 {
 1369   GEN E, G, L, g, e;
 1370   long i, k, l;
 1371 
 1372   if (typ(fa) != t_MAT || lgcols(fa) == 1) return fa;
 1373   g = gel(fa,1); l = lg(g);
 1374   e = gel(fa,2);
 1375   L = gen_indexsort(g, (void*)&cmp_universal, &cmp_nodata);
 1376   G = cgetg(l, t_COL);
 1377   E = cgetg(l, t_COL);
 1378   /* merge */
 1379   for (k=i=1; i<l; i++,k++)
 1380   {
 1381     gel(G,k) = gel(g,L[i]);
 1382     gel(E,k) = gel(e,L[i]);
 1383     if (k > 1 && gidentical(gel(G,k), gel(G,k-1)))
 1384     {
 1385       gel(E,k-1) = addii(gel(E,k), gel(E,k-1));
 1386       k--;
 1387     }
 1388   }
 1389   /* kill 0 exponents */
 1390   l = k;
 1391   for (k=i=1; i<l; i++)
 1392     if (!gequal0(gel(E,i)))
 1393     {
 1394       gel(G,k) = gel(G,i);
 1395       gel(E,k) = gel(E,i); k++;
 1396     }
 1397   setlg(G, k);
 1398   setlg(E, k); return mkmat2(G,E);
 1399 }
 1400 GEN
 1401 matreduce(GEN f)
 1402 { pari_sp av = avma;
 1403   switch(typ(f))
 1404   {
 1405     case t_VEC: case t_COL:
 1406     {
 1407       GEN e; f = vec_reduce(f, &e);
 1408       return gerepilecopy(av, mkmat2(f, zc_to_ZC(e)));
 1409     }
 1410     case t_MAT:
 1411       if (lg(f) == 3) break;
 1412     default:
 1413       pari_err_TYPE("matreduce", f);
 1414   }
 1415   if (typ(gel(f,1)) == t_VECSMALL)
 1416     f = famatsmall_reduce(f);
 1417   else
 1418   {
 1419     if (!RgV_is_ZV(gel(f,2))) pari_err_TYPE("matreduce",f);
 1420     f = famat_reduce(f);
 1421   }
 1422   return gerepilecopy(av, f);
 1423 }
 1424 
 1425 GEN
 1426 famatsmall_reduce(GEN fa)
 1427 {
 1428   GEN E, G, L, g, e;
 1429   long i, k, l;
 1430   if (lgcols(fa) == 1) return fa;
 1431   g = gel(fa,1); l = lg(g);
 1432   e = gel(fa,2);
 1433   L = vecsmall_indexsort(g);
 1434   G = cgetg(l, t_VECSMALL);
 1435   E = cgetg(l, t_VECSMALL);
 1436   /* merge */
 1437   for (k=i=1; i<l; i++,k++)
 1438   {
 1439     G[k] = g[L[i]];
 1440     E[k] = e[L[i]];
 1441     if (k > 1 && G[k] == G[k-1])
 1442     {
 1443       E[k-1] += E[k];
 1444       k--;
 1445     }
 1446   }
 1447   /* kill 0 exponents */
 1448   l = k;
 1449   for (k=i=1; i<l; i++)
 1450     if (E[i])
 1451     {
 1452       G[k] = G[i];
 1453       E[k] = E[i]; k++;
 1454     }
 1455   setlg(G, k);
 1456   setlg(E, k); return mkmat2(G,E);
 1457 }
 1458 
 1459 GEN
 1460 famat_remove_trivial(GEN fa)
 1461 {
 1462   GEN P, E, p = gel(fa,1), e = gel(fa,2);
 1463   long j, k, l = lg(p);
 1464   P = cgetg(l, t_COL);
 1465   E = cgetg(l, t_COL);
 1466   for (j = k = 1; j < l; j++)
 1467     if (signe(gel(e,j))) { gel(P,k) = gel(p,j); gel(E,k++) = gel(e,j); }
 1468   setlg(P, k); setlg(E, k); return mkmat2(P,E);
 1469 }
 1470 
 1471 GEN
 1472 famatV_factorback(GEN v, GEN e)
 1473 {
 1474   long i, l = lg(e);
 1475   GEN V;
 1476   if (l == 1) return trivial_fact();
 1477   V = signe(gel(e,1))? famat_pow_shallow(gel(v,1), gel(e,1)): trivial_fact();
 1478   for (i = 2; i < l; i++) V = famat_mulpow_shallow(V, gel(v,i), gel(e,i));
 1479   return V;
 1480 }
 1481 
 1482 GEN
 1483 famatV_zv_factorback(GEN v, GEN e)
 1484 {
 1485   long i, l = lg(e);
 1486   GEN V;
 1487   if (l == 1) return trivial_fact();
 1488   V = uel(e,1)? famat_pows_shallow(gel(v,1), uel(e,1)): trivial_fact();
 1489   for (i = 2; i < l; i++) V = famat_mulpows_shallow(V, gel(v,i), uel(e,i));
 1490   return V;
 1491 }
 1492 
 1493 GEN
 1494 ZM_famat_limit(GEN fa, GEN limit)
 1495 {
 1496   pari_sp av;
 1497   GEN E, G, g, e, r;
 1498   long i, k, l, n, lG;
 1499 
 1500   if (lgcols(fa) == 1) return fa;
 1501   g = gel(fa,1); l = lg(g);
 1502   e = gel(fa,2);
 1503   for(n=0, i=1; i<l; i++)
 1504     if (cmpii(gel(g,i),limit)<=0) n++;
 1505   lG = n<l-1 ? n+2 : n+1;
 1506   G = cgetg(lG, t_COL);
 1507   E = cgetg(lG, t_COL);
 1508   av = avma;
 1509   for (i=1, k=1, r = gen_1; i<l; i++)
 1510   {
 1511     if (cmpii(gel(g,i),limit)<=0)
 1512     {
 1513       gel(G,k) = gel(g,i);
 1514       gel(E,k) = gel(e,i);
 1515       k++;
 1516     } else r = mulii(r, powii(gel(g,i), gel(e,i)));
 1517   }
 1518   if (k<i)
 1519   {
 1520     gel(G, k) = gerepileuptoint(av, r);
 1521     gel(E, k) = gen_1;
 1522   }
 1523   return mkmat2(G,E);
 1524 }
 1525 
 1526 /* assume pr has degree 1 and coprime to Q_denom(x) */
 1527 static GEN
 1528 to_Fp_coprime(GEN nf, GEN x, GEN modpr)
 1529 {
 1530   GEN d, r, p = modpr_get_p(modpr);
 1531   x = nf_to_scalar_or_basis(nf,x);
 1532   if (typ(x) != t_COL) return Rg_to_Fp(x,p);
 1533   x = Q_remove_denom(x, &d);
 1534   r = zk_to_Fq(x, modpr);
 1535   if (d) r = Fp_div(r, d, p);
 1536   return r;
 1537 }
 1538 
 1539 /* pr coprime to all denominators occurring in x */
 1540 static GEN
 1541 famat_to_Fp_coprime(GEN nf, GEN x, GEN modpr)
 1542 {
 1543   GEN p = modpr_get_p(modpr);
 1544   GEN t = NULL, g = gel(x,1), e = gel(x,2), q = subiu(p,1);
 1545   long i, l = lg(g);
 1546   for (i = 1; i < l; i++)
 1547   {
 1548     GEN n = modii(gel(e,i), q);
 1549     if (signe(n))
 1550     {
 1551       GEN h = to_Fp_coprime(nf, gel(g,i), modpr);
 1552       h = Fp_pow(h, n, p);
 1553       t = t? Fp_mul(t, h, p): h;
 1554     }
 1555   }
 1556   return t? modii(t, p): gen_1;
 1557 }
 1558 
 1559 /* cf famat_to_nf_modideal_coprime, modpr attached to prime of degree 1 */
 1560 GEN
 1561 nf_to_Fp_coprime(GEN nf, GEN x, GEN modpr)
 1562 {
 1563   return typ(x)==t_MAT? famat_to_Fp_coprime(nf, x, modpr)
 1564                       : to_Fp_coprime(nf, x, modpr);
 1565 }
 1566 
 1567 static long
 1568 zk_pvalrem(GEN x, GEN p, GEN *py)
 1569 { return (typ(x) == t_INT)? Z_pvalrem(x, p, py): ZV_pvalrem(x, p, py); }
 1570 /* x a QC or Q. Return a ZC or Z, whose content is coprime to Z. Set v, dx
 1571  * such that x = p^v (newx / dx); dx = NULL if 1 */
 1572 static GEN
 1573 nf_remove_denom_p(GEN nf, GEN x, GEN p, GEN *pdx, long *pv)
 1574 {
 1575   long vcx;
 1576   GEN dx;
 1577   x = nf_to_scalar_or_basis(nf, x);
 1578   x = Q_remove_denom(x, &dx);
 1579   if (dx)
 1580   {
 1581     vcx = - Z_pvalrem(dx, p, &dx);
 1582     if (!vcx) vcx = zk_pvalrem(x, p, &x);
 1583     if (isint1(dx)) dx = NULL;
 1584   }
 1585   else
 1586   {
 1587     vcx = zk_pvalrem(x, p, &x);
 1588     dx = NULL;
 1589   }
 1590   *pv = vcx;
 1591   *pdx = dx; return x;
 1592 }
 1593 /* x = b^e/p^(e-1) in Z_K; x = 0 mod p/pr^e, (x,pr) = 1. Return NULL
 1594  * if p inert (instead of 1) */
 1595 static GEN
 1596 p_makecoprime(GEN pr)
 1597 {
 1598   GEN B = pr_get_tau(pr), b;
 1599   long i, e;
 1600 
 1601   if (typ(B) == t_INT) return NULL;
 1602   b = gel(B,1); /* B = multiplication table by b */
 1603   e = pr_get_e(pr);
 1604   if (e == 1) return b;
 1605   /* one could also divide (exactly) by p in each iteration */
 1606   for (i = 1; i < e; i++) b = ZM_ZC_mul(B, b);
 1607   return ZC_Z_divexact(b, powiu(pr_get_p(pr), e-1));
 1608 }
 1609 
 1610 /* Compute A = prod g[i]^e[i] mod pr^k, assuming (A, pr) = 1.
 1611  * Method: modify each g[i] so that it becomes coprime to pr,
 1612  * g[i] *= (b/p)^v_pr(g[i]), where b/p = pr^(-1) times something integral
 1613  * and prime to p; globally, we multiply by (b/p)^v_pr(A) = 1.
 1614  * Optimizations:
 1615  * 1) remove all powers of p from contents, and consider extra generator p^vp;
 1616  * modified as p * (b/p)^e = b^e / p^(e-1)
 1617  * 2) remove denominators, coprime to p, by multiplying by inverse mod prk\cap Z
 1618  *
 1619  * EX = multiple of exponent of (O_K / pr^k)^* used to reduce the product in
 1620  * case the e[i] are large */
 1621 GEN
 1622 famat_makecoprime(GEN nf, GEN g, GEN e, GEN pr, GEN prk, GEN EX)
 1623 {
 1624   GEN G, E, t, vp = NULL, p = pr_get_p(pr), prkZ = gcoeff(prk, 1,1);
 1625   long i, l = lg(g);
 1626 
 1627   G = cgetg(l+1, t_VEC);
 1628   E = cgetg(l+1, t_VEC); /* l+1: room for "modified p" */
 1629   for (i=1; i < l; i++)
 1630   {
 1631     long vcx;
 1632     GEN dx, x = nf_remove_denom_p(nf, gel(g,i), p, &dx, &vcx);
 1633     if (vcx) /* = v_p(content(g[i])) */
 1634     {
 1635       GEN a = mulsi(vcx, gel(e,i));
 1636       vp = vp? addii(vp, a): a;
 1637     }
 1638     /* x integral, content coprime to p; dx coprime to p */
 1639     if (typ(x) == t_INT)
 1640     { /* x coprime to p, hence to pr */
 1641       x = modii(x, prkZ);
 1642       if (dx) x = Fp_div(x, dx, prkZ);
 1643     }
 1644     else
 1645     {
 1646       (void)ZC_nfvalrem(x, pr, &x); /* x *= (b/p)^v_pr(x) */
 1647       x = ZC_hnfrem(FpC_red(x,prkZ), prk);
 1648       if (dx) x = FpC_Fp_mul(x, Fp_inv(dx,prkZ), prkZ);
 1649     }
 1650     gel(G,i) = x;
 1651     gel(E,i) = gel(e,i);
 1652   }
 1653 
 1654   t = vp? p_makecoprime(pr): NULL;
 1655   if (!t)
 1656   { /* no need for extra generator */
 1657     setlg(G,l);
 1658     setlg(E,l);
 1659   }
 1660   else
 1661   {
 1662     gel(G,i) = FpC_red(t, prkZ);
 1663     gel(E,i) = vp;
 1664   }
 1665   return famat_to_nf_modideal_coprime(nf, G, E, prk, EX);
 1666 }
 1667 
 1668 /* simplified version of famat_makecoprime for X = SUnits[1] */
 1669 GEN
 1670 sunits_makecoprime(GEN X, GEN pr, GEN prk)
 1671 {
 1672   GEN G, p = pr_get_p(pr), prkZ = gcoeff(prk,1,1);
 1673   long i, l = lg(X);
 1674 
 1675   G = cgetg(l, t_VEC);
 1676   for (i = 1; i < l; i++)
 1677   {
 1678     GEN x = gel(X,i);
 1679     if (typ(x) == t_INT) /* a prime */
 1680       x = equalii(x,p)? p_makecoprime(pr): modii(x, prkZ);
 1681     else
 1682     {
 1683       (void)ZC_nfvalrem(x, pr, &x); /* x *= (b/p)^v_pr(x) */
 1684       x = ZC_hnfrem(FpC_red(x,prkZ), prk);
 1685     }
 1686     gel(G,i) = x;
 1687   }
 1688   return G;
 1689 }
 1690 
 1691 /* prod g[i]^e[i] mod bid, assume (g[i], id) = 1 and 1 < lg(g) <= lg(e) */
 1692 GEN
 1693 famat_to_nf_moddivisor(GEN nf, GEN g, GEN e, GEN bid)
 1694 {
 1695   GEN t, cyc = bid_get_cyc(bid);
 1696   if (lg(cyc) == 1)
 1697     t = gen_1;
 1698   else
 1699     t = famat_to_nf_modideal_coprime(nf, g, e, bid_get_ideal(bid),
 1700                                      cyc_get_expo(cyc));
 1701   return set_sign_mod_divisor(nf, mkmat2(g,e), t, bid_get_sarch(bid));
 1702 }
 1703 
 1704 GEN
 1705 vecmul(GEN x, GEN y)
 1706 {
 1707   if (is_scalar_t(typ(x))) return gmul(x,y);
 1708   pari_APPLY_same(vecmul(gel(x,i), gel(y,i)))
 1709 }
 1710 
 1711 GEN
 1712 vecinv(GEN x)
 1713 {
 1714   if (is_scalar_t(typ(x))) return ginv(x);
 1715   pari_APPLY_same(vecinv(gel(x,i)))
 1716 }
 1717 
 1718 GEN
 1719 vecpow(GEN x, GEN n)
 1720 {
 1721   if (is_scalar_t(typ(x))) return powgi(x,n);
 1722   pari_APPLY_same(vecpow(gel(x,i), n))
 1723 }
 1724 
 1725 GEN
 1726 vecdiv(GEN x, GEN y)
 1727 {
 1728   if (is_scalar_t(typ(x))) return gdiv(x,y);
 1729   pari_APPLY_same(vecdiv(gel(x,i), gel(y,i)))
 1730 }
 1731 
 1732 /* A ideal as a square t_MAT */
 1733 static GEN
 1734 idealmulelt(GEN nf, GEN x, GEN A)
 1735 {
 1736   long i, lx;
 1737   GEN dx, dA, D;
 1738   if (lg(A) == 1) return cgetg(1, t_MAT);
 1739   x = nf_to_scalar_or_basis(nf,x);
 1740   if (typ(x) != t_COL)
 1741     return isintzero(x)? cgetg(1,t_MAT): RgM_Rg_mul(A, Q_abs_shallow(x));
 1742   x = Q_remove_denom(x, &dx);
 1743   A = Q_remove_denom(A, &dA);
 1744   x = zk_multable(nf, x);
 1745   D = mulii(zkmultable_capZ(x), gcoeff(A,1,1));
 1746   x = zkC_multable_mul(A, x);
 1747   settyp(x, t_MAT); lx = lg(x);
 1748   /* x may contain scalars (at most 1 since the ideal is nonzero)*/
 1749   for (i=1; i<lx; i++)
 1750     if (typ(gel(x,i)) == t_INT)
 1751     {
 1752       if (i > 1) swap(gel(x,1), gel(x,i)); /* help HNF */
 1753       gel(x,1) = scalarcol_shallow(gel(x,1), lx-1);
 1754       break;
 1755     }
 1756   x = ZM_hnfmodid(x, D);
 1757   dx = mul_denom(dx,dA);
 1758   return dx? gdiv(x,dx): x;
 1759 }
 1760 
 1761 /* nf a true nf, tx <= ty */
 1762 static GEN
 1763 idealmul_aux(GEN nf, GEN x, GEN y, long tx, long ty)
 1764 {
 1765   GEN z, cx, cy;
 1766   switch(tx)
 1767   {
 1768     case id_PRINCIPAL:
 1769       switch(ty)
 1770       {
 1771         case id_PRINCIPAL:
 1772           return idealhnf_principal(nf, nfmul(nf,x,y));
 1773         case id_PRIME:
 1774         {
 1775           GEN p = pr_get_p(y), pi = pr_get_gen(y), cx;
 1776           if (pr_is_inert(y)) return RgM_Rg_mul(idealhnf_principal(nf,x),p);
 1777 
 1778           x = nf_to_scalar_or_basis(nf, x);
 1779           switch(typ(x))
 1780           {
 1781             case t_INT:
 1782               if (!signe(x)) return cgetg(1,t_MAT);
 1783               return ZM_Z_mul(pr_hnf(nf,y), absi_shallow(x));
 1784             case t_FRAC:
 1785               return RgM_Rg_mul(pr_hnf(nf,y), Q_abs_shallow(x));
 1786           }
 1787           /* t_COL */
 1788           x = Q_primitive_part(x, &cx);
 1789           x = zk_multable(nf, x);
 1790           z = shallowconcat(ZM_Z_mul(x,p), ZM_ZC_mul(x,pi));
 1791           z = ZM_hnfmodid(z, mulii(p, zkmultable_capZ(x)));
 1792           return cx? ZM_Q_mul(z, cx): z;
 1793         }
 1794         default: /* id_MAT */
 1795           return idealmulelt(nf, x,y);
 1796       }
 1797     case id_PRIME:
 1798       if (ty==id_PRIME)
 1799       { y = pr_hnf(nf,y); cy = NULL; }
 1800       else
 1801         y = Q_primitive_part(y, &cy);
 1802       y = idealHNF_mul_two(nf,y,x);
 1803       return cy? ZM_Q_mul(y,cy): y;
 1804 
 1805     default: /* id_MAT */
 1806     {
 1807       long N = nf_get_degree(nf);
 1808       if (lg(x)-1 != N || lg(y)-1 != N) pari_err_DIM("idealmul");
 1809       x = Q_primitive_part(x, &cx);
 1810       y = Q_primitive_part(y, &cy); cx = mul_content(cx,cy);
 1811       y = idealHNF_mul(nf,x,y);
 1812       return cx? ZM_Q_mul(y,cx): y;
 1813     }
 1814   }
 1815 }
 1816 
 1817 /* output the ideal product ix.iy */
 1818 GEN
 1819 idealmul(GEN nf, GEN x, GEN y)
 1820 {
 1821   pari_sp av;
 1822   GEN res, ax, ay, z;
 1823   long tx = idealtyp(&x,&ax);
 1824   long ty = idealtyp(&y,&ay), f;
 1825   if (tx>ty) { swap(ax,ay); swap(x,y); lswap(tx,ty); }
 1826   f = (ax||ay); res = f? cgetg(3,t_VEC): NULL; /*product is an extended ideal*/
 1827   av = avma;
 1828   z = gerepileupto(av, idealmul_aux(checknf(nf), x,y, tx,ty));
 1829   if (!f) return z;
 1830   if (ax && ay)
 1831     ax = ext_mul(nf, ax, ay);
 1832   else
 1833     ax = gcopy(ax? ax: ay);
 1834   gel(res,1) = z; gel(res,2) = ax; return res;
 1835 }
 1836 
 1837 /* Return x, integral in 2-elt form, such that pr^2 = c * x. cf idealpowprime
 1838  * nf = true nf */
 1839 static GEN
 1840 idealsqrprime(GEN nf, GEN pr, GEN *pc)
 1841 {
 1842   GEN p = pr_get_p(pr), q, gen;
 1843   long e = pr_get_e(pr), f = pr_get_f(pr);
 1844 
 1845   q = (e == 1)? sqri(p): p;
 1846   if (e <= 2 && e * f == nf_get_degree(nf))
 1847   { /* pr^e = (p) */
 1848     *pc = q;
 1849     return mkvec2(gen_1,gen_0);
 1850   }
 1851   gen = nfsqr(nf, pr_get_gen(pr));
 1852   gen = FpC_red(gen, q);
 1853   *pc = NULL;
 1854   return mkvec2(q, gen);
 1855 }
 1856 /* cf idealpow_aux */
 1857 static GEN
 1858 idealsqr_aux(GEN nf, GEN x, long tx)
 1859 {
 1860   GEN T = nf_get_pol(nf), m, cx, a, alpha;
 1861   long N = degpol(T);
 1862   switch(tx)
 1863   {
 1864     case id_PRINCIPAL:
 1865       return idealhnf_principal(nf, nfsqr(nf,x));
 1866     case id_PRIME:
 1867       if (pr_is_inert(x)) return scalarmat(sqri(gel(x,1)), N);
 1868       x = idealsqrprime(nf, x, &cx);
 1869       x = idealhnf_two(nf,x);
 1870       return cx? ZM_Z_mul(x, cx): x;
 1871     default:
 1872       x = Q_primitive_part(x, &cx);
 1873       a = mat_ideal_two_elt(nf,x); alpha = gel(a,2); a = gel(a,1);
 1874       alpha = nfsqr(nf,alpha);
 1875       m = zk_scalar_or_multable(nf, alpha);
 1876       if (typ(m) == t_INT) {
 1877         x = gcdii(sqri(a), m);
 1878         if (cx) x = gmul(x, gsqr(cx));
 1879         x = scalarmat(x, N);
 1880       }
 1881       else
 1882       { /* could use gcdii(sqri(a), zkmultable_capZ(m)), but costly */
 1883         x = ZM_hnfmodid(m, sqri(a));
 1884         if (cx) cx = gsqr(cx);
 1885         if (cx) x = ZM_Q_mul(x, cx);
 1886       }
 1887       return x;
 1888   }
 1889 }
 1890 GEN
 1891 idealsqr(GEN nf, GEN x)
 1892 {
 1893   pari_sp av;
 1894   GEN res, ax, z;
 1895   long tx = idealtyp(&x,&ax);
 1896   res = ax? cgetg(3,t_VEC): NULL; /*product is an extended ideal*/
 1897   av = avma;
 1898   z = gerepileupto(av, idealsqr_aux(checknf(nf), x, tx));
 1899   if (!ax) return z;
 1900   gel(res,1) = z;
 1901   gel(res,2) = ext_sqr(nf, ax); return res;
 1902 }
 1903 
 1904 /* norm of an ideal */
 1905 GEN
 1906 idealnorm(GEN nf, GEN x)
 1907 {
 1908   pari_sp av;
 1909   GEN y;
 1910   long tx;
 1911 
 1912   switch(idealtyp(&x,&y))
 1913   {
 1914     case id_PRIME: return pr_norm(x);
 1915     case id_MAT: return RgM_det_triangular(x);
 1916   }
 1917   /* id_PRINCIPAL */
 1918   nf = checknf(nf); av = avma;
 1919   x = nfnorm(nf, x);
 1920   tx = typ(x);
 1921   if (tx == t_INT) return gerepileuptoint(av, absi(x));
 1922   if (tx != t_FRAC) pari_err_TYPE("idealnorm",x);
 1923   return gerepileupto(av, Q_abs(x));
 1924 }
 1925 
 1926 /* x \cap Z */
 1927 GEN
 1928 idealdown(GEN nf, GEN x)
 1929 {
 1930   pari_sp av = avma;
 1931   GEN y, c;
 1932   switch(idealtyp(&x,&y))
 1933   {
 1934     case id_PRIME: return icopy(pr_get_p(x));
 1935     case id_MAT: return gcopy(gcoeff(x,1,1));
 1936   }
 1937   /* id_PRINCIPAL */
 1938   nf = checknf(nf); av = avma;
 1939   x = nf_to_scalar_or_basis(nf, x);
 1940   if (is_rational_t(typ(x))) return Q_abs(x);
 1941   x = Q_primitive_part(x, &c);
 1942   y = zkmultable_capZ(zk_multable(nf, x));
 1943   return gerepilecopy(av, mul_content(c, y));
 1944 }
 1945 
 1946 /* true nf */
 1947 static GEN
 1948 idealismaximal_int(GEN nf, GEN p)
 1949 {
 1950   GEN L;
 1951   if (!BPSW_psp(p)) return NULL;
 1952   if (!dvdii(nf_get_index(nf), p) &&
 1953       !FpX_is_irred(FpX_red(nf_get_pol(nf),p), p)) return NULL;
 1954   L = idealprimedec(nf, p);
 1955   return lg(L) == 2? gel(L,1): NULL;
 1956 }
 1957 /* true nf */
 1958 static GEN
 1959 idealismaximal_mat(GEN nf, GEN x)
 1960 {
 1961   GEN p, c, L;
 1962   long i, l, f;
 1963   x = Q_primitive_part(x, &c);
 1964   p = gcoeff(x,1,1);
 1965   if (c)
 1966   {
 1967     if (typ(c) == t_FRAC || !equali1(p)) return NULL;
 1968     return idealismaximal_int(nf, p);
 1969   }
 1970   if (!BPSW_psp(p)) return NULL;
 1971   l = lg(x); f = 1;
 1972   for (i = 2; i < l; i++)
 1973   {
 1974     c = gcoeff(x,i,i);
 1975     if (equalii(c, p)) f++; else if (!equali1(c)) return NULL;
 1976   }
 1977   L = idealprimedec_limit_f(nf, p, f);
 1978   for (i = lg(L)-1; i; i--)
 1979   {
 1980     GEN pr = gel(L,i);
 1981     if (pr_get_f(pr) != f) break;
 1982     if (idealval(nf, x, pr) == 1) return pr;
 1983   }
 1984   return NULL;
 1985 }
 1986 /* true nf */
 1987 static GEN
 1988 idealismaximal_i(GEN nf, GEN x)
 1989 {
 1990   GEN L, p, pr, c;
 1991   long i, l;
 1992   switch(idealtyp(&x,&c))
 1993   {
 1994     case id_PRIME: return x;
 1995     case id_MAT: return idealismaximal_mat(nf, x);
 1996   }
 1997   /* id_PRINCIPAL */
 1998   nf = checknf(nf);
 1999   x = nf_to_scalar_or_basis(nf, x);
 2000   switch(typ(x))
 2001   {
 2002     case t_INT: return idealismaximal_int(nf, absi_shallow(x));
 2003     case t_FRAC: return NULL;
 2004   }
 2005   x = Q_primitive_part(x, &c);
 2006   if (c) return NULL;
 2007   p = zkmultable_capZ(zk_multable(nf, x));
 2008   L = idealprimedec(nf, p); l = lg(L); pr = NULL;
 2009   for (i = 1; i < l; i++)
 2010   {
 2011     long v = ZC_nfval(x, gel(L,i));
 2012     if (v > 1 || (v && pr)) return NULL;
 2013     pr = gel(L,i);
 2014   }
 2015   return pr;
 2016 }
 2017 GEN
 2018 idealismaximal(GEN nf, GEN x)
 2019 {
 2020   pari_sp av = avma;
 2021   x = idealismaximal_i(checknf(nf), x);
 2022   if (!x) { set_avma(av); return gen_0; }
 2023   return gerepilecopy(av, x);
 2024 }
 2025 
 2026 /* I^(-1) = { x \in K, Tr(x D^(-1) I) \in Z }, D different of K/Q
 2027  *
 2028  * nf[5][6] = pp( D^(-1) ) = pp( HNF( T^(-1) ) ), T = (Tr(wi wj))
 2029  * nf[5][7] = same in 2-elt form.
 2030  * Assume I integral. Return the integral ideal (I\cap Z) I^(-1) */
 2031 GEN
 2032 idealHNF_inv_Z(GEN nf, GEN I)
 2033 {
 2034   GEN J, dual, IZ = gcoeff(I,1,1); /* I \cap Z */
 2035   if (isint1(IZ)) return matid(lg(I)-1);
 2036   J = idealHNF_mul(nf,I, gmael(nf,5,7));
 2037  /* I in HNF, hence easily inverted; multiply by IZ to get integer coeffs
 2038   * missing content cancels while solving the linear equation */
 2039   dual = shallowtrans( hnf_divscale(J, gmael(nf,5,6), IZ) );
 2040   return ZM_hnfmodid(dual, IZ);
 2041 }
 2042 /* I HNF with rational coefficients (denominator d). */
 2043 GEN
 2044 idealHNF_inv(GEN nf, GEN I)
 2045 {
 2046   GEN J, IQ = gcoeff(I,1,1); /* I \cap Q; d IQ = dI \cap Z */
 2047   J = idealHNF_inv_Z(nf, Q_remove_denom(I, NULL)); /* = (dI)^(-1) * (d IQ) */
 2048   return equali1(IQ)? J: RgM_Rg_div(J, IQ);
 2049 }
 2050 
 2051 /* return p * P^(-1)  [integral] */
 2052 GEN
 2053 pr_inv_p(GEN pr)
 2054 {
 2055   if (pr_is_inert(pr)) return matid(pr_get_f(pr));
 2056   return ZM_hnfmodid(pr_get_tau(pr), pr_get_p(pr));
 2057 }
 2058 GEN
 2059 pr_inv(GEN pr)
 2060 {
 2061   GEN p = pr_get_p(pr);
 2062   if (pr_is_inert(pr)) return scalarmat(ginv(p), pr_get_f(pr));
 2063   return RgM_Rg_div(ZM_hnfmodid(pr_get_tau(pr),p), p);
 2064 }
 2065 
 2066 GEN
 2067 idealinv(GEN nf, GEN x)
 2068 {
 2069   GEN res, ax;
 2070   pari_sp av;
 2071   long tx = idealtyp(&x,&ax), N;
 2072 
 2073   res = ax? cgetg(3,t_VEC): NULL;
 2074   nf = checknf(nf); av = avma;
 2075   N = nf_get_degree(nf);
 2076   switch (tx)
 2077   {
 2078     case id_MAT:
 2079       if (lg(x)-1 != N) pari_err_DIM("idealinv");
 2080       x = idealHNF_inv(nf,x); break;
 2081     case id_PRINCIPAL:
 2082       x = nf_to_scalar_or_basis(nf, x);
 2083       if (typ(x) != t_COL)
 2084         x = idealhnf_principal(nf,ginv(x));
 2085       else
 2086       { /* nfinv + idealhnf where we already know (x) \cap Z */
 2087         GEN c, d;
 2088         x = Q_remove_denom(x, &c);
 2089         x = zk_inv(nf, x);
 2090         x = Q_remove_denom(x, &d); /* true inverse is c/d * x */
 2091         if (!d) /* x and x^(-1) integral => x a unit */
 2092           x = c? scalarmat(c, N): matid(N);
 2093         else
 2094         {
 2095           c = c? gdiv(c,d): ginv(d);
 2096           x = zk_multable(nf, x);
 2097           x = ZM_Q_mul(ZM_hnfmodid(x,d), c);
 2098         }
 2099       }
 2100       break;
 2101     case id_PRIME:
 2102       x = pr_inv(x); break;
 2103   }
 2104   x = gerepileupto(av,x); if (!ax) return x;
 2105   gel(res,1) = x;
 2106   gel(res,2) = ext_inv(nf, ax); return res;
 2107 }
 2108 
 2109 /* write x = A/B, A,B coprime integral ideals */
 2110 GEN
 2111 idealnumden(GEN nf, GEN x)
 2112 {
 2113   pari_sp av = avma;
 2114   GEN x0, ax, c, d, A, B, J;
 2115   long tx = idealtyp(&x,&ax);
 2116   nf = checknf(nf);
 2117   switch (tx)
 2118   {
 2119     case id_PRIME:
 2120       retmkvec2(idealhnf(nf, x), gen_1);
 2121     case id_PRINCIPAL:
 2122     {
 2123       GEN xZ, mx;
 2124       x = nf_to_scalar_or_basis(nf, x);
 2125       switch(typ(x))
 2126       {
 2127         case t_INT: return gerepilecopy(av, mkvec2(absi_shallow(x),gen_1));
 2128         case t_FRAC:return gerepilecopy(av, mkvec2(absi_shallow(gel(x,1)), gel(x,2)));
 2129       }
 2130       /* t_COL */
 2131       x = Q_remove_denom(x, &d);
 2132       if (!d) return gerepilecopy(av, mkvec2(idealhnf(nf, x), gen_1));
 2133       mx = zk_multable(nf, x);
 2134       xZ = zkmultable_capZ(mx);
 2135       x = ZM_hnfmodid(mx, xZ); /* principal ideal (x) */
 2136       x0 = mkvec2(xZ, mx); /* same, for fast multiplication */
 2137       break;
 2138     }
 2139     default: /* id_MAT */
 2140     {
 2141       long n = lg(x)-1;
 2142       if (n == 0) return mkvec2(gen_0, gen_1);
 2143       if (n != nf_get_degree(nf)) pari_err_DIM("idealnumden");
 2144       x0 = x = Q_remove_denom(x, &d);
 2145       if (!d) return gerepilecopy(av, mkvec2(x, gen_1));
 2146       break;
 2147     }
 2148   }
 2149   J = hnfmodid(x, d); /* = d/B */
 2150   c = gcoeff(J,1,1); /* (d/B) \cap Z, divides d */
 2151   B = idealHNF_inv_Z(nf, J); /* (d/B \cap Z) B/d */
 2152   if (!equalii(c,d)) B = ZM_Z_mul(B, diviiexact(d,c)); /* = B ! */
 2153   A = idealHNF_mul(nf, B, x0); /* d * (original x) * B = d A */
 2154   A = ZM_Z_divexact(A, d); /* = A ! */
 2155   return gerepilecopy(av, mkvec2(A, B));
 2156 }
 2157 
 2158 /* Return x, integral in 2-elt form, such that pr^n = c * x. Assume n != 0.
 2159  * nf = true nf */
 2160 static GEN
 2161 idealpowprime(GEN nf, GEN pr, GEN n, GEN *pc)
 2162 {
 2163   GEN p = pr_get_p(pr), q, gen;
 2164 
 2165   *pc = NULL;
 2166   if (is_pm1(n)) /* n = 1 special cased for efficiency */
 2167   {
 2168     q = p;
 2169     if (typ(pr_get_tau(pr)) == t_INT) /* inert */
 2170     {
 2171       *pc = (signe(n) >= 0)? p: ginv(p);
 2172       return mkvec2(gen_1,gen_0);
 2173     }
 2174     if (signe(n) >= 0) gen = pr_get_gen(pr);
 2175     else
 2176     {
 2177       gen = pr_get_tau(pr); /* possibly t_MAT */
 2178       *pc = ginv(p);
 2179     }
 2180   }
 2181   else if (equalis(n,2)) return idealsqrprime(nf, pr, pc);
 2182   else
 2183   {
 2184     long e = pr_get_e(pr), f = pr_get_f(pr);
 2185     GEN r, m = truedvmdis(n, e, &r);
 2186     if (e * f == nf_get_degree(nf))
 2187     { /* pr^e = (p) */
 2188       if (signe(m)) *pc = powii(p,m);
 2189       if (!signe(r)) return mkvec2(gen_1,gen_0);
 2190       q = p;
 2191       gen = nfpow(nf, pr_get_gen(pr), r);
 2192     }
 2193     else
 2194     {
 2195       m = absi_shallow(m);
 2196       if (signe(r)) m = addiu(m,1);
 2197       q = powii(p,m); /* m = ceil(|n|/e) */
 2198       if (signe(n) >= 0) gen = nfpow(nf, pr_get_gen(pr), n);
 2199       else
 2200       {
 2201         gen = pr_get_tau(pr);
 2202         if (typ(gen) == t_MAT) gen = gel(gen,1);
 2203         n = negi(n);
 2204         gen = ZC_Z_divexact(nfpow(nf, gen, n), powii(p, subii(n,m)));
 2205         *pc = ginv(q);
 2206       }
 2207     }
 2208     gen = FpC_red(gen, q);
 2209   }
 2210   return mkvec2(q, gen);
 2211 }
 2212 
 2213 /* x * pr^n. Assume x in HNF or scalar (possibly nonintegral) */
 2214 GEN
 2215 idealmulpowprime(GEN nf, GEN x, GEN pr, GEN n)
 2216 {
 2217   GEN c, cx, y;
 2218   long N;
 2219 
 2220   nf = checknf(nf);
 2221   N = nf_get_degree(nf);
 2222   if (!signe(n)) return typ(x) == t_MAT? x: scalarmat_shallow(x, N);
 2223 
 2224   /* inert, special cased for efficiency */
 2225   if (pr_is_inert(pr))
 2226   {
 2227     GEN q = powii(pr_get_p(pr), n);
 2228     return typ(x) == t_MAT? RgM_Rg_mul(x,q)
 2229                           : scalarmat_shallow(gmul(Q_abs(x),q), N);
 2230   }
 2231 
 2232   y = idealpowprime(nf, pr, n, &c);
 2233   if (typ(x) == t_MAT)
 2234   { x = Q_primitive_part(x, &cx); if (is_pm1(gcoeff(x,1,1))) x = NULL; }
 2235   else
 2236   { cx = x; x = NULL; }
 2237   cx = mul_content(c,cx);
 2238   if (x)
 2239     x = idealHNF_mul_two(nf,x,y);
 2240   else
 2241     x = idealhnf_two(nf,y);
 2242   if (cx) x = ZM_Q_mul(x,cx);
 2243   return x;
 2244 }
 2245 GEN
 2246 idealdivpowprime(GEN nf, GEN x, GEN pr, GEN n)
 2247 {
 2248   return idealmulpowprime(nf,x,pr, negi(n));
 2249 }
 2250 
 2251 /* nf = true nf */
 2252 static GEN
 2253 idealpow_aux(GEN nf, GEN x, long tx, GEN n)
 2254 {
 2255   GEN T = nf_get_pol(nf), m, cx, n1, a, alpha;
 2256   long N = degpol(T), s = signe(n);
 2257   if (!s) return matid(N);
 2258   switch(tx)
 2259   {
 2260     case id_PRINCIPAL:
 2261       return idealhnf_principal(nf, nfpow(nf,x,n));
 2262     case id_PRIME:
 2263       if (pr_is_inert(x)) return scalarmat(powii(gel(x,1), n), N);
 2264       x = idealpowprime(nf, x, n, &cx);
 2265       x = idealhnf_two(nf,x);
 2266       return cx? ZM_Q_mul(x, cx): x;
 2267     default:
 2268       if (is_pm1(n)) return (s < 0)? idealinv(nf, x): gcopy(x);
 2269       n1 = (s < 0)? negi(n): n;
 2270 
 2271       x = Q_primitive_part(x, &cx);
 2272       a = mat_ideal_two_elt(nf,x); alpha = gel(a,2); a = gel(a,1);
 2273       alpha = nfpow(nf,alpha,n1);
 2274       m = zk_scalar_or_multable(nf, alpha);
 2275       if (typ(m) == t_INT) {
 2276         x = gcdii(powii(a,n1), m);
 2277         if (s<0) x = ginv(x);
 2278         if (cx) x = gmul(x, powgi(cx,n));
 2279         x = scalarmat(x, N);
 2280       }
 2281       else
 2282       { /* could use gcdii(powii(a,n1), zkmultable_capZ(m)), but costly */
 2283         x = ZM_hnfmodid(m, powii(a,n1));
 2284         if (cx) cx = powgi(cx,n);
 2285         if (s<0) {
 2286           GEN xZ = gcoeff(x,1,1);
 2287           cx = cx ? gdiv(cx, xZ): ginv(xZ);
 2288           x = idealHNF_inv_Z(nf,x);
 2289         }
 2290         if (cx) x = ZM_Q_mul(x, cx);
 2291       }
 2292       return x;
 2293   }
 2294 }
 2295 
 2296 /* raise the ideal x to the power n (in Z) */
 2297 GEN
 2298 idealpow(GEN nf, GEN x, GEN n)
 2299 {
 2300   pari_sp av;
 2301   long tx;
 2302   GEN res, ax;
 2303 
 2304   if (typ(n) != t_INT) pari_err_TYPE("idealpow",n);
 2305   tx = idealtyp(&x,&ax);
 2306   res = ax? cgetg(3,t_VEC): NULL;
 2307   av = avma;
 2308   x = gerepileupto(av, idealpow_aux(checknf(nf), x, tx, n));
 2309   if (!ax) return x;
 2310   ax = ext_pow(nf, ax, n);
 2311   gel(res,1) = x;
 2312   gel(res,2) = ax;
 2313   return res;
 2314 }
 2315 
 2316 /* Return ideal^e in number field nf. e is a C integer. */
 2317 GEN
 2318 idealpows(GEN nf, GEN ideal, long e)
 2319 {
 2320   long court[] = {evaltyp(t_INT) | _evallg(3),0,0};
 2321   affsi(e,court); return idealpow(nf,ideal,court);
 2322 }
 2323 
 2324 static GEN
 2325 _idealmulred(GEN nf, GEN x, GEN y)
 2326 { return idealred(nf,idealmul(nf,x,y)); }
 2327 static GEN
 2328 _idealsqrred(GEN nf, GEN x)
 2329 { return idealred(nf,idealsqr(nf,x)); }
 2330 static GEN
 2331 _mul(void *data, GEN x, GEN y) { return _idealmulred((GEN)data,x,y); }
 2332 static GEN
 2333 _sqr(void *data, GEN x) { return _idealsqrred((GEN)data, x); }
 2334 
 2335 /* compute x^n (x ideal, n integer), reducing along the way */
 2336 GEN
 2337 idealpowred(GEN nf, GEN x, GEN n)
 2338 {
 2339   pari_sp av = avma, av2;
 2340   long s;
 2341   GEN y;
 2342 
 2343   if (typ(n) != t_INT) pari_err_TYPE("idealpowred",n);
 2344   s = signe(n); if (s == 0) return idealpow(nf,x,n);
 2345   y = gen_pow_i(x, n, (void*)nf, &_sqr, &_mul);
 2346   av2 = avma;
 2347   if (s < 0) y = idealinv(nf,y);
 2348   if (s < 0 || is_pm1(n)) y = idealred(nf,y);
 2349   return avma == av2? gerepilecopy(av,y): gerepileupto(av,y);
 2350 }
 2351 
 2352 GEN
 2353 idealmulred(GEN nf, GEN x, GEN y)
 2354 {
 2355   pari_sp av = avma;
 2356   return gerepileupto(av, _idealmulred(nf,x,y));
 2357 }
 2358 
 2359 long
 2360 isideal(GEN nf,GEN x)
 2361 {
 2362   long N, i, j, lx, tx = typ(x);
 2363   pari_sp av;
 2364   GEN T, xZ;
 2365 
 2366   nf = checknf(nf); T = nf_get_pol(nf); lx = lg(x);
 2367   if (tx==t_VEC && lx==3) { x = gel(x,1); tx = typ(x); lx = lg(x); }
 2368   switch(tx)
 2369   {
 2370     case t_INT: case t_FRAC: return 1;
 2371     case t_POL: return varn(x) == varn(T);
 2372     case t_POLMOD: return RgX_equal_var(T, gel(x,1));
 2373     case t_VEC: return get_prid(x)? 1 : 0;
 2374     case t_MAT: break;
 2375     default: return 0;
 2376   }
 2377   N = degpol(T);
 2378   if (lx-1 != N) return (lx == 1);
 2379   if (nbrows(x) != N) return 0;
 2380 
 2381   av = avma; x = Q_primpart(x);
 2382   if (!ZM_ishnf(x)) return 0;
 2383   xZ = gcoeff(x,1,1);
 2384   for (j=2; j<=N; j++)
 2385     if (!dvdii(xZ, gcoeff(x,j,j))) return gc_long(av,0);
 2386   for (i=2; i<=N; i++)
 2387     for (j=2; j<=N; j++)
 2388        if (! hnf_invimage(x, zk_ei_mul(nf,gel(x,i),j))) return gc_long(av,0);
 2389   return gc_long(av,1);
 2390 }
 2391 
 2392 GEN
 2393 idealdiv(GEN nf, GEN x, GEN y)
 2394 {
 2395   pari_sp av = avma, tetpil;
 2396   GEN z = idealinv(nf,y);
 2397   tetpil = avma; return gerepile(av,tetpil, idealmul(nf,x,z));
 2398 }
 2399 
 2400 /* This routine computes the quotient x/y of two ideals in the number field nf.
 2401  * It assumes that the quotient is an integral ideal.  The idea is to find an
 2402  * ideal z dividing y such that gcd(Nx/Nz, Nz) = 1.  Then
 2403  *
 2404  *   x + (Nx/Nz)    x
 2405  *   ----------- = ---
 2406  *   y + (Ny/Nz)    y
 2407  *
 2408  * Proof: we can assume x and y are integral. Let p be any prime ideal
 2409  *
 2410  * If p | Nz, then it divides neither Nx/Nz nor Ny/Nz (since Nx/Nz is the
 2411  * product of the integers N(x/y) and N(y/z)).  Both the numerator and the
 2412  * denominator on the left will be coprime to p.  So will x/y, since x/y is
 2413  * assumed integral and its norm N(x/y) is coprime to p.
 2414  *
 2415  * If instead p does not divide Nz, then v_p (Nx/Nz) = v_p (Nx) >= v_p(x).
 2416  * Hence v_p (x + Nx/Nz) = v_p(x).  Likewise for the denominators.  QED.
 2417  *
 2418  *                Peter Montgomery.  July, 1994. */
 2419 static void
 2420 err_divexact(GEN x, GEN y)
 2421 { pari_err_DOMAIN("idealdivexact","denominator(x/y)", "!=",
 2422                   gen_1,mkvec2(x,y)); }
 2423 GEN
 2424 idealdivexact(GEN nf, GEN x0, GEN y0)
 2425 {
 2426   pari_sp av = avma;
 2427   GEN x, y, xZ, yZ, Nx, Ny, Nz, cy, q, r;
 2428 
 2429   nf = checknf(nf);
 2430   x = idealhnf_shallow(nf, x0);
 2431   y = idealhnf_shallow(nf, y0);
 2432   if (lg(y) == 1) pari_err_INV("idealdivexact", y0);
 2433   if (lg(x) == 1) { set_avma(av); return cgetg(1, t_MAT); } /* numerator is zero */
 2434   y = Q_primitive_part(y, &cy);
 2435   if (cy) x = RgM_Rg_div(x,cy);
 2436   xZ = gcoeff(x,1,1); if (typ(xZ) != t_INT) err_divexact(x,y);
 2437   yZ = gcoeff(y,1,1); if (isint1(yZ)) return gerepilecopy(av, x);
 2438   Nx = idealnorm(nf,x);
 2439   Ny = idealnorm(nf,y);
 2440   if (typ(Nx) != t_INT) err_divexact(x,y);
 2441   q = dvmdii(Nx,Ny, &r);
 2442   if (signe(r)) err_divexact(x,y);
 2443   if (is_pm1(q)) { set_avma(av); return matid(nf_get_degree(nf)); }
 2444   /* Find a norm Nz | Ny such that gcd(Nx/Nz, Nz) = 1 */
 2445   for (Nz = Ny;;) /* q = Nx/Nz */
 2446   {
 2447     GEN p1 = gcdii(Nz, q);
 2448     if (is_pm1(p1)) break;
 2449     Nz = diviiexact(Nz,p1);
 2450     q = mulii(q,p1);
 2451   }
 2452   xZ = gcoeff(x,1,1); q = gcdii(q, xZ);
 2453   if (!equalii(xZ,q))
 2454   { /* Replace x/y  by  x+(Nx/Nz) / y+(Ny/Nz) */
 2455     x = ZM_hnfmodid(x, q);
 2456     /* y reduced to unit ideal ? */
 2457     if (Nz == Ny) return gerepileupto(av, x);
 2458 
 2459     yZ = gcoeff(y,1,1); q = gcdii(diviiexact(Ny,Nz), yZ);
 2460     y = ZM_hnfmodid(y, q);
 2461   }
 2462   yZ = gcoeff(y,1,1);
 2463   y = idealHNF_mul(nf,x, idealHNF_inv_Z(nf,y));
 2464   return gerepileupto(av, ZM_Z_divexact(y, yZ));
 2465 }
 2466 
 2467 GEN
 2468 idealintersect(GEN nf, GEN x, GEN y)
 2469 {
 2470   pari_sp av = avma;
 2471   long lz, lx, i;
 2472   GEN z, dx, dy, xZ, yZ;;
 2473 
 2474   nf = checknf(nf);
 2475   x = idealhnf_shallow(nf,x);
 2476   y = idealhnf_shallow(nf,y);
 2477   if (lg(x) == 1 || lg(y) == 1) { set_avma(av); return cgetg(1,t_MAT); }
 2478   x = Q_remove_denom(x, &dx);
 2479   y = Q_remove_denom(y, &dy);
 2480   if (dx) y = ZM_Z_mul(y, dx);
 2481   if (dy) x = ZM_Z_mul(x, dy);
 2482   xZ = gcoeff(x,1,1);
 2483   yZ = gcoeff(y,1,1);
 2484   dx = mul_denom(dx,dy);
 2485   z = ZM_lll(shallowconcat(x,y), 0.99, LLL_KER); lz = lg(z);
 2486   lx = lg(x);
 2487   for (i=1; i<lz; i++) setlg(z[i], lx);
 2488   z = ZM_hnfmodid(ZM_mul(x,z), lcmii(xZ, yZ));
 2489   if (dx) z = RgM_Rg_div(z,dx);
 2490   return gerepileupto(av,z);
 2491 }
 2492 
 2493 /*******************************************************************/
 2494 /*                                                                 */
 2495 /*                      T2-IDEAL REDUCTION                         */
 2496 /*                                                                 */
 2497 /*******************************************************************/
 2498 
 2499 static GEN
 2500 chk_vdir(GEN nf, GEN vdir)
 2501 {
 2502   long i, l = lg(vdir);
 2503   GEN v;
 2504   if (l != lg(nf_get_roots(nf))) pari_err_DIM("idealred");
 2505   switch(typ(vdir))
 2506   {
 2507     case t_VECSMALL: return vdir;
 2508     case t_VEC: break;
 2509     default: pari_err_TYPE("idealred",vdir);
 2510   }
 2511   v = cgetg(l, t_VECSMALL);
 2512   for (i = 1; i < l; i++) v[i] = itos(gceil(gel(vdir,i)));
 2513   return v;
 2514 }
 2515 
 2516 static void
 2517 twistG(GEN G, long r1, long i, long v)
 2518 {
 2519   long j, lG = lg(G);
 2520   if (i <= r1) {
 2521     for (j=1; j<lG; j++) gcoeff(G,i,j) = gmul2n(gcoeff(G,i,j), v);
 2522   } else {
 2523     long k = (i<<1) - r1;
 2524     for (j=1; j<lG; j++)
 2525     {
 2526       gcoeff(G,k-1,j) = gmul2n(gcoeff(G,k-1,j), v);
 2527       gcoeff(G,k  ,j) = gmul2n(gcoeff(G,k  ,j), v);
 2528     }
 2529   }
 2530 }
 2531 
 2532 GEN
 2533 nf_get_Gtwist(GEN nf, GEN vdir)
 2534 {
 2535   long i, l, v, r1;
 2536   GEN G;
 2537 
 2538   if (!vdir) return nf_get_roundG(nf);
 2539   if (typ(vdir) == t_MAT)
 2540   {
 2541     long N = nf_get_degree(nf);
 2542     if (lg(vdir) != N+1 || lgcols(vdir) != N+1) pari_err_DIM("idealred");
 2543     return vdir;
 2544   }
 2545   vdir = chk_vdir(nf, vdir);
 2546   G = RgM_shallowcopy(nf_get_G(nf));
 2547   r1 = nf_get_r1(nf);
 2548   l = lg(vdir);
 2549   for (i=1; i<l; i++)
 2550   {
 2551     v = vdir[i]; if (!v) continue;
 2552     twistG(G, r1, i, v);
 2553   }
 2554   return RM_round_maxrank(G);
 2555 }
 2556 GEN
 2557 nf_get_Gtwist1(GEN nf, long i)
 2558 {
 2559   GEN G = RgM_shallowcopy( nf_get_G(nf) );
 2560   long r1 = nf_get_r1(nf);
 2561   twistG(G, r1, i, 10);
 2562   return RM_round_maxrank(G);
 2563 }
 2564 
 2565 GEN
 2566 RM_round_maxrank(GEN G0)
 2567 {
 2568   long e, r = lg(G0)-1;
 2569   pari_sp av = avma;
 2570   for (e = 4; ; e <<= 1, set_avma(av))
 2571   {
 2572     GEN G = gmul2n(G0, e), H = ground(G);
 2573     if (ZM_rank(H) == r) return H; /* maximal rank ? */
 2574   }
 2575 }
 2576 
 2577 GEN
 2578 idealred0(GEN nf, GEN I, GEN vdir)
 2579 {
 2580   pari_sp av = avma;
 2581   GEN G, aI, IZ, J, y, yZ, my, c1 = NULL;
 2582   long N;
 2583 
 2584   nf = checknf(nf);
 2585   N = nf_get_degree(nf);
 2586   /* put first for sanity checks, unused when I obviously principal */
 2587   G = nf_get_Gtwist(nf, vdir);
 2588   switch (idealtyp(&I,&aI))
 2589   {
 2590     case id_PRIME:
 2591       if (pr_is_inert(I)) {
 2592         if (!aI) { set_avma(av); return matid(N); }
 2593         c1 = gel(I,1); I = matid(N);
 2594         goto END;
 2595       }
 2596       IZ = pr_get_p(I);
 2597       J = pr_inv_p(I);
 2598       I = idealhnf_two(nf,I);
 2599       break;
 2600     case id_MAT:
 2601       I = Q_primitive_part(I, &c1);
 2602       IZ = gcoeff(I,1,1);
 2603       if (is_pm1(IZ))
 2604       {
 2605         if (!aI) { set_avma(av); return matid(N); }
 2606         goto END;
 2607       }
 2608       J = idealHNF_inv_Z(nf, I);
 2609       break;
 2610     default: /* id_PRINCIPAL, silly case */
 2611       if (gequal0(I)) I = cgetg(1,t_MAT); else { c1 = I; I = matid(N); }
 2612       if (!aI) return I;
 2613       goto END;
 2614   }
 2615   /* now I integral, HNF; and J = (I\cap Z) I^(-1), integral */
 2616   y = idealpseudomin(J, G); /* small elt in (I\cap Z)I^(-1), integral */
 2617   if (equalii(ZV_content(y), IZ))
 2618   { /* already reduced */
 2619     if (!aI) return gerepilecopy(av, I);
 2620     goto END;
 2621   }
 2622 
 2623   my = zk_multable(nf, y);
 2624   I = ZM_Z_divexact(ZM_mul(my, I), IZ); /* y I / (I\cap Z), integral */
 2625   c1 = mul_content(c1, IZ);
 2626   my = ZM_gauss(my, col_ei(N,1)); /* y^-1 */
 2627   yZ = Q_denom(my); /* (y) \cap Z */
 2628   I = hnfmodid(I, yZ);
 2629   if (!aI) return gerepileupto(av, I);
 2630   c1 = RgC_Rg_mul(my, c1);
 2631 END:
 2632   if (c1) aI = ext_mul(nf, aI,c1);
 2633   return gerepilecopy(av, mkvec2(I, aI));
 2634 }
 2635 
 2636 GEN
 2637 idealmin(GEN nf, GEN x, GEN vdir)
 2638 {
 2639   pari_sp av = avma;
 2640   GEN y, dx;
 2641   nf = checknf(nf);
 2642   switch( idealtyp(&x,&y) )
 2643   {
 2644     case id_PRINCIPAL: return gcopy(x);
 2645     case id_PRIME: x = pr_hnf(nf,x); break;
 2646     case id_MAT: if (lg(x) == 1) return gen_0;
 2647   }
 2648   x = Q_remove_denom(x, &dx);
 2649   y = idealpseudomin(x, nf_get_Gtwist(nf,vdir));
 2650   if (dx) y = RgC_Rg_div(y, dx);
 2651   return gerepileupto(av, y);
 2652 }
 2653 
 2654 /*******************************************************************/
 2655 /*                                                                 */
 2656 /*                   APPROXIMATION THEOREM                         */
 2657 /*                                                                 */
 2658 /*******************************************************************/
 2659 /* a = ppi(a,b) ppo(a,b), where ppi regroups primes common to a and b
 2660  * and ppo(a,b) = Z_ppo(a,b) */
 2661 /* return gcd(a,b),ppi(a,b),ppo(a,b) */
 2662 GEN
 2663 Z_ppio(GEN a, GEN b)
 2664 {
 2665   GEN x, y, d = gcdii(a,b);
 2666   if (is_pm1(d)) return mkvec3(gen_1, gen_1, a);
 2667   x = d; y = diviiexact(a,d);
 2668   for(;;)
 2669   {
 2670     GEN g = gcdii(x,y);
 2671     if (is_pm1(g)) return mkvec3(d, x, y);
 2672     x = mulii(x,g); y = diviiexact(y,g);
 2673   }
 2674 }
 2675 /* a = ppg(a,b)pple(a,b), where ppg regroups primes such that v(a) > v(b)
 2676  * and pple all others */
 2677 /* return gcd(a,b),ppg(a,b),pple(a,b) */
 2678 GEN
 2679 Z_ppgle(GEN a, GEN b)
 2680 {
 2681   GEN x, y, g, d = gcdii(a,b);
 2682   if (equalii(a, d)) return mkvec3(a, gen_1, a);
 2683   x = diviiexact(a,d); y = d;
 2684   for(;;)
 2685   {
 2686     g = gcdii(x,y);
 2687     if (is_pm1(g)) return mkvec3(d, x, y);
 2688     x = mulii(x,g); y = diviiexact(y,g);
 2689   }
 2690 }
 2691 static void
 2692 Z_dcba_rec(GEN L, GEN a, GEN b)
 2693 {
 2694   GEN x, r, v, g, h, c, c0;
 2695   long n;
 2696   if (is_pm1(b)) {
 2697     if (!is_pm1(a)) vectrunc_append(L, a);
 2698     return;
 2699   }
 2700   v = Z_ppio(a,b);
 2701   a = gel(v,2);
 2702   r = gel(v,3);
 2703   if (!is_pm1(r)) vectrunc_append(L, r);
 2704   v = Z_ppgle(a,b);
 2705   g = gel(v,1);
 2706   h = gel(v,2);
 2707   x = c0 = gel(v,3);
 2708   for (n = 1; !is_pm1(h); n++)
 2709   {
 2710     GEN d, y;
 2711     long i;
 2712     v = Z_ppgle(h,sqri(g));
 2713     g = gel(v,1);
 2714     h = gel(v,2);
 2715     c = gel(v,3); if (is_pm1(c)) continue;
 2716     d = gcdii(c,b);
 2717     x = mulii(x,d);
 2718     y = d; for (i=1; i < n; i++) y = sqri(y);
 2719     Z_dcba_rec(L, diviiexact(c,y), d);
 2720   }
 2721   Z_dcba_rec(L,diviiexact(b,x), c0);
 2722 }
 2723 static GEN
 2724 Z_cba_rec(GEN L, GEN a, GEN b)
 2725 {
 2726   GEN g;
 2727   if (lg(L) > 10)
 2728   { /* a few naive steps before switching to dcba */
 2729     Z_dcba_rec(L, a, b);
 2730     return gel(L, lg(L)-1);
 2731   }
 2732   if (is_pm1(a)) return b;
 2733   g = gcdii(a,b);
 2734   if (is_pm1(g)) { vectrunc_append(L, a); return b; }
 2735   a = diviiexact(a,g);
 2736   b = diviiexact(b,g);
 2737   return Z_cba_rec(L, Z_cba_rec(L, a, g), b);
 2738 }
 2739 GEN
 2740 Z_cba(GEN a, GEN b)
 2741 {
 2742   GEN L = vectrunc_init(expi(a) + expi(b) + 2);
 2743   GEN t = Z_cba_rec(L, a, b);
 2744   if (!is_pm1(t)) vectrunc_append(L, t);
 2745   return L;
 2746 }
 2747 /* P = coprime base, extend it by b; TODO: quadratic for now */
 2748 GEN
 2749 ZV_cba_extend(GEN P, GEN b)
 2750 {
 2751   long i, l = lg(P);
 2752   GEN w = cgetg(l+1, t_VEC);
 2753   for (i = 1; i < l; i++)
 2754   {
 2755     GEN v = Z_cba(gel(P,i), b);
 2756     long nv = lg(v)-1;
 2757     gel(w,i) = vecslice(v, 1, nv-1); /* those divide P[i] but not b */
 2758     b = gel(v,nv);
 2759   }
 2760   gel(w,l) = b; return shallowconcat1(w);
 2761 }
 2762 GEN
 2763 ZV_cba(GEN v)
 2764 {
 2765   long i, l = lg(v);
 2766   GEN P;
 2767   if (l <= 2) return v;
 2768   P = Z_cba(gel(v,1), gel(v,2));
 2769   for (i = 3; i < l; i++) P = ZV_cba_extend(P, gel(v,i));
 2770   return P;
 2771 }
 2772 
 2773 /* write x = x1 x2, x2 maximal s.t. (x2,f) = 1, return x2 */
 2774 GEN
 2775 Z_ppo(GEN x, GEN f)
 2776 {
 2777   for (;;)
 2778   {
 2779     f = gcdii(x, f); if (is_pm1(f)) break;
 2780     x = diviiexact(x, f);
 2781   }
 2782   return x;
 2783 }
 2784 /* write x = x1 x2, x2 maximal s.t. (x2,f) = 1, return x2 */
 2785 ulong
 2786 u_ppo(ulong x, ulong f)
 2787 {
 2788   for (;;)
 2789   {
 2790     f = ugcd(x, f); if (f == 1) break;
 2791     x /= f;
 2792   }
 2793   return x;
 2794 }
 2795 
 2796 /* result known to be representable as an ulong */
 2797 static ulong
 2798 lcmuu(ulong a, ulong b) { ulong d = ugcd(a,b); return (a/d) * b; }
 2799 
 2800 /* assume 0 < x < N; return u in (Z/NZ)^* such that u x = gcd(x,N) (mod N);
 2801  * set *pd = gcd(x,N) */
 2802 ulong
 2803 Fl_invgen(ulong x, ulong N, ulong *pd)
 2804 {
 2805   ulong d, d0, e, v, v1;
 2806   long s;
 2807   *pd = d = xgcduu(N, x, 0, &v, &v1, &s);
 2808   if (s > 0) v = N - v;
 2809   if (d == 1) return v;
 2810   /* vx = gcd(x,N) (mod N), v coprime to N/d but need not be coprime to N */
 2811   e = N / d;
 2812   d0 = u_ppo(d, e); /* d = d0 d1, d0 coprime to N/d, rad(d1) | N/d */
 2813   if (d0 == 1) return v;
 2814   e = lcmuu(e, d / d0);
 2815   return u_chinese_coprime(v, 1, e, d0, e*d0);
 2816 }
 2817 
 2818 /* x t_INT, f ideal. Write x = x1 x2, sqf(x1) | f, (x2,f) = 1. Return x2 */
 2819 static GEN
 2820 nf_coprime_part(GEN nf, GEN x, GEN listpr)
 2821 {
 2822   long v, j, lp = lg(listpr), N = nf_get_degree(nf);
 2823   GEN x1, x2, ex;
 2824 
 2825 #if 0 /*1) via many gcds. Expensive ! */
 2826   GEN f = idealprodprime(nf, listpr);
 2827   f = ZM_hnfmodid(f, x); /* first gcd is less expensive since x in Z */
 2828   x = scalarmat(x, N);
 2829   for (;;)
 2830   {
 2831     if (gequal1(gcoeff(f,1,1))) break;
 2832     x = idealdivexact(nf, x, f);
 2833     f = ZM_hnfmodid(shallowconcat(f,x), gcoeff(x,1,1)); /* gcd(f,x) */
 2834   }
 2835   x2 = x;
 2836 #else /*2) from prime decomposition */
 2837   x1 = NULL;
 2838   for (j=1; j<lp; j++)
 2839   {
 2840     GEN pr = gel(listpr,j);
 2841     v = Z_pval(x, pr_get_p(pr)); if (!v) continue;
 2842 
 2843     ex = muluu(v, pr_get_e(pr)); /* = v_pr(x) > 0 */
 2844     x1 = x1? idealmulpowprime(nf, x1, pr, ex)
 2845            : idealpow(nf, pr, ex);
 2846   }
 2847   x = scalarmat(x, N);
 2848   x2 = x1? idealdivexact(nf, x, x1): x;
 2849 #endif
 2850   return x2;
 2851 }
 2852 
 2853 /* L0 in K^*, assume (L0,f) = 1. Return L integral, L0 = L mod f  */
 2854 GEN
 2855 make_integral(GEN nf, GEN L0, GEN f, GEN listpr)
 2856 {
 2857   GEN fZ, t, L, D2, d1, d2, d;
 2858 
 2859   L = Q_remove_denom(L0, &d);
 2860   if (!d) return L0;
 2861 
 2862   /* L0 = L / d, L integral */
 2863   fZ = gcoeff(f,1,1);
 2864   if (typ(L) == t_INT) return Fp_mul(L, Fp_inv(d, fZ), fZ);
 2865   /* Kill denom part coprime to fZ */
 2866   d2 = Z_ppo(d, fZ);
 2867   t = Fp_inv(d2, fZ); if (!is_pm1(t)) L = ZC_Z_mul(L,t);
 2868   if (equalii(d, d2)) return L;
 2869 
 2870   d1 = diviiexact(d, d2);
 2871   /* L0 = (L / d1) mod f. d1 not coprime to f
 2872    * write (d1) = D1 D2, D2 minimal, (D2,f) = 1. */
 2873   D2 = nf_coprime_part(nf, d1, listpr);
 2874   t = idealaddtoone_i(nf, D2, f); /* in D2, 1 mod f */
 2875   L = nfmuli(nf,t,L);
 2876 
 2877   /* if (L0, f) = 1, then L in D1 ==> in D1 D2 = (d1) */
 2878   return Q_div_to_int(L, d1); /* exact division */
 2879 }
 2880 
 2881 /* assume L is a list of prime ideals. Return the product */
 2882 GEN
 2883 idealprodprime(GEN nf, GEN L)
 2884 {
 2885   long l = lg(L), i;
 2886   GEN z;
 2887   if (l == 1) return matid(nf_get_degree(nf));
 2888   z = pr_hnf(nf, gel(L,1));
 2889   for (i=2; i<l; i++) z = idealHNF_mul_two(nf,z, gel(L,i));
 2890   return z;
 2891 }
 2892 
 2893 /* optimize for the frequent case I = nfhnf()[2]: lots of them are 1 */
 2894 GEN
 2895 idealprod(GEN nf, GEN I)
 2896 {
 2897   long i, l = lg(I);
 2898   GEN z;
 2899   for (i = 1; i < l; i++)
 2900     if (!equali1(gel(I,i))) break;
 2901   if (i == l) return gen_1;
 2902   z = gel(I,i);
 2903   for (i++; i<l; i++) z = idealmul(nf, z, gel(I,i));
 2904   return z;
 2905 }
 2906 
 2907 /* v_pr(idealprod(nf,I)) */
 2908 long
 2909 idealprodval(GEN nf, GEN I, GEN pr)
 2910 {
 2911   long i, l = lg(I), v = 0;
 2912   for (i = 1; i < l; i++)
 2913     if (!equali1(gel(I,i))) v += idealval(nf, gel(I,i), pr);
 2914   return v;
 2915 }
 2916 
 2917 /* assume L is a list of prime ideals. Return prod L[i]^e[i] */
 2918 GEN
 2919 factorbackprime(GEN nf, GEN L, GEN e)
 2920 {
 2921   long l = lg(L), i;
 2922   GEN z;
 2923 
 2924   if (l == 1) return matid(nf_get_degree(nf));
 2925   z = idealpow(nf, gel(L,1), gel(e,1));
 2926   for (i=2; i<l; i++)
 2927     if (signe(gel(e,i))) z = idealmulpowprime(nf,z, gel(L,i),gel(e,i));
 2928   return z;
 2929 }
 2930 
 2931 /* F in Z, divisible exactly by pr.p. Return F-uniformizer for pr, i.e.
 2932  * a t in Z_K such that v_pr(t) = 1 and (t, F/pr) = 1 */
 2933 GEN
 2934 pr_uniformizer(GEN pr, GEN F)
 2935 {
 2936   GEN p = pr_get_p(pr), t = pr_get_gen(pr);
 2937   if (!equalii(F, p))
 2938   {
 2939     long e = pr_get_e(pr);
 2940     GEN u, v, q = (e == 1)? sqri(p): p;
 2941     u = mulii(q, Fp_inv(q, diviiexact(F,p))); /* 1 mod F/p, 0 mod q */
 2942     v = subui(1UL, u); /* 0 mod F/p, 1 mod q */
 2943     if (pr_is_inert(pr))
 2944       t = addii(mulii(p, v), u);
 2945     else
 2946     {
 2947       t = ZC_Z_mul(t, v);
 2948       gel(t,1) = addii(gel(t,1), u); /* return u + vt */
 2949     }
 2950   }
 2951   return t;
 2952 }
 2953 /* L = list of prime ideals, return lcm_i (L[i] \cap \ZM) */
 2954 GEN
 2955 prV_lcm_capZ(GEN L)
 2956 {
 2957   long i, r = lg(L);
 2958   GEN F;
 2959   if (r == 1) return gen_1;
 2960   F = pr_get_p(gel(L,1));
 2961   for (i = 2; i < r; i++)
 2962   {
 2963     GEN pr = gel(L,i), p = pr_get_p(pr);
 2964     if (!dvdii(F, p)) F = mulii(F,p);
 2965   }
 2966   return F;
 2967 }
 2968 
 2969 /* Given a prime ideal factorization with possibly zero or negative
 2970  * exponents, gives b such that v_p(b) = v_p(x) for all prime ideals pr | x
 2971  * and v_pr(b) >= 0 for all other pr.
 2972  * For optimal performance, all [anti-]uniformizers should be precomputed,
 2973  * but no support for this yet.
 2974  *
 2975  * If nored, do not reduce result.
 2976  * No garbage collecting */
 2977 static GEN
 2978 idealapprfact_i(GEN nf, GEN x, int nored)
 2979 {
 2980   GEN z, d, L, e, e2, F;
 2981   long i, r;
 2982   int flagden;
 2983 
 2984   nf = checknf(nf);
 2985   L = gel(x,1);
 2986   e = gel(x,2);
 2987   F = prV_lcm_capZ(L);
 2988   flagden = 0;
 2989   z = NULL; r = lg(e);
 2990   for (i = 1; i < r; i++)
 2991   {
 2992     long s = signe(gel(e,i));
 2993     GEN pi, q;
 2994     if (!s) continue;
 2995     if (s < 0) flagden = 1;
 2996     pi = pr_uniformizer(gel(L,i), F);
 2997     q = nfpow(nf, pi, gel(e,i));
 2998     z = z? nfmul(nf, z, q): q;
 2999   }
 3000   if (!z) return gen_1;
 3001   if (nored || typ(z) != t_COL) return z;
 3002   e2 = cgetg(r, t_VEC);
 3003   for (i=1; i<r; i++) gel(e2,i) = addiu(gel(e,i), 1);
 3004   x = factorbackprime(nf, L,e2);
 3005   if (flagden) /* denominator */
 3006   {
 3007     z = Q_remove_denom(z, &d);
 3008     d = diviiexact(d, Z_ppo(d, F));
 3009     x = RgM_Rg_mul(x, d);
 3010   }
 3011   else
 3012     d = NULL;
 3013   z = ZC_reducemodlll(z, x);
 3014   return d? RgC_Rg_div(z,d): z;
 3015 }
 3016 
 3017 GEN
 3018 idealapprfact(GEN nf, GEN x) {
 3019   pari_sp av = avma;
 3020   return gerepileupto(av, idealapprfact_i(nf, x, 0));
 3021 }
 3022 GEN
 3023 idealappr(GEN nf, GEN x) {
 3024   pari_sp av = avma;
 3025   if (!is_nf_extfactor(x)) x = idealfactor(nf, x);
 3026   return gerepileupto(av, idealapprfact_i(nf, x, 0));
 3027 }
 3028 
 3029 /* OBSOLETE */
 3030 GEN
 3031 idealappr0(GEN nf, GEN x, long fl) { (void)fl; return idealappr(nf, x); }
 3032 
 3033 static GEN
 3034 mat_ideal_two_elt2(GEN nf, GEN x, GEN a)
 3035 {
 3036   GEN F = idealfactor(nf,a), P = gel(F,1), E = gel(F,2);
 3037   long i, r = lg(E);
 3038   for (i=1; i<r; i++) gel(E,i) = stoi( idealval(nf,x,gel(P,i)) );
 3039   return idealapprfact_i(nf,F,1);
 3040 }
 3041 
 3042 static void
 3043 not_in_ideal(GEN a) {
 3044   pari_err_DOMAIN("idealtwoelt2","element mod ideal", "!=", gen_0, a);
 3045 }
 3046 /* x integral in HNF, a an 'nf' */
 3047 static int
 3048 in_ideal(GEN x, GEN a)
 3049 {
 3050   switch(typ(a))
 3051   {
 3052     case t_INT: return dvdii(a, gcoeff(x,1,1));
 3053     case t_COL: return RgV_is_ZV(a) && !!hnf_invimage(x, a);
 3054     default: return 0;
 3055   }
 3056 }
 3057 
 3058 /* Given an integral ideal x and a in x, gives a b such that
 3059  * x = aZ_K + bZ_K using the approximation theorem */
 3060 GEN
 3061 idealtwoelt2(GEN nf, GEN x, GEN a)
 3062 {
 3063   pari_sp av = avma;
 3064   GEN cx, b;
 3065 
 3066   nf = checknf(nf);
 3067   a = nf_to_scalar_or_basis(nf, a);
 3068   x = idealhnf_shallow(nf,x);
 3069   if (lg(x) == 1)
 3070   {
 3071     if (!isintzero(a)) not_in_ideal(a);
 3072     set_avma(av); return gen_0;
 3073   }
 3074   x = Q_primitive_part(x, &cx);
 3075   if (cx) a = gdiv(a, cx);
 3076   if (!in_ideal(x, a)) not_in_ideal(a);
 3077   b = mat_ideal_two_elt2(nf, x, a);
 3078   if (typ(b) == t_COL)
 3079   {
 3080     GEN mod = idealhnf_principal(nf,a);
 3081     b = ZC_hnfrem(b,mod);
 3082     if (ZV_isscalar(b)) b = gel(b,1);
 3083   }
 3084   else
 3085   {
 3086     GEN aZ = typ(a) == t_COL? Q_denom(zk_inv(nf,a)): a; /* (a) \cap Z */
 3087     b = centermodii(b, aZ, shifti(aZ,-1));
 3088   }
 3089   b = cx? gmul(b,cx): gcopy(b);
 3090   return gerepileupto(av, b);
 3091 }
 3092 
 3093 /* Given 2 integral ideals x and y in nf, returns a beta in nf such that
 3094  * beta * x is an integral ideal coprime to y */
 3095 GEN
 3096 idealcoprimefact(GEN nf, GEN x, GEN fy)
 3097 {
 3098   GEN L = gel(fy,1), e;
 3099   long i, r = lg(L);
 3100 
 3101   e = cgetg(r, t_COL);
 3102   for (i=1; i<r; i++) gel(e,i) = stoi( -idealval(nf,x,gel(L,i)) );
 3103   return idealapprfact_i(nf, mkmat2(L,e), 0);
 3104 }
 3105 GEN
 3106 idealcoprime(GEN nf, GEN x, GEN y)
 3107 {
 3108   pari_sp av = avma;
 3109   return gerepileupto(av, idealcoprimefact(nf, x, idealfactor(nf,y)));
 3110 }
 3111 
 3112 GEN
 3113 nfmulmodpr(GEN nf, GEN x, GEN y, GEN modpr)
 3114 {
 3115   pari_sp av = avma;
 3116   GEN z, p, pr = modpr, T;
 3117 
 3118   nf = checknf(nf); modpr = nf_to_Fq_init(nf,&pr,&T,&p);
 3119   x = nf_to_Fq(nf,x,modpr);
 3120   y = nf_to_Fq(nf,y,modpr);
 3121   z = Fq_mul(x,y,T,p);
 3122   return gerepileupto(av, algtobasis(nf, Fq_to_nf(z,modpr)));
 3123 }
 3124 
 3125 GEN
 3126 nfdivmodpr(GEN nf, GEN x, GEN y, GEN modpr)
 3127 {
 3128   pari_sp av = avma;
 3129   nf = checknf(nf);
 3130   return gerepileupto(av, nfreducemodpr(nf, nfdiv(nf,x,y), modpr));
 3131 }
 3132 
 3133 GEN
 3134 nfpowmodpr(GEN nf, GEN x, GEN k, GEN modpr)
 3135 {
 3136   pari_sp av=avma;
 3137   GEN z, T, p, pr = modpr;
 3138 
 3139   nf = checknf(nf); modpr = nf_to_Fq_init(nf,&pr,&T,&p);
 3140   z = nf_to_Fq(nf,x,modpr);
 3141   z = Fq_pow(z,k,T,p);
 3142   return gerepileupto(av, algtobasis(nf, Fq_to_nf(z,modpr)));
 3143 }
 3144 
 3145 GEN
 3146 nfkermodpr(GEN nf, GEN x, GEN modpr)
 3147 {
 3148   pari_sp av = avma;
 3149   GEN T, p, pr = modpr;
 3150 
 3151   nf = checknf(nf); modpr = nf_to_Fq_init(nf, &pr,&T,&p);
 3152   if (typ(x)!=t_MAT) pari_err_TYPE("nfkermodpr",x);
 3153   x = nfM_to_FqM(x, nf, modpr);
 3154   return gerepilecopy(av, FqM_to_nfM(FqM_ker(x,T,p), modpr));
 3155 }
 3156 
 3157 GEN
 3158 nfsolvemodpr(GEN nf, GEN a, GEN b, GEN pr)
 3159 {
 3160   const char *f = "nfsolvemodpr";
 3161   pari_sp av = avma;
 3162   GEN T, p, modpr;
 3163 
 3164   nf = checknf(nf);
 3165   modpr = nf_to_Fq_init(nf, &pr,&T,&p);
 3166   if (typ(a)!=t_MAT) pari_err_TYPE(f,a);
 3167   a = nfM_to_FqM(a, nf, modpr);
 3168   switch(typ(b))
 3169   {
 3170     case t_MAT:
 3171       b = nfM_to_FqM(b, nf, modpr);
 3172       b = FqM_gauss(a,b,T,p);
 3173       if (!b) pari_err_INV(f,a);
 3174       a = FqM_to_nfM(b, modpr);
 3175       break;
 3176     case t_COL:
 3177       b = nfV_to_FqV(b, nf, modpr);
 3178       b = FqM_FqC_gauss(a,b,T,p);
 3179       if (!b) pari_err_INV(f,a);
 3180       a = FqV_to_nfV(b, modpr);
 3181       break;
 3182     default: pari_err_TYPE(f,b);
 3183   }
 3184   return gerepilecopy(av, a);
 3185 }