Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; 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 : /* RNF STRUCTURE AND OPERATIONS */
18 : /* */
19 : /*******************************************************************/
20 : #include "pari.h"
21 : #include "paripriv.h"
22 :
23 : #define DEBUGLEVEL DEBUGLEVEL_rnf
24 :
25 : /* eq is an rnfeq; must return a t_POL */
26 : GEN
27 87161 : eltreltoabs(GEN eq, GEN x)
28 : {
29 87161 : GEN Pabs = gel(eq,1), a = gel(eq,2), k = gel(eq,3), T = gel(eq,4), b, s;
30 87161 : long i, v = varn(Pabs);
31 87161 : pari_sp av = avma;
32 :
33 87161 : if (varncmp(gvar(x), v) > 0) x = scalarpol(x,v);
34 87161 : x = RgX_nffix("eltreltoabs", T, x, 1);
35 : /* Mod(X - k a, Pabs(X)) is a root of the relative polynomial */
36 87154 : if (signe(k))
37 15397 : x = RgXQX_translate(x, deg1pol_shallow(negi(k), gen_0, varn(T)), T);
38 87154 : b = pol_x(v);
39 87154 : s = gen_0;
40 258781 : for (i=lg(x)-1; i>1; i--)
41 : {
42 171627 : GEN c = gel(x,i);
43 171627 : if (typ(c) == t_POL) c = RgX_RgXQ_eval(c, a, Pabs);
44 171627 : s = RgX_rem(gadd(c, gmul(b,s)), Pabs);
45 : }
46 87154 : return gerepileupto(av, s);
47 : }
48 : GEN
49 134568 : rnfeltreltoabs(GEN rnf,GEN x)
50 : {
51 134568 : const char *f = "rnfeltreltoabs";
52 : GEN pol;
53 134568 : checkrnf(rnf);
54 134568 : pol = rnf_get_polabs(rnf);
55 134568 : switch(typ(x))
56 : {
57 24563 : case t_INT: return icopy(x);
58 1869 : case t_FRAC: return gcopy(x);
59 104720 : case t_POLMOD:
60 104720 : if (RgX_equal_var(gel(x,1), pol))
61 : { /* already in 'abs' form, unless possibly if nf = Q */
62 14805 : if (rnf_get_nfdegree(rnf) == 1)
63 : {
64 14784 : GEN y = gel(x,2);
65 14784 : pari_sp av = avma;
66 14784 : y = simplify_shallow(liftpol_shallow(y));
67 14784 : return gerepilecopy(av, mkpolmod(y, pol));
68 : }
69 21 : return gcopy(x);
70 : }
71 89915 : x = polmod_nffix(f,rnf,x,0);
72 89901 : if (typ(x) == t_POLMOD) return rnfeltup(rnf,x);
73 80917 : retmkpolmod(eltreltoabs(rnf_get_map(rnf), x), ZX_copy(pol));
74 3367 : case t_POL:
75 3367 : if (varn(x) == rnf_get_nfvarn(rnf)) return rnfeltup(rnf,x);
76 616 : retmkpolmod(eltreltoabs(rnf_get_map(rnf), x), ZX_copy(pol));
77 : }
78 49 : pari_err_TYPE(f,x);
79 : return NULL;/*LCOV_EXCL_LINE*/
80 : }
81 :
82 : GEN
83 28588 : eltabstorel_lift(GEN rnfeq, GEN P)
84 : {
85 28588 : GEN k, T = gel(rnfeq,4), R = gel(rnfeq,5);
86 28588 : if (is_scalar_t(typ(P))) return P;
87 27762 : k = gel(rnfeq,3);
88 27762 : P = lift_shallow(P);
89 27762 : if (signe(k))
90 6572 : P = RgXQX_translate(P, deg1pol_shallow(k, gen_0, varn(T)), T);
91 27762 : P = RgXQX_rem(P, R, T);
92 27762 : return QXQX_to_mod_shallow(P, T);
93 : }
94 : /* rnfeq = [pol,a,k,T,R], P a t_POL or scalar
95 : * Return Mod(P(x + k Mod(y, T(y))), pol(x)) */
96 : GEN
97 25837 : eltabstorel(GEN rnfeq, GEN P)
98 : {
99 25837 : GEN T = gel(rnfeq,4), R = gel(rnfeq,5);
100 25837 : return mkpolmod(eltabstorel_lift(rnfeq,P), QXQX_to_mod_shallow(R,T));
101 : }
102 : GEN
103 54838 : rnfeltabstorel(GEN rnf,GEN x)
104 : {
105 54838 : const char *f = "rnfeltabstorel";
106 54838 : pari_sp av = avma;
107 : GEN pol, T, P, NF;
108 54838 : checkrnf(rnf);
109 54838 : T = rnf_get_nfpol(rnf);
110 54838 : P = rnf_get_pol(rnf);
111 54838 : pol = rnf_get_polabs(rnf);
112 54838 : switch(typ(x))
113 : {
114 4914 : case t_INT: return icopy(x);
115 63 : case t_FRAC: return gcopy(x);
116 48335 : case t_POLMOD:
117 48335 : if (RgX_equal_var(P, gel(x,1)))
118 : {
119 14924 : x = polmod_nffix(f, rnf, x, 0);
120 14924 : P = QXQX_to_mod_shallow(P,T);
121 14924 : return gerepilecopy(av, mkpolmod(x,P));
122 : }
123 33411 : if (RgX_equal_var(T, gel(x,1))) { x = Rg_nffix(f, T, x, 0); goto END; }
124 33327 : if (!RgX_equal_var(pol, gel(x,1))) pari_err_MODULUS(f, gel(x,1),pol);
125 33285 : x = gel(x,2); break;
126 1085 : case t_POL: break;
127 441 : case t_COL:
128 441 : NF = obj_check(rnf, rnf_NFABS);
129 441 : if (!NF) pari_err_TYPE("rnfeltabstorel, apply nfinit(rnf)",x);
130 294 : x = nf_to_scalar_or_alg(NF,x); break;
131 0 : default:
132 0 : pari_err_TYPE(f,x);
133 : return NULL;/*LCOV_EXCL_LINE*/
134 : }
135 34664 : switch(typ(x))
136 : {
137 13027 : case t_INT: return icopy(x);
138 70 : case t_FRAC: return gcopy(x);
139 21567 : case t_POL: break;
140 0 : default: pari_err_TYPE(f, x);
141 : }
142 21567 : RgX_check_QX(x,f);
143 21490 : if (varn(x) != varn(pol))
144 : {
145 70 : if (varn(x) == varn(T)) { x = Rg_nffix(f,T,x,0); goto END; }
146 28 : pari_err_VAR(f, x,pol);
147 : }
148 21420 : switch(lg(x))
149 : {
150 0 : case 2: set_avma(av); return gen_0;
151 70 : case 3: return gerepilecopy(av, gel(x,2));
152 : }
153 21476 : END:
154 21476 : return gerepilecopy(av, eltabstorel(rnf_get_map(rnf), x));
155 : }
156 :
157 : /* x a t_VEC of rnf elements in 'alg' form (t_POL). Assume maximal rank or 0 */
158 : static GEN
159 2247 : modulereltoabs(GEN rnf, GEN x)
160 : {
161 2247 : GEN W=gel(x,1), I=gel(x,2), rnfeq = rnf_get_map(rnf), polabs = gel(rnfeq,1);
162 2247 : long i, j, k, m, N = lg(W)-1;
163 : GEN zknf, dzknf, M;
164 :
165 2247 : if (!N) return cgetg(1, t_VEC);
166 2184 : zknf = rnf_get_nfzk(rnf);
167 2184 : dzknf = gel(zknf,1);
168 2184 : m = rnf_get_nfdegree(rnf);
169 2184 : M = cgetg(N*m+1, t_VEC);
170 7861 : for (k=i=1; i<=N; i++)
171 : {
172 5677 : GEN c0, cid, w = gel(W,i), id = gel(I,i);
173 :
174 5677 : if (lg(id) == 1) continue; /* must be a t_MAT */
175 5628 : id = Q_primitive_part(id, &cid);
176 5628 : w = Q_primitive_part(eltreltoabs(rnfeq,w), &c0);
177 5628 : c0 = div_content(mul_content(c0,cid), dzknf);
178 5628 : if (typ(id) == t_INT)
179 11993 : for (j=1; j<=m; j++)
180 : {
181 7946 : GEN z = RgX_rem(gmul(w, gel(zknf,j)), polabs);
182 7946 : if (c0) z = RgX_Rg_mul(z, c0);
183 7946 : gel(M,k++) = z;
184 : }
185 : else
186 5493 : for (j=1; j<=m; j++)
187 : {
188 3912 : GEN c, z = Q_primitive_part(RgV_RgC_mul(zknf,gel(id,j)), &c);
189 3912 : z = RgX_rem(gmul(w, z), polabs);
190 3912 : c = mul_content(c, c0); if (c) z = RgX_Rg_mul(z, c);
191 3912 : gel(M,k++) = z;
192 : }
193 : }
194 2184 : setlg(M, k); return M;
195 : }
196 :
197 : /* Z-basis for absolute maximal order: [NF.pol, NF.zk] */
198 : GEN
199 1785 : rnf_zkabs(GEN rnf)
200 : {
201 1785 : GEN d, v, M = modulereltoabs(rnf, rnf_get_zk(rnf));
202 1785 : GEN T = rnf_get_polabs(rnf);
203 1785 : long n = degpol(T);
204 1785 : M = Q_remove_denom(M, &d); /* t_VEC of t_POL */
205 1785 : if (d)
206 : {
207 1246 : M = RgXV_to_RgM(M,n);
208 1246 : M = ZM_hnfmodall(M, d, hnf_MODID|hnf_CENTER);
209 1246 : M = RgM_Rg_div(M, d);
210 : }
211 : else
212 539 : M = matid(n);
213 1785 : v = rnf_get_ramified_primes(rnf);
214 1785 : if (lg(v) == 1)
215 : {
216 175 : GEN D = gel(rnf_get_disc(rnf),1);
217 175 : if (!isint1(D)) pari_err_TYPE("rnf_zkabs (old style rnf)", rnf);
218 : }
219 1785 : v = shallowconcat(nf_get_ramified_primes(rnf_get_nf(rnf)), v);
220 1785 : return mkvec3(T, RgM_to_RgXV(M, varn(T)), ZV_sort_uniq_shallow(v));
221 : }
222 :
223 : static GEN
224 1561 : mknfabs(GEN rnf, long prec)
225 : {
226 : GEN NF;
227 1561 : if ((NF = obj_check(rnf,rnf_NFABS)))
228 0 : { if (nf_get_prec(NF) < prec) NF = nfnewprec_shallow(NF,prec); }
229 : else
230 1561 : NF = nfinit(rnf_zkabs(rnf), prec);
231 1561 : return NF;
232 : }
233 :
234 : static GEN
235 1561 : mkupdown(GEN NF, GEN rnf)
236 : {
237 : GEN M, zknf, dzknf;
238 : long i, l;
239 1561 : zknf = rnf_get_nfzk(rnf);
240 1561 : dzknf = gel(zknf,1); if (gequal1(dzknf)) dzknf = NULL;
241 1561 : l = lg(zknf); M = cgetg(l, t_MAT);
242 1561 : gel(M,1) = vec_ei(nf_get_degree(NF), 1);
243 3052 : for (i = 2; i < l; i++)
244 : {
245 1491 : GEN c = poltobasis(NF, gel(zknf,i));
246 1491 : if (dzknf) c = gdiv(c, dzknf);
247 1491 : gel(M,i) = c;
248 : }
249 1561 : return Qevproj_init(M);
250 : }
251 : GEN
252 103921 : rnf_build_nfabs(GEN rnf, long prec)
253 : {
254 103921 : GEN NF = obj_checkbuild_prec(rnf, rnf_NFABS, &mknfabs, &nf_get_prec, prec);
255 103921 : GEN O = obj_check(rnf, rnf_MAPS);
256 103921 : if (!O)
257 1561 : { pari_sp av = avma; O = obj_insert(rnf, rnf_MAPS, mkupdown(NF, rnf)); set_avma(av); }
258 103921 : return NF;
259 : }
260 :
261 : void
262 32724 : rnfcomplete(GEN rnf)
263 32724 : { (void)rnf_build_nfabs(rnf, nf_get_prec(rnf_get_nf(rnf))); }
264 :
265 : GEN
266 2023 : nf_nfzk(GEN nf, GEN rnfeq)
267 : {
268 2023 : GEN pol = gel(rnfeq,1), a = gel(rnfeq,2);
269 2023 : return Q_primpart(QXV_QXQ_eval(nf_get_zkprimpart(nf), a, pol));
270 : }
271 :
272 : static GEN
273 3647 : rnfdisc_get_T_i(GEN P, GEN *lim)
274 : {
275 3647 : *lim = NULL;
276 3647 : if (typ(P) == t_VEC && lg(P) == 3)
277 : {
278 357 : GEN L = gel(P,2);
279 : long i, l;
280 357 : *lim = L;
281 357 : switch(typ(L))
282 : {
283 56 : case t_INT:
284 56 : if (signe(L) <= 0) return NULL;
285 56 : break;
286 301 : case t_VEC: case t_COL:
287 301 : l = lg(L);
288 917 : for (i = 1; i < l; i++)
289 : {
290 616 : GEN p = gel(L,i);
291 616 : if (typ(p) == t_INT)
292 602 : { if (signe(p) <= 0) return NULL; }
293 14 : else checkprid(p);
294 : }
295 301 : break;
296 0 : default: return NULL;
297 : }
298 357 : P = gel(P,1);
299 : }
300 3647 : return (typ(P) == t_POL)? P: NULL;
301 : }
302 : /* true nf */
303 : GEN
304 3647 : rnfdisc_get_T(GEN nf, GEN P, GEN *lim)
305 : {
306 3647 : GEN T = rnfdisc_get_T_i(P, lim);
307 3647 : if (!T) pari_err_TYPE("rnfdisc",P);
308 3647 : return RgX_nffix("rnfdisc", nf_get_pol(nf), T, 1);
309 : }
310 :
311 : GEN
312 126 : rnfpseudobasis(GEN nf, GEN pol)
313 : {
314 126 : pari_sp av = avma;
315 : GEN D, z, lim;
316 126 : nf = checknf(nf);
317 126 : pol = rnfdisc_get_T(nf, pol, &lim);
318 126 : z = rnfallbase(nf, pol, lim, NULL, &D, NULL, NULL);
319 105 : return gerepilecopy(av, shallowconcat(z,D));
320 : }
321 :
322 : GEN
323 1974 : rnfinit0(GEN nf, GEN T, long flag)
324 : {
325 1974 : pari_sp av = avma;
326 1974 : GEN lim, bas, D, f, B, DKP, rnfeq, rnf = obj_init(11, 2);
327 1974 : nf = checknf(nf);
328 1974 : T = rnfdisc_get_T(nf, T, &lim);
329 1974 : gel(rnf,11) = rnfeq = nf_rnfeq(nf,T);
330 1967 : gel(rnf,2) = nf_nfzk(nf, rnfeq);
331 1967 : bas = rnfallbase(nf, T, lim, rnf, &D, &f, &DKP);
332 1967 : B = matbasistoalg(nf,gel(bas,1));
333 1967 : gel(bas,1) = lift_if_rational( RgM_to_RgXV(B,varn(T)) );
334 1967 : gel(rnf,1) = T;
335 1967 : gel(rnf,3) = D;
336 1967 : gel(rnf,4) = f;
337 1967 : gel(rnf,5) = DKP;
338 1967 : gel(rnf,6) = cgetg(1, t_VEC); /* dummy */
339 1967 : gel(rnf,7) = bas;
340 1967 : gel(rnf,8) = lift_if_rational( RgM_inv_upper(B) );
341 1282 : gel(rnf,9) = typ(f) == t_INT? powiu(f, nf_get_degree(nf))
342 1967 : : RgM_det_triangular(f);
343 1967 : gel(rnf,10)= nf;
344 1967 : rnf = gerepilecopy(av, rnf);
345 1967 : if (flag) rnfcomplete(rnf);
346 1967 : return rnf;
347 : }
348 : GEN
349 721 : rnfinit(GEN nf, GEN T) { return rnfinit0(nf,T,0); }
350 :
351 : GEN
352 27873 : rnfeltup0(GEN rnf, GEN x, long flag)
353 : {
354 27873 : pari_sp av = avma;
355 : GEN zknf, nf, NF, POL;
356 27873 : long tx = typ(x);
357 27873 : checkrnf(rnf);
358 27873 : if (flag) rnfcomplete(rnf);
359 27873 : NF = obj_check(rnf,rnf_NFABS);
360 27873 : POL = rnf_get_polabs(rnf);
361 27873 : if (tx == t_POLMOD && RgX_equal_var(gel(x,1), POL))
362 : {
363 42 : if (flag) x = nf_to_scalar_or_basis(NF,x);
364 42 : return gerepilecopy(av, x);
365 : }
366 27831 : nf = rnf_get_nf(rnf);
367 27831 : if (NF && tx == t_COL && lg(x)-1 == degpol(POL) && nf_get_degree(rnf) > 1)
368 : {
369 0 : x = flag? nf_to_scalar_or_basis(NF,x)
370 0 : : mkpolmod(nf_to_scalar_or_alg(NF,x), POL);
371 0 : return gerepilecopy(av, x);
372 : }
373 27831 : if (NF)
374 : {
375 : GEN d, proj;
376 27614 : x = nf_to_scalar_or_basis(nf, x);
377 27614 : if (typ(x) != t_COL) return gerepilecopy(av, x);
378 26774 : proj = obj_check(rnf,rnf_MAPS);
379 26774 : x = Q_remove_denom(x,&d);
380 26774 : x = ZM_ZC_mul(gel(proj,1), x);
381 26774 : if (d) x = gdiv(x,d);
382 26774 : if (!flag) x = basistoalg(NF,x);
383 : }
384 : else
385 : {
386 217 : zknf = rnf_get_nfzk(rnf);
387 217 : x = nfeltup(nf, x, zknf);
388 112 : if (typ(x) == t_POL) x = mkpolmod(x, POL);
389 : }
390 26886 : return gerepilecopy(av, x);
391 : }
392 : GEN
393 13450 : rnfeltup(GEN rnf, GEN x) { return rnfeltup0(rnf,x,0); }
394 :
395 : GEN
396 245 : nfeltup(GEN nf, GEN x, GEN zknf)
397 : {
398 245 : GEN c, dzknf = gel(zknf,1);
399 245 : x = nf_to_scalar_or_basis(nf, x);
400 140 : if (typ(x) != t_COL) return x;
401 42 : x = Q_primitive_part(x, &c);
402 42 : if (!RgV_is_ZV(x)) pari_err_TYPE("rnfeltup", x);
403 42 : if (gequal1(dzknf)) dzknf = NULL;
404 42 : c = div_content(c, dzknf);
405 42 : x = RgV_RgC_mul(zknf, x); if (c) x = RgX_Rg_mul(x, c);
406 42 : return x;
407 : }
408 :
409 : static void
410 49 : fail(const char *f, GEN x)
411 49 : { pari_err_DOMAIN(f,"element","not in", strtoGENstr("the base field"),x); }
412 : /* x t_COL of length degabs */
413 : static GEN
414 0 : eltdown(GEN rnf, GEN x, long flag)
415 : {
416 0 : GEN y, d, proj = obj_check(rnf,rnf_MAPS);
417 0 : GEN M = gel(proj,1), iM = gel(proj,2), diM = gel(proj,3), perm = gel(proj,4);
418 0 : x = Q_remove_denom(x,&d);
419 0 : if (!RgV_is_ZV(x)) pari_err_TYPE("rnfeltdown", x);
420 0 : y = ZM_ZC_mul(iM, vecpermute(x, perm));
421 0 : if (!ZV_equal(ZM_ZC_mul(M,y),
422 0 : isint1(diM)? x: ZC_Z_mul(x,diM))) fail("rnfeltdown",x);
423 :
424 0 : d = mul_denom(d, diM);
425 0 : if (d) y = gdiv(y,d);
426 0 : if (!flag) y = basistoalg(rnf_get_nf(rnf), y);
427 0 : return y;
428 : }
429 : GEN
430 2716 : rnfeltdown0(GEN rnf, GEN x, long flag)
431 : {
432 2716 : const char *f = "rnfeltdown";
433 2716 : pari_sp av = avma;
434 : GEN z, T, NF, nf;
435 : long v;
436 :
437 2716 : checkrnf(rnf);
438 2716 : NF = obj_check(rnf,rnf_NFABS);
439 2716 : nf = rnf_get_nf(rnf);
440 2716 : T = nf_get_pol(nf);
441 2716 : v = varn(T);
442 2716 : switch(typ(x))
443 : { /* directly belonging to base field ? */
444 567 : case t_INT: return icopy(x);
445 84 : case t_FRAC:return gcopy(x);
446 1883 : case t_POLMOD:
447 1883 : if (RgX_equal_var(gel(x,1), rnf_get_polabs(rnf)))
448 : {
449 210 : if (degpol(T) == 1)
450 : {
451 182 : x = simplify_shallow(liftpol_shallow(gel(x,2)));
452 182 : if (typ(x) != t_POL) return gerepilecopy(av,x);
453 : }
454 49 : break;
455 : }
456 1673 : x = polmod_nffix(f,rnf,x,0);
457 : /* x was defined mod the relative polynomial & non constant => fail */
458 1659 : if (typ(x) == t_POL) fail(f,x);
459 1652 : if (flag) x = nf_to_scalar_or_basis(nf,x);
460 1652 : return gerepilecopy(av, x);
461 :
462 133 : case t_POL:
463 133 : if (varn(x) != v) break;
464 91 : x = Rg_nffix(f,T,x,0);
465 84 : if (flag) x = nf_to_scalar_or_basis(nf,x);
466 84 : return gerepilecopy(av, x);
467 49 : case t_COL:
468 : {
469 49 : long n = lg(x)-1;
470 49 : if (n == degpol(T) && RgV_is_QV(x))
471 : {
472 7 : if (RgV_isscalar(x)) return gcopy(gel(x,1));
473 0 : if (!flag) return gcopy(x);
474 0 : return basistoalg(nf,x);
475 : }
476 42 : if (NF) break;
477 : }
478 42 : default: pari_err_TYPE(f, x);
479 : }
480 : /* x defined mod the absolute equation */
481 91 : if (NF)
482 : {
483 0 : x = nf_to_scalar_or_basis(NF, x);
484 0 : if (typ(x) == t_COL) x = eltdown(rnf,x,flag);
485 0 : return gerepilecopy(av, x);
486 : }
487 91 : z = rnfeltabstorel(rnf,x);
488 70 : switch(typ(z))
489 : {
490 14 : case t_INT:
491 14 : case t_FRAC: return z;
492 : }
493 : /* typ(z) = t_POLMOD, varn of both components is rnf_get_varn(rnf) */
494 56 : z = gel(z,2);
495 56 : if (typ(z) == t_POL)
496 : {
497 56 : if (lg(z) != 3) fail(f,x);
498 14 : z = gel(z,2);
499 : }
500 14 : return gerepilecopy(av, z);
501 : }
502 : GEN
503 2471 : rnfeltdown(GEN rnf, GEN x) { return rnfeltdown0(rnf,x,0); }
504 :
505 : /* vector of rnf elt -> matrix of nf elts */
506 : static GEN
507 483 : rnfV_to_nfM(GEN rnf, GEN x)
508 : {
509 483 : long i, l = lg(x);
510 483 : GEN y = cgetg(l, t_MAT);
511 1463 : for (i = 1; i < l; i++) gel(y,i) = rnfalgtobasis(rnf,gel(x,i));
512 483 : return y;
513 : }
514 :
515 : static GEN
516 770 : rnfprincipaltohnf(GEN rnf,GEN x)
517 : {
518 770 : pari_sp av = avma;
519 770 : GEN bas = rnf_get_zk(rnf), nf = rnf_get_nf(rnf);
520 770 : x = rnfbasistoalg(rnf,x);
521 434 : x = gmul(x, gmodulo(gel(bas,1), rnf_get_pol(rnf)));
522 434 : return gerepileupto(av, nfhnf(nf, mkvec2(rnfV_to_nfM(rnf,x), gel(bas,2))));
523 : }
524 :
525 : /* pseudo-basis for the 0 ideal */
526 : static GEN
527 154 : rnfideal0(void) { retmkvec2(cgetg(1,t_MAT),cgetg(1,t_VEC)); }
528 :
529 : GEN
530 1330 : rnfidealhnf(GEN rnf, GEN x)
531 : {
532 : GEN z, nf, bas;
533 :
534 1330 : checkrnf(rnf); nf = rnf_get_nf(rnf);
535 1330 : switch(typ(x))
536 : {
537 182 : case t_INT: case t_FRAC:
538 182 : if (isintzero(x)) return rnfideal0();
539 126 : bas = rnf_get_zk(rnf); z = cgetg(3,t_VEC);
540 126 : gel(z,1) = matid(rnf_get_degree(rnf));
541 126 : gel(z,2) = gmul(x, gel(bas,2)); return z;
542 :
543 266 : case t_VEC:
544 266 : if (lg(x) == 3 && typ(gel(x,1)) == t_MAT) return nfhnf(nf, x);
545 : case t_MAT:
546 252 : return rnfidealabstorel(rnf, x);
547 :
548 770 : case t_POLMOD: case t_POL: case t_COL:
549 770 : return rnfprincipaltohnf(rnf,x);
550 : }
551 0 : pari_err_TYPE("rnfidealhnf",x);
552 : return NULL; /* LCOV_EXCL_LINE */
553 : }
554 :
555 : static GEN
556 105 : prodidnorm(GEN nf, GEN I)
557 : {
558 105 : long i, l = lg(I);
559 : GEN z;
560 105 : if (l == 1) return gen_1;
561 105 : z = idealnorm(nf, gel(I,1));
562 210 : for (i=2; i<l; i++) z = gmul(z, idealnorm(nf, gel(I,i)));
563 105 : return z;
564 : }
565 :
566 : GEN
567 196 : rnfidealnormrel(GEN rnf, GEN id)
568 : {
569 196 : pari_sp av = avma;
570 196 : GEN nf, z = gel(rnfidealhnf(rnf,id), 2);
571 126 : if (lg(z) == 1) return cgetg(1, t_MAT);
572 98 : nf = rnf_get_nf(rnf); z = idealprod(nf, z);
573 98 : return gerepileupto(av, idealmul(nf,z, rnf_get_index(rnf)));
574 : }
575 :
576 : GEN
577 203 : rnfidealnormabs(GEN rnf, GEN id)
578 : {
579 203 : pari_sp av = avma;
580 203 : GEN nf, z = gel(rnfidealhnf(rnf,id), 2);
581 133 : if (lg(z) == 1) return gen_0;
582 105 : nf = rnf_get_nf(rnf); z = prodidnorm(nf, z);
583 105 : return gerepileupto(av, gmul(z, gel(rnf,9)));
584 : }
585 :
586 : static GEN
587 497 : rnfidealreltoabs_i(GEN rnf, GEN x)
588 : {
589 : long i, l;
590 : GEN w;
591 497 : x = rnfidealhnf(rnf,x);
592 357 : w = gel(x,1); l = lg(w); settyp(w, t_VEC);
593 966 : for (i=1; i<l; i++) gel(w,i) = lift_shallow( rnfbasistoalg(rnf, gel(w,i)) );
594 357 : return modulereltoabs(rnf, x);
595 : }
596 : GEN
597 0 : rnfidealreltoabs(GEN rnf, GEN x)
598 : {
599 0 : pari_sp av = avma;
600 0 : return gerepilecopy(av, rnfidealreltoabs_i(rnf,x));
601 : }
602 : GEN
603 238 : rnfidealreltoabs0(GEN rnf, GEN x, long flag)
604 : {
605 238 : pari_sp av = avma;
606 : long i, l;
607 : GEN NF;
608 :
609 238 : x = rnfidealreltoabs_i(rnf, x);
610 168 : if (!flag) return gerepilecopy(av,x);
611 35 : rnfcomplete(rnf);
612 35 : NF = obj_check(rnf,rnf_NFABS);
613 35 : l = lg(x); settyp(x, t_MAT);
614 245 : for (i=1; i<l; i++) gel(x,i) = algtobasis(NF, gel(x,i));
615 35 : return gerepileupto(av, idealhnf(NF,x));
616 : }
617 :
618 : GEN
619 455 : rnfidealabstorel(GEN rnf, GEN x)
620 : {
621 455 : long n, N, j, tx = typ(x);
622 455 : pari_sp av = avma;
623 : GEN A, I, invbas;
624 :
625 455 : checkrnf(rnf);
626 455 : invbas = rnf_get_invzk(rnf);
627 455 : if (tx != t_VEC && tx != t_MAT) pari_err_TYPE("rnfidealabstorel",x);
628 315 : N = lg(x)-1;
629 315 : if (N != rnf_get_absdegree(rnf))
630 : {
631 196 : if (!N) return rnfideal0();
632 105 : pari_err_DIM("rnfidealabstorel");
633 : }
634 119 : n = rnf_get_degree(rnf);
635 119 : A = cgetg(N+1,t_MAT);
636 119 : I = cgetg(N+1,t_VEC);
637 833 : for (j=1; j<=N; j++)
638 : {
639 714 : GEN t = lift_shallow( rnfeltabstorel(rnf, gel(x,j)) );
640 714 : if (typ(t) == t_POL)
641 595 : t = RgM_RgX_mul(invbas, t);
642 : else
643 119 : t = scalarcol_shallow(t, n);
644 714 : gel(A,j) = t;
645 714 : gel(I,j) = gen_1;
646 : }
647 119 : return gerepileupto(av, nfhnf(rnf_get_nf(rnf), mkvec2(A,I)));
648 : }
649 :
650 : GEN
651 217 : rnfidealdown(GEN rnf,GEN x)
652 : {
653 217 : pari_sp av = avma;
654 : GEN I;
655 217 : if (typ(x) == t_MAT)
656 : {
657 : GEN d;
658 28 : x = Q_remove_denom(x,&d);
659 28 : if (RgM_is_ZM(x))
660 : {
661 28 : GEN NF = obj_check(rnf,rnf_NFABS);
662 28 : if (NF)
663 : {
664 28 : GEN z, proj = obj_check(rnf,rnf_MAPS), ZK = gel(proj,1);
665 : long i, lz, l;
666 28 : x = idealhnf_shallow(NF,x);
667 35 : if (lg(x) == 1) { set_avma(av); return cgetg(1,t_MAT); }
668 14 : z = ZM_lll(shallowconcat(ZK,x), 0.99, LLL_KER);
669 14 : lz = lg(z); l = lg(ZK);
670 56 : for (i = 1; i < lz; i++) setlg(gel(z,i), l);
671 14 : z = ZM_hnfmodid(z, gcoeff(x,1,1));
672 14 : if (d) z = gdiv(z,d);
673 14 : return gerepileupto(av, z);
674 : }
675 : }
676 : }
677 189 : x = rnfidealhnf(rnf,x); I = gel(x,2);
678 126 : if (lg(I) == 1) { set_avma(av); return cgetg(1,t_MAT); }
679 105 : return gerepilecopy(av, gel(I,1));
680 : }
681 :
682 : /* lift ideal x to the relative extension, returns a Z-basis */
683 : GEN
684 224 : rnfidealup(GEN rnf,GEN x)
685 : {
686 224 : pari_sp av = avma;
687 : long i, n;
688 : GEN nf, bas, bas2, I, x2, dx;
689 :
690 224 : checkrnf(rnf); nf = rnf_get_nf(rnf);
691 224 : n = rnf_get_degree(rnf);
692 224 : bas = rnf_get_zk(rnf); bas2 = gel(bas,2);
693 :
694 224 : (void)idealtyp(&x, NULL);
695 210 : x = Q_remove_denom(x, &dx);
696 210 : x2 = idealtwoelt(nf,x);
697 105 : I = cgetg(n+1,t_VEC);
698 308 : for (i=1; i<=n; i++)
699 : {
700 203 : GEN c = gel(bas2,i), d;
701 203 : if (typ(c) == t_MAT)
702 : {
703 7 : c = Q_remove_denom(c,&d);
704 7 : d = mul_denom(d, dx);
705 7 : c = idealHNF_mul(nf,c,x2);
706 : }
707 : else
708 : {
709 196 : c = idealmul(nf,c,x);
710 196 : d = dx;
711 : }
712 203 : if (d) c = gdiv(c,d);
713 203 : gel(I,i) = c;
714 : }
715 105 : return gerepilecopy(av, modulereltoabs(rnf, mkvec2(gel(bas,1), I)));
716 : }
717 : GEN
718 245 : rnfidealup0(GEN rnf,GEN x, long flag)
719 : {
720 245 : pari_sp av = avma;
721 : GEN NF, nf, proj, d, x2;
722 :
723 245 : if (!flag) return rnfidealup(rnf,x);
724 21 : checkrnf(rnf); nf = rnf_get_nf(rnf);
725 21 : rnfcomplete(rnf);
726 21 : proj = obj_check(rnf,rnf_MAPS);
727 21 : NF = obj_check(rnf,rnf_NFABS);
728 :
729 21 : (void)idealtyp(&x, NULL);
730 21 : x2 = idealtwoelt(nf,x);
731 21 : x2 = Q_remove_denom(x2,&d);
732 21 : if (typ(gel(x2,2)) == t_COL) gel(x2,2) = ZM_ZC_mul(gel(proj,1),gel(x2,2));
733 21 : x2 = idealhnf_two(NF, x2);
734 21 : if (d) x2 = gdiv(x2,d);
735 21 : return gerepileupto(av, x2);
736 : }
737 :
738 : /* x a relative HNF => vector of 2 generators (relative polmods) */
739 : GEN
740 259 : rnfidealtwoelement(GEN rnf, GEN x)
741 : {
742 259 : pari_sp av = avma;
743 : GEN y, cy, z, NF;
744 :
745 259 : y = rnfidealreltoabs_i(rnf,x);
746 189 : rnfcomplete(rnf);
747 189 : NF = obj_check(rnf,rnf_NFABS);
748 189 : y = matalgtobasis(NF, y); settyp(y, t_MAT);
749 189 : y = Q_primitive_part(y, &cy);
750 189 : y = ZM_hnf(y);
751 189 : if (lg(y) == 1) { set_avma(av); return mkvec2(gen_0, gen_0); }
752 154 : y = idealtwoelt(NF, y);
753 147 : if (cy) y = RgV_Rg_mul(y, cy);
754 147 : z = gel(y,2);
755 147 : if (typ(z) == t_COL) z = rnfeltabstorel(rnf, nf_to_scalar_or_alg(NF, z));
756 147 : return gerepilecopy(av, mkvec2(gel(y,1), z));
757 : }
758 :
759 : GEN
760 56 : rnfidealmul(GEN rnf,GEN x,GEN y)
761 : {
762 56 : pari_sp av = avma;
763 : GEN nf, z, x1, x2, p1, p2, bas;
764 :
765 56 : y = rnfidealtwoelement(rnf,y);
766 56 : if (isintzero(gel(y,1))) { set_avma(av); return rnfideal0(); }
767 49 : nf = rnf_get_nf(rnf);
768 49 : bas = rnf_get_zk(rnf);
769 49 : x = rnfidealhnf(rnf,x);
770 49 : x1 = gmodulo(gmul(gel(bas,1), matbasistoalg(nf,gel(x,1))), rnf_get_pol(rnf));
771 49 : x2 = gel(x,2);
772 49 : p1 = gmul(gel(y,1), gel(x,1));
773 49 : p2 = rnfV_to_nfM(rnf, gmul(gel(y,2), x1));
774 49 : z = mkvec2(shallowconcat(p1, p2), shallowconcat(x2, x2));
775 49 : return gerepileupto(av, nfhnf(nf,z));
776 : }
777 :
778 : /* prK wrt NF ~ Q[x]/(polabs) */
779 : static GEN
780 17384 : rnfidealprimedec_1(GEN rnf, GEN SL, GEN prK)
781 : {
782 17384 : GEN v, piL, piK = pr_get_gen(prK);
783 : long i, c, l;
784 17384 : if (pr_is_inert(prK)) return SL;
785 14192 : piL = rnfeltup0(rnf, piK, 1);
786 14192 : v = cgetg_copy(SL, &l);
787 60926 : for (i = c = 1; i < l; i++)
788 : {
789 46734 : GEN P = gel(SL,i);
790 46734 : if (ZC_prdvd(piL, P)) gel(v,c++) = P;
791 : }
792 14192 : setlg(v, c); return v;
793 : }
794 : GEN
795 17384 : rnfidealprimedec(GEN rnf, GEN pr)
796 : {
797 17384 : pari_sp av = avma;
798 : GEN p, z, NF, nf, SL;
799 17384 : checkrnf(rnf);
800 17384 : rnfcomplete(rnf);
801 17384 : NF = obj_check(rnf,rnf_NFABS);
802 17384 : nf = rnf_get_nf(rnf);
803 17384 : if (typ(pr) == t_INT) { p = pr; pr = NULL; }
804 17356 : else { checkprid(pr); p = pr_get_p(pr); }
805 17384 : SL = idealprimedec(NF, p);
806 17384 : if (pr) z = rnfidealprimedec_1(rnf, SL, pr);
807 : else
808 : {
809 28 : GEN vK = idealprimedec(nf, p), vL;
810 28 : long l = lg(vK), i;
811 28 : vL = cgetg(l, t_VEC);
812 56 : for (i = 1; i < l; i++) gel(vL,i) = rnfidealprimedec_1(rnf, SL, gel(vK,i));
813 28 : z = mkvec2(vK, vL);
814 : }
815 17384 : return gerepilecopy(av, z);
816 : }
817 :
818 : GEN
819 35 : rnfidealfactor(GEN rnf, GEN x)
820 : {
821 35 : pari_sp av = avma;
822 : GEN NF;
823 35 : checkrnf(rnf);
824 35 : rnfcomplete(rnf);
825 35 : NF = obj_check(rnf,rnf_NFABS);
826 35 : return gerepileupto(av, idealfactor(NF, rnfidealreltoabs0(rnf, x, 1)));
827 : }
828 :
829 : GEN
830 47425 : rnfequationall(GEN A, GEN B, long *pk, GEN *pLPRS)
831 : {
832 : long lA, lB;
833 : GEN nf, C;
834 :
835 47425 : A = get_nfpol(A, &nf); lA = lg(A);
836 47425 : if (!nf) {
837 9254 : if (lA<=3) pari_err_CONSTPOL("rnfequation");
838 9254 : RgX_check_ZX(A,"rnfequation");
839 : }
840 47425 : B = RgX_nffix("rnfequation", A,B,1); lB = lg(B);
841 47425 : if (lB<=3) pari_err_CONSTPOL("rnfequation");
842 47425 : B = Q_primpart(B);
843 :
844 47425 : if (!nfissquarefree(A,B))
845 7 : pari_err_DOMAIN("rnfequation","issquarefree(B)","=",gen_0,B);
846 :
847 47417 : *pk = 0; C = ZX_ZXY_resultant_all(A, B, pk, pLPRS);
848 47418 : if (signe(leading_coeff(C)) < 0) C = ZX_neg(C);
849 47418 : *pk = -*pk; return Q_primpart(C);
850 : }
851 :
852 : GEN
853 47047 : rnfequation0(GEN A, GEN B, long flall)
854 : {
855 47047 : pari_sp av = avma;
856 : GEN LPRS, C;
857 : long k;
858 :
859 47047 : C = rnfequationall(A, B, &k, flall? &LPRS: NULL);
860 47039 : if (flall)
861 : { /* a,b,c root of A,B,C = compositum, c = b + k a */
862 9737 : GEN a, mH0 = RgX_neg(gel(LPRS,1)), H1 = gel(LPRS,2);
863 9737 : a = QXQ_div(mH0, H1, C);
864 9737 : C = mkvec3(C, mkpolmod(a, C), stoi(k));
865 : }
866 47039 : return gerepilecopy(av, C);
867 : }
868 : GEN
869 34111 : rnfequation(GEN nf, GEN pol) { return rnfequation0(nf,pol,0); }
870 : GEN
871 9632 : rnfequation2(GEN nf, GEN pol) { return rnfequation0(nf,pol,1); }
872 : GEN
873 2037 : nf_rnfeq(GEN nf, GEN R)
874 : {
875 : GEN pol, a, k, junk, eq;
876 2037 : R = liftpol_shallow(R);
877 2037 : eq = rnfequation2(nf, R);
878 2030 : pol = gel(eq,1);
879 2030 : a = gel(eq,2); if (typ(a) == t_POLMOD) a = gel(a,2);
880 2030 : k = gel(eq,3);
881 2030 : return mkvec5(pol,a,k,get_nfpol(nf, &junk),R);
882 : }
883 : /* only allow abstorel */
884 : GEN
885 378 : nf_rnfeqsimple(GEN nf, GEN R)
886 : {
887 : long sa;
888 : GEN junk, pol;
889 378 : R = liftpol_shallow(R);
890 378 : pol = rnfequationall(nf, R, &sa, NULL);
891 378 : return mkvec5(pol,gen_0/*dummy*/,stoi(sa),get_nfpol(nf, &junk),R);
892 : }
893 :
894 : /*******************************************************************/
895 : /* */
896 : /* RELATIVE LLL */
897 : /* */
898 : /*******************************************************************/
899 : static GEN
900 196 : nftau(long r1, GEN x)
901 : {
902 196 : long i, l = lg(x);
903 196 : GEN s = r1? gel(x,1): gmul2n(real_i(gel(x,1)),1);
904 392 : for (i=2; i<=r1; i++) s = gadd(s, gel(x,i));
905 196 : for ( ; i < l; i++) s = gadd(s, gmul2n(real_i(gel(x,i)),1));
906 196 : return s;
907 : }
908 :
909 : static GEN
910 28 : initmat(long l)
911 : {
912 28 : GEN x = cgetg(l, t_MAT);
913 : long i;
914 196 : for (i = 1; i < l; i++) gel(x,i) = cgetg(l, t_COL);
915 28 : return x;
916 : }
917 :
918 : static GEN
919 1022 : nftocomplex(GEN nf, GEN x)
920 : {
921 1022 : GEN M = nf_get_M(nf);
922 1022 : x = nf_to_scalar_or_basis(nf,x);
923 1022 : if (typ(x) != t_COL) return const_col(nbrows(M), x);
924 161 : return RgM_RgC_mul(M, x);
925 : }
926 : /* assume x a square t_MAT, return a t_VEC of embeddings of its columns */
927 : static GEN
928 14 : mattocomplex(GEN nf, GEN x)
929 : {
930 14 : long i,j, l = lg(x);
931 14 : GEN v = cgetg(l, t_VEC);
932 98 : for (j=1; j<l; j++)
933 : {
934 84 : GEN c = gel(x,j), b = cgetg(l, t_MAT);
935 714 : for (i=1; i<l; i++) gel(b,i) = nftocomplex(nf, gel(c,i));
936 84 : b = shallowtrans(b); settyp(b, t_COL);
937 84 : gel(v,j) = b;
938 : }
939 14 : return v;
940 : }
941 :
942 : static GEN
943 14 : nf_all_roots(GEN nf, GEN x, long prec)
944 : {
945 14 : long i, j, l = lg(x), ru = lg(nf_get_roots(nf));
946 14 : GEN y = cgetg(l, t_POL), v, z;
947 :
948 14 : x = RgX_to_nfX(nf, x);
949 14 : y[1] = x[1];
950 112 : for (i=2; i<l; i++) gel(y,i) = nftocomplex(nf, gel(x,i));
951 14 : i = gprecision(y); if (i && i <= 3) return NULL;
952 :
953 14 : v = cgetg(ru, t_VEC);
954 14 : z = cgetg(l, t_POL); z[1] = x[1];
955 42 : for (i=1; i<ru; i++)
956 : {
957 224 : for (j = 2; j < l; j++) gel(z,j) = gmael(y,j,i);
958 28 : gel(v,i) = cleanroots(z, prec);
959 : }
960 14 : return v;
961 : }
962 :
963 : static GEN
964 357 : rnfscal(GEN m, GEN x, GEN y)
965 : {
966 357 : long i, l = lg(m);
967 357 : GEN z = cgetg(l, t_COL);
968 1071 : for (i = 1; i < l; i++)
969 714 : gel(z,i) = gmul(conj_i(shallowtrans(gel(x,i))), gmul(gel(m,i), gel(y,i)));
970 357 : return z;
971 : }
972 :
973 : /* x ideal in HNF */
974 : static GEN
975 364 : findmin(GEN nf, GEN x, GEN muf)
976 : {
977 364 : pari_sp av = avma;
978 : long e;
979 364 : GEN cx, y, m, M = nf_get_M(nf);
980 :
981 364 : x = Q_primitive_part(x, &cx);
982 364 : if (gequal1(gcoeff(x,1,1))) y = M;
983 : else
984 : {
985 210 : GEN G = nf_get_G(nf);
986 210 : m = lllfp(RgM_mul(G,x), 0.75, 0);
987 210 : if (typ(m) != t_MAT)
988 : {
989 0 : x = ZM_lll(x, 0.75, LLL_INPLACE);
990 0 : m = lllfp(RgM_mul(G,x), 0.75, 0);
991 0 : if (typ(m) != t_MAT) pari_err_PREC("rnflllgram");
992 : }
993 210 : x = ZM_mul(x, m);
994 210 : y = RgM_mul(M, x);
995 : }
996 364 : m = RgM_solve_realimag(y, muf);
997 364 : if (!m) return NULL; /* precision problem */
998 364 : if (cx) m = RgC_Rg_div(m, cx);
999 364 : m = grndtoi(m, &e);
1000 364 : if (e >= 0) return NULL; /* precision problem */
1001 364 : m = ZM_ZC_mul(x, m);
1002 364 : if (cx) m = ZC_Q_mul(m, cx);
1003 364 : return gerepileupto(av, m);
1004 : }
1005 :
1006 : static int
1007 364 : RED(long k, long l, GEN U, GEN mu, GEN MC, GEN nf, GEN I, GEN *Ik_inv)
1008 : {
1009 : GEN x, xc, ideal;
1010 : long i;
1011 :
1012 364 : if (!*Ik_inv) *Ik_inv = idealinv(nf, gel(I,k));
1013 364 : ideal = idealmul(nf,gel(I,l), *Ik_inv);
1014 364 : x = findmin(nf, ideal, gcoeff(mu,k,l));
1015 364 : if (!x) return 0;
1016 364 : if (gequal0(x)) return 1;
1017 :
1018 294 : xc = nftocomplex(nf,x);
1019 294 : gel(MC,k) = gsub(gel(MC,k), vecmul(xc,gel(MC,l)));
1020 294 : gel(U,k) = gsub(gel(U,k), gmul(coltoalg(nf,x), gel(U,l)));
1021 294 : gcoeff(mu,k,l) = gsub(gcoeff(mu,k,l), xc);
1022 1029 : for (i=1; i<l; i++)
1023 735 : gcoeff(mu,k,i) = gsub(gcoeff(mu,k,i), vecmul(xc,gcoeff(mu,l,i)));
1024 294 : return 1;
1025 : }
1026 :
1027 : static int
1028 84 : check_0(GEN B)
1029 : {
1030 84 : long i, l = lg(B);
1031 252 : for (i = 1; i < l; i++)
1032 168 : if (gsigne(gel(B,i)) <= 0) return 1;
1033 84 : return 0;
1034 : }
1035 :
1036 : static int
1037 98 : do_SWAP(GEN I, GEN MC, GEN MCS, GEN h, GEN mu, GEN B, long kmax, long k,
1038 : const long alpha, long r1)
1039 : {
1040 : GEN p1, p2, muf, mufc, Bf, temp;
1041 : long i, j;
1042 :
1043 98 : p1 = nftau(r1, gadd(gel(B,k),
1044 98 : gmul(gnorml2(gcoeff(mu,k,k-1)), gel(B,k-1))));
1045 98 : p2 = nftau(r1, gel(B,k-1));
1046 98 : if (gcmp(gmulsg(alpha,p1), gmulsg(alpha-1,p2)) > 0) return 0;
1047 :
1048 14 : swap(gel(MC,k-1),gel(MC,k));
1049 14 : swap(gel(h,k-1), gel(h,k));
1050 14 : swap(gel(I,k-1), gel(I,k));
1051 91 : for (j=1; j<=k-2; j++) swap(gcoeff(mu,k-1,j),gcoeff(mu,k,j));
1052 14 : muf = gcoeff(mu,k,k-1);
1053 14 : mufc = conj_i(muf);
1054 14 : Bf = gadd(gel(B,k), vecmul(real_i(vecmul(muf,mufc)), gel(B,k-1)));
1055 14 : if (check_0(Bf)) return 1; /* precision problem */
1056 :
1057 14 : p1 = vecdiv(gel(B,k-1),Bf);
1058 14 : gcoeff(mu,k,k-1) = vecmul(mufc,p1);
1059 14 : temp = gel(MCS,k-1);
1060 14 : gel(MCS,k-1) = gadd(gel(MCS,k), vecmul(muf,gel(MCS,k-1)));
1061 14 : gel(MCS,k) = gsub(vecmul(vecdiv(gel(B,k),Bf), temp),
1062 14 : vecmul(gcoeff(mu,k,k-1), gel(MCS,k)));
1063 14 : gel(B,k) = vecmul(gel(B,k),p1);
1064 14 : gel(B,k-1) = Bf;
1065 14 : for (i=k+1; i<=kmax; i++)
1066 : {
1067 0 : temp = gcoeff(mu,i,k);
1068 0 : gcoeff(mu,i,k) = gsub(gcoeff(mu,i,k-1), vecmul(muf, gcoeff(mu,i,k)));
1069 0 : gcoeff(mu,i,k-1) = gadd(temp, vecmul(gcoeff(mu,k,k-1),gcoeff(mu,i,k)));
1070 : }
1071 14 : return 1;
1072 : }
1073 :
1074 : static GEN
1075 14 : rel_T2(GEN nf, GEN pol, long lx, long prec)
1076 : {
1077 : long ru, i, j, k, l;
1078 : GEN T2, s, unro, roorder, powreorder;
1079 :
1080 14 : roorder = nf_all_roots(nf, pol, prec);
1081 14 : if (!roorder) return NULL;
1082 14 : ru = lg(roorder);
1083 98 : unro = cgetg(lx,t_COL); for (i=1; i<lx; i++) gel(unro,i) = gen_1;
1084 14 : powreorder = cgetg(lx,t_MAT); gel(powreorder,1) = unro;
1085 14 : T2 = cgetg(ru, t_VEC);
1086 42 : for (i = 1; i < ru; i++)
1087 : {
1088 28 : GEN ro = gel(roorder,i);
1089 28 : GEN m = initmat(lx);
1090 168 : for (k=2; k<lx; k++)
1091 : {
1092 140 : GEN c = cgetg(lx, t_COL); gel(powreorder,k) = c;
1093 1232 : for (j=1; j < lx; j++)
1094 1092 : gel(c,j) = gmul(gel(ro,j), gmael(powreorder,k-1,j));
1095 : }
1096 196 : for (l = 1; l < lx; l++)
1097 882 : for (k = 1; k <= l; k++)
1098 : {
1099 714 : s = gen_0;
1100 6636 : for (j = 1; j < lx; j++)
1101 5922 : s = gadd(s, gmul(conj_i(gmael(powreorder,k,j)),
1102 5922 : gmael(powreorder,l,j)));
1103 714 : if (l == k)
1104 168 : gcoeff(m, l, l) = real_i(s);
1105 : else
1106 : {
1107 546 : gcoeff(m, k, l) = s;
1108 546 : gcoeff(m, l, k) = conj_i(s);
1109 : }
1110 : }
1111 28 : gel(T2,i) = m;
1112 : }
1113 14 : return T2;
1114 : }
1115 :
1116 : /* given a base field nf (e.g main variable y), a polynomial pol with
1117 : * coefficients in nf (e.g main variable x), and an order as output
1118 : * by rnfpseudobasis, outputs a reduced order. */
1119 : GEN
1120 14 : rnflllgram(GEN nf, GEN pol, GEN order,long prec)
1121 : {
1122 14 : pari_sp av = avma;
1123 14 : long j, k, l, kmax, r1, lx, count = 0;
1124 : GEN M, I, h, H, mth, MC, MPOL, MCS, B, mu;
1125 14 : const long alpha = 10, MAX_COUNT = 4;
1126 :
1127 14 : nf = checknf(nf); r1 = nf_get_r1(nf);
1128 14 : check_ZKmodule(order, "rnflllgram");
1129 14 : M = gel(order,1);
1130 14 : I = gel(order,2); lx = lg(I);
1131 14 : if (lx < 3) return gcopy(order);
1132 14 : if (lx-1 != degpol(pol)) pari_err_DIM("rnflllgram");
1133 14 : I = leafcopy(I);
1134 14 : H = NULL;
1135 14 : MPOL = matbasistoalg(nf, M);
1136 14 : MCS = matid(lx-1); /* dummy for gerepile */
1137 14 : PRECNF:
1138 14 : if (count == MAX_COUNT)
1139 : {
1140 0 : prec = precdbl(prec); count = 0;
1141 0 : if (DEBUGLEVEL) pari_warn(warnprec,"rnflllgram",prec);
1142 0 : nf = nfnewprec_shallow(nf,prec);
1143 : }
1144 14 : mth = rel_T2(nf, pol, lx, prec);
1145 14 : if (!mth) { count = MAX_COUNT; goto PRECNF; }
1146 14 : h = NULL;
1147 14 : PRECPB:
1148 14 : if (h)
1149 : { /* precision problem, recompute. If no progress, increase nf precision */
1150 0 : if (++count == MAX_COUNT || RgM_isidentity(h)) {count = MAX_COUNT; goto PRECNF;}
1151 0 : H = H? gmul(H, h): h;
1152 0 : MPOL = gmul(MPOL, h);
1153 : }
1154 14 : h = matid(lx-1);
1155 14 : MC = mattocomplex(nf, MPOL);
1156 14 : mu = cgetg(lx,t_MAT);
1157 14 : B = cgetg(lx,t_COL);
1158 98 : for (j=1; j<lx; j++)
1159 : {
1160 84 : gel(mu,j) = zerocol(lx - 1);
1161 84 : gel(B,j) = gen_0;
1162 : }
1163 14 : if (DEBUGLEVEL) err_printf("k = ");
1164 14 : gel(B,1) = real_i(rnfscal(mth,gel(MC,1),gel(MC,1)));
1165 14 : gel(MCS,1) = gel(MC,1);
1166 14 : kmax = 1; k = 2;
1167 : do
1168 : {
1169 98 : GEN Ik_inv = NULL;
1170 98 : if (DEBUGLEVEL) err_printf("%ld ",k);
1171 98 : if (k > kmax)
1172 : { /* Incremental Gram-Schmidt */
1173 70 : kmax = k; gel(MCS,k) = gel(MC,k);
1174 343 : for (j=1; j<k; j++)
1175 : {
1176 273 : gcoeff(mu,k,j) = vecdiv(rnfscal(mth,gel(MCS,j),gel(MC,k)),
1177 273 : gel(B,j));
1178 273 : gel(MCS,k) = gsub(gel(MCS,k), vecmul(gcoeff(mu,k,j),gel(MCS,j)));
1179 : }
1180 70 : gel(B,k) = real_i(rnfscal(mth,gel(MCS,k),gel(MCS,k)));
1181 70 : if (check_0(gel(B,k))) goto PRECPB;
1182 : }
1183 98 : if (!RED(k, k-1, h, mu, MC, nf, I, &Ik_inv)) goto PRECPB;
1184 98 : if (do_SWAP(I,MC,MCS,h,mu,B,kmax,k,alpha, r1))
1185 : {
1186 14 : if (!B[k]) goto PRECPB;
1187 14 : if (k > 2) k--;
1188 : }
1189 : else
1190 : {
1191 350 : for (l=k-2; l; l--)
1192 266 : if (!RED(k, l, h, mu, MC, nf, I, &Ik_inv)) goto PRECPB;
1193 84 : k++;
1194 : }
1195 98 : if (gc_needed(av,2))
1196 : {
1197 0 : if(DEBUGMEM>1) pari_warn(warnmem,"rnflllgram");
1198 0 : gerepileall(av, H?10:9, &nf,&mth,&h,&MPOL,&B,&MC,&MCS,&mu,&I,&H);
1199 : }
1200 : }
1201 98 : while (k < lx);
1202 14 : MPOL = gmul(MPOL,h);
1203 14 : if (H) h = gmul(H, h);
1204 14 : if (DEBUGLEVEL) err_printf("\n");
1205 14 : MPOL = RgM_to_nfM(nf,MPOL);
1206 14 : h = RgM_to_nfM(nf,h);
1207 14 : return gerepilecopy(av, mkvec2(mkvec2(MPOL,I), h));
1208 : }
1209 :
1210 : GEN
1211 7 : rnfpolred(GEN nf, GEN pol, long prec)
1212 : {
1213 7 : pari_sp av = avma;
1214 7 : long i, j, n, v = varn(pol);
1215 7 : GEN id, w, I, O, nfpol, bnf = checkbnf_i(nf);
1216 :
1217 7 : if (typ(pol)!=t_POL) pari_err_TYPE("rnfpolred",pol);
1218 7 : nf = bnf? bnf_get_nf(bnf): checknf(nf);
1219 7 : if (degpol(pol) <= 1) { w = cgetg(2, t_VEC); gel(w,1) = pol_x(v); return w; }
1220 7 : nfpol = nf_get_pol(nf);
1221 :
1222 7 : id = rnfpseudobasis(nf,pol);
1223 7 : if (bnf && is_pm1( bnf_get_no(bnf) )) /* if bnf is principal */
1224 : {
1225 : GEN newI, newO;
1226 0 : O = gel(id,1);
1227 0 : I = gel(id,2); n = lg(I)-1;
1228 0 : newI = cgetg(n+1,t_VEC);
1229 0 : newO = cgetg(n+1,t_MAT);
1230 0 : for (j=1; j<=n; j++)
1231 : {
1232 0 : GEN al = gen_if_principal(bnf,gel(I,j));
1233 0 : gel(newI,j) = gen_1;
1234 0 : gel(newO,j) = nfC_nf_mul(nf, gel(O,j), al);
1235 : }
1236 0 : id = mkvec2(newO, newI);
1237 : }
1238 :
1239 7 : id = gel(rnflllgram(nf,pol,id,prec),1);
1240 7 : O = gel(id,1);
1241 7 : I = gel(id,2); n = lg(I)-1;
1242 7 : w = cgetg(n+1,t_VEC);
1243 7 : pol = lift_shallow(pol);
1244 70 : for (j=1; j<=n; j++)
1245 : {
1246 63 : GEN newpol, L, a, Ij = gel(I,j);
1247 63 : a = RgC_Rg_mul(gel(O,j), (typ(Ij) == t_MAT)? gcoeff(Ij,1,1): Ij);
1248 630 : for (i=n; i; i--) gel(a,i) = nf_to_scalar_or_alg(nf, gel(a,i));
1249 63 : a = RgV_to_RgX(a, v);
1250 63 : newpol = RgXQX_red(RgXQ_charpoly(a, pol, v), nfpol);
1251 63 : newpol = Q_primpart(newpol);
1252 :
1253 63 : (void)nfgcd_all(newpol, RgX_deriv(newpol), nfpol, nf_get_index(nf), &newpol);
1254 63 : L = leading_coeff(newpol);
1255 63 : gel(w,j) = (typ(L) == t_POL)? RgXQX_div(newpol, L, nfpol)
1256 63 : : RgX_Rg_div(newpol, L);
1257 : }
1258 7 : return gerepilecopy(av,w);
1259 : }
1260 :
1261 : /*******************************************************************/
1262 : /* */
1263 : /* LINEAR ALGEBRA OVER Z_K (HNF,SNF) */
1264 : /* */
1265 : /*******************************************************************/
1266 : /* A torsion-free module M over Z_K is given by [A,I].
1267 : * I=[a_1,...,a_k] is a row vector of k fractional ideals given in HNF.
1268 : * A is an n x k matrix (same k) such that if A_j is the j-th column of A then
1269 : * M=a_1 A_1+...+a_k A_k. We say that [A,I] is a pseudo-basis if k=n */
1270 :
1271 : /* Given an element x and an ideal I in HNF, gives an r such that x-r is in H
1272 : * and r is small */
1273 : GEN
1274 7 : nfreduce(GEN nf, GEN x, GEN I)
1275 : {
1276 7 : pari_sp av = avma;
1277 7 : x = nf_to_scalar_or_basis(checknf(nf), x);
1278 7 : if (idealtyp(&I, NULL) != id_MAT || lg(I)==1) pari_err_TYPE("nfreduce",I);
1279 7 : if (typ(x) != t_COL) x = scalarcol( gmod(x, gcoeff(I,1,1)), lg(I)-1 );
1280 7 : else x = reducemodinvertible(x, I);
1281 7 : return gerepileupto(av, x);
1282 : }
1283 : /* Given an element x and an ideal in HNF, gives an a in ideal such that
1284 : * x-a is small. No checks */
1285 : static GEN
1286 38359 : element_close(GEN nf, GEN x, GEN ideal)
1287 : {
1288 38359 : pari_sp av = avma;
1289 38359 : GEN y = gcoeff(ideal,1,1);
1290 38359 : x = nf_to_scalar_or_basis(nf, x);
1291 38359 : if (typ(y) == t_INT && is_pm1(y)) return ground(x);
1292 35776 : if (typ(x) == t_COL)
1293 14797 : x = closemodinvertible(x, ideal);
1294 : else
1295 20979 : x = gmul(y, gdivround(x,y));
1296 35776 : return gerepileupto(av, x);
1297 : }
1298 :
1299 : /* A + v B */
1300 : static GEN
1301 132855 : colcomb1(GEN nf, GEN v, GEN A, GEN B)
1302 : {
1303 132855 : if (isintzero(v)) return A;
1304 88694 : return RgC_to_nfC(nf, RgC_add(A, nfC_nf_mul(nf,B,v)));
1305 : }
1306 : /* u A + v B */
1307 : static GEN
1308 110694 : colcomb(GEN nf, GEN u, GEN v, GEN A, GEN B)
1309 : {
1310 110694 : if (isintzero(u)) return nfC_nf_mul(nf,B,v);
1311 94475 : if (u != gen_1) A = nfC_nf_mul(nf,A,u);
1312 94475 : return colcomb1(nf, v, A, B);
1313 : }
1314 :
1315 : /* return m[i,1..lim] * x */
1316 : static GEN
1317 315 : element_mulvecrow(GEN nf, GEN x, GEN m, long i, long lim)
1318 : {
1319 315 : long j, l = minss(lg(m), lim+1);
1320 315 : GEN dx, y = cgetg(l, t_VEC);
1321 315 : x = nf_to_scalar_or_basis(nf, x);
1322 315 : if (typ(x) == t_COL)
1323 : {
1324 91 : x = zk_multable(nf, Q_remove_denom(x, &dx));
1325 350 : for (j=1; j<l; j++)
1326 : {
1327 259 : GEN t = gcoeff(m,i,j);
1328 259 : if (!isintzero(t))
1329 : {
1330 112 : if (typ(t) == t_COL)
1331 28 : t = RgM_RgC_mul(x, t);
1332 : else
1333 84 : t = ZC_Q_mul(gel(x,1), t);
1334 112 : if (dx) t = gdiv(t, dx);
1335 112 : t = nf_to_scalar_or_basis(nf,t);
1336 : }
1337 259 : gel(y,j) = t;
1338 : }
1339 : }
1340 : else
1341 : {
1342 784 : for (j=1; j<l; j++) gel(y,j) = gmul(x, gcoeff(m,i,j));
1343 : }
1344 315 : return y;
1345 : }
1346 :
1347 : /* u Z[s,] + v Z[t,], limitied to the first lim entries */
1348 : static GEN
1349 196 : rowcomb(GEN nf, GEN u, GEN v, long s, long t, GEN Z, long lim)
1350 : {
1351 : GEN z;
1352 196 : if (gequal0(u))
1353 7 : z = element_mulvecrow(nf,v,Z,t, lim);
1354 : else
1355 : {
1356 189 : z = element_mulvecrow(nf,u,Z,s, lim);
1357 189 : if (!gequal0(v)) z = gadd(z, element_mulvecrow(nf,v,Z,t, lim));
1358 : }
1359 196 : return z;
1360 : }
1361 :
1362 : /* nfbezout(0,b,A,B). Either bB = NULL or b*B */
1363 : static GEN
1364 64438 : zero_nfbezout(GEN nf,GEN bB, GEN b, GEN A,GEN B,GEN *u,GEN *v,GEN *w,GEN *di)
1365 : {
1366 : GEN d;
1367 64438 : if (isint1(b))
1368 : {
1369 62920 : *v = gen_1;
1370 62920 : *w = A;
1371 62920 : d = B;
1372 62920 : *di = idealinv(nf,d);
1373 : }
1374 : else
1375 : {
1376 1518 : *v = nfinv(nf,b);
1377 1518 : *w = idealmul(nf,A,*v);
1378 1518 : d = bB? bB: idealmul(nf,b,B);
1379 1518 : *di = idealHNF_inv(nf,d);
1380 : }
1381 64438 : *u = gen_0; return d;
1382 : }
1383 :
1384 : /* Given elements a,b and ideals A, B, outputs d = a.A+b.B and gives
1385 : * di=d^-1, w=A.B.di, u, v such that au+bv=1 and u in A.di, v in B.di.
1386 : * Assume A, B nonzero, but a or b can be zero (not both) */
1387 : static GEN
1388 69695 : nfbezout(GEN nf,GEN a,GEN b, GEN A,GEN B, GEN *pu,GEN *pv,GEN *pw,GEN *pdi,
1389 : int red)
1390 : {
1391 : GEN w, u, v, d, di, aA, bB;
1392 :
1393 69695 : if (isintzero(a)) return zero_nfbezout(nf,NULL,b,A,B,pu,pv,pw,pdi);
1394 69695 : if (isintzero(b)) return zero_nfbezout(nf,NULL,a,B,A,pv,pu,pw,pdi);
1395 :
1396 69695 : if (a != gen_1) /* frequently called with a = gen_1 */
1397 : {
1398 39247 : a = nf_to_scalar_or_basis(nf,a);
1399 39247 : if (isint1(a)) a = gen_1;
1400 : }
1401 69695 : aA = (a == gen_1)? idealhnf_shallow(nf,A): idealmul(nf,a,A);
1402 69695 : bB = idealmul(nf,b,B);
1403 69695 : d = idealadd(nf,aA,bB);
1404 69695 : if (gequal(aA, d)) return zero_nfbezout(nf,d, a,B,A,pv,pu,pw,pdi);
1405 32738 : if (gequal(bB, d)) return zero_nfbezout(nf,d, b,A,B,pu,pv,pw,pdi);
1406 : /* general case is slow */
1407 5257 : di = idealHNF_inv(nf,d);
1408 5257 : aA = idealmul(nf,aA,di); /* integral */
1409 5257 : bB = idealmul(nf,bB,di); /* integral */
1410 :
1411 5257 : u = red? idealaddtoone_i(nf, aA, bB): idealaddtoone_raw(nf, aA, bB);
1412 5257 : w = idealmul(nf,aA,B);
1413 5257 : v = nfdiv(nf, nfsub(nf, gen_1, u), b);
1414 5257 : if (a != gen_1)
1415 : {
1416 1344 : GEN inva = nfinv(nf, a);
1417 1344 : u = nfmul(nf,u,inva);
1418 1344 : w = idealmul(nf, inva, w); /* AB/d */
1419 : }
1420 5257 : *pu = u; *pv = v; *pw = w; *pdi = di; return d;
1421 : }
1422 : /* v a vector of ideals, simplify in place the ones generated by elts of Q */
1423 : static void
1424 10197 : idV_simplify(GEN v)
1425 : {
1426 10197 : long i, l = lg(v);
1427 58220 : for (i = 1; i < l; i++)
1428 : {
1429 48023 : GEN M = gel(v,i);
1430 48023 : if (typ(M)==t_MAT && RgM_isscalar(M,NULL))
1431 15034 : gel(v,i) = Q_abs_shallow(gcoeff(M,1,1));
1432 : }
1433 10197 : }
1434 : /* Given a torsion-free module x outputs a pseudo-basis for x in HNF */
1435 : GEN
1436 7125 : nfhnf0(GEN nf, GEN x, long flag)
1437 : {
1438 : long i, j, def, idef, m, n;
1439 7125 : pari_sp av0 = avma, av;
1440 : GEN y, A, I, J, U;
1441 :
1442 7125 : nf = checknf(nf);
1443 7125 : check_ZKmodule(x, "nfhnf");
1444 7125 : A = gel(x,1); RgM_dimensions(A, &m, &n);
1445 7125 : I = gel(x,2);
1446 7125 : if (!n) {
1447 49 : if (!flag) return gcopy(x);
1448 0 : retmkvec2(gcopy(x), cgetg(1,t_MAT));
1449 : }
1450 7076 : U = flag? matid(n): NULL;
1451 7076 : idef = (n < m)? m-n : 0;
1452 7076 : av = avma;
1453 7076 : A = RgM_to_nfM(nf,A);
1454 7076 : I = leafcopy(I);
1455 7076 : J = zerovec(n); def = n;
1456 40758 : for (i=m; i>idef; i--)
1457 : {
1458 33682 : GEN d, di = NULL;
1459 :
1460 33927 : j=def; while (j>=1 && isintzero(gcoeff(A,i,j))) j--;
1461 33682 : if (!j)
1462 : { /* no pivot on line i */
1463 7 : if (idef) idef--;
1464 7 : continue;
1465 : }
1466 33675 : if (j==def) j--;
1467 : else {
1468 63 : swap(gel(A,j), gel(A,def));
1469 63 : swap(gel(I,j), gel(I,def));
1470 63 : if (U) swap(gel(U,j), gel(U,def));
1471 : }
1472 198356 : for ( ; j; j--)
1473 : {
1474 164681 : GEN a,b, u,v,w, S, T, S0, T0 = gel(A,j);
1475 164681 : b = gel(T0,i); if (isintzero(b)) continue;
1476 :
1477 40726 : S0 = gel(A,def); a = gel(S0,i);
1478 40726 : d = nfbezout(nf, a,b, gel(I,def),gel(I,j), &u,&v,&w,&di,1);
1479 40726 : S = colcomb(nf, u,v, S0,T0);
1480 40726 : T = colcomb(nf, a,gneg(b), T0,S0);
1481 40726 : gel(A,def) = S; gel(A,j) = T;
1482 40726 : gel(I,def) = d; gel(I,j) = w;
1483 40726 : if (U)
1484 : {
1485 42 : S0 = gel(U,def);
1486 42 : T0 = gel(U,j);
1487 42 : gel(U,def) = colcomb(nf, u,v, S0,T0);
1488 42 : gel(U,j) = colcomb(nf, a,gneg(b), T0,S0);
1489 : }
1490 : }
1491 33675 : y = gcoeff(A,i,def);
1492 33675 : if (!isint1(y))
1493 : {
1494 637 : GEN yi = nfinv(nf,y);
1495 637 : gel(A,def) = nfC_nf_mul(nf, gel(A,def), yi);
1496 637 : gel(I,def) = idealmul(nf, y, gel(I,def));
1497 637 : if (U) gel(U,def) = nfC_nf_mul(nf, gel(U,def), yi);
1498 637 : di = NULL;
1499 : }
1500 33675 : if (!di) di = idealinv(nf,gel(I,def));
1501 33675 : d = gel(I,def);
1502 33675 : gel(J,def) = di;
1503 110198 : for (j=def+1; j<=n; j++)
1504 : {
1505 76523 : GEN mc, c = gcoeff(A,i,j); if (isintzero(c)) continue;
1506 25011 : c = element_close(nf, c, idealmul(nf,d,gel(J,j)));
1507 25011 : mc = gneg(c);
1508 25011 : gel(A,j) = colcomb1(nf, mc, gel(A,j),gel(A,def));
1509 25011 : if (U) gel(U,j) = colcomb1(nf, mc, gel(U,j),gel(U,def));
1510 : }
1511 33675 : def--;
1512 33675 : if (gc_needed(av,2))
1513 : {
1514 0 : if(DEBUGMEM>1) pari_warn(warnmem,"nfhnf, i = %ld", i);
1515 0 : gerepileall(av,U?4:3, &A,&I,&J,&U);
1516 : }
1517 : }
1518 7076 : n -= def;
1519 7076 : A += def; A[0] = evaltyp(t_MAT)|evallg(n+1);
1520 7076 : I += def; I[0] = evaltyp(t_VEC)|evallg(n+1);
1521 7076 : idV_simplify(I);
1522 7076 : x = mkvec2(A,I);
1523 7076 : if (U) x = mkvec2(x,U);
1524 7076 : return gerepilecopy(av0, x);
1525 : }
1526 :
1527 : GEN
1528 7111 : nfhnf(GEN nf, GEN x) { return nfhnf0(nf, x, 0); }
1529 :
1530 : static long
1531 14 : RgV_find_denom(GEN x)
1532 : {
1533 14 : long i, l = lg(x);
1534 14 : for (i = 1; i < l; i++)
1535 14 : if (Q_denom(gel(x,i)) != gen_1) return i;
1536 0 : return 0;
1537 : }
1538 : /* A torsion module M over Z_K will be given by a row vector [A,I,J] with
1539 : * three components. I=[b_1,...,b_n] is a row vector of n fractional ideals
1540 : * given in HNF, J=[a_1,...,a_n] is a row vector of n fractional ideals in
1541 : * HNF. A is an nxn matrix (same n) such that if A_j is the j-th column of A
1542 : * and e_n is the canonical basis of K^n, then
1543 : * M=(b_1e_1+...+b_ne_n)/(a_1A_1+...a_nA_n) */
1544 :
1545 : /* x=[A,I,J] a torsion module as above. Output the
1546 : * smith normal form as K=[c_1,...,c_n] such that x = Z_K/c_1+...+Z_K/c_n */
1547 : GEN
1548 49 : nfsnf0(GEN nf, GEN x, long flag)
1549 : {
1550 : long i, j, k, l, n, m;
1551 : pari_sp av;
1552 : GEN z,u,v,w,d,dinv,A,I,J, U,V;
1553 :
1554 49 : nf = checknf(nf);
1555 49 : if (typ(x)!=t_VEC || lg(x)!=4) pari_err_TYPE("nfsnf",x);
1556 49 : A = gel(x,1);
1557 49 : I = gel(x,2);
1558 49 : J = gel(x,3);
1559 49 : if (typ(A)!=t_MAT) pari_err_TYPE("nfsnf",A);
1560 49 : n = lg(A)-1;
1561 49 : if (typ(I)!=t_VEC) pari_err_TYPE("nfsnf",I);
1562 49 : if (typ(J)!=t_VEC) pari_err_TYPE("nfsnf",J);
1563 49 : if (lg(I)!=n+1 || lg(J)!=n+1) pari_err_DIM("nfsnf");
1564 49 : RgM_dimensions(A, &m, &n);
1565 49 : if (!n || n != m) pari_err_IMPL("nfsnf for empty or non square matrices");
1566 :
1567 49 : av = avma;
1568 49 : if (!flag) U = V = NULL;
1569 : else
1570 : {
1571 21 : U = matid(m);
1572 21 : V = matid(n);
1573 : }
1574 49 : A = RgM_to_nfM(nf, A);
1575 49 : I = leafcopy(I);
1576 49 : J = leafcopy(J);
1577 168 : for (i = 1; i <= n; i++) gel(J,i) = idealinv(nf, gel(J,i));
1578 49 : z = zerovec(n);
1579 238 : for (i=n; i>=1; i--)
1580 : {
1581 : GEN Aii, a, b, db;
1582 189 : long c = 0;
1583 378 : for (j=i-1; j>=1; j--)
1584 : {
1585 189 : GEN S, T, S0, T0 = gel(A,j);
1586 189 : b = gel(T0,i); if (gequal0(b)) continue;
1587 :
1588 63 : S0 = gel(A,i); a = gel(S0,i);
1589 63 : d = nfbezout(nf, a,b, gel(J,i),gel(J,j), &u,&v,&w,&dinv,1);
1590 63 : S = colcomb(nf, u,v, S0,T0);
1591 63 : T = colcomb(nf, a,gneg(b), T0,S0);
1592 63 : gel(A,i) = S; gel(A,j) = T;
1593 63 : gel(J,i) = d; gel(J,j) = w;
1594 63 : if (V)
1595 : {
1596 28 : T0 = gel(V,j);
1597 28 : S0 = gel(V,i);
1598 28 : gel(V,i) = colcomb(nf, u,v, S0,T0);
1599 28 : gel(V,j) = colcomb(nf, a,gneg(b), T0,S0);
1600 : }
1601 : }
1602 378 : for (j=i-1; j>=1; j--)
1603 : {
1604 : GEN ri, rj;
1605 189 : b = gcoeff(A,j,i); if (gequal0(b)) continue;
1606 :
1607 70 : a = gcoeff(A,i,i);
1608 70 : d = nfbezout(nf, a,b, gel(I,i),gel(I,j), &u,&v,&w,&dinv,1);
1609 70 : ri = rowcomb(nf, u,v, i,j, A, i);
1610 70 : rj = rowcomb(nf, a,gneg(b), j,i, A, i);
1611 252 : for (k=1; k<=i; k++) {
1612 182 : gcoeff(A,j,k) = gel(rj,k);
1613 182 : gcoeff(A,i,k) = gel(ri,k);
1614 : }
1615 70 : if (U)
1616 : {
1617 28 : ri = rowcomb(nf, u,v, i,j, U, m);
1618 28 : rj = rowcomb(nf, a,gneg(b), j,i, U, m);
1619 105 : for (k=1; k<=m; k++) {
1620 77 : gcoeff(U,j,k) = gel(rj,k);
1621 77 : gcoeff(U,i,k) = gel(ri,k);
1622 : }
1623 : }
1624 70 : gel(I,i) = d; gel(I,j) = w; c = 1;
1625 : }
1626 189 : if (c) { i++; continue; }
1627 :
1628 133 : Aii = gcoeff(A,i,i); if (gequal0(Aii)) continue;
1629 133 : gel(J,i) = idealmul(nf, gel(J,i), Aii);
1630 133 : gcoeff(A,i,i) = gen_1;
1631 133 : if (V) gel(V,i) = nfC_nf_mul(nf, gel(V,i), nfinv(nf,Aii));
1632 133 : gel(z,i) = idealmul(nf,gel(J,i),gel(I,i));
1633 133 : b = Q_remove_denom(gel(z,i), &db);
1634 238 : for (k=1; k<i; k++)
1635 238 : for (l=1; l<i; l++)
1636 : {
1637 147 : GEN d, D, p1, p2, p3, Akl = gcoeff(A,k,l);
1638 : long t;
1639 147 : if (gequal0(Akl)) continue;
1640 :
1641 133 : p1 = idealmul(nf,Akl,gel(J,l));
1642 133 : p3 = idealmul(nf, p1, gel(I,k));
1643 133 : if (db) p3 = RgM_Rg_mul(p3, db);
1644 133 : if (RgM_is_ZM(p3) && hnfdivide(b, p3)) continue;
1645 :
1646 : /* find d in D = I[k]/I[i] not in J[i]/(A[k,l] J[l]) */
1647 14 : D = idealdiv(nf,gel(I,k),gel(I,i));
1648 14 : p2 = idealdiv(nf,gel(J,i), p1);
1649 14 : t = RgV_find_denom(QM_gauss(p2, D));
1650 14 : if (!t) pari_err_BUG("nfsnf");
1651 14 : d = gel(D,t);
1652 14 : p1 = element_mulvecrow(nf,d,A,k,i);
1653 42 : for (t=1; t<=i; t++) gcoeff(A,i,t) = gadd(gcoeff(A,i,t),gel(p1,t));
1654 14 : if (U)
1655 : {
1656 7 : p1 = element_mulvecrow(nf,d,U,k,i);
1657 21 : for (t=1; t<=i; t++) gcoeff(U,i,t) = gadd(gcoeff(U,i,t),gel(p1,t));
1658 : }
1659 :
1660 14 : k = i; c = 1; break;
1661 : }
1662 133 : if (gc_needed(av,1))
1663 : {
1664 0 : if(DEBUGMEM>1) pari_warn(warnmem,"nfsnf");
1665 0 : gerepileall(av,U?6:4, &A,&I,&J,&z,&U,&V);
1666 : }
1667 133 : if (c) i++; /* iterate on row/column i */
1668 : }
1669 49 : if (U) z = mkvec3(z,U,V);
1670 49 : return gerepilecopy(av, z);
1671 : }
1672 : GEN
1673 0 : nfsnf(GEN nf, GEN x) { return nfsnf0(nf,x,0); }
1674 :
1675 : /* Given a pseudo-basis x, outputs a multiple of its ideal determinant */
1676 : GEN
1677 329 : nfdetint(GEN nf, GEN x)
1678 : {
1679 : GEN pass,c,v,det1,piv,pivprec,vi,p1,A,I,id,idprod;
1680 329 : long i, j, k, rg, n, m, m1, cm=0, N;
1681 329 : pari_sp av = avma, av1;
1682 :
1683 329 : nf = checknf(nf); N = nf_get_degree(nf);
1684 329 : check_ZKmodule(x, "nfdetint");
1685 329 : A = gel(x,1);
1686 329 : I = gel(x,2);
1687 329 : n = lg(A)-1; if (!n) return gen_1;
1688 :
1689 329 : m1 = lgcols(A); m = m1-1;
1690 329 : id = matid(N);
1691 1463 : c = new_chunk(m1); for (k=1; k<=m; k++) c[k] = 0;
1692 329 : piv = pivprec = gen_1;
1693 :
1694 329 : av1 = avma;
1695 329 : det1 = idprod = gen_0; /* dummy for gerepileall */
1696 329 : pass = cgetg(m1,t_MAT);
1697 329 : v = cgetg(m1,t_COL);
1698 1463 : for (j=1; j<=m; j++)
1699 : {
1700 1134 : gel(pass,j) = zerocol(m);
1701 1134 : gel(v,j) = gen_0; /* dummy */
1702 : }
1703 1848 : for (rg=0,k=1; k<=n; k++)
1704 : {
1705 1519 : long t = 0;
1706 9555 : for (i=1; i<=m; i++)
1707 8036 : if (!c[i])
1708 : {
1709 4186 : vi=nfmul(nf,piv,gcoeff(A,i,k));
1710 39312 : for (j=1; j<=m; j++)
1711 35126 : if (c[j]) vi=gadd(vi,nfmul(nf,gcoeff(pass,i,j),gcoeff(A,j,k)));
1712 4186 : gel(v,i) = vi; if (!t && !gequal0(vi)) t=i;
1713 : }
1714 1519 : if (t)
1715 : {
1716 1519 : pivprec = piv;
1717 1519 : if (rg == m-1)
1718 : {
1719 714 : if (!cm)
1720 : {
1721 329 : cm=1; idprod = id;
1722 1463 : for (i=1; i<=m; i++)
1723 1134 : if (i!=t)
1724 805 : idprod = (idprod==id)? gel(I,c[i])
1725 805 : : idealmul(nf,idprod,gel(I,c[i]));
1726 : }
1727 714 : p1 = idealmul(nf,gel(v,t),gel(I,k)); c[t]=0;
1728 714 : det1 = (typ(det1)==t_INT)? p1: idealadd(nf,p1,det1);
1729 : }
1730 : else
1731 : {
1732 805 : rg++; piv=gel(v,t); c[t]=k;
1733 6139 : for (i=1; i<=m; i++)
1734 5334 : if (!c[i])
1735 : {
1736 29757 : for (j=1; j<=m; j++)
1737 27090 : if (c[j] && j!=t)
1738 : {
1739 7252 : p1 = gsub(nfmul(nf,piv,gcoeff(pass,i,j)),
1740 7252 : nfmul(nf,gel(v,i),gcoeff(pass,t,j)));
1741 7252 : gcoeff(pass,i,j) = rg>1? nfdiv(nf,p1,pivprec)
1742 7252 : : p1;
1743 : }
1744 2667 : gcoeff(pass,i,t) = gneg(gel(v,i));
1745 : }
1746 : }
1747 : }
1748 1519 : if (gc_needed(av1,1))
1749 : {
1750 0 : if(DEBUGMEM>1) pari_warn(warnmem,"nfdetint");
1751 0 : gerepileall(av1,6, &det1,&piv,&pivprec,&pass,&v,&idprod);
1752 : }
1753 : }
1754 329 : if (!cm) { set_avma(av); return cgetg(1,t_MAT); }
1755 329 : return gerepileupto(av, idealmul(nf,idprod,det1));
1756 : }
1757 :
1758 : /* reduce in place components of x[1..lim] mod D (destroy x). D in HNF */
1759 : static void
1760 26737 : nfcleanmod(GEN nf, GEN x, long lim, GEN D)
1761 : {
1762 : GEN DZ, DZ2, dD;
1763 : long i;
1764 26737 : D = Q_remove_denom(D, &dD);
1765 26737 : DZ = gcoeff(D,1,1); DZ2 = shifti(DZ, -1);
1766 130157 : for (i = 1; i <= lim; i++)
1767 : {
1768 103420 : GEN c = nf_to_scalar_or_basis(nf, gel(x,i));
1769 103420 : switch(typ(c)) /* c = centermod(c, D) */
1770 : {
1771 92551 : case t_INT:
1772 92551 : if (!signe(c)) break;
1773 50672 : if (dD) c = mulii(c, dD);
1774 50672 : c = centermodii(c, DZ, DZ2);
1775 50672 : if (dD) c = Qdivii(c,dD);
1776 50672 : break;
1777 217 : case t_FRAC: {
1778 217 : GEN dc = gel(c,2), nc = gel(c,1), N = mulii(DZ, dc);
1779 217 : if (dD) nc = mulii(nc, dD);
1780 217 : c = centermodii(nc, N, shifti(N,-1));
1781 217 : c = Qdivii(c, dD ? mulii(dc,dD): dc);
1782 217 : break;
1783 : }
1784 10652 : case t_COL: {
1785 : GEN dc;
1786 10652 : c = Q_remove_denom(c, &dc);
1787 10652 : if (dD) c = ZC_Z_mul(c, dD);
1788 10652 : c = ZC_hnfrem(c, dc? ZM_Z_mul(D,dc): D);
1789 10652 : dc = mul_content(dc, dD);
1790 10652 : if (ZV_isscalar(c))
1791 : {
1792 175 : c = gel(c,1);
1793 175 : if (dc) c = Qdivii(c,dc);
1794 : }
1795 : else
1796 10477 : if (dc) c = RgC_Rg_div(c, dc);
1797 10652 : break;
1798 : }
1799 : }
1800 103420 : gel(x,i) = c;
1801 : }
1802 26737 : }
1803 :
1804 : GEN
1805 3121 : nfhnfmod(GEN nf, GEN x, GEN D)
1806 : {
1807 : long li, co, i, j, def, ldef;
1808 3121 : pari_sp av0=avma, av;
1809 : GEN dA, dI, d0, w, p1, d, u, v, A, I, J, di;
1810 :
1811 3121 : nf = checknf(nf);
1812 3121 : check_ZKmodule(x, "nfhnfmod");
1813 3121 : A = gel(x,1);
1814 3121 : I = gel(x,2);
1815 3121 : co = lg(A); if (co==1) return cgetg(1,t_MAT);
1816 :
1817 3121 : li = lgcols(A);
1818 3121 : if (typ(D)!=t_MAT) D = idealhnf_shallow(nf, D);
1819 3121 : D = Q_remove_denom(D, NULL);
1820 3121 : RgM_check_ZM(D, "nfhnfmod");
1821 :
1822 3121 : av = avma;
1823 3121 : A = RgM_to_nfM(nf, A);
1824 3121 : A = Q_remove_denom(A, &dA);
1825 3121 : I = Q_remove_denom(leafcopy(I), &dI);
1826 3121 : dA = mul_denom(dA,dI);
1827 3121 : if (dA) D = ZM_Z_mul(D, powiu(dA, minss(li,co)));
1828 :
1829 3121 : def = co; ldef = (li>co)? li-co+1: 1;
1830 17469 : for (i=li-1; i>=ldef; i--)
1831 : {
1832 16679 : def--; j=def; while (j>=1 && isintzero(gcoeff(A,i,j))) j--;
1833 14348 : if (!j) continue;
1834 14348 : if (j==def) j--;
1835 : else {
1836 2086 : swap(gel(A,j), gel(A,def));
1837 2086 : swap(gel(I,j), gel(I,def));
1838 : }
1839 67006 : for ( ; j; j--)
1840 : {
1841 52658 : GEN a, b, S, T, S0, T0 = gel(A,j);
1842 52658 : b = gel(T0,i); if (isintzero(b)) continue;
1843 :
1844 14488 : S0 = gel(A,def); a = gel(S0,i);
1845 14488 : d = nfbezout(nf, a,b, gel(I,def),gel(I,j), &u,&v,&w,&di,0);
1846 14488 : S = colcomb(nf, u,v, S0,T0);
1847 14488 : T = colcomb(nf, a,gneg(b), T0,S0);
1848 14488 : if (u != gen_0 && v != gen_0) /* already reduced otherwise */
1849 1022 : nfcleanmod(nf, S, i, idealmul(nf,D,di));
1850 14488 : nfcleanmod(nf, T, i, idealdiv(nf,D,w));
1851 14488 : gel(A,def) = S; gel(A,j) = T;
1852 14488 : gel(I,def) = d; gel(I,j) = w;
1853 : }
1854 14348 : if (gc_needed(av,2))
1855 : {
1856 0 : if(DEBUGMEM>1) pari_warn(warnmem,"[1]: nfhnfmod, i = %ld", i);
1857 0 : gerepileall(av,dA? 4: 3, &A,&I,&D,&dA);
1858 : }
1859 : }
1860 3121 : def--; d0 = D;
1861 3121 : A += def; A[0] = evaltyp(t_MAT)|evallg(li);
1862 3121 : I += def; I[0] = evaltyp(t_VEC)|evallg(li);
1863 3121 : J = cgetg(li,t_VEC);
1864 17469 : for (i=li-1; i>=1; i--)
1865 : {
1866 14348 : GEN b = gcoeff(A,i,i);
1867 14348 : d = nfbezout(nf, gen_1,b, d0,gel(I,i), &u,&v,&w,&di,0);
1868 14348 : p1 = nfC_nf_mul(nf,gel(A,i),v);
1869 14348 : if (i > 1)
1870 : {
1871 11227 : d0 = idealmul(nf,d0,di);
1872 11227 : nfcleanmod(nf, p1, i, d0);
1873 : }
1874 14348 : gel(A,i) = p1; gel(p1,i) = gen_1;
1875 14348 : gel(I,i) = d;
1876 14348 : gel(J,i) = di;
1877 : }
1878 14348 : for (i=li-2; i>=1; i--)
1879 : {
1880 11227 : d = gel(I,i);
1881 46275 : for (j=i+1; j<li; j++)
1882 : {
1883 35048 : GEN c = gcoeff(A,i,j); if (isintzero(c)) continue;
1884 13348 : c = element_close(nf, c, idealmul(nf,d,gel(J,j)));
1885 13348 : gel(A,j) = colcomb1(nf, gneg(c), gel(A,j),gel(A,i));
1886 : }
1887 11227 : if (gc_needed(av,2))
1888 : {
1889 0 : if(DEBUGMEM>1) pari_warn(warnmem,"[2]: nfhnfmod, i = %ld", i);
1890 0 : gerepileall(av,dA? 4: 3, &A,&I,&J,&dA);
1891 : }
1892 : }
1893 3121 : idV_simplify(I);
1894 3121 : if (dA) I = gdiv(I,dA);
1895 3121 : return gerepilecopy(av0, mkvec2(A, I));
1896 : }
|
