"Fossies" - the Fresh Open Source Software Archive

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

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

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