Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - base4.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 19825-b77c7f8) Lines: 1302 1437 90.6 %
Date: 2016-12-05 05:49:04 Functions: 128 139 92.1 %
Legend: Lines: hit not hit

          Line data    Source code
       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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /*******************************************************************/
      15             : /*                                                                 */
      16             : /*                       BASIC NF OPERATIONS                       */
      17             : /*                           (continued)                           */
      18             : /*                                                                 */
      19             : /*******************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : /*******************************************************************/
      24             : /*                                                                 */
      25             : /*                     IDEAL OPERATIONS                            */
      26             : /*                                                                 */
      27             : /*******************************************************************/
      28             : 
      29             : /* A valid ideal is either principal (valid nf_element), or prime, or a matrix
      30             :  * on the integer basis in HNF.
      31             :  * A prime ideal is of the form [p,a,e,f,b], where the ideal is p.Z_K+a.Z_K,
      32             :  * p is a rational prime, a belongs to Z_K, e=e(P/p), f=f(P/p), and b
      33             :  * is Lenstra's constant, such that p.P^(-1)= p Z_K + b Z_K.
      34             :  *
      35             :  * An extended ideal is a couple [I,F] where I is a valid ideal and F is
      36             :  * either an algebraic number, or a factorization matrix attached to an
      37             :  * algebraic number. All routines work with either extended ideals or ideals
      38             :  * (an omitted F is assumed to be [;] <-> 1).
      39             :  * All ideals are output in HNF form. */
      40             : 
      41             : /* types and conversions */
      42             : 
      43             : long
      44     2489598 : idealtyp(GEN *ideal, GEN *arch)
      45             : {
      46     2489598 :   GEN x = *ideal;
      47     2489598 :   long t,lx,tx = typ(x);
      48             : 
      49     2489598 :   if (tx==t_VEC && lg(x)==3)
      50      337270 :   { *arch = gel(x,2); x = gel(x,1); tx = typ(x); }
      51             :   else
      52     2152328 :     *arch = NULL;
      53     2489598 :   switch(tx)
      54             :   {
      55     1140989 :     case t_MAT: lx = lg(x);
      56     1140989 :       if (lx == 1) { t = id_PRINCIPAL; x = gen_0; break; }
      57     1140912 :       if (lx != lgcols(x)) pari_err_TYPE("idealtyp [non-square t_MAT]",x);
      58     1140905 :       t = id_MAT;
      59     1140905 :       break;
      60             : 
      61     1124860 :     case t_VEC: if (lg(x)!=6) pari_err_TYPE("idealtyp",x);
      62     1124846 :       t = id_PRIME; break;
      63             : 
      64             :     case t_POL: case t_POLMOD: case t_COL:
      65             :     case t_INT: case t_FRAC:
      66      223749 :       t = id_PRINCIPAL; break;
      67             :     default:
      68           0 :       pari_err_TYPE("idealtyp",x);
      69           0 :       return 0; /*not reached*/
      70             :   }
      71     2489577 :   *ideal = x; return t;
      72             : }
      73             : 
      74             : /* nf a true nf; v = [a,x,...], a in Z. Return (a,x) */
      75             : GEN
      76     1237328 : idealhnf_two(GEN nf, GEN v)
      77             : {
      78     1237328 :   GEN p = gel(v,1), pi = gel(v,2), m = zk_scalar_or_multable(nf, pi);
      79     1237328 :   if (typ(m) == t_INT) return scalarmat(gcdii(m,p), nf_get_degree(nf));
      80      973841 :   return ZM_hnfmodid(m, p);
      81             : }
      82             : 
      83             : static GEN
      84       33499 : ZM_Q_mul(GEN x, GEN y)
      85       33499 : { return typ(y) == t_INT? ZM_Z_mul(x,y): RgM_Rg_mul(x,y); }
      86             : 
      87             : 
      88             : GEN
      89      180829 : idealhnf_principal(GEN nf, GEN x)
      90             : {
      91             :   GEN cx;
      92      180829 :   x = nf_to_scalar_or_basis(nf, x);
      93      180829 :   switch(typ(x))
      94             :   {
      95      108165 :     case t_COL: break;
      96       64230 :     case t_INT:  if (!signe(x)) return cgetg(1,t_MAT);
      97       64132 :       return scalarmat(absi(x), nf_get_degree(nf));
      98             :     case t_FRAC:
      99        8434 :       return scalarmat(Q_abs_shallow(x), nf_get_degree(nf));
     100           0 :     default: pari_err_TYPE("idealhnf",x);
     101             :   }
     102      108165 :   x = Q_primitive_part(x, &cx);
     103      108165 :   RgV_check_ZV(x, "idealhnf");
     104      108158 :   x = zk_multable(nf, x);
     105      108158 :   x = ZM_hnfmodid(x, zkmultable_capZ(x));
     106      108158 :   return cx? ZM_Q_mul(x,cx): x;
     107             : }
     108             : 
     109             : /* x integral ideal in t_MAT form, nx columns */
     110             : static GEN
     111           7 : vec_mulid(GEN nf, GEN x, long nx, long N)
     112             : {
     113           7 :   GEN m = cgetg(nx*N + 1, t_MAT);
     114             :   long i, j, k;
     115          21 :   for (i=k=1; i<=nx; i++)
     116          14 :     for (j=1; j<=N; j++) gel(m, k++) = zk_ei_mul(nf, gel(x,i),j);
     117           7 :   return m;
     118             : }
     119             : GEN
     120      285556 : idealhnf_shallow(GEN nf, GEN x)
     121             : {
     122      285556 :   long tx = typ(x), lx = lg(x), N;
     123             : 
     124             :   /* cannot use idealtyp because here we allow non-square matrices */
     125      285556 :   if (tx == t_VEC && lx == 3) { x = gel(x,1); tx = typ(x); lx = lg(x); }
     126      285556 :   if (tx == t_VEC && lx == 6) return idealhnf_two(nf,x); /* PRIME */
     127      144144 :   switch(tx)
     128             :   {
     129             :     case t_MAT:
     130             :     {
     131             :       GEN cx;
     132        7098 :       long nx = lx-1;
     133        7098 :       N = nf_get_degree(nf);
     134        7098 :       if (nx == 0) return cgetg(1, t_MAT);
     135        7077 :       if (nbrows(x) != N) pari_err_TYPE("idealhnf [wrong dimension]",x);
     136        7070 :       if (nx == 1) return idealhnf_principal(nf, gel(x,1));
     137             : 
     138        5775 :       if (nx == N && RgM_is_ZM(x) && ZM_ishnf(x)) return x;
     139         980 :       x = Q_primitive_part(x, &cx);
     140         980 :       if (nx < N) x = vec_mulid(nf, x, nx, N);
     141         980 :       x = ZM_hnfmod(x, ZM_detmult(x));
     142         980 :       return cx? ZM_Q_mul(x,cx): x;
     143             :     }
     144             :     case t_QFI:
     145             :     case t_QFR:
     146             :     {
     147          14 :       pari_sp av = avma;
     148          14 :       GEN u, D = nf_get_disc(nf), T = nf_get_pol(nf), f = nf_get_index(nf);
     149          14 :       GEN A = gel(x,1), B = gel(x,2);
     150          14 :       N = nf_get_degree(nf);
     151          14 :       if (N != 2)
     152           0 :         pari_err_TYPE("idealhnf [Qfb for non-quadratic fields]", x);
     153          14 :       if (!equalii(qfb_disc(x), D))
     154           7 :         pari_err_DOMAIN("idealhnf [Qfb]", "disc(q)", "!=", D, x);
     155             :       /* x -> A Z + (-B + sqrt(D)) / 2 Z
     156             :          K = Q[t]/T(t), t^2 + ut + v = 0,  u^2 - 4v = Df^2
     157             :          => t = (-u + sqrt(D) f)/2
     158             :          => sqrt(D)/2 = (t + u/2)/f */
     159           7 :       u = gel(T,3);
     160           7 :       B = deg1pol_shallow(ginv(f),
     161             :                           gsub(gdiv(u, shifti(f,1)), gdiv(B,gen_2)),
     162           7 :                           varn(T));
     163           7 :       return gerepileupto(av, idealhnf_two(nf, mkvec2(A,B)));
     164             :     }
     165      137032 :     default: return idealhnf_principal(nf, x); /* PRINCIPAL */
     166             :   }
     167             : }
     168             : GEN
     169        2590 : idealhnf(GEN nf, GEN x)
     170             : {
     171        2590 :   pari_sp av = avma;
     172        2590 :   GEN y = idealhnf_shallow(checknf(nf), x);
     173        2576 :   return (avma == av)? gcopy(y): gerepileupto(av, y);
     174             : }
     175             : 
     176             : /* GP functions */
     177             : 
     178             : GEN
     179          63 : idealtwoelt0(GEN nf, GEN x, GEN a)
     180             : {
     181          63 :   if (!a) return idealtwoelt(nf,x);
     182          42 :   return idealtwoelt2(nf,x,a);
     183             : }
     184             : 
     185             : GEN
     186          42 : idealpow0(GEN nf, GEN x, GEN n, long flag)
     187             : {
     188          42 :   if (flag) return idealpowred(nf,x,n);
     189          35 :   return idealpow(nf,x,n);
     190             : }
     191             : 
     192             : GEN
     193          56 : idealmul0(GEN nf, GEN x, GEN y, long flag)
     194             : {
     195          56 :   if (flag) return idealmulred(nf,x,y);
     196          49 :   return idealmul(nf,x,y);
     197             : }
     198             : 
     199             : GEN
     200          42 : idealdiv0(GEN nf, GEN x, GEN y, long flag)
     201             : {
     202          42 :   switch(flag)
     203             :   {
     204          21 :     case 0: return idealdiv(nf,x,y);
     205          21 :     case 1: return idealdivexact(nf,x,y);
     206           0 :     default: pari_err_FLAG("idealdiv");
     207             :   }
     208           0 :   return NULL; /* not reached */
     209             : }
     210             : 
     211             : GEN
     212          70 : idealaddtoone0(GEN nf, GEN arg1, GEN arg2)
     213             : {
     214          70 :   if (!arg2) return idealaddmultoone(nf,arg1);
     215          35 :   return idealaddtoone(nf,arg1,arg2);
     216             : }
     217             : 
     218             : /* b not a scalar */
     219             : static GEN
     220          28 : hnf_Z_ZC(GEN nf, GEN a, GEN b) { return hnfmodid(zk_multable(nf,b), a); }
     221             : /* b not a scalar */
     222             : static GEN
     223          21 : hnf_Z_QC(GEN nf, GEN a, GEN b)
     224             : {
     225             :   GEN db;
     226          21 :   b = Q_remove_denom(b, &db);
     227          21 :   if (db) a = mulii(a, db);
     228          21 :   b = hnf_Z_ZC(nf,a,b);
     229          21 :   return db? RgM_Rg_div(b, db): b;
     230             : }
     231             : /* b not a scalar (not point in trying to optimize for this case) */
     232             : static GEN
     233          28 : hnf_Q_QC(GEN nf, GEN a, GEN b)
     234             : {
     235             :   GEN da, db;
     236          28 :   if (typ(a) == t_INT) return hnf_Z_QC(nf, a, b);
     237           7 :   da = gel(a,2);
     238           7 :   a = gel(a,1);
     239           7 :   b = Q_remove_denom(b, &db);
     240             :   /* write da = d*A, db = d*B, gcd(A,B) = 1
     241             :    * gcd(a/(d A), b/(d B)) = gcd(a B, A b) / A B d = gcd(a B, b) / A B d */
     242           7 :   if (db)
     243             :   {
     244           7 :     GEN d = gcdii(da,db);
     245           7 :     if (!is_pm1(d)) db = diviiexact(db,d); /* B */
     246           7 :     if (!is_pm1(db))
     247             :     {
     248           7 :       a = mulii(a, db); /* a B */
     249           7 :       da = mulii(da, db); /* A B d = lcm(denom(a),denom(b)) */
     250             :     }
     251             :   }
     252           7 :   return RgM_Rg_div(hnf_Z_ZC(nf,a,b), da);
     253             : }
     254             : static GEN
     255           7 : hnf_QC_QC(GEN nf, GEN a, GEN b)
     256             : {
     257             :   GEN da, db, d, x;
     258           7 :   a = Q_remove_denom(a, &da);
     259           7 :   b = Q_remove_denom(b, &db);
     260           7 :   if (da) b = ZC_Z_mul(b, da);
     261           7 :   if (db) a = ZC_Z_mul(a, db);
     262           7 :   d = mul_denom(da, db);
     263           7 :   a = zk_multable(nf,a); da = zkmultable_capZ(a);
     264           7 :   b = zk_multable(nf,b); db = zkmultable_capZ(b);
     265           7 :   x = ZM_hnfmodid(shallowconcat(a,b), gcdii(da,db));
     266           7 :   return d? RgM_Rg_div(x, d): x;
     267             : }
     268             : static GEN
     269          21 : hnf_Q_Q(GEN nf, GEN a, GEN b) {return scalarmat(Q_gcd(a,b), nf_get_degree(nf));}
     270             : GEN
     271         119 : idealhnf0(GEN nf, GEN a, GEN b)
     272             : {
     273             :   long ta, tb;
     274             :   pari_sp av;
     275             :   GEN x;
     276         119 :   if (!b) return idealhnf(nf,a);
     277             : 
     278             :   /* HNF of aZ_K+bZ_K */
     279          56 :   av = avma; nf = checknf(nf);
     280          56 :   a = nf_to_scalar_or_basis(nf,a); ta = typ(a);
     281          56 :   b = nf_to_scalar_or_basis(nf,b); tb = typ(b);
     282          56 :   if (ta == t_COL)
     283          14 :     x = (tb==t_COL)? hnf_QC_QC(nf, a,b): hnf_Q_QC(nf, b,a);
     284             :   else
     285          42 :     x = (tb==t_COL)? hnf_Q_QC(nf, a,b): hnf_Q_Q(nf, a,b);
     286          56 :   return gerepileupto(av, x);
     287             : }
     288             : 
     289             : /*******************************************************************/
     290             : /*                                                                 */
     291             : /*                       TWO-ELEMENT FORM                          */
     292             : /*                                                                 */
     293             : /*******************************************************************/
     294             : static GEN idealapprfact_i(GEN nf, GEN x, int nored);
     295             : 
     296             : static int
     297      185156 : ok_elt(GEN x, GEN xZ, GEN y)
     298             : {
     299      185156 :   pari_sp av = avma;
     300      185156 :   int r = ZM_equal(x, ZM_hnfmodid(y, xZ));
     301      185156 :   avma = av; return r;
     302             : }
     303             : 
     304             : static GEN
     305       58298 : addmul_col(GEN a, long s, GEN b)
     306             : {
     307             :   long i,l;
     308       58298 :   if (!s) return a? leafcopy(a): a;
     309       58150 :   if (!a) return gmulsg(s,b);
     310       54903 :   l = lg(a);
     311      292165 :   for (i=1; i<l; i++)
     312      237262 :     if (signe(gel(b,i))) gel(a,i) = addii(gel(a,i), mulsi(s, gel(b,i)));
     313       54903 :   return a;
     314             : }
     315             : 
     316             : /* a <-- a + s * b, all coeffs integers */
     317             : static GEN
     318       25003 : addmul_mat(GEN a, long s, GEN b)
     319             : {
     320             :   long j,l;
     321             :   /* copy otherwise next call corrupts a */
     322       25003 :   if (!s) return a? RgM_shallowcopy(a): a;
     323       23404 :   if (!a) return gmulsg(s,b);
     324       12931 :   l = lg(a);
     325       62943 :   for (j=1; j<l; j++)
     326       50012 :     (void)addmul_col(gel(a,j), s, gel(b,j));
     327       12931 :   return a;
     328             : }
     329             : 
     330             : static GEN
     331      117063 : get_random_a(GEN nf, GEN x, GEN xZ)
     332             : {
     333             :   pari_sp av;
     334      117063 :   long i, lm, l = lg(x);
     335             :   GEN a, z, beta, mul;
     336             : 
     337      117063 :   beta= cgetg(l, t_VEC);
     338      117063 :   mul = cgetg(l, t_VEC); lm = 1; /* = lg(mul) */
     339             :   /* look for a in x such that a O/xZ = x O/xZ */
     340      229082 :   for (i = 2; i < l; i++)
     341             :   {
     342      225835 :     GEN xi = gel(x,i);
     343      225835 :     GEN t = FpM_red(zk_multable(nf,xi), xZ); /* ZM, cannot be a scalar */
     344      225835 :     if (gequal0(t)) continue;
     345      174683 :     if (ok_elt(x,xZ, t)) return xi;
     346       60867 :     gel(beta,lm) = xi;
     347             :     /* mul[i] = { canonical generators for x[i] O/xZ as Z-module } */
     348       60867 :     gel(mul,lm) = t; lm++;
     349             :   }
     350        3247 :   setlg(mul, lm);
     351        3247 :   setlg(beta,lm);
     352        3247 :   z = cgetg(lm, t_VECSMALL);
     353       10515 :   for(av = avma;; avma = av)
     354             :   {
     355       35518 :     for (a=NULL,i=1; i<lm; i++)
     356             :     {
     357       25003 :       long t = random_bits(4) - 7; /* in [-7,8] */
     358       25003 :       z[i] = t;
     359       25003 :       a = addmul_mat(a, t, gel(mul,i));
     360             :     }
     361             :     /* a = matrix (NOT HNF) of ideal generated by beta.z in O/xZ */
     362       10515 :     if (a && ok_elt(x,xZ, a)) break;
     363        7268 :   }
     364       11533 :   for (a=NULL,i=1; i<lm; i++)
     365        8286 :     a = addmul_col(a, z[i], gel(beta,i));
     366        3247 :   return a;
     367             : }
     368             : 
     369             : /* x square matrix, assume it is HNF */
     370             : static GEN
     371      321064 : mat_ideal_two_elt(GEN nf, GEN x)
     372             : {
     373             :   GEN y, a, cx, xZ;
     374      321064 :   long N = nf_get_degree(nf);
     375             :   pari_sp av, tetpil;
     376             : 
     377      321064 :   if (lg(x)-1 != N) pari_err_DIM("idealtwoelt");
     378      321050 :   if (N == 2) return mkvec2copy(gcoeff(x,1,1), gel(x,2));
     379             : 
     380      126572 :   y = cgetg(3,t_VEC); av = avma;
     381      126572 :   cx = Q_content(x);
     382      126572 :   xZ = gcoeff(x,1,1);
     383      126572 :   if (gequal(xZ, cx)) /* x = (cx) */
     384             :   {
     385        3479 :     gel(y,1) = cx;
     386        3479 :     gel(y,2) = gen_0; return y;
     387             :   }
     388      123093 :   if (equali1(cx)) cx = NULL;
     389             :   else
     390             :   {
     391         413 :     x = Q_div_to_int(x, cx);
     392         413 :     xZ = gcoeff(x,1,1);
     393             :   }
     394      123093 :   if (N < 6)
     395      112664 :     a = get_random_a(nf, x, xZ);
     396             :   else
     397             :   {
     398       10429 :     const long FB[] = { _evallg(15+1) | evaltyp(t_VECSMALL),
     399             :       2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47
     400             :     };
     401       10429 :     GEN P, E, a1 = Z_smoothen(xZ, (GEN)FB, &P, &E);
     402       10429 :     if (!a1) /* factors completely */
     403        6030 :       a = idealapprfact_i(nf, idealfactor(nf,x), 1);
     404        4399 :     else if (lg(P) == 1) /* no small factors */
     405        3184 :       a = get_random_a(nf, x, xZ);
     406             :     else /* general case */
     407             :     {
     408             :       GEN A0, A1, a0, u0, u1, v0, v1, pi0, pi1, t, u;
     409        1215 :       a0 = diviiexact(xZ, a1);
     410        1215 :       A0 = ZM_hnfmodid(x, a0); /* smooth part of x */
     411        1215 :       A1 = ZM_hnfmodid(x, a1); /* cofactor */
     412        1215 :       pi0 = idealapprfact_i(nf, idealfactor(nf,A0), 1);
     413        1215 :       pi1 = get_random_a(nf, A1, a1);
     414        1215 :       (void)bezout(a0, a1, &v0,&v1);
     415        1215 :       u0 = mulii(a0, v0);
     416        1215 :       u1 = mulii(a1, v1);
     417        1215 :       if (typ(pi0) != t_COL) t = addmulii(u0, pi0, u1);
     418             :       else
     419        1215 :       { t = ZC_Z_mul(pi0, u1); gel(t,1) = addii(gel(t,1), u0); }
     420        1215 :       u = ZC_Z_mul(pi1, u0); gel(u,1) = addii(gel(u,1), u1);
     421        1215 :       a = nfmuli(nf, centermod(u, xZ), centermod(t, xZ));
     422             :     }
     423             :   }
     424      123093 :   if (cx)
     425             :   {
     426         413 :     a = centermod(a, xZ);
     427         413 :     tetpil = avma;
     428         413 :     if (typ(cx) == t_INT)
     429             :     {
     430         371 :       gel(y,1) = mulii(xZ, cx);
     431         371 :       gel(y,2) = ZC_Z_mul(a, cx);
     432             :     }
     433             :     else
     434             :     {
     435          42 :       gel(y,1) = gmul(xZ, cx);
     436          42 :       gel(y,2) = RgC_Rg_mul(a, cx);
     437             :     }
     438             :   }
     439             :   else
     440             :   {
     441      122680 :     tetpil = avma;
     442      122680 :     gel(y,1) = icopy(xZ);
     443      122680 :     gel(y,2) = centermod(a, xZ);
     444             :   }
     445      123093 :   gerepilecoeffssp(av,tetpil,y+1,2); return y;
     446             : }
     447             : 
     448             : /* Given an ideal x, returns [a,alpha] such that a is in Q,
     449             :  * x = a Z_K + alpha Z_K, alpha in K^*
     450             :  * a = 0 or alpha = 0 are possible, but do not try to determine whether
     451             :  * x is principal. */
     452             : GEN
     453       24395 : idealtwoelt(GEN nf, GEN x)
     454             : {
     455             :   pari_sp av;
     456             :   GEN z;
     457       24395 :   long tx = idealtyp(&x,&z);
     458       24388 :   nf = checknf(nf);
     459       24388 :   if (tx == id_MAT) return mat_ideal_two_elt(nf,x);
     460        1218 :   if (tx == id_PRIME) return mkvec2copy(gel(x,1), gel(x,2));
     461             :   /* id_PRINCIPAL */
     462         455 :   av = avma; x = nf_to_scalar_or_basis(nf, x);
     463         721 :   return gerepilecopy(av, typ(x)==t_COL? mkvec2(gen_0,x):
     464         350 :                                          mkvec2(Q_abs_shallow(x),gen_0));
     465             : }
     466             : 
     467             : /*******************************************************************/
     468             : /*                                                                 */
     469             : /*                         FACTORIZATION                           */
     470             : /*                                                                 */
     471             : /*******************************************************************/
     472             : /* x integral ideal in HNF, Zval = v_p(x \cap Z) > 0; return v_p(Nx) */
     473             : static long
     474       99016 : idealHNF_norm_pval(GEN x, GEN p, long Zval)
     475             : {
     476       99016 :   long i, v = Zval, l = lg(x);
     477       99016 :   for (i = 2; i < l; i++) v += Z_pval(gcoeff(x,i,i), p);
     478       99016 :   return v;
     479             : }
     480             : 
     481             : /* return P, primes dividing Nx and xZ = x\cap Z, set v_p(Nx), v_p(xZ);
     482             :  * x integral in HNF */
     483             : GEN
     484       15246 : idealHNF_Z_factor(GEN x, GEN *pvN, GEN *pvZ)
     485             : {
     486       15246 :   GEN xZ = gcoeff(x,1,1), f, P, E, vN, vZ;
     487             :   long i, l;
     488       15246 :   if (typ(xZ) != t_INT) pari_err_TYPE("idealfactor",x);
     489       15246 :   f = Z_factor(xZ);
     490       15246 :   P = gel(f,1); l = lg(P);
     491       15246 :   E = gel(f,2);
     492       15246 :   *pvN = vN = cgetg(l, t_VECSMALL);
     493       15246 :   *pvZ = vZ = cgetg(l, t_VECSMALL);
     494       28342 :   for (i = 1; i < l; i++)
     495             :   {
     496       13096 :     vZ[i] = itou(gel(E,i));
     497       13096 :     vN[i] = idealHNF_norm_pval(x,gel(P,i), vZ[i]);
     498             :   }
     499       15246 :   return P;
     500             : }
     501             : 
     502             : /* v_P(A)*f(P) <= Nval [e.g. Nval = v_p(Norm A)], Zval = v_p(A \cap Z).
     503             :  * Return v_P(A) */
     504             : static long
     505      113218 : idealHNF_val(GEN A, GEN P, long Nval, long Zval)
     506             : {
     507      113218 :   long f = pr_get_f(P), vmax, v, e, i, j, k, l;
     508             :   GEN mul, B, a, y, r, p, pk, cx, vals;
     509             :   pari_sp av;
     510             : 
     511      113218 :   if (Nval < f) return 0;
     512      113148 :   p = pr_get_p(P);
     513      113148 :   e = pr_get_e(P);
     514             :   /* v_P(A) <= max [ e * v_p(A \cap Z), floor[v_p(Nix) / f ] */
     515      113148 :   vmax = minss(Zval * e, Nval / f);
     516      113148 :   mul = pr_get_tau(P);
     517      113148 :   l = lg(mul);
     518      113148 :   B = cgetg(l,t_MAT);
     519             :   /* B[1] not needed: v_pr(A[1]) = v_pr(A \cap Z) is known already */
     520      113148 :   gel(B,1) = gen_0; /* dummy */
     521      463800 :   for (j = 2; j < l; j++)
     522             :   {
     523      393976 :     GEN x = gel(A,j);
     524      393976 :     gel(B,j) = y = cgetg(l, t_COL);
     525     3628790 :     for (i = 1; i < l; i++)
     526             :     { /* compute a = (x.t0)_i, A in HNF ==> x[j+1..l-1] = 0 */
     527     3278138 :       a = mulii(gel(x,1), gcoeff(mul,i,1));
     528     3278138 :       for (k = 2; k <= j; k++) a = addii(a, mulii(gel(x,k), gcoeff(mul,i,k)));
     529             :       /* p | a ? */
     530     3278138 :       gel(y,i) = dvmdii(a,p,&r); if (signe(r)) return 0;
     531             :     }
     532             :   }
     533       69824 :   vals = cgetg(l, t_VECSMALL);
     534             :   /* vals[1] not needed */
     535      356803 :   for (j = 2; j < l; j++)
     536             :   {
     537      286979 :     gel(B,j) = Q_primitive_part(gel(B,j), &cx);
     538      286979 :     vals[j] = cx? 1 + e * Q_pval(cx, p): 1;
     539             :   }
     540       69824 :   pk = powiu(p, ceildivuu(vmax, e));
     541       69824 :   av = avma; y = cgetg(l,t_COL);
     542             :   /* can compute mod p^ceil((vmax-v)/e) */
     543      119901 :   for (v = 1; v < vmax; v++)
     544             :   { /* we know v_pr(Bj) >= v for all j */
     545       52335 :     if (e == 1 || (vmax - v) % e == 0) pk = diviiexact(pk, p);
     546      420912 :     for (j = 2; j < l; j++)
     547             :     {
     548      370835 :       GEN x = gel(B,j); if (v < vals[j]) continue;
     549     3746069 :       for (i = 1; i < l; i++)
     550             :       {
     551     3473922 :         pari_sp av2 = avma;
     552     3473922 :         a = mulii(gel(x,1), gcoeff(mul,i,1));
     553     3473922 :         for (k = 2; k < l; k++) a = addii(a, mulii(gel(x,k), gcoeff(mul,i,k)));
     554             :         /* a = (x.t_0)_i; p | a ? */
     555     3473922 :         a = dvmdii(a,p,&r); if (signe(r)) return v;
     556     3471664 :         if (lgefint(a) > lgefint(pk)) a = remii(a, pk);
     557     3471664 :         gel(y,i) = gerepileuptoint(av2, a);
     558             :       }
     559      272147 :       gel(B,j) = y; y = x;
     560      272147 :       if (gc_needed(av,3))
     561             :       {
     562           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"idealval");
     563           0 :         gerepileall(av,3, &y,&B,&pk);
     564             :       }
     565             :     }
     566             :   }
     567       67566 :   return v;
     568             : }
     569             : /* true nf, x integral ideal */
     570             : static GEN
     571       15246 : idealHNF_factor(GEN nf, GEN x)
     572             : {
     573       15246 :   const long N = lg(x)-1;
     574             :   long i, j, k, l, v;
     575             :   GEN vp, vN, vZ, vP, vE, cx;
     576             : 
     577       15246 :   x = Q_primitive_part(x, &cx);
     578       15246 :   vp = idealHNF_Z_factor(x, &vN,&vZ);
     579       15246 :   l = lg(vp);
     580       15246 :   i = cx? expi(cx)+1: 1;
     581       15246 :   vP = cgetg((l+i-2)*N+1, t_COL);
     582       15246 :   vE = cgetg((l+i-2)*N+1, t_COL);
     583       28342 :   for (i = k = 1; i < l; i++)
     584             :   {
     585       13096 :     GEN L, p = gel(vp,i);
     586       13096 :     long Nval = vN[i], Zval = vZ[i], vc = cx? Z_pvalrem(cx,p,&cx): 0;
     587       13096 :     if (vc)
     588             :     {
     589         385 :       L = idealprimedec(nf,p);
     590         385 :       if (is_pm1(cx)) cx = NULL;
     591             :     }
     592             :     else
     593       12711 :       L = idealprimedec_limit_f(nf,p,Nval);
     594       27305 :     for (j = 1; j < lg(L); j++)
     595             :     {
     596       27298 :       GEN P = gel(L,j);
     597       27298 :       pari_sp av = avma;
     598       27298 :       v = idealHNF_val(x, P, Nval, Zval);
     599       27298 :       avma = av;
     600       27298 :       Nval -= v*pr_get_f(P);
     601       27298 :       v += vc * pr_get_e(P); if (!v) continue;
     602       16400 :       gel(vP,k) = P;
     603       16400 :       gel(vE,k) = utoipos(v); k++;
     604       16400 :       if (!Nval) break; /* now only the content contributes */
     605             :     }
     606       13118 :     if (vc) for (j++; j<lg(L); j++)
     607             :     {
     608          22 :       GEN P = gel(L,j);
     609          22 :       gel(vP,k) = P;
     610          22 :       gel(vE,k) = utoipos(vc * pr_get_e(P)); k++;
     611             :     }
     612             :   }
     613       15246 :   if (cx)
     614             :   {
     615        2744 :     GEN f = Z_factor(cx), cP = gel(f,1), cE = gel(f,2);
     616        2744 :     long lc = lg(cP);
     617        5775 :     for (i=1; i<lc; i++)
     618             :     {
     619        3031 :       GEN p = gel(cP,i), L = idealprimedec(nf,p);
     620        3031 :       long vc = itos(gel(cE,i));
     621        6825 :       for (j=1; j<lg(L); j++)
     622             :       {
     623        3794 :         GEN P = gel(L,j);
     624        3794 :         gel(vP,k) = P;
     625        3794 :         gel(vE,k) = utoipos(vc * pr_get_e(P)); k++;
     626             :       }
     627             :     }
     628             :   }
     629       15246 :   setlg(vP, k);
     630       15246 :   setlg(vE, k); return mkmat2(vP, vE);
     631             : }
     632             : /* c * vector(#L,i,L[i].e), assume results fit in ulong */
     633             : static GEN
     634        2898 : prV_e_muls(GEN L, long c)
     635             : {
     636        2898 :   long j, l = lg(L);
     637        2898 :   GEN z = cgetg(l, t_COL);
     638        2898 :   for (j = 1; j < l; j++) gel(z,j) = stoi(c * pr_get_e(gel(L,j)));
     639        2898 :   return z;
     640             : }
     641             : /* true nf, y in Q */
     642             : static GEN
     643        2898 : Q_nffactor(GEN nf, GEN y)
     644             : {
     645             :   GEN f, P, E;
     646             :   long lfa, i;
     647        2898 :   if (typ(y) == t_INT)
     648             :   {
     649        2884 :     if (!signe(y)) pari_err_DOMAIN("idealfactor", "ideal", "=",gen_0,y);
     650        2863 :     if (is_pm1(y)) return trivial_fact();
     651             :   }
     652        2163 :   f = factor(Q_abs_shallow(y));
     653        2163 :   P = gel(f,1); lfa = lg(P);
     654        2163 :   E = gel(f,2);
     655        5061 :   for (i = 1; i < lfa; i++)
     656             :   {
     657        2898 :     gel(P,i) = idealprimedec(nf, gel(P,i));
     658        2898 :     gel(E,i) = prV_e_muls(gel(P,i), itos(gel(E,i)));
     659             :   }
     660        2163 :   settyp(P,t_VEC); P = shallowconcat1(P);
     661        2163 :   settyp(E,t_VEC); E = shallowconcat1(E);
     662        2163 :   gel(f,1) = P; settyp(P, t_COL);
     663        2163 :   gel(f,2) = E; return f;
     664             : }
     665             : 
     666             : GEN
     667       18165 : idealfactor(GEN nf, GEN x)
     668             : {
     669       18165 :   pari_sp av = avma;
     670             :   GEN fa, y;
     671       18165 :   long tx = idealtyp(&x,&y);
     672             : 
     673       18165 :   nf = checknf(nf);
     674       18165 :   if (tx == id_PRIME) retmkmat2(mkcolcopy(x), mkcol(gen_1));
     675       18137 :   if (tx == id_PRINCIPAL)
     676             :   {
     677        4382 :     y = nf_to_scalar_or_basis(nf, x);
     678        4382 :     if (typ(y) != t_COL) return gerepilecopy(av, Q_nffactor(nf, y));
     679             :   }
     680       15239 :   y = idealnumden(nf, x);
     681       15239 :   fa = idealHNF_factor(nf, gel(y,1));
     682       15239 :   if (!isint1(gel(y,2)))
     683             :   {
     684           7 :     GEN F = idealHNF_factor(nf, gel(y,2));
     685           7 :     fa = famat_mul_shallow(fa, famat_inv_shallow(F));
     686             :   }
     687       15239 :   fa = gerepilecopy(av, fa);
     688       15239 :   return sort_factor(fa, (void*)&cmp_prime_ideal, &cmp_nodata);
     689             : }
     690             : 
     691             : /* P prime ideal in idealprimedec format. Return valuation(A) at P */
     692             : long
     693      310236 : idealval(GEN nf, GEN A, GEN P)
     694             : {
     695      310236 :   pari_sp av = avma;
     696             :   GEN a, p, cA;
     697      310236 :   long vcA, v, Zval, tx = idealtyp(&A,&a);
     698             : 
     699      310236 :   if (tx == id_PRINCIPAL) return nfval(nf,A,P);
     700      309074 :   checkprid(P);
     701      309074 :   if (tx == id_PRIME) return pr_equal(P, A)? 1: 0;
     702             :   /* id_MAT */
     703      309046 :   nf = checknf(nf);
     704      309046 :   A = Q_primitive_part(A, &cA);
     705      309046 :   p = pr_get_p(P);
     706      309046 :   vcA = cA? Q_pval(cA,p): 0;
     707      309046 :   if (pr_is_inert(P)) { avma = av; return vcA; }
     708      307660 :   Zval = Z_pval(gcoeff(A,1,1), p);
     709      307660 :   if (!Zval) v = 0;
     710             :   else
     711             :   {
     712       85920 :     long Nval = idealHNF_norm_pval(A, p, Zval);
     713       85920 :     v = idealHNF_val(A, P, Nval, Zval);
     714             :   }
     715      307660 :   avma = av; return vcA? v + vcA*pr_get_e(P): v;
     716             : }
     717             : GEN
     718          42 : gpidealval(GEN nf, GEN ix, GEN P)
     719             : {
     720          42 :   long v = idealval(nf,ix,P);
     721          42 :   return v == LONG_MAX? mkoo(): stoi(v);
     722             : }
     723             : 
     724             : /* gcd and generalized Bezout */
     725             : 
     726             : GEN
     727       28224 : idealadd(GEN nf, GEN x, GEN y)
     728             : {
     729       28224 :   pari_sp av = avma;
     730             :   long tx, ty;
     731             :   GEN z, a, dx, dy, dz;
     732             : 
     733       28224 :   tx = idealtyp(&x,&z);
     734       28224 :   ty = idealtyp(&y,&z); nf = checknf(nf);
     735       28224 :   if (tx != id_MAT) x = idealhnf_shallow(nf,x);
     736       28224 :   if (ty != id_MAT) y = idealhnf_shallow(nf,y);
     737       28224 :   if (lg(x) == 1) return gerepilecopy(av,y);
     738       28224 :   if (lg(y) == 1) return gerepilecopy(av,x); /* check for 0 ideal */
     739       28224 :   dx = Q_denom(x);
     740       28224 :   dy = Q_denom(y); dz = lcmii(dx,dy);
     741       28224 :   if (is_pm1(dz)) dz = NULL; else {
     742        5670 :     x = Q_muli_to_int(x, dz);
     743        5670 :     y = Q_muli_to_int(y, dz);
     744             :   }
     745       28224 :   a = gcdii(gcoeff(x,1,1), gcoeff(y,1,1));
     746       28224 :   if (is_pm1(a))
     747             :   {
     748       12578 :     long N = lg(x)-1;
     749       12578 :     if (!dz) { avma = av; return matid(N); }
     750         692 :     return gerepileupto(av, scalarmat(ginv(dz), N));
     751             :   }
     752       15646 :   z = ZM_hnfmodid(shallowconcat(x,y), a);
     753       15646 :   if (dz) z = RgM_Rg_div(z,dz);
     754       15646 :   return gerepileupto(av,z);
     755             : }
     756             : 
     757             : static GEN
     758          28 : trivial_merge(GEN x)
     759          28 : { return (lg(x) == 1 || !is_pm1(gcoeff(x,1,1)))? NULL: gen_1; }
     760             : GEN
     761      152607 : idealaddtoone_i(GEN nf, GEN x, GEN y)
     762             : {
     763             :   GEN a;
     764      152607 :   long tx = idealtyp(&x, &a/*junk*/);
     765      152607 :   long ty = idealtyp(&y, &a/*junk*/);
     766      152607 :   if (tx != id_MAT) x = idealhnf_shallow(nf, x);
     767      152607 :   if (ty != id_MAT) y = idealhnf_shallow(nf, y);
     768      152607 :   if (lg(x) == 1)
     769          14 :     a = trivial_merge(y);
     770      152593 :   else if (lg(y) == 1)
     771          14 :     a = trivial_merge(x);
     772             :   else {
     773      152579 :     a = hnfmerge_get_1(x, y);
     774      152579 :     if (a && typ(a) == t_COL) a = ZC_reducemodlll(a, idealHNF_mul(nf,x,y));
     775             :   }
     776      152607 :   if (!a) pari_err_COPRIME("idealaddtoone",x,y);
     777      152593 :   return a;
     778             : }
     779             : 
     780             : GEN
     781        2842 : idealaddtoone(GEN nf, GEN x, GEN y)
     782             : {
     783        2842 :   GEN z = cgetg(3,t_VEC), a;
     784        2842 :   pari_sp av = avma;
     785        2842 :   nf = checknf(nf);
     786        2842 :   a = gerepileupto(av, idealaddtoone_i(nf,x,y));
     787        2828 :   gel(z,1) = a;
     788        2828 :   gel(z,2) = typ(a) == t_COL? Z_ZC_sub(gen_1,a): subui(1,a);
     789        2828 :   return z;
     790             : }
     791             : 
     792             : /* assume elements of list are integral ideals */
     793             : GEN
     794          35 : idealaddmultoone(GEN nf, GEN list)
     795             : {
     796          35 :   pari_sp av = avma;
     797          35 :   long N, i, l, nz, tx = typ(list);
     798             :   GEN H, U, perm, L;
     799             : 
     800          35 :   nf = checknf(nf); N = nf_get_degree(nf);
     801          35 :   if (!is_vec_t(tx)) pari_err_TYPE("idealaddmultoone",list);
     802          35 :   l = lg(list);
     803          35 :   L = cgetg(l, t_VEC);
     804          35 :   if (l == 1)
     805           0 :     pari_err_DOMAIN("idealaddmultoone", "sum(ideals)", "!=", gen_1, L);
     806          35 :   nz = 0; /* number of non-zero ideals in L */
     807          98 :   for (i=1; i<l; i++)
     808             :   {
     809          70 :     GEN I = gel(list,i);
     810          70 :     if (typ(I) != t_MAT) I = idealhnf_shallow(nf,I);
     811          70 :     if (lg(I) != 1)
     812             :     {
     813          42 :       nz++; RgM_check_ZM(I,"idealaddmultoone");
     814          35 :       if (lgcols(I) != N+1) pari_err_TYPE("idealaddmultoone [not an ideal]", I);
     815             :     }
     816          63 :     gel(L,i) = I;
     817             :   }
     818          28 :   H = ZM_hnfperm(shallowconcat1(L), &U, &perm);
     819          28 :   if (lg(H) == 1 || !equali1(gcoeff(H,1,1)))
     820           7 :     pari_err_DOMAIN("idealaddmultoone", "sum(ideals)", "!=", gen_1, L);
     821          49 :   for (i=1; i<=N; i++)
     822          49 :     if (perm[i] == 1) break;
     823          21 :   U = gel(U,(nz-1)*N + i); /* (L[1]|...|L[nz]) U = 1 */
     824          21 :   nz = 0;
     825          63 :   for (i=1; i<l; i++)
     826             :   {
     827          42 :     GEN c = gel(L,i);
     828          42 :     if (lg(c) == 1)
     829          14 :       c = gen_0;
     830             :     else {
     831          28 :       c = ZM_ZC_mul(c, vecslice(U, nz*N + 1, (nz+1)*N));
     832          28 :       nz++;
     833             :     }
     834          42 :     gel(L,i) = c;
     835             :   }
     836          21 :   return gerepilecopy(av, L);
     837             : }
     838             : 
     839             : /* multiplication */
     840             : 
     841             : /* x integral ideal (without archimedean component) in HNF form
     842             :  * y = [a,alpha] corresponds to the integral ideal aZ_K+alpha Z_K, a in Z,
     843             :  * alpha a ZV or a ZM (multiplication table). Multiply them */
     844             : static GEN
     845      741449 : idealHNF_mul_two(GEN nf, GEN x, GEN y)
     846             : {
     847      741449 :   GEN m, a = gel(y,1), alpha = gel(y,2);
     848             :   long i, N;
     849             : 
     850      741449 :   if (typ(alpha) != t_MAT)
     851             :   {
     852      563134 :     alpha = zk_scalar_or_multable(nf, alpha);
     853      563134 :     if (typ(alpha) == t_INT) /* e.g. y inert ? 0 should not (but may) occur */
     854        3220 :       return signe(a)? ZM_Z_mul(x, gcdii(a, alpha)): cgetg(1,t_MAT);
     855             :   }
     856      738229 :   N = lg(x)-1; m = cgetg((N<<1)+1,t_MAT);
     857      738229 :   for (i=1; i<=N; i++) gel(m,i)   = ZM_ZC_mul(alpha,gel(x,i));
     858      738229 :   for (i=1; i<=N; i++) gel(m,i+N) = ZC_Z_mul(gel(x,i), a);
     859      738229 :   return ZM_hnfmodid(m, mulii(a, gcoeff(x,1,1)));
     860             : }
     861             : 
     862             : /* Assume ix and iy are integral in HNF form [NOT extended]. Not memory clean.
     863             :  * HACK: ideal in iy can be of the form [a,b], a in Z, b in Z_K */
     864             : GEN
     865      417716 : idealHNF_mul(GEN nf, GEN x, GEN y)
     866             : {
     867             :   GEN z;
     868      417716 :   if (typ(y) == t_VEC)
     869      192295 :     z = idealHNF_mul_two(nf,x,y);
     870             :   else
     871             :   { /* reduce one ideal to two-elt form. The smallest */
     872      225421 :     GEN xZ = gcoeff(x,1,1), yZ = gcoeff(y,1,1);
     873      225421 :     if (cmpii(xZ, yZ) < 0)
     874             :     {
     875       30531 :       if (is_pm1(xZ)) return gcopy(y);
     876       23551 :       z = idealHNF_mul_two(nf, y, mat_ideal_two_elt(nf,x));
     877             :     }
     878             :     else
     879             :     {
     880      194890 :       if (is_pm1(yZ)) return gcopy(x);
     881      180442 :       z = idealHNF_mul_two(nf, x, mat_ideal_two_elt(nf,y));
     882             :     }
     883             :   }
     884      396288 :   return z;
     885             : }
     886             : 
     887             : /* operations on elements in factored form */
     888             : 
     889             : GEN
     890       86300 : famat_mul_shallow(GEN f, GEN g)
     891             : {
     892       86300 :   if (typ(f) != t_MAT) f = to_famat_shallow(f,gen_1);
     893       86300 :   if (typ(g) != t_MAT) g = to_famat_shallow(g,gen_1);
     894       86300 :   if (lg(f) == 1) return g;
     895       72662 :   if (lg(g) == 1) return f;
     896      142148 :   return mkmat2(shallowconcat(gel(f,1), gel(g,1)),
     897      142148 :                 shallowconcat(gel(f,2), gel(g,2)));
     898             : }
     899             : GEN
     900       58975 : famat_mulpow_shallow(GEN f, GEN g, GEN e)
     901             : {
     902       58975 :   if (!signe(e)) return f;
     903       57190 :   return famat_mul_shallow(f, famat_pow_shallow(g, e));
     904             : }
     905             : 
     906             : GEN
     907           0 : to_famat(GEN x, GEN y) { retmkmat2(mkcolcopy(x), mkcolcopy(y)); }
     908             : GEN
     909      783056 : to_famat_shallow(GEN x, GEN y) { return mkmat2(mkcol(x), mkcol(y)); }
     910             : 
     911             : /* concat the single elt x; not gconcat since x may be a t_COL */
     912             : static GEN
     913       62906 : append(GEN v, GEN x)
     914             : {
     915       62906 :   long i, l = lg(v);
     916       62906 :   GEN w = cgetg(l+1, typ(v));
     917       62906 :   for (i=1; i<l; i++) gel(w,i) = gcopy(gel(v,i));
     918       62906 :   gel(w,i) = gcopy(x); return w;
     919             : }
     920             : /* add x^1 to famat f */
     921             : static GEN
     922       87261 : famat_add(GEN f, GEN x)
     923             : {
     924       87261 :   GEN h = cgetg(3,t_MAT);
     925       87261 :   if (lg(f) == 1)
     926             :   {
     927       24355 :     gel(h,1) = mkcolcopy(x);
     928       24355 :     gel(h,2) = mkcol(gen_1);
     929             :   }
     930             :   else
     931             :   {
     932       62906 :     gel(h,1) = append(gel(f,1), x);
     933       62906 :     gel(h,2) = gconcat(gel(f,2), gen_1);
     934             :   }
     935       87261 :   return h;
     936             : }
     937             : 
     938             : GEN
     939      113729 : famat_mul(GEN f, GEN g)
     940             : {
     941             :   GEN h;
     942      113729 :   if (typ(g) != t_MAT) {
     943       87261 :     if (typ(f) == t_MAT) return famat_add(f, g);
     944           0 :     h = cgetg(3, t_MAT);
     945           0 :     gel(h,1) = mkcol2(gcopy(f), gcopy(g));
     946           0 :     gel(h,2) = mkcol2(gen_1, gen_1);
     947             :   }
     948       26468 :   if (typ(f) != t_MAT) return famat_add(g, f);
     949       26468 :   if (lg(f) == 1) return gcopy(g);
     950        4214 :   if (lg(g) == 1) return gcopy(f);
     951        1890 :   h = cgetg(3,t_MAT);
     952        1890 :   gel(h,1) = gconcat(gel(f,1), gel(g,1));
     953        1890 :   gel(h,2) = gconcat(gel(f,2), gel(g,2));
     954        1890 :   return h;
     955             : }
     956             : 
     957             : GEN
     958       55486 : famat_sqr(GEN f)
     959             : {
     960             :   GEN h;
     961       55486 :   if (lg(f) == 1) return cgetg(1,t_MAT);
     962       26452 :   if (typ(f) != t_MAT) return to_famat(f,gen_2);
     963       26452 :   h = cgetg(3,t_MAT);
     964       26452 :   gel(h,1) = gcopy(gel(f,1));
     965       26452 :   gel(h,2) = gmul2n(gel(f,2),1);
     966       26452 :   return h;
     967             : }
     968             : 
     969             : GEN
     970       25109 : famat_inv_shallow(GEN f)
     971             : {
     972       25109 :   if (lg(f) == 1) return f;
     973       25109 :   if (typ(f) != t_MAT) return to_famat_shallow(f,gen_m1);
     974          14 :   return mkmat2(gel(f,1), ZC_neg(gel(f,2)));
     975             : }
     976             : GEN
     977        4616 : famat_inv(GEN f)
     978             : {
     979        4616 :   if (lg(f) == 1) return cgetg(1,t_MAT);
     980        2532 :   if (typ(f) != t_MAT) return to_famat(f,gen_m1);
     981        2532 :   retmkmat2(gcopy(gel(f,1)), ZC_neg(gel(f,2)));
     982             : }
     983             : GEN
     984         421 : famat_pow(GEN f, GEN n)
     985             : {
     986         421 :   if (lg(f) == 1) return cgetg(1,t_MAT);
     987           0 :   if (typ(f) != t_MAT) return to_famat(f,n);
     988           0 :   retmkmat2(gcopy(gel(f,1)), ZC_Z_mul(gel(f,2),n));
     989             : }
     990             : GEN
     991       57190 : famat_pow_shallow(GEN f, GEN n)
     992             : {
     993       57190 :   if (is_pm1(n)) return signe(n) > 0? f: famat_inv_shallow(f);
     994       29064 :   if (lg(f) == 1) return f;
     995       29064 :   if (typ(f) != t_MAT) return to_famat_shallow(f,n);
     996        1260 :   return mkmat2(gel(f,1), ZC_Z_mul(gel(f,2),n));
     997             : }
     998             : 
     999             : GEN
    1000           0 : famat_Z_gcd(GEN M, GEN n)
    1001             : {
    1002           0 :   pari_sp av=avma;
    1003           0 :   long i, j, l=lgcols(M);
    1004           0 :   GEN F=cgetg(3,t_MAT);
    1005           0 :   gel(F,1)=cgetg(l,t_COL);
    1006           0 :   gel(F,2)=cgetg(l,t_COL);
    1007           0 :   for (i=1, j=1; i<l; i++)
    1008             :   {
    1009           0 :     GEN p = gcoeff(M,i,1);
    1010           0 :     GEN e = gminsg(Z_pval(n,p),gcoeff(M,i,2));
    1011           0 :     if (signe(e))
    1012             :     {
    1013           0 :       gcoeff(F,j,1)=p;
    1014           0 :       gcoeff(F,j,2)=e;
    1015           0 :       j++;
    1016             :     }
    1017             :   }
    1018           0 :   setlg(gel(F,1),j); setlg(gel(F,2),j);
    1019           0 :   return gerepilecopy(av,F);
    1020             : }
    1021             : 
    1022             : /* x assumed to be a t_MATs (factorization matrix), or compatible with
    1023             :  * the element_* functions. */
    1024             : static GEN
    1025       66049 : ext_sqr(GEN nf, GEN x)
    1026       66049 : { return (typ(x)==t_MAT)? famat_sqr(x): nfsqr(nf, x); }
    1027             : static GEN
    1028      166866 : ext_mul(GEN nf, GEN x, GEN y)
    1029      166866 : { return (typ(x)==t_MAT)? famat_mul(x,y): nfmul(nf, x, y); }
    1030             : static GEN
    1031        4476 : ext_inv(GEN nf, GEN x)
    1032        4476 : { return (typ(x)==t_MAT)? famat_inv(x): nfinv(nf, x); }
    1033             : static GEN
    1034         421 : ext_pow(GEN nf, GEN x, GEN n)
    1035         421 : { return (typ(x)==t_MAT)? famat_pow(x,n): nfpow(nf, x, n); }
    1036             : 
    1037             : GEN
    1038           0 : famat_to_nf(GEN nf, GEN f)
    1039             : {
    1040             :   GEN t, x, e;
    1041             :   long i;
    1042           0 :   if (lg(f) == 1) return gen_1;
    1043             : 
    1044           0 :   x = gel(f,1);
    1045           0 :   e = gel(f,2);
    1046           0 :   t = nfpow(nf, gel(x,1), gel(e,1));
    1047           0 :   for (i=lg(x)-1; i>1; i--)
    1048           0 :     t = nfmul(nf, t, nfpow(nf, gel(x,i), gel(e,i)));
    1049           0 :   return t;
    1050             : }
    1051             : 
    1052             : GEN
    1053        2569 : famat_reduce(GEN fa)
    1054             : {
    1055             :   GEN E, G, L, g, e;
    1056             :   long i, k, l;
    1057             : 
    1058        2569 :   if (lg(fa) == 1) return fa;
    1059        2513 :   g = gel(fa,1); l = lg(g);
    1060        2513 :   e = gel(fa,2);
    1061        2513 :   L = gen_indexsort(g, (void*)&cmp_universal, &cmp_nodata);
    1062        2513 :   G = cgetg(l, t_COL);
    1063        2513 :   E = cgetg(l, t_COL);
    1064             :   /* merge */
    1065       10976 :   for (k=i=1; i<l; i++,k++)
    1066             :   {
    1067        8463 :     gel(G,k) = gel(g,L[i]);
    1068        8463 :     gel(E,k) = gel(e,L[i]);
    1069        8463 :     if (k > 1 && gidentical(gel(G,k), gel(G,k-1)))
    1070             :     {
    1071         798 :       gel(E,k-1) = addii(gel(E,k), gel(E,k-1));
    1072         798 :       k--;
    1073             :     }
    1074             :   }
    1075             :   /* kill 0 exponents */
    1076        2513 :   l = k;
    1077       10178 :   for (k=i=1; i<l; i++)
    1078        7665 :     if (!gequal0(gel(E,i)))
    1079             :     {
    1080        7070 :       gel(G,k) = gel(G,i);
    1081        7070 :       gel(E,k) = gel(E,i); k++;
    1082             :     }
    1083        2513 :   setlg(G, k);
    1084        2513 :   setlg(E, k); return mkmat2(G,E);
    1085             : }
    1086             : 
    1087             : GEN
    1088       12690 : famatsmall_reduce(GEN fa)
    1089             : {
    1090             :   GEN E, G, L, g, e;
    1091             :   long i, k, l;
    1092       12690 :   if (lg(fa) == 1) return fa;
    1093       12690 :   g = gel(fa,1); l = lg(g);
    1094       12690 :   e = gel(fa,2);
    1095       12690 :   L = vecsmall_indexsort(g);
    1096       12691 :   G = cgetg(l, t_VECSMALL);
    1097       12691 :   E = cgetg(l, t_VECSMALL);
    1098             :   /* merge */
    1099      113620 :   for (k=i=1; i<l; i++,k++)
    1100             :   {
    1101      100929 :     G[k] = g[L[i]];
    1102      100929 :     E[k] = e[L[i]];
    1103      100929 :     if (k > 1 && G[k] == G[k-1])
    1104             :     {
    1105        5898 :       E[k-1] += E[k];
    1106        5898 :       k--;
    1107             :     }
    1108             :   }
    1109             :   /* kill 0 exponents */
    1110       12691 :   l = k;
    1111      107722 :   for (k=i=1; i<l; i++)
    1112       95031 :     if (E[i])
    1113             :     {
    1114       92100 :       G[k] = G[i];
    1115       92100 :       E[k] = E[i]; k++;
    1116             :     }
    1117       12691 :   setlg(G, k);
    1118       12691 :   setlg(E, k); return mkmat2(G,E);
    1119             : }
    1120             : 
    1121             : GEN
    1122       48307 : ZM_famat_limit(GEN fa, GEN limit)
    1123             : {
    1124             :   pari_sp av;
    1125             :   GEN E, G, g, e, r;
    1126             :   long i, k, l, n, lG;
    1127             : 
    1128       48307 :   if (lg(fa) == 1) return fa;
    1129       48307 :   g = gel(fa,1); l = lg(g);
    1130       48307 :   e = gel(fa,2);
    1131      108899 :   for(n=0, i=1; i<l; i++)
    1132       60592 :     if (cmpii(gel(g,i),limit)<=0) n++;
    1133       48307 :   lG = n<l-1 ? n+2 : n+1;
    1134       48307 :   G = cgetg(lG, t_COL);
    1135       48307 :   E = cgetg(lG, t_COL);
    1136       48307 :   av = avma;
    1137      108899 :   for (i=1, k=1, r = gen_1; i<l; i++)
    1138             :   {
    1139       60592 :     if (cmpii(gel(g,i),limit)<=0)
    1140             :     {
    1141       60508 :       gel(G,k) = gel(g,i);
    1142       60508 :       gel(E,k) = gel(e,i);
    1143       60508 :       k++;
    1144          84 :     } else r = mulii(r, powii(gel(g,i), gel(e,i)));
    1145             :   }
    1146       48307 :   if (k<i)
    1147             :   {
    1148          84 :     gel(G, k) = gerepileuptoint(av, r);
    1149          84 :     gel(E, k) = gen_1;
    1150             :   }
    1151       48307 :   return mkmat2(G,E);
    1152             : }
    1153             : 
    1154             : /* assume pr has degree 1 and coprime to Q_denom(x) */
    1155             : static GEN
    1156        5118 : to_Fp_coprime(GEN nf, GEN x, GEN modpr)
    1157             : {
    1158        5118 :   GEN d, r, p = modpr_get_p(modpr);
    1159        5118 :   x = nf_to_scalar_or_basis(nf,x);
    1160        5118 :   if (typ(x) != t_COL) return Rg_to_Fp(x,p);
    1161        4754 :   x = Q_remove_denom(x, &d);
    1162        4754 :   r = zk_to_Fq(x, modpr);
    1163        4754 :   if (d) r = Fp_div(r, d, p);
    1164        4754 :   return r;
    1165             : }
    1166             : 
    1167             : /* pr coprime to all denominators occurring in x */
    1168             : static GEN
    1169         789 : famat_to_Fp_coprime(GEN nf, GEN x, GEN modpr)
    1170             : {
    1171         789 :   GEN p = modpr_get_p(modpr);
    1172         789 :   GEN t = NULL, g = gel(x,1), e = gel(x,2), q = subiu(p,1);
    1173         789 :   long i, l = lg(g);
    1174        2428 :   for (i = 1; i < l; i++)
    1175             :   {
    1176        1639 :     GEN n = modii(gel(e,i), q);
    1177        1639 :     if (signe(n))
    1178             :     {
    1179        1639 :       GEN h = to_Fp_coprime(nf, gel(g,i), modpr);
    1180        1639 :       h = Fp_pow(h, n, p);
    1181        1639 :       t = t? Fp_mul(t, h, p): h;
    1182             :     }
    1183             :   }
    1184         789 :   return t? modii(t, p): gen_1;
    1185             : }
    1186             : 
    1187             : /* cf famat_to_nf_modideal_coprime, modpr attached to prime of degree 1 */
    1188             : GEN
    1189        4268 : nf_to_Fp_coprime(GEN nf, GEN x, GEN modpr)
    1190             : {
    1191        8536 :   return typ(x)==t_MAT? famat_to_Fp_coprime(nf, x, modpr)
    1192        4268 :                       : to_Fp_coprime(nf, x, modpr);
    1193             : }
    1194             : 
    1195             : static long
    1196      110405 : zk_pvalrem(GEN x, GEN p, GEN *py)
    1197      110405 : { return (typ(x) == t_INT)? Z_pvalrem(x, p, py): ZV_pvalrem(x, p, py); }
    1198             : /* x a QC or Q. Return a ZC or Z, whose content is coprime to Z. Set v, dx
    1199             :  * such that x = p^v (newx / dx); dx = NULL if 1 */
    1200             : static GEN
    1201      232548 : nf_remove_denom_p(GEN nf, GEN x, GEN p, GEN *pdx, long *pv)
    1202             : {
    1203             :   long vcx;
    1204             :   GEN dx;
    1205      232548 :   x = nf_to_scalar_or_basis(nf, x);
    1206      232548 :   x = Q_remove_denom(x, &dx);
    1207      232548 :   if (dx)
    1208             :   {
    1209      162764 :     vcx = - Z_pvalrem(dx, p, &dx);
    1210      162764 :     if (!vcx) vcx = zk_pvalrem(x, p, &x);
    1211      162764 :     if (isint1(dx)) dx = NULL;
    1212             :   }
    1213             :   else
    1214             :   {
    1215       69784 :     vcx = zk_pvalrem(x, p, &x);
    1216       69784 :     dx = NULL;
    1217             :   }
    1218      232548 :   *pv = vcx;
    1219      232548 :   *pdx = dx; return x;
    1220             : }
    1221             : /* x = b^e/p^(e-1) in Z_K; x = 0 mod p/pr^e, (x,pr) = 1. Return NULL
    1222             :  * if p inert (instead of 1) */
    1223             : static GEN
    1224       58233 : p_makecoprime(GEN pr)
    1225             : {
    1226       58233 :   GEN B = pr_get_tau(pr), b;
    1227             :   long i, e;
    1228             : 
    1229       58233 :   if (typ(B) == t_INT) return NULL;
    1230       58093 :   b = gel(B,1); /* B = multiplication table by b */
    1231       58093 :   e = pr_get_e(pr);
    1232       58093 :   if (e == 1) return b;
    1233             :   /* one could also divide (exactly) by p in each iteration */
    1234       16058 :   for (i = 1; i < e; i++) b = ZM_ZC_mul(B, b);
    1235       16058 :   return ZC_Z_divexact(b, powiu(pr_get_p(pr), e-1));
    1236             : }
    1237             : 
    1238             : /* Compute A = prod g[i]^e[i] mod pr^k, assuming (A, pr) = 1.
    1239             :  * Method: modify each g[i] so that it becomes coprime to pr,
    1240             :  * g[i] *= (b/p)^v_pr(g[i]), where b/p = pr^(-1) times something integral
    1241             :  * and prime to p; globally, we multiply by (b/p)^v_pr(A) = 1.
    1242             :  * Optimizations:
    1243             :  * 1) remove all powers of p from contents, and consider extra generator p^vp;
    1244             :  * modified as p * (b/p)^e = b^e / p^(e-1)
    1245             :  * 2) remove denominators, coprime to p, by multiplying by inverse mod prk\cap Z
    1246             :  *
    1247             :  * EX = multiple of exponent of (O_K / pr^k)^* used to reduce the product in
    1248             :  * case the e[i] are large */
    1249             : GEN
    1250       92502 : famat_makecoprime(GEN nf, GEN g, GEN e, GEN pr, GEN prk, GEN EX)
    1251             : {
    1252       92502 :   GEN G, E, t, vp = NULL, p = pr_get_p(pr), prkZ = gcoeff(prk, 1,1);
    1253       92502 :   long i, l = lg(g);
    1254             : 
    1255       92502 :   G = cgetg(l+1, t_VEC);
    1256       92502 :   E = cgetg(l+1, t_VEC); /* l+1: room for "modified p" */
    1257      325050 :   for (i=1; i < l; i++)
    1258             :   {
    1259             :     long vcx;
    1260      232548 :     GEN dx, x = nf_remove_denom_p(nf, gel(g,i), p, &dx, &vcx);
    1261      232548 :     if (vcx) /* = v_p(content(g[i])) */
    1262             :     {
    1263      122934 :       GEN a = mulsi(vcx, gel(e,i));
    1264      122934 :       vp = vp? addii(vp, a): a;
    1265             :     }
    1266             :     /* x integral, content coprime to p; dx coprime to p */
    1267      232548 :     if (typ(x) == t_INT)
    1268             :     { /* x coprime to p, hence to pr */
    1269       30216 :       x = modii(x, prkZ);
    1270       30216 :       if (dx) x = Fp_div(x, dx, prkZ);
    1271             :     }
    1272             :     else
    1273             :     {
    1274      202332 :       (void)ZC_nfvalrem(x, pr, &x); /* x *= (b/p)^v_pr(x) */
    1275      202332 :       x = ZC_hnfrem(FpC_red(x,prkZ), prk);
    1276      202332 :       if (dx) x = FpC_Fp_mul(x, Fp_inv(dx,prkZ), prkZ);
    1277             :     }
    1278      232548 :     gel(G,i) = x;
    1279      232548 :     gel(E,i) = gel(e,i);
    1280             :   }
    1281             : 
    1282       92502 :   t = vp? p_makecoprime(pr): NULL;
    1283       92502 :   if (!t)
    1284             :   { /* no need for extra generator */
    1285       34409 :     setlg(G,l);
    1286       34409 :     setlg(E,l);
    1287             :   }
    1288             :   else
    1289             :   {
    1290       58093 :     gel(G,i) = FpC_red(t, prkZ);
    1291       58093 :     gel(E,i) = vp;
    1292             :   }
    1293       92502 :   return famat_to_nf_modideal_coprime(nf, G, E, prk, EX);
    1294             : }
    1295             : 
    1296             : /* prod g[i]^e[i] mod bid, assume (g[i], id) = 1 */
    1297             : GEN
    1298        1008 : famat_to_nf_moddivisor(GEN nf, GEN g, GEN e, GEN bid)
    1299             : {
    1300             :   GEN t, cyc;
    1301        1008 :   if (lg(g) == 1) return gen_1;
    1302        1008 :   cyc = bid_get_cyc(bid);
    1303        1008 :   if (lg(cyc) == 1)
    1304           0 :     t = gen_1;
    1305             :   else
    1306        1008 :     t = famat_to_nf_modideal_coprime(nf, g, e, bid_get_ideal(bid), gel(cyc,1));
    1307        1008 :   return set_sign_mod_divisor(nf, mkmat2(g,e), t, bid_get_sarch(bid));
    1308             : }
    1309             : 
    1310             : GEN
    1311      178024 : vecmul(GEN x, GEN y)
    1312             : {
    1313      178024 :   long i,lx, tx = typ(x);
    1314             :   GEN z;
    1315      178024 :   if (is_scalar_t(tx)) return gmul(x,y);
    1316       15358 :   z = cgetg_copy(x, &lx);
    1317       15358 :   for (i=1; i<lx; i++) gel(z,i) = vecmul(gel(x,i), gel(y,i));
    1318       15358 :   return z;
    1319             : }
    1320             : 
    1321             : GEN
    1322           0 : vecinv(GEN x)
    1323             : {
    1324           0 :   long i,lx, tx = typ(x);
    1325             :   GEN z;
    1326           0 :   if (is_scalar_t(tx)) return ginv(x);
    1327           0 :   z = cgetg_copy(x, &lx);
    1328           0 :   for (i=1; i<lx; i++) gel(z,i) = vecinv(gel(x,i));
    1329           0 :   return z;
    1330             : }
    1331             : 
    1332             : GEN
    1333       15729 : vecpow(GEN x, GEN n)
    1334             : {
    1335       15729 :   long i,lx, tx = typ(x);
    1336             :   GEN z;
    1337       15729 :   if (is_scalar_t(tx)) return powgi(x,n);
    1338        4270 :   z = cgetg_copy(x, &lx);
    1339        4270 :   for (i=1; i<lx; i++) gel(z,i) = vecpow(gel(x,i), n);
    1340        4270 :   return z;
    1341             : }
    1342             : 
    1343             : GEN
    1344         903 : vecdiv(GEN x, GEN y)
    1345             : {
    1346         903 :   long i,lx, tx = typ(x);
    1347             :   GEN z;
    1348         903 :   if (is_scalar_t(tx)) return gdiv(x,y);
    1349         301 :   z = cgetg_copy(x, &lx);
    1350         301 :   for (i=1; i<lx; i++) gel(z,i) = vecdiv(gel(x,i), gel(y,i));
    1351         301 :   return z;
    1352             : }
    1353             : 
    1354             : /* A ideal as a square t_MAT */
    1355             : static GEN
    1356      129140 : idealmulelt(GEN nf, GEN x, GEN A)
    1357             : {
    1358             :   long i, lx;
    1359             :   GEN dx, dA, D;
    1360      129140 :   if (lg(A) == 1) return cgetg(1, t_MAT);
    1361      129140 :   x = nf_to_scalar_or_basis(nf,x);
    1362      129140 :   if (typ(x) != t_COL)
    1363       18038 :     return isintzero(x)? cgetg(1,t_MAT): RgM_Rg_mul(A, Q_abs_shallow(x));
    1364      111102 :   x = Q_remove_denom(x, &dx);
    1365      111102 :   A = Q_remove_denom(A, &dA);
    1366      111102 :   x = zk_multable(nf, x);
    1367      111102 :   D = mulii(zkmultable_capZ(x), gcoeff(A,1,1));
    1368      111102 :   x = zkC_multable_mul(A, x);
    1369      111102 :   settyp(x, t_MAT); lx = lg(x);
    1370             :   /* x may contain scalars (at most 1 since the ideal is non-0)*/
    1371      399926 :   for (i=1; i<lx; i++)
    1372      290839 :     if (typ(gel(x,i)) == t_INT)
    1373             :     {
    1374        2015 :       if (i > 1) swap(gel(x,1), gel(x,i)); /* help HNF */
    1375        2015 :       gel(x,1) = scalarcol_shallow(gel(x,1), lx-1);
    1376        2015 :       break;
    1377             :     }
    1378      111102 :   x = ZM_hnfmodid(x, D);
    1379      111102 :   dx = mul_denom(dx,dA);
    1380      111102 :   return dx? gdiv(x,dx): x;
    1381             : }
    1382             : 
    1383             : /* nf a true nf, tx <= ty */
    1384             : static GEN
    1385      533048 : idealmul_aux(GEN nf, GEN x, GEN y, long tx, long ty)
    1386             : {
    1387             :   GEN z, cx, cy;
    1388      533048 :   switch(tx)
    1389             :   {
    1390             :     case id_PRINCIPAL:
    1391      153779 :       switch(ty)
    1392             :       {
    1393             :         case id_PRINCIPAL:
    1394       24513 :           return idealhnf_principal(nf, nfmul(nf,x,y));
    1395             :         case id_PRIME:
    1396             :         {
    1397         126 :           GEN p = gel(y,1), pi = gel(y,2), cx;
    1398         126 :           if (pr_is_inert(y)) return RgM_Rg_mul(idealhnf_principal(nf,x),p);
    1399             : 
    1400          42 :           x = nf_to_scalar_or_basis(nf, x);
    1401          42 :           switch(typ(x))
    1402             :           {
    1403             :             case t_INT:
    1404          28 :               if (!signe(x)) return cgetg(1,t_MAT);
    1405          28 :               return ZM_Z_mul(idealhnf_two(nf,y), absi(x));
    1406             :             case t_FRAC:
    1407           7 :               return RgM_Rg_mul(idealhnf_two(nf,y), Q_abs_shallow(x));
    1408             :           }
    1409             :           /* t_COL */
    1410           7 :           x = Q_primitive_part(x, &cx);
    1411           7 :           x = zk_multable(nf, x);
    1412           7 :           z = shallowconcat(ZM_Z_mul(x,p), ZM_ZC_mul(x,pi));
    1413           7 :           z = ZM_hnfmodid(z, mulii(p, zkmultable_capZ(x)));
    1414           7 :           return cx? ZM_Q_mul(z, cx): z;
    1415             :         }
    1416             :         default: /* id_MAT */
    1417      129140 :           return idealmulelt(nf, x,y);
    1418             :       }
    1419             :     case id_PRIME:
    1420      325190 :       if (ty==id_PRIME)
    1421      300421 :       { y = idealhnf_two(nf,y); cy = NULL; }
    1422             :       else
    1423       24769 :         y = Q_primitive_part(y, &cy);
    1424      325190 :       y = idealHNF_mul_two(nf,y,x);
    1425      325190 :       return cy? RgM_Rg_mul(y,cy): y;
    1426             : 
    1427             :     default: /* id_MAT */
    1428             :     {
    1429       54079 :       long N = nf_get_degree(nf);
    1430       54079 :       if (lg(x)-1 != N || lg(y)-1 != N) pari_err_DIM("idealmul");
    1431       54065 :       x = Q_primitive_part(x, &cx);
    1432       54065 :       y = Q_primitive_part(y, &cy); cx = mul_content(cx,cy);
    1433       54065 :       y = idealHNF_mul(nf,x,y);
    1434       54065 :       return cx? ZM_Q_mul(y,cx): y;
    1435             :     }
    1436             :   }
    1437             : }
    1438             : 
    1439             : /* output the ideal product ix.iy */
    1440             : GEN
    1441      533048 : idealmul(GEN nf, GEN x, GEN y)
    1442             : {
    1443             :   pari_sp av;
    1444             :   GEN res, ax, ay, z;
    1445      533048 :   long tx = idealtyp(&x,&ax);
    1446      533048 :   long ty = idealtyp(&y,&ay), f;
    1447      533048 :   if (tx>ty) { swap(ax,ay); swap(x,y); lswap(tx,ty); }
    1448      533048 :   f = (ax||ay); res = f? cgetg(3,t_VEC): NULL; /*product is an extended ideal*/
    1449      533048 :   av = avma;
    1450      533048 :   z = gerepileupto(av, idealmul_aux(checknf(nf), x,y, tx,ty));
    1451      533027 :   if (!f) return z;
    1452       45293 :   if (ax && ay)
    1453       43366 :     ax = ext_mul(nf, ax, ay);
    1454             :   else
    1455        1927 :     ax = gcopy(ax? ax: ay);
    1456       45293 :   gel(res,1) = z; gel(res,2) = ax; return res;
    1457             : }
    1458             : 
    1459             : /* Return x, integral in 2-elt form, such that pr^2 = c * x. cf idealpowprime
    1460             :  * nf = true nf */
    1461             : static GEN
    1462       37795 : idealsqrprime(GEN nf, GEN pr, GEN *pc)
    1463             : {
    1464       37795 :   GEN p = pr_get_p(pr), q, gen;
    1465       37795 :   long e = pr_get_e(pr), f = pr_get_f(pr);
    1466             : 
    1467       37795 :   q = (e == 1)? sqri(p): p;
    1468       37795 :   if (e <= 2 && e * f == nf_get_degree(nf))
    1469             :   { /* pr^e = (p) */
    1470        7189 :     *pc = q;
    1471        7189 :     return mkvec2(gen_1,gen_0);
    1472             :   }
    1473       30606 :   gen = nfsqr(nf, pr_get_gen(pr));
    1474       30606 :   gen = FpC_red(gen, q);
    1475       30606 :   *pc = NULL;
    1476       30606 :   return mkvec2(q, gen);
    1477             : }
    1478             : /* cf idealpow_aux */
    1479             : static GEN
    1480       66476 : idealsqr_aux(GEN nf, GEN x, long tx)
    1481             : {
    1482       66476 :   GEN T = nf_get_pol(nf), m, cx, a, alpha;
    1483       66476 :   long N = degpol(T);
    1484       66476 :   switch(tx)
    1485             :   {
    1486             :     case id_PRINCIPAL:
    1487          63 :       return idealhnf_principal(nf, nfsqr(nf,x));
    1488             :     case id_PRIME:
    1489       25364 :       if (pr_is_inert(x)) return scalarmat(sqri(gel(x,1)), N);
    1490       25196 :       x = idealsqrprime(nf, x, &cx);
    1491       25196 :       x = idealhnf_two(nf,x);
    1492       25196 :       return cx? ZM_Z_mul(x, cx): x;
    1493             :     default:
    1494       41049 :       x = Q_primitive_part(x, &cx);
    1495       41049 :       a = mat_ideal_two_elt(nf,x); alpha = gel(a,2); a = gel(a,1);
    1496       41049 :       alpha = nfsqr(nf,alpha);
    1497       41049 :       m = zk_scalar_or_multable(nf, alpha);
    1498       41049 :       if (typ(m) == t_INT) {
    1499        1449 :         x = gcdii(sqri(a), m);
    1500        1449 :         if (cx) x = gmul(x, gsqr(cx));
    1501        1449 :         x = scalarmat(x, N);
    1502             :       }
    1503             :       else
    1504             :       {
    1505       39600 :         x = ZM_hnfmodid(m, gcdii(sqri(a), zkmultable_capZ(m)));
    1506       39600 :         if (cx) cx = gsqr(cx);
    1507       39600 :         if (cx) x = RgM_Rg_mul(x, cx);
    1508             :       }
    1509       41049 :       return x;
    1510             :   }
    1511             : }
    1512             : GEN
    1513       66476 : idealsqr(GEN nf, GEN x)
    1514             : {
    1515             :   pari_sp av;
    1516             :   GEN res, ax, z;
    1517       66476 :   long tx = idealtyp(&x,&ax);
    1518       66476 :   res = ax? cgetg(3,t_VEC): NULL; /*product is an extended ideal*/
    1519       66476 :   av = avma;
    1520       66476 :   z = gerepileupto(av, idealsqr_aux(checknf(nf), x, tx));
    1521       66476 :   if (!ax) return z;
    1522       66049 :   gel(res,1) = z;
    1523       66049 :   gel(res,2) = ext_sqr(nf, ax); return res;
    1524             : }
    1525             : 
    1526             : /* norm of an ideal */
    1527             : GEN
    1528        5439 : idealnorm(GEN nf, GEN x)
    1529             : {
    1530             :   pari_sp av;
    1531             :   GEN y, T;
    1532             :   long tx;
    1533             : 
    1534        5439 :   switch(idealtyp(&x,&y))
    1535             :   {
    1536         175 :     case id_PRIME: return pr_norm(x);
    1537        3514 :     case id_MAT: return RgM_det_triangular(x);
    1538             :   }
    1539             :   /* id_PRINCIPAL */
    1540        1750 :   nf = checknf(nf); T = nf_get_pol(nf); av = avma;
    1541        1750 :   x = nf_to_scalar_or_alg(nf, x);
    1542        1750 :   x = (typ(x) == t_POL)? RgXQ_norm(x, T): gpowgs(x, degpol(T));
    1543        1750 :   tx = typ(x);
    1544        1750 :   if (tx == t_INT) return gerepileuptoint(av, absi(x));
    1545         385 :   if (tx != t_FRAC) pari_err_TYPE("idealnorm",x);
    1546         385 :   return gerepileupto(av, Q_abs(x));
    1547             : }
    1548             : 
    1549             : /* I^(-1) = { x \in K, Tr(x D^(-1) I) \in Z }, D different of K/Q
    1550             :  *
    1551             :  * nf[5][6] = pp( D^(-1) ) = pp( HNF( T^(-1) ) ), T = (Tr(wi wj))
    1552             :  * nf[5][7] = same in 2-elt form.
    1553             :  * Assume I integral. Return the integral ideal (I\cap Z) I^(-1) */
    1554             : GEN
    1555      162618 : idealHNF_inv_Z(GEN nf, GEN I)
    1556             : {
    1557      162618 :   GEN J, dual, IZ = gcoeff(I,1,1); /* I \cap Z */
    1558      162618 :   if (isint1(IZ)) return matid(lg(I)-1);
    1559      158985 :   J = idealHNF_mul(nf,I, gmael(nf,5,7));
    1560             :  /* I in HNF, hence easily inverted; multiply by IZ to get integer coeffs
    1561             :   * missing content cancels while solving the linear equation */
    1562      158985 :   dual = shallowtrans( hnf_divscale(J, gmael(nf,5,6), IZ) );
    1563      158985 :   return ZM_hnfmodid(dual, IZ);
    1564             : }
    1565             : /* I HNF with rational coefficients (denominator d). */
    1566             : GEN
    1567       31140 : idealHNF_inv(GEN nf, GEN I)
    1568             : {
    1569       31140 :   GEN J, IQ = gcoeff(I,1,1); /* I \cap Q; d IQ = dI \cap Z */
    1570       31140 :   J = idealHNF_inv_Z(nf, Q_remove_denom(I, NULL)); /* = (dI)^(-1) * (d IQ) */
    1571       31140 :   return equali1(IQ)? J: RgM_Rg_div(J, IQ);
    1572             : }
    1573             : 
    1574             : /* return p * P^(-1)  [integral] */
    1575             : GEN
    1576       39355 : pr_inv_p(GEN pr)
    1577             : {
    1578       39355 :   if (pr_is_inert(pr)) return matid(pr_get_f(pr));
    1579       39075 :   return ZM_hnfmodid(pr_get_tau(pr), pr_get_p(pr));
    1580             : }
    1581             : GEN
    1582        1379 : pr_inv(GEN pr)
    1583             : {
    1584        1379 :   GEN p = pr_get_p(pr);
    1585        1379 :   if (pr_is_inert(pr)) return scalarmat(ginv(p), pr_get_f(pr));
    1586        1204 :   return RgM_Rg_div(ZM_hnfmodid(pr_get_tau(pr),p), p);
    1587             : }
    1588             : 
    1589             : GEN
    1590       46245 : idealinv(GEN nf, GEN x)
    1591             : {
    1592             :   GEN res, ax;
    1593             :   pari_sp av;
    1594       46245 :   long tx = idealtyp(&x,&ax), N;
    1595             : 
    1596       46245 :   res = ax? cgetg(3,t_VEC): NULL;
    1597       46245 :   nf = checknf(nf); av = avma;
    1598       46245 :   N = nf_get_degree(nf);
    1599       46245 :   switch (tx)
    1600             :   {
    1601             :     case id_MAT:
    1602       27374 :       if (lg(x)-1 != N) pari_err_DIM("idealinv");
    1603       27374 :       x = idealHNF_inv(nf,x); break;
    1604             :     case id_PRINCIPAL:
    1605       17814 :       x = nf_to_scalar_or_basis(nf, x);
    1606       17814 :       if (typ(x) != t_COL)
    1607       17772 :         x = idealhnf_principal(nf,ginv(x));
    1608             :       else
    1609             :       { /* nfinv + idealhnf where we already know (x) \cap Z */
    1610             :         GEN c, d;
    1611          42 :         x = Q_remove_denom(x, &c);
    1612          42 :         x = zk_inv(nf, x);
    1613          42 :         x = Q_remove_denom(x, &d); /* true inverse is c/d * x */
    1614          42 :         if (!d) /* x and x^(-1) integral => x a unit */
    1615           7 :           x = scalarmat_shallow(c? c: gen_1, N);
    1616             :         else
    1617             :         {
    1618          35 :           c = c? gdiv(c,d): ginv(d);
    1619          35 :           x = zk_multable(nf, x);
    1620          35 :           x = ZM_Q_mul(ZM_hnfmodid(x,d), c);
    1621             :         }
    1622             :       }
    1623       17814 :       break;
    1624             :     case id_PRIME:
    1625        1057 :       x = pr_inv(x); break;
    1626             :   }
    1627       46245 :   x = gerepileupto(av,x); if (!ax) return x;
    1628        4476 :   gel(res,1) = x;
    1629        4476 :   gel(res,2) = ext_inv(nf, ax); return res;
    1630             : }
    1631             : 
    1632             : /* write x = A/B, A,B coprime integral ideals */
    1633             : GEN
    1634       15372 : idealnumden(GEN nf, GEN x)
    1635             : {
    1636       15372 :   pari_sp av = avma;
    1637             :   GEN x0, ax, c, d, A, B, J;
    1638       15372 :   long tx = idealtyp(&x,&ax);
    1639       15372 :   nf = checknf(nf);
    1640       15372 :   switch (tx)
    1641             :   {
    1642             :     case id_PRIME:
    1643           7 :       retmkvec2(idealhnf(nf, x), gen_1);
    1644             :     case id_PRINCIPAL:
    1645             :     {
    1646             :       GEN xZ, mx;
    1647        1603 :       x = nf_to_scalar_or_basis(nf, x);
    1648        1603 :       switch(typ(x))
    1649             :       {
    1650          56 :         case t_INT: return gerepilecopy(av, mkvec2(absi(x),gen_1));
    1651          14 :         case t_FRAC:return gerepilecopy(av, mkvec2(absi(gel(x,1)), gel(x,2)));
    1652             :       }
    1653             :       /* t_COL */
    1654        1533 :       x = Q_remove_denom(x, &d);
    1655        1533 :       if (!d) return gerepilecopy(av, mkvec2(idealhnf(nf, x), gen_1));
    1656          14 :       mx = zk_multable(nf, x);
    1657          14 :       xZ = zkmultable_capZ(mx);
    1658          14 :       x = ZM_hnfmodid(mx, xZ); /* principal ideal (x) */
    1659          14 :       x0 = mkvec2(xZ, mx); /* same, for fast multiplication */
    1660          14 :       break;
    1661             :     }
    1662             :     default: /* id_MAT */
    1663             :     {
    1664       13762 :       long n = lg(x)-1;
    1665       13762 :       if (n == 0) return mkvec2(gen_0, gen_1);
    1666       13762 :       if (n != nf_get_degree(nf)) pari_err_DIM("idealnumden");
    1667       13762 :       x0 = x = Q_remove_denom(x, &d);
    1668       13762 :       if (!d) return gerepilecopy(av, mkvec2(x, gen_1));
    1669          14 :       break;
    1670             :     }
    1671             :   }
    1672          28 :   J = hnfmodid(x, d); /* = d/B */
    1673          28 :   c = gcoeff(J,1,1); /* (d/B) \cap Z, divides d */
    1674          28 :   B = idealHNF_inv_Z(nf, J); /* (d/B \cap Z) B/d */
    1675          28 :   if (!equalii(c,d)) B = ZM_Z_mul(B, diviiexact(d,c)); /* = B ! */
    1676          28 :   A = idealHNF_mul(nf, B, x0); /* d * (original x) * B = d A */
    1677          28 :   A = ZM_Z_divexact(A, d); /* = A ! */
    1678          28 :   return gerepilecopy(av, mkvec2(A, B));
    1679             : }
    1680             : 
    1681             : /* Return x, integral in 2-elt form, such that pr^n = c * x. Assume n != 0.
    1682             :  * nf = true nf */
    1683             : static GEN
    1684       49264 : idealpowprime(GEN nf, GEN pr, GEN n, GEN *pc)
    1685             : {
    1686       49264 :   GEN p = pr_get_p(pr), q, gen;
    1687             : 
    1688       49264 :   *pc = NULL;
    1689       49264 :   if (is_pm1(n)) /* n = 1 special cased for efficiency */
    1690             :   {
    1691       25027 :     q = p;
    1692       25027 :     if (typ(pr_get_tau(pr)) == t_INT) /* inert */
    1693             :     {
    1694           0 :       *pc = (signe(n) >= 0)? p: ginv(p);
    1695           0 :       return mkvec2(gen_1,gen_0);
    1696             :     }
    1697       25027 :     if (signe(n) >= 0) gen = pr_get_gen(pr);
    1698             :     else
    1699             :     {
    1700        2002 :       gen = pr_get_tau(pr); /* possibly t_MAT */
    1701        2002 :       *pc = ginv(p);
    1702             :     }
    1703             :   }
    1704       24237 :   else if (equalis(n,2)) return idealsqrprime(nf, pr, pc);
    1705             :   else
    1706             :   {
    1707       11638 :     long e = pr_get_e(pr), f = pr_get_f(pr);
    1708       11638 :     GEN r, m = truedvmdis(n, e, &r);
    1709       11638 :     if (e * f == nf_get_degree(nf))
    1710             :     { /* pr^e = (p) */
    1711        5180 :       if (signe(m)) *pc = powii(p,m);
    1712        5180 :       if (!signe(r)) return mkvec2(gen_1,gen_0);
    1713        1862 :       q = p;
    1714        1862 :       gen = nfpow(nf, pr_get_gen(pr), r);
    1715             :     }
    1716             :     else
    1717             :     {
    1718        6458 :       m = absi(m);
    1719        6458 :       if (signe(r)) m = addiu(m,1);
    1720        6458 :       q = powii(p,m); /* m = ceil(|n|/e) */
    1721        6458 :       if (signe(n) >= 0) gen = nfpow(nf, pr_get_gen(pr), n);
    1722             :       else
    1723             :       {
    1724         917 :         gen = pr_get_tau(pr);
    1725         917 :         if (typ(gen) == t_MAT) gen = gel(gen,1);
    1726         917 :         n = negi(n);
    1727         917 :         gen = ZC_Z_divexact(nfpow(nf, gen, n), powii(p, subii(n,m)));
    1728         917 :         *pc = ginv(q);
    1729             :       }
    1730             :     }
    1731        8320 :     gen = FpC_red(gen, q);
    1732             :   }
    1733       33347 :   return mkvec2(q, gen);
    1734             : }
    1735             : 
    1736             : /* x * pr^n. Assume x in HNF or scalar (possibly non-integral) */
    1737             : GEN
    1738       34783 : idealmulpowprime(GEN nf, GEN x, GEN pr, GEN n)
    1739             : {
    1740             :   GEN c, cx, y;
    1741             :   long N;
    1742             : 
    1743       34783 :   nf = checknf(nf);
    1744       34783 :   N = nf_get_degree(nf);
    1745       34783 :   if (!signe(n)) return typ(x) == t_MAT? x: scalarmat_shallow(x, N);
    1746             : 
    1747             :   /* inert, special cased for efficiency */
    1748       34671 :   if (pr_is_inert(pr))
    1749             :   {
    1750        2877 :     GEN q = powii(pr_get_p(pr), n);
    1751        2877 :     return typ(x) == t_MAT? RgM_Rg_mul(x,q): scalarmat_shallow(gmul(x,q), N);
    1752             :   }
    1753             : 
    1754       31794 :   y = idealpowprime(nf, pr, n, &c);
    1755       31794 :   if (typ(x) == t_MAT)
    1756       31318 :   { x = Q_primitive_part(x, &cx); if (is_pm1(gcoeff(x,1,1))) x = NULL; }
    1757             :   else
    1758         476 :   { cx = x; x = NULL; }
    1759       31794 :   cx = mul_content(c,cx);
    1760       31794 :   if (x)
    1761       19950 :     x = idealHNF_mul_two(nf,x,y);
    1762             :   else
    1763       11844 :     x = idealhnf_two(nf,y);
    1764       31794 :   if (cx) x = RgM_Rg_mul(x,cx);
    1765       31794 :   return x;
    1766             : }
    1767             : GEN
    1768        4795 : idealdivpowprime(GEN nf, GEN x, GEN pr, GEN n)
    1769             : {
    1770        4795 :   return idealmulpowprime(nf,x,pr, negi(n));
    1771             : }
    1772             : 
    1773             : /* nf = true nf */
    1774             : static GEN
    1775      157349 : idealpow_aux(GEN nf, GEN x, long tx, GEN n)
    1776             : {
    1777      157349 :   GEN T = nf_get_pol(nf), m, cx, n1, a, alpha;
    1778      157349 :   long N = degpol(T), s = signe(n);
    1779      157349 :   if (!s) return matid(N);
    1780      154499 :   switch(tx)
    1781             :   {
    1782             :     case id_PRINCIPAL:
    1783           0 :       return idealhnf_principal(nf, nfpow(nf,x,n));
    1784             :     case id_PRIME:
    1785       60779 :       if (pr_is_inert(x)) return scalarmat(powii(gel(x,1), n), N);
    1786       17470 :       x = idealpowprime(nf, x, n, &cx);
    1787       17470 :       x = idealhnf_two(nf,x);
    1788       17470 :       return cx? RgM_Rg_mul(x, cx): x;
    1789             :     default:
    1790       93720 :       if (is_pm1(n)) return (s < 0)? idealinv(nf, x): gcopy(x);
    1791       52852 :       n1 = (s < 0)? negi(n): n;
    1792             : 
    1793       52852 :       x = Q_primitive_part(x, &cx);
    1794       52852 :       a = mat_ideal_two_elt(nf,x); alpha = gel(a,2); a = gel(a,1);
    1795       52852 :       alpha = nfpow(nf,alpha,n1);
    1796       52852 :       m = zk_scalar_or_multable(nf, alpha);
    1797       52852 :       if (typ(m) == t_INT) {
    1798          91 :         x = gcdii(powii(a,n1), m);
    1799          91 :         if (s<0) x = ginv(x);
    1800          91 :         if (cx) x = gmul(x, powgi(cx,n));
    1801          91 :         x = scalarmat(x, N);
    1802             :       }
    1803             :       else
    1804             :       {
    1805       52761 :         x = ZM_hnfmodid(m, gcdii(powii(a,n1), zkmultable_capZ(m)));
    1806       52761 :         if (cx) cx = powgi(cx,n);
    1807       52761 :         if (s<0) {
    1808           7 :           GEN xZ = gcoeff(x,1,1);
    1809           7 :           cx = cx ? gdiv(cx, xZ): ginv(xZ);
    1810           7 :           x = idealHNF_inv_Z(nf,x);
    1811             :         }
    1812       52761 :         if (cx) x = RgM_Rg_mul(x, cx);
    1813             :       }
    1814       52852 :       return x;
    1815             :   }
    1816             : }
    1817             : 
    1818             : /* raise the ideal x to the power n (in Z) */
    1819             : GEN
    1820      157349 : idealpow(GEN nf, GEN x, GEN n)
    1821             : {
    1822             :   pari_sp av;
    1823             :   long tx;
    1824             :   GEN res, ax;
    1825             : 
    1826      157349 :   if (typ(n) != t_INT) pari_err_TYPE("idealpow",n);
    1827      157349 :   tx = idealtyp(&x,&ax);
    1828      157349 :   res = ax? cgetg(3,t_VEC): NULL;
    1829      157349 :   av = avma;
    1830      157349 :   x = gerepileupto(av, idealpow_aux(checknf(nf), x, tx, n));
    1831      157349 :   if (!ax) return x;
    1832         421 :   ax = ext_pow(nf, ax, n);
    1833         421 :   gel(res,1) = x;
    1834         421 :   gel(res,2) = ax;
    1835         421 :   return res;
    1836             : }
    1837             : 
    1838             : /* Return ideal^e in number field nf. e is a C integer. */
    1839             : GEN
    1840       14268 : idealpows(GEN nf, GEN ideal, long e)
    1841             : {
    1842       14268 :   long court[] = {evaltyp(t_INT) | _evallg(3),0,0};
    1843       14268 :   affsi(e,court); return idealpow(nf,ideal,court);
    1844             : }
    1845             : 
    1846             : static GEN
    1847       45342 : _idealmulred(GEN nf, GEN x, GEN y)
    1848       45342 : { return idealred(nf,idealmul(nf,x,y)); }
    1849             : static GEN
    1850       66105 : _idealsqrred(GEN nf, GEN x)
    1851       66105 : { return idealred(nf,idealsqr(nf,x)); }
    1852             : static GEN
    1853       28664 : _mul(void *data, GEN x, GEN y) { return _idealmulred((GEN)data,x,y); }
    1854             : static GEN
    1855       66105 : _sqr(void *data, GEN x) { return _idealsqrred((GEN)data, x); }
    1856             : 
    1857             : /* compute x^n (x ideal, n integer), reducing along the way */
    1858             : GEN
    1859       69200 : idealpowred(GEN nf, GEN x, GEN n)
    1860             : {
    1861       69200 :   pari_sp av = avma;
    1862             :   long s;
    1863             :   GEN y;
    1864             : 
    1865       69200 :   if (typ(n) != t_INT) pari_err_TYPE("idealpowred",n);
    1866       69200 :   s = signe(n); if (s == 0) return idealpow(nf,x,n);
    1867       68779 :   y = gen_pow(x, n, (void*)nf, &_sqr, &_mul);
    1868             : 
    1869       68779 :   if (s < 0) y = idealinv(nf,y);
    1870       68779 :   if (s < 0 || is_pm1(n)) y = idealred(nf,y);
    1871       68779 :   return gerepileupto(av,y);
    1872             : }
    1873             : 
    1874             : GEN
    1875       16678 : idealmulred(GEN nf, GEN x, GEN y)
    1876             : {
    1877       16678 :   pari_sp av = avma;
    1878       16678 :   return gerepileupto(av, _idealmulred(nf,x,y));
    1879             : }
    1880             : 
    1881             : long
    1882          91 : isideal(GEN nf,GEN x)
    1883             : {
    1884          91 :   long N, i, j, lx, tx = typ(x);
    1885             :   pari_sp av;
    1886             :   GEN T, xZ;
    1887             : 
    1888          91 :   nf = checknf(nf); T = nf_get_pol(nf); lx = lg(x);
    1889          91 :   if (tx==t_VEC && lx==3) { x = gel(x,1); tx = typ(x); lx = lg(x); }
    1890          91 :   switch(tx)
    1891             :   {
    1892          14 :     case t_INT: case t_FRAC: return 1;
    1893           7 :     case t_POL: return varn(x) == varn(T);
    1894           7 :     case t_POLMOD: return RgX_equal_var(T, gel(x,1));
    1895          14 :     case t_VEC: return get_prid(x)? 1 : 0;
    1896          42 :     case t_MAT: break;
    1897           7 :     default: return 0;
    1898             :   }
    1899          42 :   N = degpol(T);
    1900          42 :   if (lx-1 != N) return (lx == 1);
    1901          28 :   if (nbrows(x) != N) return 0;
    1902             : 
    1903          28 :   av = avma; x = Q_primpart(x);
    1904          28 :   if (!ZM_ishnf(x)) return 0;
    1905          14 :   xZ = gcoeff(x,1,1);
    1906          21 :   for (j=2; j<=N; j++)
    1907          14 :     if (!dvdii(xZ, gcoeff(x,j,j))) { avma = av; return 0; }
    1908          14 :   for (i=2; i<=N; i++)
    1909          14 :     for (j=2; j<=N; j++)
    1910           7 :       if (! hnf_invimage(x, zk_ei_mul(nf,gel(x,i),j))) { avma = av; return 0; }
    1911           7 :   avma=av; return 1;
    1912             : }
    1913             : 
    1914             : GEN
    1915       13748 : idealdiv(GEN nf, GEN x, GEN y)
    1916             : {
    1917       13748 :   pari_sp av = avma, tetpil;
    1918       13748 :   GEN z = idealinv(nf,y);
    1919       13748 :   tetpil = avma; return gerepile(av,tetpil, idealmul(nf,x,z));
    1920             : }
    1921             : 
    1922             : /* This routine computes the quotient x/y of two ideals in the number field nf.
    1923             :  * It assumes that the quotient is an integral ideal.  The idea is to find an
    1924             :  * ideal z dividing y such that gcd(Nx/Nz, Nz) = 1.  Then
    1925             :  *
    1926             :  *   x + (Nx/Nz)    x
    1927             :  *   ----------- = ---
    1928             :  *   y + (Ny/Nz)    y
    1929             :  *
    1930             :  * Proof: we can assume x and y are integral. Let p be any prime ideal
    1931             :  *
    1932             :  * If p | Nz, then it divides neither Nx/Nz nor Ny/Nz (since Nx/Nz is the
    1933             :  * product of the integers N(x/y) and N(y/z)).  Both the numerator and the
    1934             :  * denominator on the left will be coprime to p.  So will x/y, since x/y is
    1935             :  * assumed integral and its norm N(x/y) is coprime to p.
    1936             :  *
    1937             :  * If instead p does not divide Nz, then v_p (Nx/Nz) = v_p (Nx) >= v_p(x).
    1938             :  * Hence v_p (x + Nx/Nz) = v_p(x).  Likewise for the denominators.  QED.
    1939             :  *
    1940             :  *                Peter Montgomery.  July, 1994. */
    1941             : static void
    1942           7 : err_divexact(GEN x, GEN y)
    1943           7 : { pari_err_DOMAIN("idealdivexact","denominator(x/y)", "!=",
    1944           0 :                   gen_1,mkvec2(x,y)); }
    1945             : GEN
    1946         574 : idealdivexact(GEN nf, GEN x0, GEN y0)
    1947             : {
    1948         574 :   pari_sp av = avma;
    1949             :   GEN x, y, yZ, Nx, Ny, Nz, cy, q, r;
    1950             : 
    1951         574 :   nf = checknf(nf);
    1952         574 :   x = idealhnf_shallow(nf, x0);
    1953         574 :   y = idealhnf_shallow(nf, y0);
    1954         574 :   if (lg(y) == 1) pari_err_INV("idealdivexact", y0);
    1955         567 :   if (lg(x) == 1) { avma = av; return cgetg(1, t_MAT); } /* numerator is zero */
    1956         567 :   y = Q_primitive_part(y, &cy);
    1957         567 :   if (cy) x = RgM_Rg_div(x,cy);
    1958         567 :   Nx = idealnorm(nf,x);
    1959         567 :   Ny = idealnorm(nf,y);
    1960         567 :   if (typ(Nx) != t_INT) err_divexact(x,y);
    1961         560 :   q = dvmdii(Nx,Ny, &r);
    1962         560 :   if (signe(r)) err_divexact(x,y);
    1963         560 :   if (is_pm1(q)) { avma = av; return matid(nf_get_degree(nf)); }
    1964             :   /* Find a norm Nz | Ny such that gcd(Nx/Nz, Nz) = 1 */
    1965         343 :   for (Nz = Ny;;) /* q = Nx/Nz */
    1966             :   {
    1967         364 :     GEN p1 = gcdii(Nz, q);
    1968         364 :     if (is_pm1(p1)) break;
    1969          21 :     Nz = diviiexact(Nz,p1);
    1970          21 :     q = mulii(q,p1);
    1971          21 :   }
    1972             :   /* Replace x/y  by  x+(Nx/Nz) / y+(Ny/Nz) */
    1973         343 :   x = ZM_hnfmodid(x, q);
    1974             :   /* y reduced to unit ideal ? */
    1975         343 :   if (Nz == Ny) return gerepileupto(av, x);
    1976             : 
    1977          21 :   y = ZM_hnfmodid(y, diviiexact(Ny,Nz));
    1978          21 :   yZ = gcoeff(y,1,1);
    1979          21 :   y = idealHNF_mul(nf,x, idealHNF_inv_Z(nf,y));
    1980          21 :   return gerepileupto(av, RgM_Rg_div(y, yZ));
    1981             : }
    1982             : 
    1983             : GEN
    1984          21 : idealintersect(GEN nf, GEN x, GEN y)
    1985             : {
    1986          21 :   pari_sp av = avma;
    1987             :   long lz, lx, i;
    1988             :   GEN z, dx, dy, xZ, yZ;;
    1989             : 
    1990          21 :   nf = checknf(nf);
    1991          21 :   x = idealhnf_shallow(nf,x);
    1992          21 :   y = idealhnf_shallow(nf,y);
    1993          21 :   if (lg(x) == 1 || lg(y) == 1) { avma = av; return cgetg(1,t_MAT); }
    1994          14 :   x = Q_remove_denom(x, &dx);
    1995          14 :   y = Q_remove_denom(y, &dy);
    1996          14 :   if (dx) y = ZM_Z_mul(y, dx);
    1997          14 :   if (dy) x = ZM_Z_mul(x, dy);
    1998          14 :   xZ = gcoeff(x,1,1);
    1999          14 :   yZ = gcoeff(y,1,1);
    2000          14 :   dx = mul_denom(dx,dy);
    2001          14 :   z = ZM_lll(shallowconcat(x,y), 0.99, LLL_KER); lz = lg(z);
    2002          14 :   lx = lg(x);
    2003          14 :   for (i=1; i<lz; i++) setlg(z[i], lx);
    2004          14 :   z = ZM_hnfmodid(ZM_mul(x,z), lcmii(xZ, yZ));
    2005          14 :   if (dx) z = RgM_Rg_div(z,dx);
    2006          14 :   return gerepileupto(av,z);
    2007             : }
    2008             : 
    2009             : /*******************************************************************/
    2010             : /*                                                                 */
    2011             : /*                      T2-IDEAL REDUCTION                         */
    2012             : /*                                                                 */
    2013             : /*******************************************************************/
    2014             : 
    2015             : static GEN
    2016          21 : chk_vdir(GEN nf, GEN vdir)
    2017             : {
    2018          21 :   long i, l = lg(vdir);
    2019             :   GEN v;
    2020          21 :   if (l != lg(nf_get_roots(nf))) pari_err_DIM("idealred");
    2021          14 :   switch(typ(vdir))
    2022             :   {
    2023           0 :     case t_VECSMALL: return vdir;
    2024          14 :     case t_VEC: break;
    2025           0 :     default: pari_err_TYPE("idealred",vdir);
    2026             :   }
    2027          14 :   v = cgetg(l, t_VECSMALL);
    2028          14 :   for (i = 1; i < l; i++) v[i] = itos(gceil(gel(vdir,i)));
    2029          14 :   return v;
    2030             : }
    2031             : 
    2032             : static void
    2033       26983 : twistG(GEN G, long r1, long i, long v)
    2034             : {
    2035       26983 :   long j, lG = lg(G);
    2036       26983 :   if (i <= r1) {
    2037       24582 :     for (j=1; j<lG; j++) gcoeff(G,i,j) = gmul2n(gcoeff(G,i,j), v);
    2038             :   } else {
    2039        2401 :     long k = (i<<1) - r1;
    2040       13251 :     for (j=1; j<lG; j++)
    2041             :     {
    2042       10850 :       gcoeff(G,k-1,j) = gmul2n(gcoeff(G,k-1,j), v);
    2043       10850 :       gcoeff(G,k  ,j) = gmul2n(gcoeff(G,k  ,j), v);
    2044             :     }
    2045             :   }
    2046       26983 : }
    2047             : 
    2048             : GEN
    2049      178134 : nf_get_Gtwist(GEN nf, GEN vdir)
    2050             : {
    2051             :   long i, l, v, r1;
    2052             :   GEN G;
    2053             : 
    2054      178134 :   if (!vdir) return nf_get_roundG(nf);
    2055        3092 :   if (typ(vdir) == t_MAT)
    2056             :   {
    2057        3071 :     long N = nf_get_degree(nf);
    2058        3071 :     if (lg(vdir) != N+1 || lgcols(vdir) != N+1) pari_err_DIM("idealred");
    2059        3071 :     return vdir;
    2060             :   }
    2061          21 :   vdir = chk_vdir(nf, vdir);
    2062          14 :   G = RgM_shallowcopy(nf_get_G(nf));
    2063          14 :   r1 = nf_get_r1(nf);
    2064          14 :   l = lg(vdir);
    2065          56 :   for (i=1; i<l; i++)
    2066             :   {
    2067          42 :     v = vdir[i]; if (!v) continue;
    2068          42 :     twistG(G, r1, i, v);
    2069             :   }
    2070          14 :   return RM_round_maxrank(G);
    2071             : }
    2072             : GEN
    2073       26941 : nf_get_Gtwist1(GEN nf, long i)
    2074             : {
    2075       26941 :   GEN G = RgM_shallowcopy( nf_get_G(nf) );
    2076       26941 :   long r1 = nf_get_r1(nf);
    2077       26941 :   twistG(G, r1, i, 10);
    2078       26941 :   return RM_round_maxrank(G);
    2079             : }
    2080             : 
    2081             : GEN
    2082       32324 : RM_round_maxrank(GEN G0)
    2083             : {
    2084       32324 :   long e, r = lg(G0)-1;
    2085       32324 :   pari_sp av = avma;
    2086       32324 :   GEN G = G0;
    2087       32324 :   for (e = 4; ; e <<= 1)
    2088             :   {
    2089       32324 :     GEN H = ground(G);
    2090       64648 :     if (ZM_rank(H) == r) return H; /* maximal rank ? */
    2091           0 :     avma = av;
    2092           0 :     G = gmul2n(G0, e);
    2093           0 :   }
    2094             : }
    2095             : 
    2096             : GEN
    2097      178127 : idealred0(GEN nf, GEN I, GEN vdir)
    2098             : {
    2099      178127 :   pari_sp av = avma;
    2100      178127 :   GEN G, aI, IZ, J, y, yZ, my, c1 = NULL;
    2101             :   long N;
    2102             : 
    2103      178127 :   nf = checknf(nf);
    2104      178127 :   N = nf_get_degree(nf);
    2105             :   /* put first for sanity checks, unused when I obviously principal */
    2106      178127 :   G = nf_get_Gtwist(nf, vdir);
    2107      178120 :   switch (idealtyp(&I,&aI))
    2108             :   {
    2109             :     case id_PRIME:
    2110       39180 :       if (pr_is_inert(I)) {
    2111         581 :         if (!aI) { avma = av; return matid(N); }
    2112         581 :         c1 = gel(I,1); I = matid(N);
    2113         581 :         goto END;
    2114             :       }
    2115       38599 :       IZ = pr_get_p(I);
    2116       38599 :       J = pr_inv_p(I);
    2117       38599 :       I = idealhnf_two(nf,I);
    2118       38599 :       break;
    2119             :     case id_MAT:
    2120      138926 :       I = Q_primitive_part(I, &c1);
    2121      138926 :       IZ = gcoeff(I,1,1);
    2122      138926 :       if (is_pm1(IZ))
    2123             :       {
    2124        7504 :         if (!aI) { avma = av; return matid(N); }
    2125        7448 :         goto END;
    2126             :       }
    2127      131422 :       J = idealHNF_inv_Z(nf, I);
    2128      131422 :       break;
    2129             :     default: /* id_PRINCIPAL, silly case */
    2130          14 :       if (gequal0(I)) I = cgetg(1,t_MAT); else { c1 = I; I = matid(N); }
    2131          14 :       if (!aI) return I;
    2132           7 :       goto END;
    2133             :   }
    2134             :   /* now I integral, HNF; and J = (I\cap Z) I^(-1), integral */
    2135      170021 :   y = idealpseudomin(J, G); /* small elt in (I\cap Z)I^(-1), integral */
    2136      170021 :   if (ZV_isscalar(y))
    2137             :   { /* already reduced */
    2138       60497 :     if (!aI) return gerepilecopy(av, I);
    2139       60098 :     goto END;
    2140             :   }
    2141             : 
    2142      109524 :   my = zk_multable(nf, y);
    2143      109524 :   I = ZM_Z_divexact(ZM_mul(my, I), IZ); /* y I / (I\cap Z), integral */
    2144      109524 :   c1 = mul_content(c1, IZ);
    2145      109524 :   my = ZM_gauss(my, col_ei(N,1)); /* y^-1 */
    2146      109524 :   yZ = Q_denom(my); /* (y) \cap Z */
    2147      109524 :   I = hnfmodid(I, yZ);
    2148      109524 :   if (!aI) return gerepileupto(av, I);
    2149      109258 :   c1 = RgC_Rg_mul(my, c1);
    2150             : END:
    2151      177392 :   if (c1) aI = ext_mul(nf, aI,c1);
    2152      177392 :   return gerepilecopy(av, mkvec2(I, aI));
    2153             : }
    2154             : 
    2155             : GEN
    2156           7 : idealmin(GEN nf, GEN x, GEN vdir)
    2157             : {
    2158           7 :   pari_sp av = avma;
    2159             :   GEN y, dx;
    2160           7 :   nf = checknf(nf);
    2161           7 :   switch( idealtyp(&x,&y) )
    2162             :   {
    2163           0 :     case id_PRINCIPAL: return gcopy(x);
    2164           0 :     case id_PRIME: x = idealhnf_two(nf,x); break;
    2165           7 :     case id_MAT: if (lg(x) == 1) return gen_0;
    2166             :   }
    2167           7 :   x = Q_remove_denom(x, &dx);
    2168           7 :   y = idealpseudomin(x, nf_get_Gtwist(nf,vdir));
    2169           7 :   if (dx) y = RgC_Rg_div(y, dx);
    2170           7 :   return gerepileupto(av, y);
    2171             : }
    2172             : 
    2173             : /*******************************************************************/
    2174             : /*                                                                 */
    2175             : /*                   APPROXIMATION THEOREM                         */
    2176             : /*                                                                 */
    2177             : /*******************************************************************/
    2178             : /* a = ppi(a,b) ppo(a,b), where ppi regroups primes common to a and b
    2179             :  * and ppo(a,b) = Z_ppo(a,b) */
    2180             : /* return gcd(a,b),ppi(a,b),ppo(a,b) */
    2181             : GEN
    2182      452130 : Z_ppio(GEN a, GEN b)
    2183             : {
    2184      452130 :   GEN x, y, d = gcdii(a,b);
    2185      452130 :   if (is_pm1(d)) return mkvec3(gen_1, gen_1, a);
    2186      343959 :   x = d; y = diviiexact(a,d);
    2187             :   for(;;)
    2188             :   {
    2189      406434 :     GEN g = gcdii(x,y);
    2190      406434 :     if (is_pm1(g)) return mkvec3(d, x, y);
    2191       62475 :     x = mulii(x,g); y = diviiexact(y,g);
    2192       62475 :   }
    2193             : }
    2194             : /* a = ppg(a,b)pple(a,b), where ppg regroups primes such that v(a) > v(b)
    2195             :  * and pple all others */
    2196             : /* return gcd(a,b),ppg(a,b),pple(a,b) */
    2197             : GEN
    2198           0 : Z_ppgle(GEN a, GEN b)
    2199             : {
    2200           0 :   GEN x, y, g, d = gcdii(a,b);
    2201           0 :   if (equalii(a, d)) return mkvec3(a, gen_1, a);
    2202           0 :   x = diviiexact(a,d); y = d;
    2203             :   for(;;)
    2204             :   {
    2205           0 :     g = gcdii(x,y);
    2206           0 :     if (is_pm1(g)) return mkvec3(d, x, y);
    2207           0 :     x = mulii(x,g); y = diviiexact(y,g);
    2208           0 :   }
    2209             : }
    2210             : static void
    2211           0 : Z_dcba_rec(GEN L, GEN a, GEN b)
    2212             : {
    2213             :   GEN x, r, v, g, h, c, c0;
    2214             :   long n;
    2215           0 :   if (is_pm1(b)) {
    2216           0 :     if (!is_pm1(a)) vectrunc_append(L, a);
    2217           0 :     return;
    2218             :   }
    2219           0 :   v = Z_ppio(a,b);
    2220           0 :   a = gel(v,2);
    2221           0 :   r = gel(v,3);
    2222           0 :   if (!is_pm1(r)) vectrunc_append(L, r);
    2223           0 :   v = Z_ppgle(a,b);
    2224           0 :   g = gel(v,1);
    2225           0 :   h = gel(v,2);
    2226           0 :   x = c0 = gel(v,3);
    2227           0 :   for (n = 1; !is_pm1(h); n++)
    2228             :   {
    2229             :     GEN d, y;
    2230             :     long i;
    2231           0 :     v = Z_ppgle(h,sqri(g));
    2232           0 :     g = gel(v,1);
    2233           0 :     h = gel(v,2);
    2234           0 :     c = gel(v,3); if (is_pm1(c)) continue;
    2235           0 :     d = gcdii(c,b);
    2236           0 :     x = mulii(x,d);
    2237           0 :     y = d; for (i=1; i < n; i++) y = sqri(y);
    2238           0 :     Z_dcba_rec(L, diviiexact(c,y), d);
    2239             :   }
    2240           0 :   Z_dcba_rec(L,diviiexact(b,x), c0);
    2241             : }
    2242             : static GEN
    2243     3062871 : Z_cba_rec(GEN L, GEN a, GEN b)
    2244             : {
    2245             :   GEN g;
    2246     3062871 :   if (lg(L) > 10)
    2247             :   { /* a few naive steps before switching to dcba */
    2248           0 :     Z_dcba_rec(L, a, b);
    2249           0 :     return gel(L, lg(L)-1);
    2250             :   }
    2251     3062871 :   if (is_pm1(a)) return b;
    2252     1819839 :   g = gcdii(a,b);
    2253     1819839 :   if (is_pm1(g)) { vectrunc_append(L, a); return b; }
    2254     1359456 :   a = diviiexact(a,g);
    2255     1359456 :   b = diviiexact(b,g);
    2256     1359456 :   return Z_cba_rec(L, Z_cba_rec(L, a, g), b);
    2257             : }
    2258             : GEN
    2259      343959 : Z_cba(GEN a, GEN b)
    2260             : {
    2261      343959 :   GEN L = vectrunc_init(expi(a) + expi(b) + 2);
    2262      343959 :   GEN t = Z_cba_rec(L, a, b);
    2263      343959 :   if (!is_pm1(t)) vectrunc_append(L, t);
    2264      343959 :   return L;
    2265             : }
    2266             : 
    2267             : /* write x = x1 x2, x2 maximal s.t. (x2,f) = 1, return x2 */
    2268             : GEN
    2269     1104451 : Z_ppo(GEN x, GEN f)
    2270             : {
    2271             :   for (;;)
    2272             :   {
    2273     1104451 :     f = gcdii(x, f); if (is_pm1(f)) break;
    2274      752261 :     x = diviiexact(x, f);
    2275      752261 :   }
    2276      352190 :   return x;
    2277             : }
    2278             : /* write x = x1 x2, x2 maximal s.t. (x2,f) = 1, return x2 */
    2279             : ulong
    2280      234178 : u_ppo(ulong x, ulong f)
    2281             : {
    2282             :   for (;;)
    2283             :   {
    2284      234178 :     f = ugcd(x, f); if (f == 1) break;
    2285      104132 :     x /= f;
    2286      104132 :   }
    2287      130046 :   return x;
    2288             : }
    2289             : 
    2290             : /* x t_INT, f ideal. Write x = x1 x2, sqf(x1) | f, (x2,f) = 1. Return x2 */
    2291             : static GEN
    2292           7 : nf_coprime_part(GEN nf, GEN x, GEN listpr)
    2293             : {
    2294           7 :   long v, j, lp = lg(listpr), N = nf_get_degree(nf);
    2295             :   GEN x1, x2, ex;
    2296             : 
    2297             : #if 0 /*1) via many gcds. Expensive ! */
    2298             :   GEN f = idealprodprime(nf, listpr);
    2299             :   f = ZM_hnfmodid(f, x); /* first gcd is less expensive since x in Z */
    2300             :   x = scalarmat(x, N);
    2301             :   for (;;)
    2302             :   {
    2303             :     if (gequal1(gcoeff(f,1,1))) break;
    2304             :     x = idealdivexact(nf, x, f);
    2305             :     f = ZM_hnfmodid(shallowconcat(f,x), gcoeff(x,1,1)); /* gcd(f,x) */
    2306             :   }
    2307             :   x2 = x;
    2308             : #else /*2) from prime decomposition */
    2309           7 :   x1 = NULL;
    2310          21 :   for (j=1; j<lp; j++)
    2311             :   {
    2312          14 :     GEN pr = gel(listpr,j);
    2313          14 :     v = Z_pval(x, pr_get_p(pr)); if (!v) continue;
    2314             : 
    2315           7 :     ex = muluu(v, pr_get_e(pr)); /* = v_pr(x) > 0 */
    2316           7 :     x1 = x1? idealmulpowprime(nf, x1, pr, ex)
    2317           7 :            : idealpow(nf, pr, ex);
    2318             :   }
    2319           7 :   x = scalarmat(x, N);
    2320           7 :   x2 = x1? idealdivexact(nf, x, x1): x;
    2321             : #endif
    2322           7 :   return x2;
    2323             : }
    2324             : 
    2325             : /* L0 in K^*, assume (L0,f) = 1. Return L integral, L0 = L mod f  */
    2326             : GEN
    2327        2002 : make_integral(GEN nf, GEN L0, GEN f, GEN listpr)
    2328             : {
    2329             :   GEN fZ, t, L, D2, d1, d2, d;
    2330             : 
    2331        2002 :   L = Q_remove_denom(L0, &d);
    2332        2002 :   if (!d) return L0;
    2333             : 
    2334             :   /* L0 = L / d, L integral */
    2335        1365 :   fZ = gcoeff(f,1,1);
    2336        1365 :   if (typ(L) == t_INT) return Fp_mul(L, Fp_inv(d, fZ), fZ);
    2337             :   /* Kill denom part coprime to fZ */
    2338        1351 :   d2 = Z_ppo(d, fZ);
    2339        1351 :   t = Fp_inv(d2, fZ); if (!is_pm1(t)) L = ZC_Z_mul(L,t);
    2340        1351 :   if (equalii(d, d2)) return L;
    2341             : 
    2342           7 :   d1 = diviiexact(d, d2);
    2343             :   /* L0 = (L / d1) mod f. d1 not coprime to f
    2344             :    * write (d1) = D1 D2, D2 minimal, (D2,f) = 1. */
    2345           7 :   D2 = nf_coprime_part(nf, d1, listpr);
    2346           7 :   t = idealaddtoone_i(nf, D2, f); /* in D2, 1 mod f */
    2347           7 :   L = nfmuli(nf,t,L);
    2348             : 
    2349             :   /* if (L0, f) = 1, then L in D1 ==> in D1 D2 = (d1) */
    2350           7 :   return Q_div_to_int(L, d1); /* exact division */
    2351             : }
    2352             : 
    2353             : /* assume L is a list of prime ideals. Return the product */
    2354             : GEN
    2355         126 : idealprodprime(GEN nf, GEN L)
    2356             : {
    2357         126 :   long l = lg(L), i;
    2358             :   GEN z;
    2359         126 :   if (l == 1) return matid(nf_get_degree(nf));
    2360         126 :   z = idealhnf_two(nf, gel(L,1));
    2361         126 :   for (i=2; i<l; i++) z = idealHNF_mul_two(nf,z, gel(L,i));
    2362         126 :   return z;
    2363             : }
    2364             : 
    2365             : /* optimize for the frequent case I = nfhnf()[2]: lots of them are 1 */
    2366             : GEN
    2367        1134 : idealprod(GEN nf, GEN I)
    2368             : {
    2369        1134 :   long i, l = lg(I);
    2370             :   GEN z;
    2371        1904 :   for (i = 1; i < l; i++)
    2372        1876 :     if (!equali1(gel(I,i))) break;
    2373        1134 :   if (i == l) return gen_1;
    2374        1106 :   z = gel(I,i);
    2375        1106 :   for (i++; i<l; i++) z = idealmul(nf, z, gel(I,i));
    2376        1106 :   return z;
    2377             : }
    2378             : 
    2379             : /* assume L is a list of prime ideals. Return prod L[i]^e[i] */
    2380             : GEN
    2381        5005 : factorbackprime(GEN nf, GEN L, GEN e)
    2382             : {
    2383        5005 :   long l = lg(L), i;
    2384             :   GEN z;
    2385             : 
    2386        5005 :   if (l == 1) return matid(nf_get_degree(nf));
    2387        4991 :   z = idealpow(nf, gel(L,1), gel(e,1));
    2388        7133 :   for (i=2; i<l; i++)
    2389        2142 :     if (signe(gel(e,i))) z = idealmulpowprime(nf,z, gel(L,i),gel(e,i));
    2390        4991 :   return z;
    2391             : }
    2392             : 
    2393             : /* F in Z, divisible exactly by pr.p. Return F-uniformizer for pr, i.e.
    2394             :  * a t in Z_K such that v_pr(t) = 1 and (t, F/pr) = 1 */
    2395             : GEN
    2396       11578 : pr_uniformizer(GEN pr, GEN F)
    2397             : {
    2398       11578 :   GEN p = pr_get_p(pr), t = pr_get_gen(pr);
    2399       11578 :   if (!equalii(F, p))
    2400             :   {
    2401        3843 :     long e = pr_get_e(pr);
    2402        3843 :     GEN u, v, q = (e == 1)? sqri(p): p;
    2403        3843 :     u = mulii(q, Fp_inv(q, diviiexact(F,p))); /* 1 mod F/p, 0 mod q */
    2404        3843 :     v = subui(1UL, u); /* 0 mod F/p, 1 mod q */
    2405        3843 :     if (pr_is_inert(pr))
    2406           0 :       t = addii(mulii(p, v), u);
    2407             :     else
    2408             :     {
    2409        3843 :       t = ZC_Z_mul(t, v);
    2410        3843 :       gel(t,1) = addii(gel(t,1), u); /* return u + vt */
    2411             :     }
    2412             :   }
    2413       11578 :   return t;
    2414             : }
    2415             : /* L = list of prime ideals, return lcm_i (L[i] \cap \ZM) */
    2416             : GEN
    2417       11578 : prV_lcm_capZ(GEN L)
    2418             : {
    2419       11578 :   long i, r = lg(L);
    2420             :   GEN F;
    2421       11578 :   if (r == 1) return gen_1;
    2422       10493 :   F = pr_get_p(gel(L,1));
    2423       14770 :   for (i = 2; i < r; i++)
    2424             :   {
    2425        4277 :     GEN pr = gel(L,i), p = pr_get_p(pr);
    2426        4277 :     if (!dvdii(F, p)) F = mulii(F,p);
    2427             :   }
    2428       10493 :   return F;
    2429             : }
    2430             : 
    2431             : /* Given a prime ideal factorization with possibly zero or negative
    2432             :  * exponents, gives b such that v_p(b) = v_p(x) for all prime ideals pr | x
    2433             :  * and v_pr(b)> = 0 for all other pr.
    2434             :  * For optimal performance, all [anti-]uniformizers should be precomputed,
    2435             :  * but no support for this yet.
    2436             :  *
    2437             :  * If nored, do not reduce result.
    2438             :  * No garbage collecting */
    2439             : static GEN
    2440       10724 : idealapprfact_i(GEN nf, GEN x, int nored)
    2441             : {
    2442             :   GEN z, d, L, e, e2, F;
    2443             :   long i, r;
    2444             :   int flagden;
    2445             : 
    2446       10724 :   nf = checknf(nf);
    2447       10724 :   L = gel(x,1);
    2448       10724 :   e = gel(x,2);
    2449       10724 :   F = prV_lcm_capZ(L);
    2450       10724 :   flagden = 0;
    2451       10724 :   z = NULL; r = lg(e);
    2452       24444 :   for (i = 1; i < r; i++)
    2453             :   {
    2454       13720 :     long s = signe(gel(e,i));
    2455             :     GEN pi, q;
    2456       13720 :     if (!s) continue;
    2457       11508 :     if (s < 0) flagden = 1;
    2458       11508 :     pi = pr_uniformizer(gel(L,i), F);
    2459       11508 :     q = nfpow(nf, pi, gel(e,i));
    2460       11508 :     z = z? nfmul(nf, z, q): q;
    2461             :   }
    2462       10724 :   if (!z) return gen_1;
    2463        7952 :   if (nored || typ(z) != t_COL) return z;
    2464         693 :   e2 = cgetg(r, t_VEC);
    2465         693 :   for (i=1; i<r; i++) gel(e2,i) = addis(gel(e,i), 1);
    2466         693 :   x = factorbackprime(nf, L,e2);
    2467         693 :   if (flagden) /* denominator */
    2468             :   {
    2469         679 :     z = Q_remove_denom(z, &d);
    2470         679 :     d = diviiexact(d, Z_ppo(d, F));
    2471         679 :     x = RgM_Rg_mul(x, d);
    2472             :   }
    2473             :   else
    2474          14 :     d = NULL;
    2475         693 :   z = ZC_reducemodlll(z, x);
    2476         693 :   return d? RgC_Rg_div(z,d): z;
    2477             : }
    2478             : 
    2479             : GEN
    2480           0 : idealapprfact(GEN nf, GEN x) {
    2481           0 :   pari_sp av = avma;
    2482           0 :   return gerepileupto(av, idealapprfact_i(nf, x, 0));
    2483             : }
    2484             : GEN
    2485          14 : idealappr(GEN nf, GEN x) {
    2486          14 :   pari_sp av = avma;
    2487          14 :   if (!is_nf_extfactor(x)) x = idealfactor(nf, x);
    2488          14 :   return gerepileupto(av, idealapprfact_i(nf, x, 0));
    2489             : }
    2490             : 
    2491             : /* OBSOLETE */
    2492             : GEN
    2493          14 : idealappr0(GEN nf, GEN x, long fl) { (void)fl; return idealappr(nf, x); }
    2494             : 
    2495             : static GEN
    2496          21 : mat_ideal_two_elt2(GEN nf, GEN x, GEN a)
    2497             : {
    2498          21 :   GEN F = idealfactor(nf,a), P = gel(F,1), E = gel(F,2);
    2499          21 :   long i, r = lg(E);
    2500          21 :   for (i=1; i<r; i++) gel(E,i) = stoi( idealval(nf,x,gel(P,i)) );
    2501          21 :   return idealapprfact_i(nf,F,1);
    2502             : }
    2503             : 
    2504             : static void
    2505          14 : not_in_ideal(GEN a) {
    2506          14 :   pari_err_DOMAIN("idealtwoelt2","element mod ideal", "!=", gen_0, a);
    2507           0 : }
    2508             : /* x integral in HNF, a an 'nf' */
    2509             : static int
    2510          28 : in_ideal(GEN x, GEN a)
    2511             : {
    2512          28 :   switch(typ(a))
    2513             :   {
    2514          14 :     case t_INT: return dvdii(a, gcoeff(x,1,1));
    2515           7 :     case t_COL: return RgV_is_ZV(a) && !!hnf_invimage(x, a);
    2516           7 :     default: return 0;
    2517             :   }
    2518             : }
    2519             : 
    2520             : /* Given an integral ideal x and a in x, gives a b such that
    2521             :  * x = aZ_K + bZ_K using the approximation theorem */
    2522             : GEN
    2523          42 : idealtwoelt2(GEN nf, GEN x, GEN a)
    2524             : {
    2525          42 :   pari_sp av = avma;
    2526             :   GEN cx, b;
    2527             : 
    2528          42 :   nf = checknf(nf);
    2529          42 :   a = nf_to_scalar_or_basis(nf, a);
    2530          42 :   x = idealhnf_shallow(nf,x);
    2531          42 :   if (lg(x) == 1)
    2532             :   {
    2533          14 :     if (!isintzero(a)) not_in_ideal(a);
    2534           7 :     avma = av; return gen_0;
    2535             :   }
    2536          28 :   x = Q_primitive_part(x, &cx);
    2537          28 :   if (cx) a = gdiv(a, cx);
    2538          28 :   if (!in_ideal(x, a)) not_in_ideal(a);
    2539          21 :   b = mat_ideal_two_elt2(nf, x, a);
    2540          21 :   if (typ(b) == t_COL)
    2541             :   {
    2542          14 :     GEN mod = idealhnf_principal(nf,a);
    2543          14 :     b = ZC_hnfrem(b,mod);
    2544          14 :     if (ZV_isscalar(b)) b = gel(b,1);
    2545             :   }
    2546             :   else
    2547             :   {
    2548           7 :     GEN aZ = typ(a) == t_COL? Q_denom(zk_inv(nf,a)): a; /* (a) \cap Z */
    2549           7 :     b = centermodii(b, aZ, shifti(aZ,-1));
    2550             :   }
    2551          21 :   b = cx? gmul(b,cx): gcopy(b);
    2552          21 :   return gerepileupto(av, b);
    2553             : }
    2554             : 
    2555             : /* Given 2 integral ideals x and y in nf, returns a beta in nf such that
    2556             :  * beta * x is an integral ideal coprime to y */
    2557             : GEN
    2558        3444 : idealcoprimefact(GEN nf, GEN x, GEN fy)
    2559             : {
    2560        3444 :   GEN L = gel(fy,1), e;
    2561        3444 :   long i, r = lg(L);
    2562             : 
    2563        3444 :   e = cgetg(r, t_COL);
    2564        3444 :   for (i=1; i<r; i++) gel(e,i) = stoi( -idealval(nf,x,gel(L,i)) );
    2565        3444 :   return idealapprfact_i(nf, mkmat2(L,e), 0);
    2566             : }
    2567             : GEN
    2568          70 : idealcoprime(GEN nf, GEN x, GEN y)
    2569             : {
    2570          70 :   pari_sp av = avma;
    2571          70 :   return gerepileupto(av, idealcoprimefact(nf, x, idealfactor(nf,y)));
    2572             : }
    2573             : 
    2574             : GEN
    2575           7 : nfmulmodpr(GEN nf, GEN x, GEN y, GEN modpr)
    2576             : {
    2577           7 :   pari_sp av = avma;
    2578           7 :   GEN z, p, pr = modpr, T;
    2579             : 
    2580           7 :   nf = checknf(nf); modpr = nf_to_Fq_init(nf,&pr,&T,&p);
    2581           0 :   x = nf_to_Fq(nf,x,modpr);
    2582           0 :   y = nf_to_Fq(nf,y,modpr);
    2583           0 :   z = Fq_mul(x,y,T,p);
    2584           0 :   return gerepileupto(av, algtobasis(nf, Fq_to_nf(z,modpr)));
    2585             : }
    2586             : 
    2587             : GEN
    2588           0 : nfdivmodpr(GEN nf, GEN x, GEN y, GEN modpr)
    2589             : {
    2590           0 :   pari_sp av = avma;
    2591           0 :   nf = checknf(nf);
    2592           0 :   return gerepileupto(av, nfreducemodpr(nf, nfdiv(nf,x,y), modpr));
    2593             : }
    2594             : 
    2595             : GEN
    2596           0 : nfpowmodpr(GEN nf, GEN x, GEN k, GEN modpr)
    2597             : {
    2598           0 :   pari_sp av=avma;
    2599           0 :   GEN z, T, p, pr = modpr;
    2600             : 
    2601           0 :   nf = checknf(nf); modpr = nf_to_Fq_init(nf,&pr,&T,&p);
    2602           0 :   z = nf_to_Fq(nf,x,modpr);
    2603           0 :   z = Fq_pow(z,k,T,p);
    2604           0 :   return gerepileupto(av, algtobasis(nf, Fq_to_nf(z,modpr)));
    2605             : }
    2606             : 
    2607             : GEN
    2608           0 : nfkermodpr(GEN nf, GEN x, GEN modpr)
    2609             : {
    2610           0 :   pari_sp av = avma;
    2611           0 :   GEN T, p, pr = modpr;
    2612             : 
    2613           0 :   nf = checknf(nf); modpr = nf_to_Fq_init(nf, &pr,&T,&p);
    2614           0 :   if (typ(x)!=t_MAT) pari_err_TYPE("nfkermodpr",x);
    2615           0 :   x = nfM_to_FqM(x, nf, modpr);
    2616           0 :   return gerepilecopy(av, FqM_to_nfM(FqM_ker(x,T,p), modpr));
    2617             : }
    2618             : 
    2619             : GEN
    2620           0 : nfsolvemodpr(GEN nf, GEN a, GEN b, GEN pr)
    2621             : {
    2622           0 :   const char *f = "nfsolvemodpr";
    2623           0 :   pari_sp av = avma;
    2624             :   GEN T, p, modpr;
    2625             : 
    2626           0 :   nf = checknf(nf);
    2627           0 :   modpr = nf_to_Fq_init(nf, &pr,&T,&p);
    2628           0 :   if (typ(a)!=t_MAT) pari_err_TYPE(f,a);
    2629           0 :   a = nfM_to_FqM(a, nf, modpr);
    2630           0 :   switch(typ(b))
    2631             :   {
    2632             :     case t_MAT:
    2633           0 :       b = nfM_to_FqM(b, nf, modpr);
    2634           0 :       b = FqM_gauss(a,b,T,p);
    2635           0 :       if (!b) pari_err_INV(f,a);
    2636           0 :       a = FqM_to_nfM(b, modpr);
    2637           0 :       break;
    2638             :     case t_COL:
    2639           0 :       b = nfV_to_FqV(b, nf, modpr);
    2640           0 :       b = FqM_FqC_gauss(a,b,T,p);
    2641           0 :       if (!b) pari_err_INV(f,a);
    2642           0 :       a = FqV_to_nfV(b, modpr);
    2643           0 :       break;
    2644           0 :     default: pari_err_TYPE(f,b);
    2645             :   }
    2646           0 :   return gerepilecopy(av, a);
    2647             : }

Generated by: LCOV version 1.11