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