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 : /* NUMBER FIELDS */
18 : /* */
19 : /**************************************************************/
20 : #include "pari.h"
21 : #include "paripriv.h"
22 :
23 : #define DEBUGLEVEL DEBUGLEVEL_nf
24 :
25 : int new_galois_format = 0;
26 :
27 : /* v a t_VEC, lg(v) = 13, sanity check for true rnf */
28 : static int
29 289004 : v13checkrnf(GEN v)
30 289004 : { return typ(gel(v,6)) == t_VEC; }
31 : static int
32 46410 : rawcheckbnf(GEN v) { return typ(v)==t_VEC && lg(v)==11; }
33 : static int
34 45430 : rawchecknf(GEN v) { return typ(v)==t_VEC && lg(v)==10; }
35 : /* v a t_VEC, lg(v) = 11, sanity check for true bnf */
36 : static int
37 37247 : v11checkbnf(GEN v) { return rawchecknf(bnf_get_nf(v)); }
38 : /* v a t_VEC, lg(v) = 13, sanity check for true nf and true bnf */
39 : static int
40 798 : v13checkgchar(GEN v) { return rawchecknf(gchar_get_nf(v)) && rawcheckbnf(gchar_get_bnf(v)); }
41 : /* v a t_VEC, lg(v) = 10, sanity check for true nf */
42 : static int
43 39480 : v10checknf(GEN v) { return typ(gel(v,1))==t_POL; }
44 : /* v a t_VEC, lg(v) = 9, sanity check for true gal */
45 : static int
46 658 : v9checkgal(GEN v)
47 658 : { GEN x = gel(v,2); return typ(x) == t_VEC && lg(x) == 4; }
48 :
49 : int
50 298359 : checkrnf_i(GEN rnf)
51 298359 : { return (typ(rnf)==t_VEC && lg(rnf)==13 && v13checkrnf(rnf)); }
52 :
53 : void
54 288416 : checkrnf(GEN rnf)
55 288416 : { if (!checkrnf_i(rnf)) pari_err_TYPE("checkrnf",rnf); }
56 :
57 : GEN
58 5238068 : checkbnf_i(GEN X)
59 : {
60 5238068 : if (typ(X) == t_VEC)
61 5237902 : switch (lg(X))
62 : {
63 5234192 : case 11:
64 5234192 : if (typ(gel(X,6)) != t_INT) return NULL; /* pre-2.2.4 format */
65 5234192 : if (lg(gel(X,10)) != 4) return NULL; /* pre-2.8.1 format */
66 5234192 : return X;
67 3240 : case 7: return checkbnf_i(bnr_get_bnf(X));
68 : }
69 636 : return NULL;
70 : }
71 :
72 : GEN
73 78195691 : checknf_i(GEN X)
74 : {
75 78195691 : if (typ(X)==t_VEC)
76 77546648 : switch(lg(X))
77 : {
78 77065822 : case 10: return X;
79 473942 : case 11: return checknf_i(bnf_get_nf(X));
80 2884 : case 7: return checknf_i(bnr_get_bnf(X));
81 2212 : case 3: if (typ(gel(X,2)) == t_POLMOD) return checknf_i(gel(X,1));
82 : }
83 653036 : return NULL;
84 : }
85 :
86 : GEN
87 3230123 : checkbnf(GEN x)
88 : {
89 3230123 : GEN bnf = checkbnf_i(x);
90 3230111 : if (!bnf) pari_err_TYPE("checkbnf [please apply bnfinit()]",x);
91 3230090 : return bnf;
92 : }
93 :
94 : GEN
95 76312837 : checknf(GEN x)
96 : {
97 76312837 : GEN nf = checknf_i(x);
98 76312623 : if (!nf) pari_err_TYPE("checknf [please apply nfinit()]",x);
99 76312594 : return nf;
100 : }
101 :
102 : GEN
103 2003046 : checkbnr_i(GEN bnr)
104 : {
105 2003046 : if (typ(bnr)!=t_VEC || lg(bnr)!=7 || !checkbnf_i(bnr_get_bnf(bnr)))
106 84 : return NULL;
107 2002963 : return bnr;
108 : }
109 : void
110 1819959 : checkbnr(GEN bnr)
111 : {
112 1819959 : if (!checkbnr_i(bnr))
113 14 : pari_err_TYPE("checkbnr [please apply bnrinit()]",bnr);
114 1819945 : }
115 :
116 : void
117 0 : checksqmat(GEN x, long N)
118 : {
119 0 : if (typ(x)!=t_MAT) pari_err_TYPE("checksqmat",x);
120 0 : if (lg(x) == 1 || lgcols(x) != N+1) pari_err_DIM("checksqmat");
121 0 : }
122 :
123 : GEN
124 1636413 : checkbid_i(GEN bid)
125 : {
126 : GEN f;
127 1636413 : if (typ(bid)!=t_VEC || lg(bid)!=6 || typ(bid_get_U(bid)) != t_VEC)
128 254112 : return NULL;
129 1382301 : f = bid_get_mod(bid);
130 1382297 : if (typ(f)!=t_VEC || lg(f)!=3) return NULL;
131 1382297 : return bid;
132 : }
133 : void
134 1379757 : checkbid(GEN bid)
135 : {
136 1379757 : if (!checkbid_i(bid)) pari_err_TYPE("checkbid",bid);
137 1379750 : }
138 : void
139 19810 : checkabgrp(GEN v)
140 : {
141 19810 : if (typ(v) == t_VEC) switch(lg(v))
142 : {
143 17423 : case 4: if (typ(gel(v,3)) != t_VEC) break;
144 19810 : case 3: if (typ(gel(v,2)) != t_VEC) break;
145 19782 : if (typ(gel(v,1)) != t_INT) break;
146 19782 : return;/*OK*/
147 0 : default: break;
148 : }
149 28 : pari_err_TYPE("checkabgrp",v);
150 : }
151 :
152 : GEN
153 267800 : checknfelt_mod(GEN nf, GEN x, const char *s)
154 : {
155 267800 : GEN T = gel(x,1), a = gel(x,2), Tnf = nf_get_pol(nf);
156 267797 : if (!RgX_equal_var(T, Tnf)) pari_err_MODULUS(s, T, Tnf);
157 267754 : return a;
158 : }
159 :
160 : int
161 10645 : check_ZKmodule_i(GEN M)
162 : {
163 10645 : return (typ(M) ==t_VEC && lg(M) >= 3
164 10645 : && typ(gel(M,1)) == t_MAT
165 10645 : && typ(gel(M,2)) == t_VEC
166 21290 : && lgcols(M) == lg(gel(M,2)));
167 : }
168 : void
169 10589 : check_ZKmodule(GEN M, const char *s)
170 10589 : { if (!check_ZKmodule_i(M)) pari_err_TYPE(s, M); }
171 :
172 : static long
173 110810 : typv6(GEN x)
174 : {
175 110810 : if (typ(gel(x,1)) == t_VEC && lg(gel(x,3)) == 3)
176 : {
177 12705 : GEN t = gel(x,3);
178 12705 : if (typ(t) != t_VEC) return typ_NULL;
179 12705 : t = gel(x,5);
180 12705 : switch(typ(gel(x,5)))
181 : {
182 448 : case t_VEC: return typ_BID;
183 12257 : case t_MAT: return typ_BIDZ;
184 0 : default: return typ_NULL;
185 : }
186 : }
187 98105 : if (typ(gel(x,2)) == t_COL && typ(gel(x,3)) == t_INT) return typ_PRID;
188 196 : return typ_NULL;
189 : }
190 :
191 : GEN
192 22071 : get_bnf(GEN x, long *t)
193 : {
194 22071 : switch(typ(x))
195 : {
196 56 : case t_POL: *t = typ_POL; return NULL;
197 56 : case t_QUAD: *t = typ_Q ; return NULL;
198 21448 : case t_VEC:
199 21448 : switch(lg(x))
200 : {
201 4445 : case 5: if (typ(gel(x,1)) != t_INT) break;
202 4389 : *t = typ_QUA; return NULL;
203 6265 : case 6: *t = typv6(x); return NULL;
204 2373 : case 7: *t = typ_BNR;
205 2373 : x = bnr_get_bnf(x);
206 2373 : if (!rawcheckbnf(x)) break;
207 2317 : return x;
208 84 : case 9:
209 84 : if (!v9checkgal(x)) break;
210 84 : *t = typ_GAL; return NULL;
211 448 : case 10:
212 448 : if (!v10checknf(x)) break;
213 448 : *t = typ_NF; return NULL;
214 469 : case 11:
215 469 : if (!v11checkbnf(x)) break;
216 469 : *t = typ_BNF; return x;
217 238 : case 13:
218 238 : if (v13checkgchar(x)) { *t = typ_GCHAR; return gchar_get_bnf(x); }
219 56 : if (!v13checkrnf(x)) break;
220 56 : *t = typ_RNF; return NULL;
221 266 : case 17: *t = typ_ELL; return NULL;
222 : }
223 6972 : break;
224 112 : case t_COL:
225 112 : if (get_prid(x)) { *t = typ_MODPR; return NULL; }
226 56 : break;
227 : }
228 7427 : *t = typ_NULL; return NULL;
229 : }
230 :
231 : GEN
232 115507 : get_nf(GEN x, long *t)
233 : {
234 115507 : switch(typ(x))
235 : {
236 133 : case t_POL : *t = typ_POL; return NULL;
237 133 : case t_QUAD: *t = typ_Q ; return NULL;
238 112560 : case t_VEC:
239 112560 : switch(lg(x))
240 : {
241 140 : case 3:
242 140 : if (typ(gel(x,2)) != t_POLMOD) break;
243 140 : return get_nf(gel(x,1),t);
244 266 : case 5:
245 266 : if (typ(gel(x,1)) != t_INT) break;
246 133 : *t = typ_QUA; return NULL;
247 98945 : case 6: *t = typv6(x); return NULL;
248 7301 : case 7:
249 7301 : x = bnr_get_bnf(x);
250 7301 : if (!rawcheckbnf(x) || !rawchecknf(x = bnf_get_nf(x))) break;
251 7168 : *t = typ_BNR; return x;
252 567 : case 9:
253 567 : if (!v9checkgal(x)) break;
254 567 : *t = typ_GAL; return NULL;
255 1120 : case 10:
256 1120 : if (!v10checknf(x)) break;
257 1120 : *t = typ_NF; return x;
258 217 : case 11:
259 217 : if (!rawchecknf(x = bnf_get_nf(x))) break;
260 217 : *t = typ_BNF; return x;
261 385 : case 13:
262 385 : if (v13checkgchar(x)) { *t = typ_GCHAR; return gchar_get_nf(x); }
263 357 : if (!v13checkrnf(x)) break;
264 357 : *t = typ_RNF; return NULL;
265 3486 : case 17: *t = typ_ELL; return NULL;
266 : }
267 399 : break;
268 133 : case t_QFB: *t = typ_QFB; return NULL;
269 266 : case t_COL:
270 266 : if (get_prid(x)) { *t = typ_MODPR; return NULL; }
271 133 : break;
272 : }
273 2814 : *t = typ_NULL; return NULL;
274 : }
275 :
276 : long
277 80605 : nftyp(GEN x)
278 : {
279 80605 : switch(typ(x))
280 : {
281 28 : case t_POL : return typ_POL;
282 7 : case t_QUAD: return typ_Q;
283 80563 : case t_VEC:
284 80563 : switch(lg(x))
285 : {
286 175 : case 13:
287 175 : if (v13checkgchar(x)) return typ_GCHAR;
288 168 : if (!v13checkrnf(x)) break;
289 168 : return typ_RNF;
290 37912 : case 10:
291 37912 : if (!v10checknf(x)) break;
292 37905 : return typ_NF;
293 273 : case 11:
294 273 : if (!v11checkbnf(x)) break;
295 273 : return typ_BNF;
296 36519 : case 7:
297 36519 : x = bnr_get_bnf(x);
298 36519 : if (!rawcheckbnf(x) || !v11checkbnf(x)) break;
299 36505 : return typ_BNR;
300 5600 : case 6:
301 5600 : return typv6(x);
302 7 : case 9:
303 7 : if (!v9checkgal(x)) break;
304 0 : return typ_GAL;
305 7 : case 17: return typ_ELL;
306 : }
307 14 : }
308 105 : return typ_NULL;
309 : }
310 :
311 : /*************************************************************************/
312 : /** **/
313 : /** GALOIS GROUP **/
314 : /** **/
315 : /*************************************************************************/
316 :
317 : GEN
318 7597 : tschirnhaus(GEN x)
319 : {
320 7597 : pari_sp av = avma, av2;
321 7597 : long a, v = varn(x);
322 7597 : GEN u, y = cgetg(5,t_POL);
323 :
324 7597 : if (typ(x)!=t_POL) pari_err_TYPE("tschirnhaus",x);
325 7597 : if (lg(x) < 4) pari_err_CONSTPOL("tschirnhaus");
326 7597 : if (v) { u = leafcopy(x); setvarn(u,0); x=u; }
327 7597 : y[1] = evalsigne(1)|evalvarn(0);
328 : do
329 : {
330 8303 : a = random_bits(2); if (a==0) a = 1; gel(y,4) = stoi(a);
331 8303 : a = random_bits(3); if (a>=4) a -= 8; gel(y,3) = stoi(a);
332 8303 : a = random_bits(3); if (a>=4) a -= 8; gel(y,2) = stoi(a);
333 8303 : u = RgXQ_charpoly(y,x,v); av2 = avma;
334 : }
335 8303 : while (degpol(RgX_gcd(u,RgX_deriv(u)))); /* while u not separable */
336 7597 : if (DEBUGLEVEL>1)
337 0 : err_printf("Tschirnhaus transform. New pol: %Ps",u);
338 7597 : set_avma(av2); return gerepileupto(av,u);
339 : }
340 :
341 : /* Assume pol in Z[X], monic of degree n. Find L in Z such that
342 : * POL = L^(-n) pol(L x) is monic in Z[X]. Return POL and set *ptk = L.
343 : * No GC. */
344 : GEN
345 843686 : ZX_Z_normalize(GEN pol, GEN *ptk)
346 : {
347 843686 : long i,j, sk, n = degpol(pol); /* > 0 */
348 : GEN k, fa, P, E, a, POL;
349 :
350 843698 : if (ptk) *ptk = gen_1;
351 843698 : if (!n) return pol;
352 843691 : a = pol + 2; k = gel(a,n-1); /* a[i] = coeff of degree i */
353 1623710 : for (i = n-2; i >= 0; i--)
354 : {
355 1368079 : k = gcdii(k, gel(a,i));
356 1367862 : if (is_pm1(k)) return pol;
357 : }
358 255631 : sk = signe(k);
359 255631 : if (!sk) return pol; /* monomial! */
360 253706 : fa = absZ_factor_limit(k, 0); k = gen_1;
361 253745 : P = gel(fa,1);
362 253745 : E = gel(fa,2);
363 253745 : POL = leafcopy(pol); a = POL+2;
364 544907 : for (i = lg(P)-1; i > 0; i--)
365 : {
366 291164 : GEN p = gel(P,i), pv, pvj;
367 291164 : long vmin = itos(gel(E,i));
368 : /* find v_p(k) = min floor( v_p(a[i]) / (n-i)) */
369 1204948 : for (j=n-1; j>=0; j--)
370 : {
371 : long v;
372 913787 : if (!signe(gel(a,j))) continue;
373 785521 : v = Z_pval(gel(a,j), p) / (n - j);
374 785518 : if (v < vmin) vmin = v;
375 : }
376 291161 : if (!vmin) continue;
377 18188 : pvj = pv = powiu(p,vmin); k = mulii(k, pv);
378 : /* a[j] /= p^(v*(n-j)) */
379 73972 : for (j=n-1; j>=0; j--)
380 : {
381 55791 : if (j < n-1) pvj = mulii(pvj, pv);
382 55790 : gel(a,j) = diviiexact(gel(a,j), pvj);
383 : }
384 : }
385 253743 : if (ptk) *ptk = k;
386 253743 : return POL;
387 : }
388 :
389 : /* Assume pol != 0 in Z[X]. Find C in Q, L in Z such that POL = C pol(x/L) monic
390 : * in Z[X]. Return POL and set *pL = L. Wasteful (but correct) if pol is not
391 : * primitive: better if caller used Q_primpart already. No GC. */
392 : GEN
393 841078 : ZX_primitive_to_monic(GEN pol, GEN *pL)
394 : {
395 841078 : long i,j, n = degpol(pol);
396 841078 : GEN lc = leading_coeff(pol), L, fa, P, E, a, POL;
397 :
398 841075 : if (is_pm1(lc))
399 : {
400 838194 : if (pL) *pL = gen_1;
401 838194 : return signe(lc) < 0? ZX_neg(pol): pol;
402 : }
403 2877 : if (signe(lc) < 0)
404 35 : POL = ZX_neg(pol);
405 : else
406 2842 : POL = leafcopy(pol);
407 2877 : a = POL+2; lc = gel(a,n);
408 2877 : fa = absZ_factor_limit(lc,0); L = gen_1;
409 2877 : P = gel(fa,1);
410 2877 : E = gel(fa,2);
411 5901 : for (i = lg(P)-1; i > 0; i--)
412 : {
413 3024 : GEN p = gel(P,i), pk, pku;
414 3024 : long v, j0, e = itos(gel(E,i)), k = e/n, d = k*n - e;
415 :
416 3024 : if (d < 0) { k++; d += n; }
417 : /* k = ceil(e[i] / n); find d, k such that p^d pol(x / p^k) monic */
418 17115 : for (j=n-1; j>0; j--)
419 : {
420 14091 : if (!signe(gel(a,j))) continue;
421 8694 : v = Z_pval(gel(a,j), p);
422 8764 : while (v + d < k * j) { k++; d += n; }
423 : }
424 3024 : pk = powiu(p,k); j0 = d/k;
425 3024 : L = mulii(L, pk);
426 :
427 3024 : pku = powiu(p,d - k*j0);
428 : /* a[j] *= p^(d - kj) */
429 7091 : for (j=j0; j>=0; j--)
430 : {
431 4067 : if (j < j0) pku = mulii(pku, pk);
432 4067 : gel(a,j) = mulii(gel(a,j), pku);
433 : }
434 3024 : j0++;
435 3024 : pku = powiu(p,k*j0 - d);
436 : /* a[j] /= p^(kj - d) */
437 19096 : for (j=j0; j<=n; j++)
438 : {
439 16072 : if (j > j0) pku = mulii(pku, pk);
440 16072 : gel(a,j) = diviiexact(gel(a,j), pku);
441 : }
442 : }
443 2877 : if (pL) *pL = L;
444 2877 : return POL;
445 : }
446 : /* Assume pol != 0 in Z[X]. Find C,L in Q such that POL = C pol(x/L)
447 : * monic in Z[X]. Return POL and set *pL = L.
448 : * Wasteful (but correct) if pol is not primitive: better if caller used
449 : * Q_primpart already. No GC. */
450 : GEN
451 840760 : ZX_Q_normalize(GEN pol, GEN *pL)
452 : {
453 840760 : GEN lc, POL = ZX_primitive_to_monic(pol, &lc);
454 840757 : POL = ZX_Z_normalize(POL, pL);
455 840651 : if (pL) *pL = gdiv(lc, *pL);
456 840767 : return POL;
457 : }
458 :
459 : GEN
460 5751439 : ZX_Q_mul(GEN A, GEN z)
461 : {
462 5751439 : pari_sp av = avma;
463 5751439 : long i, l = lg(A);
464 : GEN d, n, Ad, B, u;
465 5751439 : if (typ(z)==t_INT) return ZX_Z_mul(A,z);
466 4828323 : n = gel(z, 1); d = gel(z, 2);
467 4828323 : Ad = RgX_to_RgC(FpX_red(A, d), l-2);
468 4828323 : u = gcdii(d, FpV_factorback(Ad, NULL, d));
469 4828317 : B = cgetg(l, t_POL);
470 4828321 : B[1] = A[1];
471 4828321 : if (equali1(u))
472 : {
473 4659313 : for(i=2; i<l; i++)
474 2884216 : gel(B, i) = mkfrac(mulii(n, gel(A,i)), d);
475 : } else
476 : {
477 18456130 : for(i=2; i<l; i++)
478 : {
479 15402907 : GEN di = gcdii(gel(Ad, i-1), u);
480 15402874 : GEN ni = mulii(n, diviiexact(gel(A,i), di));
481 15402862 : if (equalii(d, di))
482 3630004 : gel(B, i) = ni;
483 : else
484 11772896 : gel(B, i) = mkfrac(ni, diviiexact(d, di));
485 : }
486 : }
487 4828320 : return gerepilecopy(av, B);
488 : }
489 :
490 : /* T != 0 in Z[x], returns a monic polynomial U in Z[x] generating the
491 : * same field: there exist C in Q, L in Z such that U(x) = C T(x/L).
492 : * Set *L = NULL if L = 1, and to L otherwise. No garbage collecting. */
493 : GEN
494 0 : ZX_to_monic(GEN T, GEN *L)
495 : {
496 0 : GEN lc = leading_coeff(T);
497 0 : if (is_pm1(lc)) { *L = gen_1; return signe(lc) > 0? T: ZX_neg(T); }
498 0 : return ZX_primitive_to_monic(Q_primpart(T), L);
499 : }
500 :
501 : GEN
502 42 : poltomonic(GEN T, GEN *L)
503 : {
504 42 : pari_sp av = avma;
505 42 : if (typ(T) != t_POL || !RgX_is_QX(T)) pari_err_TYPE("poltomonic", T);
506 42 : if (degpol(T) < 0) pari_err_ROOTS0("poltomonic");
507 35 : T = ZX_Q_normalize(Q_primpart(T), L); return gc_all(av, L? 2: 1, &T, L);
508 : }
509 :
510 : GEN
511 8897 : ZXX_Q_mul(GEN A, GEN z)
512 : {
513 : long i, l;
514 : GEN B;
515 8897 : if (typ(z)==t_INT) return ZXX_Z_mul(A,z);
516 8365 : B = cgetg_copy(A, &l);
517 8365 : B[1] = A[1];
518 116767 : for (i=2; i<l; i++)
519 : {
520 108402 : GEN Ai = gel(A,i);
521 108402 : gel(B,i) = typ(Ai)==t_POL ? ZX_Q_mul(Ai, z): gmul(Ai, z);
522 : }
523 8365 : return B;
524 : }
525 :
526 : /* Evaluate pol in s using nfelt arithmetic and Horner rule */
527 : GEN
528 11151 : nfpoleval(GEN nf, GEN pol, GEN s)
529 : {
530 11151 : pari_sp av=avma;
531 11151 : long i=lg(pol)-1;
532 : GEN res;
533 11151 : if (i==1) return gen_0;
534 11151 : res = nf_to_scalar_or_basis(nf, gel(pol,i));
535 27944 : for (i-- ; i>=2; i--)
536 16793 : res = nfadd(nf, nfmul(nf, s, res), gel(pol,i));
537 11151 : return gerepileupto(av, res);
538 : }
539 :
540 : static GEN
541 94972 : QX_table_nfpoleval(GEN nf, GEN pol, GEN m)
542 : {
543 94972 : pari_sp av = avma;
544 94972 : long i = lg(pol)-1;
545 : GEN res, den;
546 94972 : if (i==1) return gen_0;
547 94972 : pol = Q_remove_denom(pol, &den);
548 94972 : res = scalarcol_shallow(gel(pol,i), nf_get_degree(nf));
549 393376 : for (i-- ; i>=2; i--)
550 298403 : res = ZC_Z_add(ZM_ZC_mul(m, res), gel(pol,i));
551 94973 : if (den) res = RgC_Rg_div(res, den);
552 94972 : return gerepileupto(av, res);
553 : }
554 :
555 : GEN
556 8722 : FpX_FpC_nfpoleval(GEN nf, GEN pol, GEN a, GEN p)
557 : {
558 8722 : pari_sp av=avma;
559 8722 : long i=lg(pol)-1, n=nf_get_degree(nf);
560 : GEN res, Ma;
561 8722 : if (i==1) return zerocol(n);
562 8722 : Ma = FpM_red(zk_multable(nf, a), p);
563 8722 : res = scalarcol(gel(pol,i),n);
564 67900 : for (i-- ; i>=2; i--)
565 : {
566 59178 : res = FpM_FpC_mul(Ma, res, p);
567 59178 : gel(res,1) = Fp_add(gel(res,1), gel(pol,i), p);
568 : }
569 8722 : return gerepileupto(av, res);
570 : }
571 :
572 : /* compute s(x), not stack clean */
573 : static GEN
574 95159 : ZC_galoisapply(GEN nf, GEN s, GEN x)
575 : {
576 95159 : x = nf_to_scalar_or_alg(nf, x);
577 95158 : if (typ(x) != t_POL) return scalarcol(x, nf_get_degree(nf));
578 94969 : return QX_table_nfpoleval(nf, x, zk_multable(nf, s));
579 : }
580 :
581 : /* true nf; S = automorphism in basis form, return an FpC = S(z) mod p */
582 : GEN
583 6153 : zk_galoisapplymod(GEN nf, GEN z, GEN S, GEN p)
584 : {
585 : GEN den, pe, pe1, denpe, R;
586 :
587 6153 : z = nf_to_scalar_or_alg(nf, z);
588 6153 : if (typ(z) != t_POL) return z;
589 6153 : if (gequalX(z)) return FpC_red(S, p); /* common, e.g. modpr_genFq */
590 5775 : z = Q_remove_denom(z,&den);
591 5775 : denpe = pe = NULL;
592 5775 : pe1 = p;
593 5775 : if (den)
594 : {
595 5110 : ulong e = Z_pvalrem(den, p, &den);
596 5110 : if (e) { pe = powiu(p, e); pe1 = mulii(pe, p); }
597 5110 : denpe = Zp_inv(den, p, e+1);
598 : }
599 5775 : R = FpX_FpC_nfpoleval(nf, FpX_red(z, pe1), FpC_red(S, pe1), pe1);
600 5775 : if (denpe) R = FpC_Fp_mul(R, denpe, pe1);
601 5775 : if (pe) R = gdivexact(R, pe);
602 5775 : return R;
603 : }
604 :
605 : /* true nf */
606 : static GEN
607 7 : pr_make(GEN nf, GEN p, GEN u, GEN e, GEN f)
608 : {
609 7 : GEN t = FpM_deplin(zk_multable(nf, u), p);
610 7 : t = zk_scalar_or_multable(nf, t);
611 7 : return mkvec5(p, u, e, f, t);
612 : }
613 : static GEN
614 7 : pr_galoisapply(GEN nf, GEN pr, GEN aut)
615 : {
616 7 : GEN p = pr_get_p(pr), u = zk_galoisapplymod(nf, pr_get_gen(pr), aut, p);
617 7 : return pr_make(nf, p, u, gel(pr,3), gel(pr,4));
618 : }
619 : static GEN
620 0 : pr_galoismatrixapply(GEN nf, GEN pr, GEN M)
621 : {
622 0 : GEN p = pr_get_p(pr), u = FpC_red(ZM_ZC_mul(M, pr_get_gen(pr)), p);
623 0 : return pr_make(nf, p, u, gel(pr,3), gel(pr,4));
624 : }
625 :
626 : static GEN
627 7 : vecgaloisapply(GEN nf, GEN aut, GEN x)
628 21 : { pari_APPLY_same(galoisapply(nf, aut, gel(x,i))); }
629 : static GEN
630 0 : vecgaloismatrixapply(GEN nf, GEN aut, GEN x)
631 0 : { pari_APPLY_same(nfgaloismatrixapply(nf, aut, gel(x,i))); }
632 :
633 : /* x: famat or standard algebraic number, aut automorphism in ZC form
634 : * simplified from general galoisapply */
635 : static GEN
636 49 : elt_galoisapply(GEN nf, GEN aut, GEN x)
637 : {
638 49 : pari_sp av = avma;
639 49 : switch(typ(x))
640 : {
641 7 : case t_INT: return icopy(x);
642 7 : case t_FRAC: return gcopy(x);
643 7 : case t_POLMOD: x = gel(x,2); /* fall through */
644 14 : case t_POL: {
645 14 : GEN y = basistoalg(nf, ZC_galoisapply(nf, aut, x));
646 14 : return gerepileupto(av,y);
647 : }
648 7 : case t_COL:
649 7 : return gerepileupto(av, ZC_galoisapply(nf, aut, x));
650 14 : case t_MAT:
651 14 : switch(lg(x)) {
652 7 : case 1: return cgetg(1, t_MAT);
653 7 : case 3: retmkmat2(vecgaloisapply(nf,aut,gel(x,1)), ZC_copy(gel(x,2)));
654 : }
655 : }
656 0 : pari_err_TYPE("galoisapply",x);
657 : return NULL; /* LCOV_EXCL_LINE */
658 : }
659 : /* M automorphism in matrix form */
660 : static GEN
661 0 : elt_galoismatrixapply(GEN nf, GEN M, GEN x)
662 : {
663 0 : if (typ(x) == t_MAT)
664 0 : switch(lg(x)) {
665 0 : case 1: return cgetg(1, t_MAT);
666 0 : case 3: retmkmat2(vecgaloismatrixapply(nf,M,gel(x,1)), ZC_copy(gel(x,2)));
667 : }
668 0 : return nfgaloismatrixapply(nf, M, x);
669 : }
670 :
671 : GEN
672 126422 : galoisapply(GEN nf, GEN aut, GEN x)
673 : {
674 126422 : pari_sp av = avma;
675 : long lx;
676 : GEN y;
677 :
678 126422 : nf = checknf(nf);
679 126422 : switch(typ(x))
680 : {
681 378 : case t_INT: return icopy(x);
682 7 : case t_FRAC: return gcopy(x);
683 :
684 35 : case t_POLMOD: x = gel(x,2); /* fall through */
685 1302 : case t_POL:
686 1302 : aut = algtobasis(nf, aut);
687 1302 : y = basistoalg(nf, ZC_galoisapply(nf, aut, x));
688 1302 : return gerepileupto(av,y);
689 :
690 56 : case t_VEC:
691 56 : aut = algtobasis(nf, aut);
692 56 : switch(lg(x))
693 : {
694 7 : case 6:
695 7 : if (pr_is_inert(x)) { set_avma(av); return gcopy(x); }
696 7 : return gerepilecopy(av, pr_galoisapply(nf, x, aut));
697 49 : case 3: y = cgetg(3,t_VEC);
698 49 : gel(y,1) = galoisapply(nf, aut, gel(x,1));
699 49 : gel(y,2) = elt_galoisapply(nf, aut, gel(x,2));
700 49 : return gerepileupto(av, y);
701 : }
702 0 : break;
703 :
704 88601 : case t_COL:
705 88601 : aut = algtobasis(nf, aut);
706 88601 : return gerepileupto(av, ZC_galoisapply(nf, aut, x));
707 :
708 36078 : case t_MAT: /* ideal */
709 36078 : lx = lg(x); if (lx==1) return cgetg(1,t_MAT);
710 36078 : if (nbrows(x) != nf_get_degree(nf)) break;
711 36078 : y = RgM_mul(nfgaloismatrix(nf,aut), x);
712 36078 : return gerepileupto(av, idealhnf_shallow(nf,y));
713 : }
714 0 : pari_err_TYPE("galoisapply",x);
715 : return NULL; /* LCOV_EXCL_LINE */
716 : }
717 :
718 : /* M automorphism in galoismatrix form */
719 : GEN
720 13993 : nfgaloismatrixapply(GEN nf, GEN M, GEN x)
721 : {
722 13993 : pari_sp av = avma;
723 : long lx;
724 : GEN y;
725 :
726 13993 : nf = checknf(nf);
727 13993 : switch(typ(x))
728 : {
729 735 : case t_INT: return icopy(x);
730 0 : case t_FRAC: return gcopy(x);
731 :
732 0 : case t_POLMOD: x = gel(x,2); /* fall through */
733 0 : case t_POL:
734 0 : x = algtobasis(nf, x);
735 0 : return gerepileupto(av, basistoalg(nf, RgM_RgC_mul(M, x)));
736 :
737 0 : case t_VEC:
738 0 : switch(lg(x))
739 : {
740 0 : case 6:
741 0 : if (pr_is_inert(x)) { set_avma(av); return gcopy(x); }
742 0 : return gerepilecopy(av, pr_galoismatrixapply(nf, x, M));
743 0 : case 3: y = cgetg(3,t_VEC);
744 0 : gel(y,1) = nfgaloismatrixapply(nf, M, gel(x,1));
745 0 : gel(y,2) = elt_galoismatrixapply(nf, M, gel(x,2));
746 0 : return gerepileupto(av, y);
747 : }
748 0 : break;
749 :
750 6111 : case t_COL: return RgM_RgC_mul(M, x);
751 :
752 7147 : case t_MAT: /* ideal */
753 7147 : lx = lg(x); if (lx==1) return cgetg(1,t_MAT);
754 7147 : if (nbrows(x) != nf_get_degree(nf)) break;
755 7147 : return gerepileupto(av, idealhnf_shallow(nf,RgM_mul(M, x)));
756 : }
757 0 : pari_err_TYPE("galoisapply",x);
758 : return NULL; /* LCOV_EXCL_LINE */
759 : }
760 :
761 : /* compute action of automorphism s on nf.zk */
762 : GEN
763 94514 : nfgaloismatrix(GEN nf, GEN s)
764 : {
765 94514 : pari_sp av2, av = avma;
766 : GEN zk, D, M, H, m;
767 : long k, n;
768 :
769 94514 : nf = checknf(nf);
770 94514 : zk = nf_get_zkprimpart(nf); n = lg(zk)-1;
771 94514 : M = cgetg(n+1, t_MAT);
772 94513 : gel(M,1) = col_ei(n, 1); /* s(1) = 1 */
773 94514 : if (n == 1) return M;
774 94514 : av2 = avma;
775 94514 : if (typ(s) != t_COL) s = algtobasis(nf, s);
776 94514 : D = nf_get_zkden(nf);
777 94514 : H = RgV_to_RgM(zk, n);
778 94514 : if (n == 2)
779 : {
780 62958 : GEN t = gel(H,2); /* D * s(w_2) */
781 62958 : t = ZC_Z_add(ZC_Z_mul(s, gel(t,2)), gel(t,1));
782 62957 : gel(M,2) = gerepileupto(av2, gdiv(t, D));
783 62957 : return M;
784 : }
785 31556 : m = zk_multable(nf, s);
786 31555 : gel(M,2) = s; /* M[,k] = s(x^(k-1)) */
787 138205 : for (k = 3; k <= n; k++) gel(M,k) = ZM_ZC_mul(m, gel(M,k-1));
788 31554 : M = ZM_mul(M, H);
789 31556 : if (!equali1(D)) M = ZM_Z_divexact(M, D);
790 31555 : return gerepileupto(av, M);
791 : }
792 :
793 : static GEN
794 8140 : get_aut(GEN nf, GEN gal, GEN aut, GEN g)
795 : {
796 8140 : return aut ? gel(aut, g[1]): poltobasis(nf, galoispermtopol(gal, g));
797 : }
798 :
799 : static GEN
800 1435 : idealquasifrob(GEN nf, GEN gal, GEN grp, GEN pr, GEN subg, GEN *S, GEN aut)
801 : {
802 1435 : pari_sp av = avma;
803 1435 : long i, n = nf_get_degree(nf), f = pr_get_f(pr);
804 1435 : GEN pi = pr_get_gen(pr);
805 5551 : for (i=1; i<=n; i++)
806 : {
807 5551 : GEN g = gel(grp,i);
808 5551 : if ((!subg && perm_orderu(g) == (ulong)f)
809 5089 : || (subg && perm_relorder(g, subg)==f))
810 : {
811 1939 : *S = get_aut(nf, gal, aut, g);
812 1939 : if (ZC_prdvd(ZC_galoisapply(nf, *S, pi), pr)) return g;
813 504 : set_avma(av);
814 : }
815 : }
816 0 : pari_err_BUG("idealquasifrob [Frobenius not found]");
817 : return NULL;/*LCOV_EXCL_LINE*/
818 : }
819 :
820 : GEN
821 1379 : nfgaloispermtobasis(GEN nf, GEN gal)
822 : {
823 1379 : GEN grp = gal_get_group(gal);
824 1379 : long i, n = lg(grp)-1;
825 1379 : GEN aut = cgetg(n+1, t_VEC);
826 15323 : for(i=1; i<=n; i++)
827 : {
828 13944 : pari_sp av = avma;
829 13944 : GEN g = gel(grp, i);
830 13944 : GEN vec = poltobasis(nf, galoispermtopol(gal, g));
831 13944 : gel(aut, g[1]) = gerepileupto(av, vec);
832 : }
833 1379 : return aut;
834 : }
835 :
836 : static void
837 2457 : gal_check_pol(const char *f, GEN x, GEN y)
838 2457 : { if (!RgX_equal_var(x,y)) pari_err_MODULUS(f,x,y); }
839 :
840 : /* true nf */
841 : GEN
842 56 : idealfrobenius_aut(GEN nf, GEN gal, GEN pr, GEN aut)
843 : {
844 56 : pari_sp av = avma;
845 56 : GEN S=NULL, g=NULL; /*-Wall*/
846 : GEN T, p, a, b, modpr;
847 : long f, n, s;
848 56 : f = pr_get_f(pr); n = nf_get_degree(nf);
849 56 : if (f==1) { set_avma(av); return identity_perm(n); }
850 56 : g = idealquasifrob(nf, gal, gal_get_group(gal), pr, NULL, &S, aut);
851 56 : if (f==2) return gerepileuptoleaf(av, g);
852 21 : modpr = zk_to_Fq_init(nf,&pr,&T,&p);
853 21 : a = pol_x(nf_get_varn(nf));
854 21 : b = nf_to_Fq(nf, zk_galoisapplymod(nf, modpr_genFq(modpr), S, p), modpr);
855 42 : for (s = 1; s < f-1; s++)
856 : {
857 21 : a = Fq_pow(a, p, T, p);
858 21 : if (ZX_equal(a, b)) break;
859 : }
860 21 : g = perm_powu(g, Fl_inv(s, f));
861 21 : return gerepileupto(av, g);
862 : }
863 :
864 : GEN
865 63 : idealfrobenius(GEN nf, GEN gal, GEN pr)
866 : {
867 63 : nf = checknf(nf);
868 63 : checkgal(gal);
869 63 : checkprid(pr);
870 63 : gal_check_pol("idealfrobenius",nf_get_pol(nf),gal_get_pol(gal));
871 63 : if (pr_get_e(pr)>1) pari_err_DOMAIN("idealfrobenius","pr.e", ">", gen_1,pr);
872 56 : return idealfrobenius_aut(nf, gal, pr, NULL);
873 : }
874 :
875 : /* true nf */
876 : GEN
877 616 : idealramfrobenius_aut(GEN nf, GEN gal, GEN pr, GEN ram, GEN aut)
878 : {
879 616 : pari_sp av = avma;
880 616 : GEN S=NULL, g=NULL; /*-Wall*/
881 : GEN T, p, a, b, modpr;
882 : GEN isog, deco;
883 : long f, n, s;
884 616 : f = pr_get_f(pr); n = nf_get_degree(nf);
885 616 : if (f==1) { set_avma(av); return identity_perm(n); }
886 399 : modpr = zk_to_Fq_init(nf,&pr,&T,&p);
887 399 : deco = group_elts(gel(ram,1), nf_get_degree(nf));
888 399 : isog = group_set(gel(ram,2), nf_get_degree(nf));
889 399 : g = idealquasifrob(nf, gal, deco, pr, isog, &S, aut);
890 399 : a = pol_x(nf_get_varn(nf));
891 399 : b = nf_to_Fq(nf, zk_galoisapplymod(nf, modpr_genFq(modpr), S, p), modpr);
892 854 : for (s=0; !ZX_equal(a, b); s++)
893 455 : a = Fq_pow(a, p, T, p);
894 399 : g = perm_powu(g, Fl_inv(s, f));
895 399 : return gerepileupto(av, g);
896 : }
897 :
898 : GEN
899 0 : idealramfrobenius(GEN nf, GEN gal, GEN pr, GEN ram)
900 : {
901 0 : return idealramfrobenius_aut(nf, gal, pr, ram, NULL);
902 : }
903 :
904 : static GEN
905 1834 : idealinertiagroup(GEN nf, GEN gal, GEN aut, GEN pr)
906 : {
907 1834 : long i, n = nf_get_degree(nf);
908 1834 : GEN p, T, modpr = zk_to_Fq_init(nf,&pr,&T,&p);
909 1834 : GEN b = modpr_genFq(modpr);
910 1834 : long e = pr_get_e(pr), coprime = ugcd(e, pr_get_f(pr)) == 1;
911 1834 : GEN grp = gal_get_group(gal), pi = pr_get_gen(pr);
912 1834 : pari_sp ltop = avma;
913 7860 : for (i=1; i<=n; i++)
914 : {
915 7860 : GEN iso = gel(grp,i);
916 7860 : if (perm_orderu(iso) == (ulong)e)
917 : {
918 3296 : GEN S = get_aut(nf, gal, aut, iso);
919 3296 : if (ZC_prdvd(ZC_galoisapply(nf, S, pi), pr)
920 2352 : && (coprime || gequalX(nf_to_Fq(nf, galoisapply(nf,S,b), modpr))))
921 1834 : return iso;
922 1462 : set_avma(ltop);
923 : }
924 : }
925 0 : pari_err_BUG("idealinertiagroup [no isotropic element]");
926 : return NULL;/*LCOV_EXCL_LINE*/
927 : }
928 :
929 : static GEN
930 1897 : idealramgroupstame(GEN nf, GEN gal, GEN aut, GEN pr)
931 : {
932 1897 : pari_sp av = avma;
933 : GEN iso, frob, giso, isog, S, res;
934 1897 : long e = pr_get_e(pr), f = pr_get_f(pr);
935 1897 : GEN grp = gal_get_group(gal);
936 1897 : if (e == 1)
937 : {
938 63 : if (f==1)
939 0 : return cgetg(1,t_VEC);
940 63 : frob = idealquasifrob(nf, gal, grp, pr, NULL, &S, aut);
941 63 : set_avma(av);
942 63 : res = cgetg(2, t_VEC);
943 63 : gel(res, 1) = cyclicgroup(frob, f);
944 63 : return res;
945 : }
946 1834 : res = cgetg(3, t_VEC);
947 1834 : av = avma;
948 1834 : iso = idealinertiagroup(nf, gal, aut, pr);
949 1834 : set_avma(av);
950 1834 : giso = cyclicgroup(iso, e);
951 1834 : gel(res, 2) = giso;
952 1834 : if (f==1)
953 : {
954 917 : gel(res, 1) = giso;
955 917 : return res;
956 : }
957 917 : av = avma;
958 917 : isog = group_set(giso, nf_get_degree(nf));
959 917 : frob = idealquasifrob(nf, gal, grp, pr, isog, &S, aut);
960 917 : set_avma(av);
961 917 : gel(res, 1) = dicyclicgroup(iso,frob,e,f);
962 917 : return res;
963 : }
964 :
965 : /* true nf, p | e */
966 : static GEN
967 497 : idealramgroupswild(GEN nf, GEN gal, GEN aut, GEN pr)
968 : {
969 497 : pari_sp av2, av = avma;
970 497 : GEN p, T, idx, g, gbas, pi, pibas, Dpi, modpr = zk_to_Fq_init(nf,&pr,&T,&p);
971 497 : long bound, i, vDpi, vDg, n = nf_get_degree(nf);
972 497 : long e = pr_get_e(pr);
973 497 : long f = pr_get_f(pr);
974 : ulong nt,rorder;
975 497 : GEN pg, ppi, grp = gal_get_group(gal);
976 :
977 : /* G_i = {s: v(s(pi) - pi) > i} trivial for i > bound;
978 : * v_pr(Diff) = sum_{i = 0}^{bound} (#G_i - 1) >= e-1 + bound*(p-1)*/
979 497 : bound = (idealval(nf, nf_get_diff(nf), pr) - (e-1)) / (itou(p)-1);
980 497 : (void) u_pvalrem(n,p,&nt);
981 497 : rorder = e*f*(n/nt);
982 497 : idx = const_vecsmall(n,-1);
983 497 : pg = NULL;
984 497 : vDg = 0;
985 497 : if (f == 1)
986 154 : g = gbas = NULL;
987 : else
988 : {
989 : GEN Dg;
990 343 : g = nf_to_scalar_or_alg(nf, modpr_genFq(modpr));
991 343 : if (!gequalX(g)) /* p | nf.index */
992 : {
993 7 : g = Q_remove_denom(g, &Dg);
994 7 : vDg = Z_pval(Dg,p);
995 7 : pg = powiu(p, vDg + 1);
996 7 : g = FpX_red(g, pg);
997 : }
998 343 : gbas = nf_to_scalar_or_basis(nf, g);
999 : }
1000 497 : pi = nf_to_scalar_or_alg(nf, pr_get_gen(pr));
1001 497 : pi = Q_remove_denom(pi, &Dpi);
1002 497 : vDpi = Dpi ? Z_pval(Dpi, p): 0;
1003 497 : ppi = powiu(p, vDpi + (bound + e)/e);
1004 497 : pi = FpX_red(pi, ppi);
1005 497 : pibas = nf_to_scalar_or_basis(nf, pi);
1006 497 : av2 = avma;
1007 4704 : for (i = 2; i <= n; i++)
1008 : {
1009 4207 : GEN S, Spi, piso, iso = gel(grp, i);
1010 4207 : long j, o, ix = iso[1];
1011 4207 : if (idx[ix] >= 0 || rorder % (o = (long)perm_orderu(iso))) continue;
1012 :
1013 2905 : piso = iso;
1014 2905 : S = get_aut(nf, gal, aut, iso);
1015 2905 : Spi = FpX_FpC_nfpoleval(nf, pi, FpC_red(S, ppi), ppi);
1016 : /* Computation made knowing that the valuation is <= bound + 1. Correct
1017 : * to maximal value if reduction mod ppi altered this */
1018 2905 : idx[ix] = minss(bound+1, idealval(nf, gsub(Spi,pibas), pr) - e*vDpi);
1019 2905 : if (idx[ix] == 0) idx[ix] = -1;
1020 2457 : else if (g)
1021 : {
1022 1848 : GEN Sg = pg? FpX_FpC_nfpoleval(nf, g, FpC_red(S, pg), pg): S;
1023 1848 : if (vDg)
1024 42 : { if (nfval(nf, gsub(Sg, gbas), pr) - e*vDg <= 0) idx[ix] = 0; }
1025 : else /* same, more efficient */
1026 1806 : { if (!ZC_prdvd(gsub(Sg, gbas), pr)) idx[ix] = 0; }
1027 : }
1028 5488 : for (j = 2; j < o; j++)
1029 : {
1030 2583 : piso = perm_mul(piso,iso);
1031 2583 : if (ugcd(j,o)==1) idx[ piso[1] ] = idx[ix];
1032 : }
1033 2905 : set_avma(av2);
1034 : }
1035 497 : return gerepileuptoleaf(av, idx);
1036 : }
1037 :
1038 : GEN
1039 2394 : idealramgroups_aut(GEN nf, GEN gal, GEN pr, GEN aut)
1040 : {
1041 2394 : pari_sp av = avma;
1042 : GEN tbl, idx, res, set, sub;
1043 : long i, j, e, n, maxm, p;
1044 : ulong et;
1045 2394 : nf = checknf(nf);
1046 2394 : checkgal(gal);
1047 2394 : checkprid(pr);
1048 2394 : gal_check_pol("idealramgroups",nf_get_pol(nf),gal_get_pol(gal));
1049 2394 : e = pr_get_e(pr); n = nf_get_degree(nf);
1050 2394 : p = itos(pr_get_p(pr));
1051 2394 : if (e%p) return idealramgroupstame(nf, gal, aut, pr);
1052 497 : (void) u_lvalrem(e,p,&et);
1053 497 : idx = idealramgroupswild(nf, gal, aut, pr);
1054 497 : sub = group_subgroups(galois_group(gal));
1055 497 : tbl = subgroups_tableset(sub, n);
1056 497 : maxm = vecsmall_max(idx)+1;
1057 497 : res = cgetg(maxm+1,t_VEC);
1058 497 : set = zero_F2v(n); F2v_set(set,1);
1059 2499 : for(i=maxm; i>0; i--)
1060 : {
1061 : long ix;
1062 20468 : for(j=1;j<=n;j++)
1063 18466 : if (idx[j]==i-1)
1064 3521 : F2v_set(set,j);
1065 2002 : ix = tableset_find_index(tbl, set);
1066 2002 : if (ix==0) pari_err_BUG("idealramgroups");
1067 2002 : gel(res,i) = gel(sub, ix);
1068 : }
1069 497 : return gerepilecopy(av, res);
1070 : }
1071 :
1072 : GEN
1073 112 : idealramgroups(GEN nf, GEN gal, GEN pr)
1074 : {
1075 112 : return idealramgroups_aut(nf, gal, pr, NULL);
1076 : }
1077 :
1078 : /* x = relative polynomial nf = absolute nf, bnf = absolute bnf */
1079 : GEN
1080 112 : get_bnfpol(GEN x, GEN *bnf, GEN *nf)
1081 : {
1082 112 : *bnf = checkbnf_i(x);
1083 112 : *nf = checknf_i(x);
1084 112 : if (*nf) x = nf_get_pol(*nf);
1085 112 : if (typ(x) != t_POL) pari_err_TYPE("get_bnfpol",x);
1086 112 : return x;
1087 : }
1088 :
1089 : GEN
1090 667331 : get_nfpol(GEN x, GEN *nf)
1091 : {
1092 667331 : if (typ(x) == t_POL) { *nf = NULL; return x; }
1093 380616 : *nf = checknf(x); return nf_get_pol(*nf);
1094 : }
1095 :
1096 : static GEN
1097 531 : incl_disc(GEN nfa, GEN a, int nolocal)
1098 : {
1099 : GEN d;
1100 531 : if (nfa) return nf_get_disc(nfa);
1101 461 : if (nolocal) return NULL;
1102 454 : d = ZX_disc(a);
1103 454 : if (!signe(d)) pari_err_IRREDPOL("nfisincl",a);
1104 447 : return d;
1105 : }
1106 :
1107 : static int
1108 49 : badp(GEN fa, GEN db, long q)
1109 : {
1110 49 : GEN P = gel(fa,1), E = gel(fa,2);
1111 49 : long i, l = lg(P);
1112 105 : for (i = 1; i < l; i++)
1113 56 : if (mod2(gel(E,i)) && !dvdii(db, powiu(gel(P,i),q))) return 1;
1114 49 : return 0;
1115 : }
1116 :
1117 : /* is isomorphism / inclusion (a \subset b) compatible with what we know about
1118 : * basic invariants ? (degree, signature, discriminant); test for isomorphism
1119 : * if fliso is set and for inclusion otherwise */
1120 : static int
1121 283 : tests_OK(GEN a, GEN nfa, GEN b, GEN nfb, long fliso)
1122 : {
1123 : GEN da2, da, db, fa, P, U;
1124 283 : long i, l, q, m = degpol(a), n = degpol(b);
1125 :
1126 283 : if (m <= 0) pari_err_IRREDPOL("nfisincl",a);
1127 283 : if (n <= 0) pari_err_IRREDPOL("nfisincl",b);
1128 276 : q = n / m; /* relative degree */
1129 276 : if (fliso) { if (n != m) return 0; } else { if (n % m) return 0; }
1130 276 : if (m == 1) return 1;
1131 :
1132 : /*local test expensive if n^2 >> m^4 <=> q = n/m >> m */
1133 269 : db = incl_disc(nfb, b, q > m);
1134 269 : da = db? incl_disc(nfa, a, 0): NULL;
1135 262 : if (nfa && nfb) /* both nf structures available */
1136 : {
1137 7 : long r1a = nf_get_r1(nfa), r1b = nf_get_r1(nfb);
1138 0 : return fliso ? (r1a == r1b && equalii(da, db))
1139 7 : : (r1b <= r1a * q && dvdii(db, powiu(da, q)));
1140 : }
1141 255 : if (!db) return 1;
1142 248 : if (fliso) return issquare(gdiv(da,db));
1143 :
1144 203 : if (nfa)
1145 : {
1146 7 : P = nf_get_ramified_primes(nfa); l = lg(P);
1147 14 : for (i = 1; i < l; i++)
1148 7 : if (Z_pval(db, gel(P,i)) < q * Z_pval(da, gel(P,i))) return 0;
1149 7 : return 1;
1150 : }
1151 196 : else if (nfb)
1152 : {
1153 28 : P = nf_get_ramified_primes(nfb); l = lg(P);
1154 56 : for (i = 1; i < l; i++)
1155 : {
1156 28 : GEN p = gel(P,i);
1157 28 : long va = Z_pval(nfdisc(mkvec2(a, mkvec(p))), p);
1158 28 : if (va && Z_pval(db, gel(P,i)) < va * q) return 0;
1159 : }
1160 28 : return 1;
1161 : }
1162 : /* da = dK A^2, db = dL B^2, dL = dK^q * N(D)
1163 : * da = da1 * da2, da2 maximal s.t. (da2, db) = 1: let p a prime divisor of
1164 : * da2 then p \nmid da1 * dK and p | A => v_p(da) = v_p(da2) is even */
1165 168 : da2 = Z_ppo(da, db);
1166 168 : if (!is_pm1(da2))
1167 : { /* replace da by da1 all of whose prime divisors divide db */
1168 126 : da2 = absi_shallow(da2);
1169 126 : if (!Z_issquare(da2)) return 0;
1170 14 : da = diviiexact(da, da2);
1171 : }
1172 56 : if (is_pm1(da)) return 1;
1173 49 : fa = absZ_factor_limit_strict(da, 0, &U);
1174 49 : if (badp(fa, db, q)) return 0;
1175 49 : if (U && mod2(gel(U,2)) && expi(gel(U,1)) < 150)
1176 : { /* cofactor is small, finish */
1177 0 : fa = absZ_factor(gel(U,1));
1178 0 : if (badp(fa, db, q)) return 0;
1179 : }
1180 49 : return 1;
1181 : }
1182 :
1183 : GEN
1184 77 : nfisisom(GEN a, GEN b)
1185 : {
1186 77 : pari_sp av = avma;
1187 : long i, va, vb, lx;
1188 : GEN nfa, nfb, y, la, lb;
1189 77 : int newvar, sw = 0;
1190 :
1191 77 : a = get_nfpol(a, &nfa);
1192 77 : b = get_nfpol(b, &nfb);
1193 77 : if (!nfa) { a = Q_primpart(a); RgX_check_ZX(a, "nfisisom"); }
1194 77 : if (!nfb) { b = Q_primpart(b); RgX_check_ZX(b, "nfisisom"); }
1195 77 : if (ZX_equal(a, b))
1196 : {
1197 25 : y = galoisconj(nfb? nfb: b, NULL); settyp(y, t_VEC);
1198 25 : return gerepilecopy(av,y);
1199 : }
1200 52 : if (nfa && !nfb) { swap(a,b); nfb = nfa; nfa = NULL; sw = 1; }
1201 52 : if (!tests_OK(a, nfa, b, nfb, 1)) { set_avma(av); return gen_0; }
1202 :
1203 45 : if (nfb) lb = gen_1; else nfb = b = ZX_Q_normalize(b,&lb);
1204 45 : if (nfa) la = gen_1; else nfa = a = ZX_Q_normalize(a,&la);
1205 45 : va = varn(a); vb = varn(b); newvar = (varncmp(vb,va) <= 0);
1206 45 : if (newvar) { a = leafcopy(a); setvarn(a, fetch_var_higher()); }
1207 45 : y = lift_shallow(nfroots(nfb,a));
1208 45 : if (newvar) (void)delete_var();
1209 45 : lx = lg(y); if (lx==1) { set_avma(av); return gen_0; }
1210 45 : if (sw) { vb = va; b = leafcopy(b); setvarn(b, vb); }
1211 272 : for (i=1; i<lx; i++)
1212 : {
1213 227 : GEN t = gel(y,i);
1214 227 : if (typ(t) == t_POL) setvarn(t, vb); else t = scalarpol(t, vb);
1215 227 : if (lb != gen_1) t = RgX_unscale(t, lb);
1216 227 : if (la != gen_1) t = RgX_Rg_div(t, la);
1217 227 : gel(y,i) = sw? RgXQ_reverse(t, b): t;
1218 : }
1219 45 : settyp(y, t_VEC); return gerepilecopy(av,y);
1220 : }
1221 :
1222 : static GEN
1223 7314 : partmap_reverse(GEN a, GEN b, GEN t, GEN la, GEN lb, long v)
1224 : {
1225 7314 : pari_sp av = avma;
1226 7314 : GEN rnf = rnfequation2(a, t), z;
1227 7315 : if (!RgX_equal(b, gel(rnf,1)))
1228 7 : { setvarn(b,v); pari_err_IRREDPOL("nfisincl", b); }
1229 7308 : z = liftpol_shallow(gel(rnf, 2));
1230 7308 : setvarn(z, v);
1231 7308 : if (!isint1(lb)) z = RgX_unscale(z, lb);
1232 7308 : if (!isint1(la)) z = RgX_Rg_div(z, la);
1233 7308 : return gerepilecopy(av, z);
1234 : }
1235 :
1236 : static GEN
1237 8085 : partmap_reverse_frac(GEN a, GEN b, GEN t, GEN la, GEN lb, long v)
1238 : {
1239 8085 : pari_sp av = avma;
1240 8085 : long k = 0;
1241 : GEN N, D, G, L, de;
1242 8085 : GEN C = ZX_ZXY_resultant_all(a, Q_remove_denom(t,&de), &k, &L);
1243 8085 : if (k || degpol(b) != degpol(C))
1244 0 : { setvarn(b,v); pari_err_IRREDPOL("nfisincl", b); }
1245 8085 : L = Q_primpart(L);
1246 8085 : N = gel(L,1); if (!signe(N)) { set_avma(av); return pol_0(v); }
1247 8078 : D = gel(L,2);
1248 8078 : N = RgX_neg(N); setvarn(N, v); setvarn(D, v);
1249 8078 : G = QX_gcd(N,D);
1250 8078 : if (degpol(G)) { N = RgX_div(N,G); D = RgX_div(D,G); }
1251 8078 : if (!isint1(lb)) { N = RgX_unscale(N, lb); D = RgX_unscale(D, lb); }
1252 8078 : if (!isint1(la)) D = RgX_Rg_mul(D, la);
1253 8078 : return gerepilecopy(av, mkrfrac(N,D));
1254 : }
1255 :
1256 : GEN
1257 8085 : partmap_reverse_frac_worker(GEN t, GEN a, GEN b, GEN la, GEN lb, long v)
1258 8085 : { return partmap_reverse_frac(a, b, t, la, lb, v); }
1259 :
1260 : static GEN
1261 7195 : nfisincl_from_fact(GEN a, long da, GEN b, GEN la, GEN lb, long vb, GEN y,
1262 : long flag)
1263 : {
1264 7195 : long i, k, l = lg(y), db = degpol(b), d = db / da;
1265 7195 : GEN x = cgetg(l, t_VEC);
1266 7594 : for (i= k = 1; i < l; i++)
1267 : {
1268 7412 : GEN t = gel(y,i);
1269 7412 : if (degpol(t) != d) continue;
1270 7315 : gel(x, k++) = partmap_reverse(a, b, t, la, lb, vb);
1271 7308 : if (flag) return gel(x,1);
1272 : }
1273 182 : if (k==1) return gen_0;
1274 91 : setlg(x, k);
1275 91 : gen_sort_inplace(x, (void*)&cmp_RgX, &cmp_nodata, NULL);
1276 91 : return x;
1277 : }
1278 :
1279 : static GEN
1280 1232 : nfisincl_from_fact_frac(GEN a, GEN b, GEN la, GEN lb, long vb, GEN y)
1281 : {
1282 1232 : long i, k, l = lg(y), d = degpol(b) / degpol(a);
1283 1232 : GEN worker, x = cgetg(l, t_VEC);
1284 9317 : for (i = k = 1; i < l; i++)
1285 : {
1286 8085 : GEN t = gel(y,i);
1287 8085 : if (degpol(t) != d) continue;
1288 8085 : gel(x, k++) = t;
1289 : }
1290 1232 : if (k==1) return gen_0;
1291 1232 : worker = snm_closure(is_entry("_partmap_reverse_frac_worker"),
1292 : mkvec5(a,b,la,lb,stoi(vb)));
1293 1232 : setlg(x, k); return gen_parapply(worker, x);
1294 : }
1295 :
1296 : GEN
1297 7308 : nfisincl0(GEN fa, GEN fb, long flag)
1298 : {
1299 7308 : pari_sp av = avma;
1300 : long vb;
1301 : GEN a, b, nfa, nfb, x, y, la, lb;
1302 : int newvar;
1303 7308 : if (flag < 0 || flag > 3) pari_err_FLAG("nfisincl");
1304 :
1305 7308 : a = get_nfpol(fa, &nfa);
1306 7308 : b = get_nfpol(fb, &nfb);
1307 7308 : if (!nfa) { a = Q_primpart(a); RgX_check_ZX(a, "nsisincl"); }
1308 7307 : if (!nfb) { b = Q_primpart(b); RgX_check_ZX(b, "nsisincl"); }
1309 7308 : if (ZX_equal(a, b) && flag<=1)
1310 : {
1311 21 : if (flag==1)
1312 : {
1313 7 : x = pol_x(varn(b));
1314 7 : return degpol(b) > 1 ? x: RgX_rem(x,b);
1315 : }
1316 14 : x = galoisconj(fb, NULL); settyp(x, t_VEC);
1317 14 : return gerepilecopy(av,x);
1318 : }
1319 7287 : if (flag==0 && !tests_OK(a, nfa, b, nfb, 0)) { set_avma(av); return gen_0; }
1320 :
1321 7168 : if (nfb) lb = gen_1; else nfb = b = ZX_Q_normalize(b,&lb);
1322 7168 : if (nfa) la = gen_1; else nfa = a = ZX_Q_normalize(a,&la);
1323 7168 : vb = varn(b); newvar = (varncmp(varn(a),vb) <= 0);
1324 7168 : if (newvar) { b = leafcopy(b); setvarn(b, fetch_var_higher()); }
1325 7168 : y = lift_shallow(gel(nffactor(nfa,b),1));
1326 7168 : if (flag==2)
1327 0 : x = nfisincl_from_fact_frac(a, b, la, lb, vb, y);
1328 : else
1329 7168 : x = nfisincl_from_fact(nfa, degpol(a), b, la, lb, vb, y, flag);
1330 7161 : if (newvar) (void)delete_var();
1331 7161 : return gerepilecopy(av,x);
1332 : }
1333 :
1334 : GEN
1335 112 : nfisincl(GEN fa, GEN fb) { return nfisincl0(fa, fb, 0); }
1336 :
1337 : static GEN
1338 8092 : RgF_to_Flxq(GEN F, GEN T, ulong p)
1339 : {
1340 : GEN N, D, iD;
1341 8092 : if (typ(F)==t_POL) return RgX_to_Flx(F, p);
1342 8085 : N = RgX_to_Flx(gel(F,1), p); D = RgX_to_Flx(gel(F,2), p);
1343 8085 : iD = Flxq_invsafe(D, T, p);
1344 8085 : if (!iD) return NULL;
1345 8078 : return Flxq_mul(N, iD, T, p);
1346 : }
1347 :
1348 : #define pari_APPLY_abort(EXPR)\
1349 : { \
1350 : long i, _l; \
1351 : GEN _y = cgetg_copy(x, &_l);\
1352 : for (i=1; i<_l; i++) \
1353 : { GEN _z = EXPR;\
1354 : if (!_z) return _z;\
1355 : gel(_y,i) = _z;\
1356 : } return _y;\
1357 : }
1358 :
1359 : static GEN
1360 1239 : RgFV_to_FlxqV(GEN x, GEN T, ulong p)
1361 9324 : { pari_APPLY_abort(RgF_to_Flxq(gel(x,i), T, p)) }
1362 :
1363 : static GEN
1364 1232 : nfsplitting_auto(GEN g, GEN R)
1365 : {
1366 : pari_sp av;
1367 : forprime_t T;
1368 1232 : long i, d = degpol(g);
1369 : ulong p;
1370 1232 : GEN P, K, N, G, q, den = Q_denom(R), Rp, Gp;
1371 1232 : u_forprime_init(&T, d*maxss(expu(d)-3, 2), ULONG_MAX);
1372 1232 : av = avma;
1373 79408 : for(;; set_avma(av))
1374 : {
1375 80640 : p = u_forprime_next(&T);
1376 80640 : if (dvdiu(den,p)) continue;
1377 80640 : Gp = ZX_to_Flx(g, p);
1378 80640 : if (!Flx_is_totally_split(Gp, p)) continue;
1379 1239 : P = Flx_roots(Gp, p);
1380 1239 : Rp = RgFV_to_FlxqV(R, Gp, p);
1381 1239 : if (Rp) break;
1382 7 : if (DEBUGLEVEL) err_printf("nfsplitting_auto: bad p : %lu\n",p);
1383 : }
1384 1232 : if (d == 1) return mkvec3(g, mkcol(gel(Rp,1)), utoi(p));
1385 1225 : K = Flm_Flc_invimage(FlxV_to_Flm(Rp, d), vecsmall_ei(d, 2), p);
1386 1225 : if (!K) pari_err_BUG("nfsplitting_auto");
1387 1225 : N = Flm_transpose(FlxV_Flv_multieval(Rp, P, p));
1388 1225 : q = perm_inv(vecsmall_indexsort(gel(N,1)));
1389 1225 : G = cgetg(d+1, t_COL);
1390 35707 : for (i=1; i<=d; i++)
1391 : {
1392 34482 : GEN r = perm_mul(vecsmall_indexsort(gel(N,i)), q);
1393 34482 : gel(G,i) = FlxV_Flc_mul(Rp, vecpermute(K, r), p);
1394 : }
1395 1225 : return mkvec3(g, G, utoi(p));
1396 : }
1397 :
1398 : static GEN
1399 1421 : nfsplitting_composite(GEN P)
1400 : {
1401 1421 : GEN F = gel(ZX_factor(P), 1), Q = NULL;
1402 1421 : long i, n = lg(F)-1;
1403 2849 : for (i = 1; i <= n; i++)
1404 : {
1405 1428 : GEN Fi = gel(F, i);
1406 1428 : if (degpol(Fi) == 1) continue;
1407 1400 : Q = Q ? veclast(compositum(Q, Fi)): Fi;
1408 : }
1409 1421 : return Q ? Q: pol_x(varn(P));
1410 : }
1411 : GEN
1412 1435 : nfsplitting0(GEN T0, GEN D, long flag)
1413 : {
1414 1435 : pari_sp av = avma;
1415 : long d, Ds, v;
1416 1435 : GEN T, F, K, N = NULL, lT = NULL;
1417 1435 : if (flag < 0 || flag > 3) pari_err_FLAG("nfsplitting");
1418 1435 : T = T0 = get_nfpol(T0, &K);
1419 1428 : if (!K)
1420 : {
1421 : GEN c;
1422 1407 : if (typ(T) != t_POL) pari_err_TYPE("nfsplitting",T);
1423 1407 : T = Q_primitive_part(T, &c);
1424 1407 : lT = leading_coeff(T); if (isint1(lT)) lT = NULL;
1425 1407 : if (flag && (c || lT)) pari_err_TYPE("nfsplitting", T0);
1426 1400 : RgX_check_ZX(T,"nfsplitting");
1427 : }
1428 1421 : T = nfsplitting_composite(T);
1429 1421 : if (flag && !ZX_equal(T, T0)) pari_err_IRREDPOL("nfsplitting", T0);
1430 1407 : d = degpol(T); v = varn(T);
1431 1407 : if (d <= 1 && !flag) { set_avma(av); return pol_x(v); }
1432 1379 : if (!K) {
1433 1358 : if (lT) T = polredbest(T,0);
1434 1358 : K = T;
1435 : }
1436 1379 : if (D)
1437 1239 : { if (typ(D) != t_INT || signe(D) < 1) pari_err_TYPE("nfsplitting",D); }
1438 140 : else if (d <= 7 ||
1439 35 : (d <= 11 && pari_is_dir(stack_strcat(pari_datadir, "/galdata"))))
1440 126 : D = gel(polgalois(T,DEFAULTPREC), 1);
1441 : else
1442 14 : D = mpfact(d);
1443 1379 : Ds = itos_or_0(D);
1444 1379 : T = leafcopy(T); setvarn(T, fetch_var_higher());
1445 1379 : for(F = T;;)
1446 1659 : {
1447 3038 : GEN P = gel(nffactor(K, F), 1), Q = veclast(P);
1448 3038 : if (degpol(gel(P,1)) == degpol(Q))
1449 : {
1450 1288 : if (!flag) break;
1451 1260 : P = liftall_shallow(P);
1452 1260 : if (flag==1)
1453 28 : N = nfisincl_from_fact(K, d, F, gen_1, gen_1, v, P, flag);
1454 : else
1455 1232 : N = nfisincl_from_fact_frac(T0, F, gen_1, gen_1, v, P);
1456 1260 : break;
1457 : }
1458 1750 : F = rnfequation(K,Q);
1459 1750 : if (degpol(F) == Ds && !flag) break;
1460 : }
1461 1379 : if (umodiu(D,degpol(F)))
1462 : {
1463 14 : char *sD = itostr(D);
1464 14 : pari_warn(warner,stack_strcat("ignoring incorrect degree bound ",sD));
1465 : }
1466 1379 : setvarn(F, v); (void)delete_var();
1467 1379 : if (flag) F = flag == 3? nfsplitting_auto(F, N): mkvec2(F, N);
1468 1379 : return gerepilecopy(av, F);
1469 : }
1470 :
1471 : GEN
1472 0 : nfsplitting(GEN T, GEN D) { return nfsplitting0(T, D, 0); }
1473 :
1474 : /*************************************************************************/
1475 : /** **/
1476 : /** INITALG **/
1477 : /** **/
1478 : /*************************************************************************/
1479 : typedef struct {
1480 : GEN T;
1481 : GEN ro; /* roots of T */
1482 : long r1;
1483 : GEN basden;
1484 : long prec;
1485 : long extraprec; /* possibly -1 = irrelevant or not computed */
1486 : GEN M, G; /* possibly NULL = irrelevant or not computed */
1487 : } nffp_t;
1488 :
1489 : static GEN
1490 206021 : get_roots(GEN x, long r1, long prec)
1491 : {
1492 : long i, ru;
1493 : GEN z;
1494 206021 : if (typ(x) != t_POL)
1495 : {
1496 0 : z = leafcopy(x);
1497 0 : ru = (lg(z)-1 + r1) >> 1;
1498 : }
1499 : else
1500 : {
1501 206021 : long n = degpol(x);
1502 206021 : z = (r1 == n)? ZX_realroots_irred(x, prec): QX_complex_roots(x,prec);
1503 206023 : ru = (n+r1)>>1;
1504 : }
1505 449193 : for (i=r1+1; i<=ru; i++) gel(z,i) = gel(z, (i<<1)-r1);
1506 206023 : z[0]=evaltyp(t_VEC)|evallg(ru+1); return z;
1507 : }
1508 :
1509 : GEN
1510 0 : nf_get_allroots(GEN nf)
1511 : {
1512 0 : return embed_roots(nf_get_roots(nf), nf_get_r1(nf));
1513 : }
1514 :
1515 : /* For internal use. compute trace(x mod pol), sym=polsym(pol,deg(pol)-1) */
1516 : static GEN
1517 384830 : quicktrace(GEN x, GEN sym)
1518 : {
1519 384830 : GEN p1 = gen_0;
1520 : long i;
1521 :
1522 384830 : if (typ(x) != t_POL) return gmul(x, gel(sym,1));
1523 384830 : if (signe(x))
1524 : {
1525 384830 : sym--;
1526 2894745 : for (i=lg(x)-1; i>1; i--)
1527 2509991 : p1 = gadd(p1, gmul(gel(x,i),gel(sym,i)));
1528 : }
1529 384754 : return p1;
1530 : }
1531 :
1532 : static GEN
1533 83354 : get_Tr(GEN mul, GEN x, GEN basden)
1534 : {
1535 83354 : GEN t, bas = gel(basden,1), den = gel(basden,2);
1536 83354 : long i, j, n = lg(bas)-1;
1537 83354 : GEN T = cgetg(n+1,t_MAT), TW = cgetg(n+1,t_COL), sym = polsym(x, n-1);
1538 :
1539 83354 : gel(TW,1) = utoipos(n);
1540 262290 : for (i=2; i<=n; i++)
1541 : {
1542 178942 : t = quicktrace(gel(bas,i), sym);
1543 178935 : if (den && gel(den,i)) t = diviiexact(t,gel(den,i));
1544 178932 : gel(TW,i) = t; /* tr(w[i]) */
1545 : }
1546 83348 : gel(T,1) = TW;
1547 262290 : for (i=2; i<=n; i++)
1548 : {
1549 178938 : gel(T,i) = cgetg(n+1,t_COL); gcoeff(T,1,i) = gel(TW,i);
1550 706767 : for (j=2; j<=i; j++) /* Tr(W[i]W[j]) */
1551 527825 : gcoeff(T,i,j) = gcoeff(T,j,i) = ZV_dotproduct(gel(mul,j+(i-1)*n), TW);
1552 : }
1553 83352 : return T;
1554 : }
1555 :
1556 : /* return [bas[i]*denom(bas[i]), denom(bas[i])], denom 1 is given as NULL */
1557 : static GEN
1558 202097 : get_bas_den(GEN bas)
1559 : {
1560 202097 : GEN b,d,den, dbas = leafcopy(bas);
1561 202099 : long i, l = lg(bas);
1562 202099 : int power = 1;
1563 202099 : den = cgetg(l,t_VEC);
1564 937514 : for (i=1; i<l; i++)
1565 : {
1566 735415 : b = Q_remove_denom(gel(bas,i), &d);
1567 735410 : gel(dbas,i) = b;
1568 735410 : gel(den,i) = d; if (d) power = 0;
1569 : }
1570 202099 : if (power) den = NULL; /* power basis */
1571 202099 : return mkvec2(dbas, den);
1572 : }
1573 :
1574 : /* return multiplication table for S->basis */
1575 : static GEN
1576 83353 : nf_multable(nfmaxord_t *S, GEN invbas)
1577 : {
1578 83353 : GEN T = S->T, w = gel(S->basden,1), den = gel(S->basden,2);
1579 83353 : long i,j, n = degpol(T);
1580 83353 : GEN mul = cgetg(n*n+1,t_MAT);
1581 :
1582 : /* i = 1 split for efficiency, assume w[1] = 1 */
1583 345653 : for (j=1; j<=n; j++)
1584 262298 : gel(mul,j) = gel(mul,1+(j-1)*n) = col_ei(n, j);
1585 262296 : for (i=2; i<=n; i++)
1586 706770 : for (j=i; j<=n; j++)
1587 : {
1588 527829 : pari_sp av = avma;
1589 527829 : GEN z = (i == j)? ZXQ_sqr(gel(w,i), T): ZXQ_mul(gel(w,i),gel(w,j), T);
1590 527833 : z = ZM_ZX_mul(invbas, z); /* integral column */
1591 527799 : if (den)
1592 : {
1593 345464 : GEN d = mul_denom(gel(den,i), gel(den,j));
1594 345453 : if (d) z = ZC_Z_divexact(z, d);
1595 : }
1596 527797 : gel(mul,j+(i-1)*n) = gel(mul,i+(j-1)*n) = gerepileupto(av,z);
1597 : }
1598 83351 : return mul;
1599 : }
1600 :
1601 : /* as get_Tr, mul_table not precomputed */
1602 : static GEN
1603 27870 : make_Tr(nfmaxord_t *S)
1604 : {
1605 27870 : GEN T = S->T, w = gel(S->basden,1), den = gel(S->basden,2);
1606 27870 : long i,j, n = degpol(T);
1607 27870 : GEN c, t, d, M = cgetg(n+1,t_MAT), sym = polsym(T, n-1);
1608 :
1609 : /* W[i] = w[i]/den[i]; assume W[1] = 1, case i = 1 split for efficiency */
1610 27870 : c = cgetg(n+1,t_COL); gel(M,1) = c;
1611 27870 : gel(c, 1) = utoipos(n);
1612 83449 : for (j=2; j<=n; j++)
1613 : {
1614 55580 : pari_sp av = avma;
1615 55580 : t = quicktrace(gel(w,j), sym);
1616 55571 : if (den)
1617 : {
1618 35409 : d = gel(den,j);
1619 35409 : if (d) t = diviiexact(t, d);
1620 : }
1621 55572 : gel(c,j) = gerepileuptoint(av, t);
1622 : }
1623 83448 : for (i=2; i<=n; i++)
1624 : {
1625 55578 : c = cgetg(n+1,t_COL); gel(M,i) = c;
1626 205896 : for (j=1; j<i ; j++) gel(c,j) = gcoeff(M,i,j);
1627 205894 : for ( ; j<=n; j++)
1628 : {
1629 150315 : pari_sp av = avma;
1630 150315 : t = (i == j)? ZXQ_sqr(gel(w,i), T): ZXQ_mul(gel(w,i),gel(w,j), T);
1631 150315 : t = quicktrace(t, sym);
1632 150299 : if (den)
1633 : {
1634 118050 : d = mul_denom(gel(den,i),gel(den,j));
1635 118052 : if (d) t = diviiexact(t, d);
1636 : }
1637 150301 : gel(c,j) = gerepileuptoint(av, t); /* Tr (W[i]W[j]) */
1638 : }
1639 : }
1640 27870 : return M;
1641 : }
1642 :
1643 : /* [bas[i]/den[i]]= integer basis. roo = real part of the roots */
1644 : static void
1645 237207 : make_M(nffp_t *F, int trunc)
1646 : {
1647 237207 : GEN bas = gel(F->basden,1), den = gel(F->basden,2), ro = F->ro;
1648 : GEN m, d, M;
1649 237207 : long i, j, l = lg(ro), n = lg(bas);
1650 237207 : M = cgetg(n,t_MAT);
1651 237204 : gel(M,1) = const_col(l-1, gen_1); /* bas[1] = 1 */
1652 839383 : for (j=2; j<n; j++) gel(M,j) = cgetg(l,t_COL);
1653 773001 : for (i=1; i<l; i++)
1654 : {
1655 535811 : GEN r = gel(ro,i), ri;
1656 535811 : ri = (gexpo(r) > 1)? ginv(r): NULL;
1657 2674971 : for (j=2; j<n; j++) gcoeff(M,i,j) = RgX_cxeval(gel(bas,j), r, ri);
1658 : }
1659 237190 : if (den)
1660 524131 : for (j=2; j<n; j++)
1661 : {
1662 406763 : d = gel(den,j); if (!d) continue;
1663 344013 : m = gel(M,j);
1664 1721070 : for (i=1; i<l; i++) gel(m,i) = gdiv(gel(m,i), d);
1665 : }
1666 :
1667 237176 : if (trunc && gprecision(M) > F->prec)
1668 : {
1669 18136 : M = gprec_w(M, F->prec);
1670 18136 : F->ro = gprec_w(ro,F->prec);
1671 : }
1672 237176 : F->M = M;
1673 237176 : }
1674 :
1675 : /* return G real such that G~ * G = T_2 */
1676 : static void
1677 237209 : make_G(nffp_t *F)
1678 : {
1679 237209 : GEN G, M = F->M;
1680 237209 : long i, j, k, r1 = F->r1, l = lg(M);
1681 :
1682 237209 : if (r1 == l-1) { F->G = M; return; }
1683 189934 : G = cgetg(l, t_MAT);
1684 894217 : for (j = 1; j < l; j++)
1685 : {
1686 704341 : GEN g, m = gel(M,j);
1687 704341 : gel(G,j) = g = cgetg(l, t_COL);
1688 1145575 : for (k = i = 1; i <= r1; i++) gel(g,k++) = gel(m,i);
1689 2403570 : for ( ; k < l; i++)
1690 : {
1691 1699289 : GEN r = gel(m,i);
1692 1699289 : if (typ(r) == t_COMPLEX)
1693 : {
1694 1395701 : GEN a = gel(r,1), b = gel(r,2);
1695 1395701 : gel(g,k++) = mpadd(a, b);
1696 1395676 : gel(g,k++) = mpsub(a, b);
1697 : }
1698 : else
1699 : {
1700 303588 : gel(g,k++) = r;
1701 303588 : gel(g,k++) = r;
1702 : }
1703 : }
1704 : }
1705 189876 : F->G = G;
1706 : }
1707 :
1708 : static long
1709 269421 : prec_fix(long prec)
1710 : {
1711 : #ifndef LONG_IS_64BIT
1712 : /* make sure that default accuracy is the same on 32/64bit */
1713 38540 : if (odd(prec)) prec++;
1714 : #endif
1715 269421 : return prec;
1716 : }
1717 : static void
1718 237208 : make_M_G(nffp_t *F, int trunc)
1719 : {
1720 : long n, eBD, prec;
1721 237208 : if (F->extraprec < 0)
1722 : { /* not initialized yet; compute roots so that absolute accuracy
1723 : * of M & G >= prec */
1724 : double er;
1725 237188 : n = degpol(F->T);
1726 237189 : eBD = 1 + gexpo(gel(F->basden,1));
1727 237190 : er = F->ro? (1+gexpo(F->ro)): fujiwara_bound(F->T);
1728 237189 : if (er < 0) er = 0;
1729 237189 : F->extraprec = nbits2extraprec(n*er + eBD + log2(n));
1730 : }
1731 237207 : prec = prec_fix(F->prec + F->extraprec);
1732 237207 : if (!F->ro || gprecision(gel(F->ro,1)) < prec)
1733 206022 : F->ro = get_roots(F->T, F->r1, prec);
1734 :
1735 237208 : make_M(F, trunc);
1736 237208 : make_G(F);
1737 237200 : }
1738 :
1739 : static void
1740 174291 : nffp_init(nffp_t *F, nfmaxord_t *S, long prec)
1741 : {
1742 174291 : F->T = S->T;
1743 174291 : F->r1 = S->r1;
1744 174291 : F->basden = S->basden;
1745 174291 : F->ro = NULL;
1746 174291 : F->extraprec = -1;
1747 174291 : F->prec = prec;
1748 174291 : }
1749 :
1750 : /* let bas a t_VEC of QX giving a Z-basis of O_K. Return the index of the
1751 : * basis. Assume bas[1] = 1 and that the leading coefficient of elements
1752 : * of bas are of the form 1/b for a t_INT b */
1753 : static GEN
1754 1589 : get_nfindex(GEN bas)
1755 : {
1756 1589 : pari_sp av = avma;
1757 1589 : long n = lg(bas)-1, i;
1758 : GEN D, d, mat;
1759 :
1760 : /* assume bas[1] = 1 */
1761 1589 : D = gel(bas,1);
1762 1589 : if (! is_pm1(simplify_shallow(D))) pari_err_TYPE("get_nfindex", D);
1763 1589 : D = gen_1;
1764 8589 : for (i = 2; i <= n; i++)
1765 : { /* after nfbasis, basis is upper triangular! */
1766 7007 : GEN B = gel(bas,i), lc;
1767 7007 : if (degpol(B) != i-1) break;
1768 7000 : lc = gel(B, i+1);
1769 7000 : switch (typ(lc))
1770 : {
1771 2688 : case t_INT: continue;
1772 4312 : case t_FRAC: if (is_pm1(gel(lc,1)) ) {D = mulii(D, gel(lc,2)); continue;}
1773 0 : default: pari_err_TYPE("get_nfindex", B);
1774 : }
1775 : }
1776 1589 : if (i <= n)
1777 : { /* not triangular after all */
1778 7 : bas = vecslice(bas,i,n);
1779 7 : bas = Q_remove_denom(bas, &d);
1780 7 : if (!d) return D;
1781 7 : mat = RgV_to_RgM(bas, n);
1782 7 : mat = rowslice(mat, i,n);
1783 7 : D = mulii(D, diviiexact(powiu(d, n-i+1), absi_shallow(ZM_det(mat))));
1784 : }
1785 1589 : return gerepileuptoint(av, D);
1786 : }
1787 : /* make sure all components of S are initialized */
1788 : static void
1789 166700 : nfmaxord_complete(nfmaxord_t *S)
1790 : {
1791 166700 : if (!S->dT) S->dT = ZX_disc(S->T);
1792 166700 : if (!S->index)
1793 : {
1794 1589 : if (S->dK) /* fast */
1795 0 : S->index = sqrti( diviiexact(S->dT, S->dK) );
1796 : else
1797 1589 : S->index = get_nfindex(S->basis);
1798 : }
1799 166700 : if (!S->dK) S->dK = diviiexact(S->dT, sqri(S->index));
1800 166700 : if (S->r1 < 0) S->r1 = ZX_sturm_irred(S->T);
1801 166703 : if (!S->basden) S->basden = get_bas_den(S->basis);
1802 166705 : }
1803 :
1804 : GEN
1805 83351 : nfmaxord_to_nf(nfmaxord_t *S, GEN ro, long prec)
1806 : {
1807 83351 : GEN nf = cgetg(10,t_VEC);
1808 83351 : GEN T = S->T, Tr, D, w, A, dA, MDI, mat = cgetg(9,t_VEC);
1809 83354 : long n = degpol(T);
1810 : nffp_t F;
1811 83351 : nfmaxord_complete(S);
1812 83351 : nffp_init(&F,S,prec);
1813 83352 : F.ro = ro;
1814 83352 : make_M_G(&F, 0);
1815 :
1816 83352 : gel(nf,1) = S->T;
1817 83352 : gel(nf,2) = mkvec2s(S->r1, (n - S->r1)>>1);
1818 83355 : gel(nf,3) = S->dK;
1819 83355 : gel(nf,4) = S->index;
1820 83355 : gel(nf,5) = mat;
1821 83355 : if (gprecision(gel(F.ro,1)) > prec) F.ro = gprec_wtrunc(F.ro, prec);
1822 83353 : gel(nf,6) = F.ro;
1823 83353 : w = S->basis;
1824 83353 : if (!is_pm1(S->index)) w = Q_remove_denom(w, NULL);
1825 83354 : gel(nf,7) = w;
1826 83354 : gel(nf,8) = ZM_inv(RgV_to_RgM(w,n), NULL);
1827 83355 : gel(nf,9) = nf_multable(S, nf_get_invzk(nf));
1828 83354 : gel(mat,1) = F.M;
1829 83354 : gel(mat,2) = F.G;
1830 :
1831 83354 : Tr = get_Tr(gel(nf,9), T, S->basden);
1832 83353 : gel(mat,6) = A = ZM_inv(Tr, &dA); /* dA T^-1, primitive */
1833 83352 : A = ZM_hnfmodid(A, dA);
1834 : /* CAVEAT: nf is not complete yet, but the fields needed for
1835 : * idealtwoelt, zk_scalar_or_multable and idealinv are present ! */
1836 83350 : MDI = idealtwoelt(nf, A);
1837 83353 : gel(MDI,2) = zk_scalar_or_multable(nf, gel(MDI,2));
1838 83355 : gel(mat,7) = MDI;
1839 83355 : if (is_pm1(S->index))
1840 : { /* principal ideal (T'), whose norm is |dK| */
1841 48512 : D = zk_scalar_or_multable(nf, ZX_deriv(T));
1842 48513 : if (typ(D) == t_MAT) D = ZM_hnfmod(D, absi_shallow(S->dK));
1843 : }
1844 : else
1845 : {
1846 34842 : GEN c = diviiexact(dA, gcoeff(A,1,1));
1847 34839 : D = idealHNF_inv_Z(nf, A); /* (A\cap Z) / A */
1848 34842 : if (!is_pm1(c)) D = ZM_Z_mul(D, c);
1849 : }
1850 83355 : gel(mat,3) = RM_round_maxrank(F.G);
1851 83354 : gel(mat,4) = Tr;
1852 83354 : gel(mat,5) = D;
1853 83354 : w = S->dKP; if (!w) { w = gel(absZ_factor(S->dK), 1); settyp(w, t_VEC); }
1854 83354 : gel(mat,8) = w; return nf;
1855 : }
1856 :
1857 : static GEN
1858 3101 : primes_certify(GEN dK, GEN dKP)
1859 : {
1860 3101 : long i, l = lg(dKP);
1861 3101 : GEN v, w, D = dK;
1862 3101 : v = vectrunc_init(l);
1863 3101 : w = vectrunc_init(l);
1864 9695 : for (i = 1; i < l; i++)
1865 : {
1866 6594 : GEN p = gel(dKP,i);
1867 6594 : vectrunc_append(isprime(p)? w: v, p);
1868 6594 : (void)Z_pvalrem(D, p, &D);
1869 : }
1870 3101 : if (!is_pm1(D))
1871 : {
1872 0 : if (signe(D) < 0) D = negi(D);
1873 0 : vectrunc_append(isprime(D)? w: v, D);
1874 : }
1875 3101 : return mkvec2(v,w);
1876 : }
1877 : GEN
1878 3101 : nfcertify(GEN nf)
1879 : {
1880 3101 : pari_sp av = avma;
1881 : GEN vw;
1882 3101 : nf = checknf(nf);
1883 3101 : vw = primes_certify(nf_get_disc(nf), nf_get_ramified_primes(nf));
1884 3101 : return gerepilecopy(av, gel(vw,1));
1885 : }
1886 :
1887 : /* set *pro to roots of S->T */
1888 : static GEN
1889 72789 : get_red_G(nfmaxord_t *S, GEN *pro)
1890 : {
1891 72789 : pari_sp av = avma;
1892 72789 : GEN G, u, u0 = NULL;
1893 72789 : long prec, n = degpol(S->T);
1894 : nffp_t F;
1895 :
1896 72788 : prec = nbits2prec(n+32);
1897 72788 : nffp_init(&F, S, prec);
1898 : for (;;)
1899 : {
1900 72788 : F.prec = prec; make_M_G(&F, 0); G = F.G;
1901 72786 : if (u0) G = RgM_mul(G, u0);
1902 72786 : if (DEBUGLEVEL)
1903 0 : err_printf("get_red_G: starting LLL, prec = %ld (%ld + %ld)\n",
1904 0 : prec + F.extraprec, prec, F.extraprec);
1905 72786 : if ((u = lllfp(G, 0.99, LLL_KEEP_FIRST)))
1906 : {
1907 72788 : if (lg(u)-1 == n) break;
1908 : /* singular ==> loss of accuracy */
1909 0 : if (u0) u0 = gerepileupto(av, RgM_mul(u0,u));
1910 0 : else u0 = gerepilecopy(av, u);
1911 : }
1912 0 : prec = precdbl(prec) + nbits2extraprec(gexpo(u0));
1913 0 : F.ro = NULL;
1914 0 : if (DEBUGLEVEL) pari_warn(warnprec,"get_red_G", prec);
1915 : }
1916 72788 : if (u0) u = RgM_mul(u0,u);
1917 72787 : *pro = F.ro; return u;
1918 : }
1919 :
1920 : /* Compute an LLL-reduced basis for the integer basis of nf(T).
1921 : * set *pro = roots of x if computed [NULL if not computed] */
1922 : static void
1923 101239 : set_LLL_basis(nfmaxord_t *S, GEN *pro, long flag, double DELTA)
1924 : {
1925 101239 : GEN B = S->basis;
1926 101239 : long N = degpol(S->T);
1927 101239 : if (S->r1 < 0)
1928 : {
1929 17983 : S->r1 = ZX_sturm_irred(S->T);
1930 17983 : if (odd(N - S->r1)) pari_err_IRREDPOL("set_LLL_basis", S->T);
1931 : }
1932 101233 : if (!S->basden) S->basden = get_bas_den(B);
1933 101233 : *pro = NULL; if (flag & nf_NOLLL) return;
1934 100659 : if (S->r1 == N) {
1935 27870 : pari_sp av = avma;
1936 27870 : GEN u = ZM_lll(make_Tr(S), DELTA, LLL_GRAM|LLL_KEEP_FIRST|LLL_IM);
1937 27870 : B = gerepileupto(av, RgV_RgM_mul(B, u));
1938 : }
1939 : else
1940 72789 : B = RgV_RgM_mul(B, get_red_G(S, pro));
1941 100657 : S->basis = B;
1942 100657 : S->basden = get_bas_den(B);
1943 : }
1944 :
1945 : /* = 1 iff |a| > |b| or equality and a > 0 */
1946 : static int
1947 140737 : cmpii_polred(GEN a, GEN b)
1948 : {
1949 140737 : int fl = abscmpii(a, b);
1950 : long sa, sb;
1951 140737 : if (fl) return fl;
1952 117375 : sa = signe(a);
1953 117375 : sb = signe(b);
1954 117375 : if (sa == sb) return 0;
1955 755 : return sa == 1? 1: -1;
1956 : }
1957 : static int
1958 32943 : ZX_cmp(GEN x, GEN y)
1959 32943 : { return gen_cmp_RgX((void*)cmpii_polred, x, y); }
1960 : /* current best: ZX x of discriminant *dx, is ZX y better than x ?
1961 : * (if so update *dx); both x and y are monic */
1962 : static int
1963 54812 : ZX_is_better(GEN y, GEN x, GEN *dx)
1964 : {
1965 : pari_sp av;
1966 : int cmp;
1967 : GEN d;
1968 54812 : if (!*dx) *dx = ZX_disc(x);
1969 54812 : av = avma; d = ZX_disc(y);
1970 54812 : cmp = abscmpii(d, *dx);
1971 54812 : if (cmp < 0) { *dx = d; return 1; }
1972 50173 : return gc_bool(av, cmp? 0: (ZX_cmp(y, x) < 0));
1973 : }
1974 :
1975 : static void polredbest_aux(nfmaxord_t *S, GEN *pro, GEN *px, GEN *pdx, GEN *pa);
1976 : /* Seek a simpler, polynomial pol defining the same number field as
1977 : * x (assumed to be monic at this point) */
1978 : static GEN
1979 273 : nfpolred(nfmaxord_t *S, GEN *pro)
1980 : {
1981 273 : GEN x = S->T, dx, b, rev;
1982 273 : long n = degpol(x), v = varn(x);
1983 :
1984 273 : if (n == 1) {
1985 98 : S->T = pol_x(v);
1986 98 : *pro = NULL;
1987 98 : return scalarpol_shallow(negi(gel(x,2)), v);
1988 : }
1989 175 : polredbest_aux(S, pro, &x, &dx, &b);
1990 175 : if (x == S->T) return NULL; /* no improvement */
1991 140 : if (DEBUGLEVEL>1) err_printf("xbest = %Ps\n",x);
1992 :
1993 : /* update T */
1994 140 : rev = QXQ_reverse(b, S->T);
1995 140 : S->basis = QXV_QXQ_eval(S->basis, rev, x);
1996 140 : S->index = sqrti( diviiexact(dx,S->dK) );
1997 140 : S->basden = get_bas_den(S->basis);
1998 140 : S->dT = dx;
1999 140 : S->T = x;
2000 140 : *pro = NULL; /* reset */
2001 140 : return rev;
2002 : }
2003 :
2004 : /* Either nf type or ZX or [monic ZX, data], where data is either an integral
2005 : * basis (deprecated), or listP data (nfbasis input format) to specify
2006 : * a set of primes at with the basis order must be maximal.
2007 : * 1) nf type (or unrecognized): return t_VEC
2008 : * 2) ZX or [ZX, listP]: return t_POL
2009 : * 3) [ZX, order basis]: return 0 (deprecated)
2010 : * incorrect: return -1 */
2011 : static long
2012 85222 : nf_input_type(GEN x)
2013 : {
2014 85222 : GEN T, V, DKP = NULL;
2015 : long i, d, v;
2016 85222 : switch(typ(x))
2017 : {
2018 79664 : case t_POL: return t_POL;
2019 5558 : case t_VEC:
2020 5558 : switch(lg(x))
2021 : {
2022 1792 : case 4: DKP = gel(x,3);
2023 5530 : case 3: break;
2024 28 : default: return t_VEC; /* nf or incorrect */
2025 : }
2026 5530 : T = gel(x,1); V = gel(x,2);
2027 5530 : if (typ(T) != t_POL) return -1;
2028 5530 : switch(typ(V))
2029 : {
2030 224 : case t_INT: case t_MAT: return t_POL;
2031 5306 : case t_VEC: case t_COL:
2032 5306 : if (RgV_is_ZV(V)) return t_POL;
2033 2457 : break;
2034 0 : default: return -1;
2035 : }
2036 2457 : d = degpol(T); v = varn(T);
2037 2457 : if (d<1 || !RgX_is_ZX(T) || !isint1(gel(T,d+2)) || lg(V)-1!=d) return -1;
2038 16002 : for (i = 1; i <= d; i++)
2039 : { /* check integer basis */
2040 13566 : GEN c = gel(V,i);
2041 13566 : switch(typ(c))
2042 : {
2043 35 : case t_INT: break;
2044 13531 : case t_POL: if (varn(c) == v && RgX_is_QX(c) && degpol(c) < d) break;
2045 : /* fall through */
2046 14 : default: return -1;
2047 : }
2048 : }
2049 2436 : if (DKP && (typ(DKP) != t_VEC || !RgV_is_ZV(DKP))) return -1;
2050 2436 : return 0;
2051 : }
2052 0 : return t_VEC; /* nf or incorrect */
2053 : }
2054 :
2055 : /* cater for obsolete nf_PARTIALFACT flag */
2056 : static void
2057 3906 : nfinit_basic_partial(nfmaxord_t *S, GEN T)
2058 : {
2059 3906 : if (typ(T) == t_POL) { nfmaxord(S, mkvec2(T,utoipos(500000)), 0); }
2060 14 : else nfinit_basic(S, T);
2061 3906 : }
2062 : static void
2063 14147 : nfinit_basic_flag(nfmaxord_t *S, GEN x, long flag)
2064 : {
2065 14147 : if (flag & nf_PARTIALFACT)
2066 35 : nfinit_basic_partial(S, x);
2067 : else
2068 14112 : nfinit_basic(S, x);
2069 14140 : }
2070 :
2071 : /* true nf */
2072 : static GEN
2073 62925 : nf_basden(GEN nf)
2074 : {
2075 62925 : GEN zkD = nf_get_zkprimpart(nf), D = nf_get_zkden(nf);
2076 62925 : D = equali1(D)? NULL: const_vec(lg(zkD)-1, D);
2077 62925 : return mkvec2(zkD, D);
2078 : }
2079 : void
2080 85222 : nfinit_basic(nfmaxord_t *S, GEN T)
2081 : {
2082 85222 : switch (nf_input_type(T))
2083 : {
2084 82733 : case t_POL: nfmaxord(S, T, 0); return;
2085 28 : case t_VEC:
2086 : { /* nf, bnf, bnr */
2087 28 : GEN nf = checknf(T);
2088 28 : S->T = S->T0 = nf_get_pol(nf);
2089 28 : S->basis = nf_get_zk(nf); /* probably useless */
2090 28 : S->basden = nf_basden(nf);
2091 28 : S->index = nf_get_index(nf);
2092 28 : S->dK = nf_get_disc(nf);
2093 28 : S->dKP = nf_get_ramified_primes(nf);
2094 28 : S->dT = mulii(S->dK, sqri(S->index));
2095 28 : S->r1 = nf_get_r1(nf); break;
2096 : }
2097 2436 : case 0: /* monic integral polynomial + integer basis (+ ramified primes)*/
2098 2436 : S->T = S->T0 = gel(T,1);
2099 2436 : S->basis = gel(T,2);
2100 2436 : S->basden = NULL;
2101 2436 : S->index = NULL;
2102 2436 : S->dK = NULL;
2103 2436 : S->dKP = NULL;
2104 2436 : if (lg(T) == 4)
2105 : {
2106 1792 : GEN v = gel(T,3); if (typ(v) == t_COL) v = shallowtrans(v);
2107 1792 : S->dKP = v;
2108 : }
2109 2436 : S->dT = NULL;
2110 2436 : S->r1 = -1; break;
2111 21 : default: /* -1 */
2112 21 : pari_err_TYPE("nfinit_basic", T);
2113 : }
2114 2464 : S->dTP = S->dTE = S->dKE = NULL;
2115 2464 : S->unscale = gen_1;
2116 : }
2117 :
2118 : GEN
2119 83353 : nfinit_complete(nfmaxord_t *S, long flag, long prec)
2120 : {
2121 83353 : GEN nf, unscale = S->unscale, rev = NULL;
2122 :
2123 83353 : if (!ZX_is_irred(S->T)) pari_err_IRREDPOL("nfinit",S->T);
2124 83350 : if (!(flag & nf_RED) && !ZX_is_monic(S->T0))
2125 : {
2126 42 : pari_warn(warner,"nonmonic polynomial. Result of the form [nf,c]");
2127 42 : flag |= nf_RED | nf_ORIG;
2128 : }
2129 83350 : if (!(flag & nf_RED) && !isint1(unscale))
2130 : { /* implies lc(x0) = 1 and L := 1/unscale is integral */
2131 4382 : long d = degpol(S->T0);
2132 4382 : GEN L = ginv(unscale); /* x = L^(-deg(x)) x0(L X) */
2133 4382 : GEN f = powiu(L, (d*(d-1)) >> 1);
2134 4382 : S->T = S->T0; /* restore original user-supplied x0, unscale data */
2135 4382 : S->unscale = gen_1;
2136 4382 : S->dT = gmul(S->dT, sqri(f));
2137 4382 : S->basis = RgXV_unscale(S->basis, unscale);
2138 4382 : S->index = gmul(S->index, f);
2139 : }
2140 83350 : nfmaxord_complete(S); /* more expensive after set_LLL_basis */
2141 83354 : if (flag & nf_RED)
2142 : {
2143 : GEN ro;
2144 : /* lie to polred: more efficient to update *after* modreverse, than to
2145 : * unscale in the polred subsystem */
2146 273 : S->unscale = gen_1;
2147 273 : rev = nfpolred(S, &ro);
2148 273 : nf = nfmaxord_to_nf(S, ro, prec);
2149 273 : S->unscale = unscale; /* restore */
2150 : }
2151 : else
2152 : {
2153 83081 : GEN ro; set_LLL_basis(S, &ro, flag, 0.99);
2154 83078 : nf = nfmaxord_to_nf(S, ro, prec);
2155 : }
2156 83354 : if (flag & nf_ORIG)
2157 : {
2158 77 : if (!rev)
2159 : { /* no improvement */
2160 28 : long v = varn(S->T);
2161 28 : rev = degpol(S->T) == 1? pol_0(v): pol_x(v);
2162 : }
2163 77 : if (!isint1(unscale)) rev = RgX_Rg_div(rev, unscale);
2164 77 : nf = mkvec2(nf, mkpolmod(rev, S->T));
2165 : }
2166 83354 : return nf;
2167 : }
2168 : /* Initialize the number field defined by the polynomial x (in variable v)
2169 : * flag & nf_RED: try a polred first.
2170 : * flag & nf_ORIG: return [nfinit(x), Mod(a,red)], where
2171 : * Mod(a,red) = Mod(v,x) (i.e return the base change). */
2172 : GEN
2173 9943 : nfinit0(GEN x, long flag,long prec)
2174 : {
2175 9943 : const pari_sp av = avma;
2176 : nfmaxord_t S;
2177 9943 : if (flag < 0 || flag > 7) pari_err_FLAG("nfinit");
2178 9943 : if (checkrnf_i(x)) return rnf_build_nfabs(x, prec);
2179 9922 : nfinit_basic(&S, x);
2180 9901 : return gerepilecopy(av, nfinit_complete(&S, flag, prec));
2181 : }
2182 : GEN
2183 182 : nfinitred(GEN x, long prec) { return nfinit0(x, nf_RED, prec); }
2184 : GEN
2185 0 : nfinitred2(GEN x, long prec) { return nfinit0(x, nf_RED|nf_ORIG, prec); }
2186 : GEN
2187 6128 : nfinit(GEN x, long prec) { return nfinit0(x, 0, prec); }
2188 :
2189 : /* assume x a bnr/bnf/nf */
2190 : long
2191 886758 : nf_get_prec(GEN x)
2192 : {
2193 886758 : GEN nf = checknf(x), ro = nf_get_roots(nf);
2194 886770 : return (typ(ro)==t_VEC)? precision(gel(ro,1)): DEFAULTPREC;
2195 : }
2196 :
2197 : /* true nf */
2198 : GEN
2199 62897 : nfnewprec_shallow(GEN nf, long prec)
2200 : {
2201 62897 : GEN m, NF = leafcopy(nf);
2202 : nffp_t F;
2203 :
2204 62897 : F.T = nf_get_pol(nf);
2205 62897 : F.ro = NULL;
2206 62897 : F.r1 = nf_get_r1(nf);
2207 62897 : F.basden = nf_basden(nf);
2208 62897 : F.extraprec = -1;
2209 62897 : F.prec = prec; make_M_G(&F, 0);
2210 62895 : gel(NF,5) = m = leafcopy(gel(NF,5));
2211 62897 : gel(m,1) = F.M;
2212 62897 : gel(m,2) = F.G;
2213 62897 : gel(NF,6) = F.ro; return NF;
2214 : }
2215 :
2216 : GEN
2217 154 : nfnewprec(GEN nf, long prec)
2218 : {
2219 : GEN z;
2220 154 : switch(nftyp(nf))
2221 : {
2222 49 : default: pari_err_TYPE("nfnewprec", nf);
2223 14 : case typ_BNF: z = bnfnewprec(nf,prec); break;
2224 7 : case typ_BNR: z = bnrnewprec(nf,prec); break;
2225 84 : case typ_NF: {
2226 84 : pari_sp av = avma;
2227 84 : z = gerepilecopy(av, nfnewprec_shallow(checknf(nf), prec));
2228 84 : break;
2229 : }
2230 : }
2231 105 : return z;
2232 : }
2233 :
2234 : /********************************************************************/
2235 : /** **/
2236 : /** POLRED **/
2237 : /** **/
2238 : /********************************************************************/
2239 : GEN
2240 0 : embednorm_T2(GEN x, long r1)
2241 : {
2242 0 : pari_sp av = avma;
2243 0 : GEN p = RgV_sumpart(x, r1);
2244 0 : GEN q = RgV_sumpart2(x,r1+1, lg(x)-1);
2245 0 : if (q != gen_0) p = gadd(p, gmul2n(q,1));
2246 0 : return avma == av? gcopy(p): gerepileupto(av, p);
2247 : }
2248 :
2249 : /* simplified version of gnorm for scalar, noncomplex inputs, without GC */
2250 : static GEN
2251 164454 : real_norm(GEN x)
2252 : {
2253 164454 : switch(typ(x))
2254 : {
2255 0 : case t_INT: return sqri(x);
2256 164454 : case t_REAL: return sqrr(x);
2257 0 : case t_FRAC: return sqrfrac(x);
2258 : }
2259 0 : pari_err_TYPE("real_norm", x);
2260 : return NULL;/*LCOV_EXCL_LINE*/
2261 : }
2262 : /* simplified version of gnorm, without GC */
2263 : static GEN
2264 35865192 : complex_norm(GEN x)
2265 : {
2266 35865192 : return typ(x) == t_COMPLEX? cxnorm(x): real_norm(x);
2267 : }
2268 : /* return T2(x), argument r1 needed in case x has components whose type
2269 : * is unexpected, e.g. all of them t_INT for embed(gen_1) */
2270 : GEN
2271 71177 : embed_T2(GEN x, long r1)
2272 : {
2273 71177 : pari_sp av = avma;
2274 71177 : long i, l = lg(x);
2275 71177 : GEN c, s = NULL, t = NULL;
2276 71177 : if (typ(gel(x,1)) == t_INT) return muliu(gel(x,1), 2*(l-1)-r1);
2277 235632 : for (i = 1; i <= r1; i++)
2278 : {
2279 164454 : c = real_norm(gel(x,i));
2280 164457 : s = s? gadd(s, c): c;
2281 : }
2282 203335 : for (; i < l; i++)
2283 : {
2284 132160 : c = complex_norm(gel(x,i));
2285 132160 : t = t? gadd(t, c): c;
2286 : }
2287 71175 : if (t) { t = gmul2n(t,1); s = s? gadd(s,t): t; }
2288 71175 : return gerepileupto(av, s);
2289 : }
2290 : /* return N(x) */
2291 : GEN
2292 21355555 : embed_norm(GEN x, long r1)
2293 : {
2294 21355555 : pari_sp av = avma;
2295 21355555 : long i, l = lg(x);
2296 21355555 : GEN c, s = NULL, t = NULL;
2297 21355555 : if (typ(gel(x,1)) == t_INT) return powiu(gel(x,1), 2*(l-1)-r1);
2298 52779909 : for (i = 1; i <= r1; i++)
2299 : {
2300 31550457 : c = gel(x,i);
2301 31550457 : s = s? gmul(s, c): c;
2302 : }
2303 56961305 : for (; i < l; i++)
2304 : {
2305 35733078 : c = complex_norm(gel(x,i));
2306 35730833 : t = t? gmul(t, c): c;
2307 : }
2308 21228227 : if (t) s = s? gmul(s,t): t;
2309 21228257 : return gerepileupto(av, s);
2310 : }
2311 :
2312 : typedef struct {
2313 : long r1, v, prec;
2314 : GEN ZKembed; /* embeddings of fincke-pohst-reduced Zk basis */
2315 : GEN u; /* matrix giving fincke-pohst-reduced Zk basis */
2316 : GEN M; /* embeddings of initial (LLL-reduced) Zk basis */
2317 : GEN bound; /* T2 norm of the polynomial defining nf */
2318 : long expo_best_disc; /* expo(disc(x)), best generator so far */
2319 : } CG_data;
2320 :
2321 : /* characteristic pol of x (given by embeddings) */
2322 : static GEN
2323 253062 : get_pol(CG_data *d, GEN x)
2324 : {
2325 : long e;
2326 253062 : GEN g = grndtoi(roots_to_pol_r1(x, d->v, d->r1), &e);
2327 253064 : return (e > -5)? NULL: g;
2328 : }
2329 :
2330 : /* characteristic pol of x (given as vector on (w_i)) */
2331 : static GEN
2332 93154 : get_polchar(CG_data *d, GEN x)
2333 93154 : { return get_pol(d, RgM_RgC_mul(d->ZKembed,x)); }
2334 :
2335 : /* Choose a canonical polynomial in the pair { Pmin_a, Pmin_{-a} }, i.e.
2336 : * { z(X), (-1)^(deg z) z(-Z) } and keeping the smallest wrt cmpii_polred
2337 : * Either leave z alone (return 1) or set z <- (-1)^n z(-X). In place. */
2338 : int
2339 96407 : ZX_canon_neg(GEN z)
2340 : {
2341 : long i, s;
2342 202858 : for (i = lg(z)-2; i >= 2; i -= 2)
2343 : { /* examine the odd (resp. even) part of z if deg(z) even (resp. odd). */
2344 195043 : s = signe(gel(z,i));
2345 195043 : if (!s) continue;
2346 : /* non trivial */
2347 88592 : if (s < 0) break; /* z(X) < (-1)^n z(-X) */
2348 :
2349 216044 : for (; i>=2; i-=2) gel(z,i) = negi(gel(z,i));
2350 37514 : return 1;
2351 : }
2352 58893 : return 0;
2353 : }
2354 : /* return a defining polynomial for Q(alpha), v = embeddings of alpha.
2355 : * Return NULL on failure: discriminant too large or non primitive */
2356 : static GEN
2357 136059 : try_polmin(CG_data *d, nfmaxord_t *S, GEN v, long flag, GEN *ai)
2358 : {
2359 136059 : const long best = flag & nf_ABSOLUTE;
2360 : long ed;
2361 136059 : pari_sp av = avma;
2362 : GEN g;
2363 136059 : if (best)
2364 : {
2365 135185 : ed = expo(embed_disc(v, d->r1, LOWDEFAULTPREC));
2366 135184 : set_avma(av); if (d->expo_best_disc < ed) return NULL;
2367 : }
2368 : else
2369 874 : ed = 0;
2370 79621 : g = get_pol(d, v);
2371 : /* accuracy too low, compute algebraically */
2372 79629 : if (!g) { set_avma(av); g = ZXQ_charpoly(*ai, S->T, varn(S->T)); }
2373 79629 : g = ZX_radical(g);
2374 79629 : if (best && degpol(g) != degpol(S->T)) return gc_NULL(av);
2375 33769 : g = gerepilecopy(av, g);
2376 33769 : d->expo_best_disc = ed;
2377 33769 : if (flag & nf_ORIG)
2378 : {
2379 29858 : if (ZX_canon_neg(g)) *ai = RgX_neg(*ai);
2380 29858 : if (!isint1(S->unscale)) *ai = RgX_unscale(*ai, S->unscale);
2381 : }
2382 : else
2383 3911 : (void)ZX_canon_neg(g);
2384 33769 : if (DEBUGLEVEL>3) err_printf("polred: generator %Ps\n", g);
2385 33769 : return g;
2386 : }
2387 :
2388 : /* does x generate the correct field ? */
2389 : static GEN
2390 93154 : chk_gen(void *data, GEN x)
2391 : {
2392 93154 : pari_sp av = avma, av1;
2393 93154 : GEN h, g = get_polchar((CG_data*)data,x);
2394 93152 : if (!g) pari_err_PREC("chk_gen");
2395 93152 : av1 = avma;
2396 93152 : h = ZX_gcd(g, ZX_deriv(g));
2397 93154 : if (degpol(h)) return gc_NULL(av);
2398 59612 : if (DEBUGLEVEL>3) err_printf(" generator: %Ps\n",g);
2399 59612 : set_avma(av1); return gerepileupto(av, g);
2400 : }
2401 :
2402 : static long
2403 32214 : chk_gen_prec(long N, long bit)
2404 32214 : { return prec_fix(nbits2prec(10 + (long)log2((double)N) + bit)); }
2405 :
2406 : /* v = [P,A] two vectors (of ZX and ZV resp.) of same length; remove duplicate
2407 : * polynomials in P, updating A, in place. Among elements having the same
2408 : * characteristic pol, choose the smallest according to ZV_abscmp */
2409 : static void
2410 13853 : remove_duplicates(GEN v)
2411 : {
2412 13853 : GEN x, a, P = gel(v,1), A = gel(v,2);
2413 13853 : long k, i, l = lg(P);
2414 13853 : pari_sp av = avma;
2415 :
2416 13853 : if (l < 2) return;
2417 13853 : (void)sort_factor_pol(mkvec2(P, A), cmpii);
2418 13853 : x = gel(P,1); a = gel(A,1);
2419 58116 : for (k=1,i=2; i<l; i++)
2420 44263 : if (ZX_equal(gel(P,i), x))
2421 : {
2422 21743 : if (ZV_abscmp(gel(A,i), a) < 0) a = gel(A,i);
2423 : }
2424 : else
2425 : {
2426 22520 : gel(A,k) = a;
2427 22520 : gel(P,k) = x;
2428 22520 : k++;
2429 22520 : x = gel(P,i); a = gel(A,i);
2430 : }
2431 13853 : l = k+1;
2432 13853 : gel(A,k) = a; setlg(A,l);
2433 13853 : gel(P,k) = x; setlg(P,l); set_avma(av);
2434 : }
2435 :
2436 : static void
2437 18158 : polred_init(nfmaxord_t *S, nffp_t *F, CG_data *d)
2438 : {
2439 18158 : long e, prec, n = degpol(S->T);
2440 : double log2rho;
2441 : GEN ro;
2442 18158 : set_LLL_basis(S, &ro, 0, 0.9999);
2443 : /* || polchar ||_oo < 2^e ~ 2 (n * rho)^n, rho = max modulus of root */
2444 18151 : log2rho = ro ? (double)gexpo(ro): fujiwara_bound(S->T);
2445 18151 : e = n * (long)(log2rho + log2((double)n)) + 1;
2446 18151 : if (e < 0) e = 0; /* can occur if n = 1 */
2447 18151 : prec = chk_gen_prec(n, e);
2448 18151 : nffp_init(F,S,prec);
2449 18151 : F->ro = ro;
2450 18151 : make_M_G(F, 1);
2451 :
2452 18151 : d->v = varn(S->T);
2453 18151 : d->expo_best_disc = -1;
2454 18151 : d->ZKembed = NULL;
2455 18151 : d->M = NULL;
2456 18151 : d->u = NULL;
2457 18151 : d->r1= S->r1;
2458 18151 : }
2459 : static GEN
2460 14070 : findmindisc(GEN y)
2461 : {
2462 14070 : GEN x = gel(y,1), dx = NULL;
2463 14070 : long i, l = lg(y);
2464 35281 : for (i = 2; i < l; i++)
2465 : {
2466 21211 : GEN yi = gel(y,i);
2467 21211 : if (ZX_is_better(yi,x,&dx)) x = yi;
2468 : }
2469 14070 : return x;
2470 : }
2471 : /* filter [y,b] from polred_aux: keep a single polynomial of degree n in y
2472 : * [ the best wrt discriminant ordering ], but keep all imprimitive
2473 : * polynomials */
2474 : static void
2475 4088 : filter(GEN y, GEN b, long n)
2476 : {
2477 : GEN x, a, dx;
2478 4088 : long i, k = 1, l = lg(y);
2479 4088 : a = x = dx = NULL;
2480 37920 : for (i = 1; i < l; i++)
2481 : {
2482 33832 : GEN yi = gel(y,i), ai = gel(b,i);
2483 33832 : if (degpol(yi) == n)
2484 : {
2485 33664 : if (!dx) dx = ZX_disc(yi); else if (!ZX_is_better(yi,x,&dx)) continue;
2486 5874 : x = yi; a = ai; continue;
2487 : }
2488 168 : gel(y,k) = yi;
2489 168 : gel(b,k) = ai; k++;
2490 : }
2491 4088 : if (dx)
2492 : {
2493 3915 : gel(y,k) = x;
2494 3915 : gel(b,k) = a; k++;
2495 : }
2496 4088 : setlg(y, k);
2497 4088 : setlg(b, k);
2498 4088 : }
2499 :
2500 : static GEN
2501 4123 : polred_aux(nfmaxord_t *S, GEN *pro, long flag)
2502 : { /* only keep polynomials of max degree and best discriminant */
2503 4123 : const long best = flag & nf_ABSOLUTE;
2504 4123 : const long orig = flag & nf_ORIG;
2505 4123 : GEN M, b, y, x = S->T;
2506 4123 : long maxi, i, j, k, v = varn(x), n = lg(S->basis)-1;
2507 : nffp_t F;
2508 : CG_data d;
2509 :
2510 4123 : if (n == 1)
2511 : {
2512 28 : if (!best)
2513 : {
2514 14 : GEN X = pol_x(v);
2515 14 : return orig? mkmat2(mkcol(X),mkcol(gen_1)): mkvec(X);
2516 : }
2517 : else
2518 14 : return orig? trivial_fact(): cgetg(1,t_VEC);
2519 : }
2520 :
2521 4095 : polred_init(S, &F, &d);
2522 4088 : if (pro) *pro = F.ro;
2523 4088 : M = F.M;
2524 4088 : if (best)
2525 : {
2526 4025 : if (!S->dT) S->dT = ZX_disc(S->T);
2527 4025 : d.expo_best_disc = expi(S->dT);
2528 : }
2529 :
2530 : /* n + 2 sum_{1 <= i <= n} n-i = n + n(n-1) = n*n */
2531 4088 : y = cgetg(n*n + 1, t_VEC);
2532 4088 : b = cgetg(n*n + 1, t_COL);
2533 4088 : k = 1;
2534 4088 : if (!best) { gel(y,1) = pol_x(v); gel(b,1) = gen_0; k++; }
2535 26922 : for (i = 2; i <= n; i++)
2536 : {
2537 : GEN ch, ai;
2538 22834 : ai = gel(S->basis,i);
2539 22834 : ch = try_polmin(&d, S, gel(M,i), flag, &ai);
2540 22834 : if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
2541 : }
2542 4088 : maxi = minss(n, 3);
2543 15974 : for (i = 1; i <= maxi; i++)
2544 68502 : for (j = i+1; j <= n; j++)
2545 : {
2546 : GEN ch, ai, v;
2547 56616 : ai = gadd(gel(S->basis,i), gel(S->basis,j));
2548 56616 : v = RgV_add(gel(M,i), gel(M,j));
2549 : /* defining polynomial for Q(w_i+w_j) */
2550 56615 : ch = try_polmin(&d, S, v, flag, &ai);
2551 56616 : if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
2552 :
2553 56616 : ai = gsub(gel(S->basis,i), gel(S->basis,j));
2554 56616 : v = RgV_sub(gel(M,i), gel(M,j));
2555 : /* defining polynomial for Q(w_i-w_j) */
2556 56615 : ch = try_polmin(&d, S, v, flag, &ai);
2557 56616 : if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
2558 : }
2559 4088 : setlg(y, k);
2560 4088 : setlg(b, k); filter(y, b, n);
2561 4088 : if (!orig) return gen_sort_uniq(y, (void*)cmpii, &gen_cmp_RgX);
2562 3409 : settyp(y, t_COL);
2563 3409 : (void)sort_factor_pol(mkmat2(y, b), cmpii);
2564 3409 : return mkmat2(b, y);
2565 : }
2566 :
2567 : /* FIXME: obsolete */
2568 : static GEN
2569 84 : Polred(GEN x, long flag, GEN fa)
2570 : {
2571 84 : pari_sp av = avma;
2572 : nfmaxord_t S;
2573 84 : if (fa)
2574 14 : nfinit_basic(&S, mkvec2(x,fa));
2575 : else
2576 70 : nfinit_basic_flag(&S, x, flag);
2577 77 : return gerepilecopy(av, polred_aux(&S, NULL, flag));
2578 : }
2579 :
2580 : /* finds "best" polynomial in polred_aux list, defaulting to S->T if none of
2581 : * them is primitive. *px a ZX, characteristic polynomial of Mod(*pb,S->T),
2582 : * *pdx its discriminant if pdx != NULL. Set *pro = polroots(S->T) */
2583 : static void
2584 4046 : polredbest_aux(nfmaxord_t *S, GEN *pro, GEN *px, GEN *pdx, GEN *pb)
2585 : {
2586 4046 : GEN y, dx, x = S->T; /* default value */
2587 : long i, l;
2588 4046 : y = polred_aux(S, pro, pb? nf_ORIG|nf_ABSOLUTE: nf_ABSOLUTE);
2589 4039 : dx = S->dT;
2590 4039 : if (pb)
2591 : {
2592 3395 : GEN a, b = deg1pol_shallow(S->unscale, gen_0, varn(x));
2593 3395 : a = gel(y,1); l = lg(a);
2594 3395 : y = gel(y,2);
2595 6610 : for (i=1; i<l; i++)
2596 : {
2597 3215 : GEN yi = gel(y,i);
2598 3215 : if (ZX_is_better(yi,x,&dx)) { x = yi; b = gel(a,i); }
2599 : }
2600 3395 : *pb = b;
2601 : }
2602 : else
2603 : {
2604 644 : l = lg(y);
2605 1281 : for (i=1; i<l; i++)
2606 : {
2607 637 : GEN yi = gel(y,i);
2608 637 : if (ZX_is_better(yi,x,&dx)) x = yi;
2609 : }
2610 : }
2611 4039 : if (pdx) { if (!dx) dx = ZX_disc(x); *pdx = dx; }
2612 4039 : *px = x;
2613 4039 : }
2614 : static GEN
2615 3871 : polredbest_i(GEN T, long flag)
2616 : {
2617 : nfmaxord_t S;
2618 : GEN a;
2619 3871 : nfinit_basic_partial(&S, T);
2620 3871 : polredbest_aux(&S, NULL, &T, NULL, flag? &a: NULL);
2621 3864 : if (flag == 2)
2622 350 : T = mkvec2(T, a);
2623 3514 : else if (flag == 1)
2624 : {
2625 2870 : GEN b = (S.T0 == T)? pol_x(varn(T)): QXQ_reverse(a, S.T0);
2626 : /* charpoly(Mod(a,T0)) = T; charpoly(Mod(b,T)) = S.x */
2627 2870 : if (degpol(T) == 1) b = grem(b,T);
2628 2870 : T = mkvec2(T, mkpolmod(b,T));
2629 : }
2630 3864 : return T;
2631 : }
2632 : GEN
2633 3514 : polredbest(GEN T, long flag)
2634 : {
2635 3514 : pari_sp av = avma;
2636 3514 : if (flag < 0 || flag > 1) pari_err_FLAG("polredbest");
2637 3514 : return gerepilecopy(av, polredbest_i(T, flag));
2638 : }
2639 : /* DEPRECATED: backward compatibility */
2640 : GEN
2641 70 : polred0(GEN x, long flag, GEN fa)
2642 : {
2643 70 : long fl = 0;
2644 70 : if (flag & 1) fl |= nf_PARTIALFACT;
2645 70 : if (flag & 2) fl |= nf_ORIG;
2646 70 : return Polred(x, fl, fa);
2647 : }
2648 :
2649 : GEN
2650 21 : polredord(GEN x)
2651 : {
2652 21 : pari_sp av = avma;
2653 : GEN v, lt;
2654 : long i, n, vx;
2655 :
2656 21 : if (typ(x) != t_POL) pari_err_TYPE("polredord",x);
2657 21 : x = Q_primpart(x); RgX_check_ZX(x,"polredord");
2658 21 : n = degpol(x); if (n <= 0) pari_err_CONSTPOL("polredord");
2659 21 : if (n == 1) return gerepilecopy(av, mkvec(x));
2660 14 : lt = leading_coeff(x); vx = varn(x);
2661 14 : if (is_pm1(lt))
2662 : {
2663 7 : if (signe(lt) < 0) x = ZX_neg(x);
2664 7 : v = pol_x_powers(n, vx);
2665 : }
2666 : else
2667 : { GEN L;
2668 : /* basis for Dedekind order */
2669 7 : v = cgetg(n+1, t_VEC);
2670 7 : gel(v,1) = scalarpol_shallow(lt, vx);
2671 14 : for (i = 2; i <= n; i++)
2672 7 : gel(v,i) = RgX_Rg_add(RgX_mulXn(gel(v,i-1), 1), gel(x,n+3-i));
2673 7 : gel(v,1) = pol_1(vx);
2674 7 : x = ZX_Q_normalize(x, &L);
2675 7 : v = gsubst(v, vx, monomial(ginv(L),1,vx));
2676 14 : for (i=2; i <= n; i++)
2677 7 : if (Q_denom(gel(v,i)) == gen_1) gel(v,i) = pol_xn(i-1, vx);
2678 : }
2679 14 : return gerepileupto(av, polred(mkvec2(x, v)));
2680 : }
2681 :
2682 : GEN
2683 14 : polred(GEN x) { return Polred(x, 0, NULL); }
2684 : GEN
2685 0 : smallpolred(GEN x) { return Polred(x, nf_PARTIALFACT, NULL); }
2686 : GEN
2687 0 : factoredpolred(GEN x, GEN fa) { return Polred(x, 0, fa); }
2688 : GEN
2689 0 : polred2(GEN x) { return Polred(x, nf_ORIG, NULL); }
2690 : GEN
2691 0 : smallpolred2(GEN x) { return Polred(x, nf_PARTIALFACT|nf_ORIG, NULL); }
2692 : GEN
2693 0 : factoredpolred2(GEN x, GEN fa) { return Polred(x, nf_PARTIALFACT, fa); }
2694 :
2695 : /********************************************************************/
2696 : /** **/
2697 : /** POLREDABS **/
2698 : /** **/
2699 : /********************************************************************/
2700 : /* set V[k] := matrix of multiplication by nk.zk[k] */
2701 : static GEN
2702 22458 : set_mulid(GEN V, GEN M, GEN Mi, long r1, long r2, long N, long k)
2703 : {
2704 22458 : GEN v, Mk = cgetg(N+1, t_MAT);
2705 : long i, e;
2706 35617 : for (i = 1; i < k; i++) gel(Mk,i) = gmael(V, i, k);
2707 147914 : for ( ; i <=N; i++)
2708 : {
2709 125456 : v = vecmul(gel(M,k), gel(M,i));
2710 125456 : v = RgM_RgC_mul(Mi, split_realimag(v, r1, r2));
2711 125457 : gel(Mk,i) = grndtoi(v, &e);
2712 125457 : if (e > -5) return NULL;
2713 : }
2714 22458 : gel(V,k) = Mk; return Mk;
2715 : }
2716 :
2717 : static GEN
2718 12706 : ZM_image_shallow(GEN M, long *pr)
2719 : {
2720 : long j, k, r;
2721 12706 : GEN y, d = ZM_pivots(M, &k);
2722 12706 : r = lg(M)-1 - k;
2723 12706 : y = cgetg(r+1,t_MAT);
2724 57063 : for (j=k=1; j<=r; k++)
2725 44357 : if (d[k]) gel(y,j++) = gel(M,k);
2726 12706 : *pr = r; return y;
2727 : }
2728 :
2729 : /* U = base change matrix, R = Cholesky form of the quadratic form [matrix
2730 : * Q from algo 2.7.6] */
2731 : static GEN
2732 14071 : chk_gen_init(FP_chk_fun *chk, GEN R, GEN U)
2733 : {
2734 14071 : CG_data *d = (CG_data*)chk->data;
2735 : GEN P, V, D, inv, bound, S, M;
2736 14071 : long N = lg(U)-1, r1 = d->r1, r2 = (N-r1)>>1;
2737 14071 : long i, j, prec, firstprim = 0, skipfirst = 0;
2738 : pari_sp av;
2739 :
2740 14071 : d->u = U;
2741 14071 : d->ZKembed = M = RgM_mul(d->M, U);
2742 :
2743 14071 : av = avma; bound = d->bound;
2744 14071 : D = cgetg(N+1, t_VECSMALL);
2745 94314 : for (i = 1; i <= N; i++)
2746 : {
2747 80251 : pari_sp av2 = avma;
2748 80251 : P = get_pol(d, gel(M,i));
2749 80251 : if (!P) pari_err_PREC("chk_gen_init");
2750 80243 : P = gerepilecopy(av2, ZX_radical(P));
2751 80243 : D[i] = degpol(P);
2752 80243 : if (D[i] == N)
2753 : { /* primitive element */
2754 56834 : GEN B = embed_T2(gel(M,i), r1);
2755 56834 : if (!firstprim) firstprim = i; /* index of first primitive element */
2756 56834 : if (DEBUGLEVEL>2) err_printf("chk_gen_init: generator %Ps\n",P);
2757 56834 : if (gcmp(B,bound) < 0) bound = gerepileuptoleaf(av2, B);
2758 : }
2759 : else
2760 : {
2761 23409 : if (DEBUGLEVEL>2) err_printf("chk_gen_init: subfield %Ps\n",P);
2762 23409 : if (firstprim)
2763 : { /* cycle basis vectors so that primitive elements come last */
2764 2716 : GEN u = d->u, e = M;
2765 2716 : GEN te = gel(e,i), tu = gel(u,i), tR = gel(R,i);
2766 2716 : long tS = D[i];
2767 8364 : for (j = i; j > firstprim; j--)
2768 : {
2769 5648 : u[j] = u[j-1];
2770 5648 : e[j] = e[j-1];
2771 5648 : R[j] = R[j-1];
2772 5648 : D[j] = D[j-1];
2773 : }
2774 2716 : gel(u,firstprim) = tu;
2775 2716 : gel(e,firstprim) = te;
2776 2716 : gel(R,firstprim) = tR;
2777 2716 : D[firstprim] = tS; firstprim++;
2778 : }
2779 : }
2780 : }
2781 14063 : if (!firstprim)
2782 : { /* try (a little) to find primitive elements to improve bound */
2783 28 : GEN x = cgetg(N+1, t_VECSMALL);
2784 28 : if (DEBUGLEVEL>1)
2785 0 : err_printf("chk_gen_init: difficult field, trying random elements\n");
2786 308 : for (i = 0; i < 10; i++)
2787 : {
2788 : GEN e, B;
2789 3010 : for (j = 1; j <= N; j++) x[j] = (long)random_Fl(7) - 3;
2790 280 : e = RgM_zc_mul(M, x);
2791 280 : B = embed_T2(e, r1);
2792 280 : if (gcmp(B,bound) >= 0) continue;
2793 33 : P = get_pol(d, e); if (!P) pari_err_PREC( "chk_gen_init");
2794 33 : if (!ZX_is_squarefree(P)) continue;
2795 33 : if (DEBUGLEVEL>2) err_printf("chk_gen_init: generator %Ps\n",P);
2796 33 : bound = B ;
2797 : }
2798 : }
2799 :
2800 14063 : if (firstprim != 1)
2801 : {
2802 13877 : inv = ginv( split_realimag(M, r1, r2) ); /*TODO: use QR?*/
2803 13877 : V = gel(inv,1);
2804 53874 : for (i = 2; i <= r1+r2; i++) V = gadd(V, gel(inv,i));
2805 : /* V corresponds to 1_Z */
2806 13877 : V = grndtoi(V, &j);
2807 13877 : if (j > -5) pari_err_BUG("precision too low in chk_gen_init");
2808 13877 : S = mkmat(V); /* 1 */
2809 :
2810 13877 : V = cgetg(N+1, t_VEC);
2811 35641 : for (i = 1; i <= N; i++,skipfirst++)
2812 : { /* S = Q-basis of subfield generated by nf.zk[1..i-1] */
2813 : GEN Mx, M2;
2814 35641 : long j, k, h, rkM, dP = D[i];
2815 :
2816 35641 : if (dP == N) break; /* primitive */
2817 22458 : Mx = set_mulid(V, M, inv, r1, r2, N, i);
2818 22458 : if (!Mx) break; /* prec. problem. Stop */
2819 22458 : if (dP == 1) continue;
2820 9302 : rkM = lg(S)-1;
2821 9302 : M2 = cgetg(N+1, t_MAT); /* we will add to S the elts of M2 */
2822 9302 : gel(M2,1) = col_ei(N, i); /* nf.zk[i] */
2823 9302 : k = 2;
2824 18620 : for (h = 1; h < dP; h++)
2825 : {
2826 : long r; /* add to M2 the elts of S * nf.zk[i] */
2827 39490 : for (j = 1; j <= rkM; j++) gel(M2,k++) = ZM_ZC_mul(Mx, gel(S,j));
2828 12706 : setlg(M2, k); k = 1;
2829 12706 : S = ZM_image_shallow(shallowconcat(S,M2), &r);
2830 13400 : if (r == rkM) break;
2831 10012 : if (r > rkM)
2832 : {
2833 10012 : rkM = r;
2834 10012 : if (rkM == N) break;
2835 : }
2836 : }
2837 9302 : if (rkM == N) break;
2838 : /* Q(w[1],...,w[i-1]) is a strict subfield of nf */
2839 : }
2840 : }
2841 : /* x_1,...,x_skipfirst generate a strict subfield [unless N=skipfirst=1] */
2842 14063 : chk->skipfirst = skipfirst;
2843 14063 : if (DEBUGLEVEL>2) err_printf("chk_gen_init: skipfirst = %ld\n",skipfirst);
2844 :
2845 : /* should be DEF + gexpo( max_k C^n_k (bound/k)^(k/2) ) */
2846 14063 : bound = gerepileuptoleaf(av, bound);
2847 14063 : prec = chk_gen_prec(N, (gexpo(bound)*N)/2);
2848 14063 : if (DEBUGLEVEL)
2849 0 : err_printf("chk_gen_init: new prec = %ld (initially %ld)\n", prec, d->prec);
2850 14063 : if (prec > d->prec) pari_err_BUG("polredabs (precision problem)");
2851 14063 : if (prec < d->prec) d->ZKembed = gprec_w(M, prec);
2852 14063 : return bound;
2853 : }
2854 :
2855 : static GEN
2856 14077 : polredabs_i(GEN x, nfmaxord_t *S, GEN *u, long flag)
2857 : {
2858 14077 : FP_chk_fun chk = { &chk_gen, &chk_gen_init, NULL, NULL, 0 };
2859 : nffp_t F;
2860 : CG_data d;
2861 : GEN v, y, a;
2862 : long i, l;
2863 :
2864 14077 : nfinit_basic_flag(S, x, flag);
2865 14077 : x = S->T0;
2866 14077 : if (degpol(x) == 1)
2867 : {
2868 14 : long vx = varn(x);
2869 14 : *u = NULL;
2870 14 : return mkvec2(mkvec( pol_x(vx) ),
2871 14 : mkvec( deg1pol_shallow(gen_1, negi(gel(S->T,2)), vx) ));
2872 : }
2873 14063 : chk.data = (void*)&d;
2874 14063 : polred_init(S, &F, &d);
2875 14063 : d.bound = embed_T2(F.ro, d.r1);
2876 14063 : if (realprec(d.bound) > F.prec) d.bound = rtor(d.bound, F.prec);
2877 : for (;;)
2878 20 : {
2879 14083 : GEN R = R_from_QR(F.G, F.prec);
2880 14083 : if (R)
2881 : {
2882 14071 : d.prec = F.prec;
2883 14071 : d.M = F.M;
2884 14071 : v = fincke_pohst(mkvec(R),NULL,-1, 0, &chk);
2885 14071 : if (v) break;
2886 : }
2887 20 : F.prec = precdbl(F.prec);
2888 20 : F.ro = NULL;
2889 20 : make_M_G(&F, 1);
2890 20 : if (DEBUGLEVEL) pari_warn(warnprec,"polredabs0",F.prec);
2891 : }
2892 14063 : y = gel(v,1);
2893 14063 : a = gel(v,2); l = lg(a);
2894 73201 : for (i = 1; i < l; i++) /* normalize wrt z -> -z */
2895 59138 : if (ZX_canon_neg(gel(y,i)) && (flag & (nf_ORIG|nf_RAW)))
2896 546 : gel(a,i) = ZC_neg(gel(a,i));
2897 14063 : *u = d.u; return v;
2898 : }
2899 :
2900 : GEN
2901 13853 : polredabs0(GEN x, long flag)
2902 : {
2903 13853 : pari_sp av = avma;
2904 : GEN Y, A, u, v;
2905 : nfmaxord_t S;
2906 : long i, l;
2907 :
2908 13853 : v = polredabs_i(x, &S, &u, flag);
2909 13853 : remove_duplicates(v);
2910 13853 : Y = gel(v,1);
2911 13853 : A = gel(v,2);
2912 13853 : l = lg(A); if (l == 1) pari_err_BUG("polredabs (missing vector)");
2913 13853 : if (DEBUGLEVEL) err_printf("Found %ld minimal polynomials.\n",l-1);
2914 13853 : if (!(flag & nf_ALL))
2915 : {
2916 13846 : GEN y = findmindisc(Y);
2917 28127 : for (i = 1; i < l; i++)
2918 28127 : if (ZX_equal(gel(Y,i), y)) break;
2919 13846 : Y = mkvec(gel(Y,i));
2920 13846 : A = mkvec(gel(A,i)); l = 2;
2921 : }
2922 14119 : if (flag & (nf_RAW|nf_ORIG)) for (i = 1; i < l; i++)
2923 : {
2924 266 : GEN y = gel(Y,i), a = gel(A,i);
2925 266 : if (u) a = RgV_RgC_mul(S.basis, ZM_ZC_mul(u, a));
2926 266 : if (flag & nf_ORIG)
2927 : {
2928 259 : a = QXQ_reverse(a, S.T);
2929 259 : if (!isint1(S.unscale)) a = gdiv(a, S.unscale); /* not RgX_Rg_div */
2930 259 : a = mkpolmod(a,y);
2931 : }
2932 266 : gel(Y,i) = mkvec2(y, a);
2933 : }
2934 13853 : return gerepilecopy(av, (flag & nf_ALL)? Y: gel(Y,1));
2935 : }
2936 :
2937 : GEN
2938 0 : polredabsall(GEN x, long flun) { return polredabs0(x, flun | nf_ALL); }
2939 : GEN
2940 13398 : polredabs(GEN x) { return polredabs0(x,0); }
2941 : GEN
2942 0 : polredabs2(GEN x) { return polredabs0(x,nf_ORIG); }
2943 :
2944 : /* relative polredabs/best. Returns relative polynomial by default (flag = 0)
2945 : * flag & nf_ORIG: + element (base change)
2946 : * flag & nf_ABSOLUTE: absolute polynomial */
2947 : static GEN
2948 588 : rnfpolred_i(GEN nf, GEN R, long flag, long best)
2949 : {
2950 588 : const char *f = best? "rnfpolredbest": "rnfpolredabs";
2951 588 : const long abs = ((flag & nf_ORIG) && (flag & nf_ABSOLUTE));
2952 588 : GEN listP = NULL, red, pol, A, P, T, rnfeq;
2953 588 : pari_sp av = avma;
2954 :
2955 588 : if (typ(R) == t_VEC) {
2956 14 : if (lg(R) != 3) pari_err_TYPE(f,R);
2957 14 : listP = gel(R,2);
2958 14 : R = gel(R,1);
2959 : }
2960 588 : if (typ(R) != t_POL) pari_err_TYPE(f,R);
2961 588 : nf = checknf(nf);
2962 588 : T = nf_get_pol(nf);
2963 588 : R = RgX_nffix(f, T, R, 0);
2964 588 : if (best || (flag & nf_PARTIALFACT))
2965 : {
2966 364 : rnfeq = abs? nf_rnfeq(nf, R): nf_rnfeqsimple(nf, R);
2967 364 : pol = gel(rnfeq,1);
2968 364 : if (listP) pol = mkvec2(pol, listP);
2969 357 : red = best? polredbest_i(pol, abs? 1: 2)
2970 364 : : polredabs0(pol, (abs? nf_ORIG: nf_RAW)|nf_PARTIALFACT);
2971 364 : P = gel(red,1);
2972 364 : A = gel(red,2);
2973 : }
2974 : else
2975 : {
2976 : nfmaxord_t S;
2977 : GEN rnf, u, v, y, a;
2978 : long i, j, l;
2979 : pari_timer ti;
2980 224 : if (DEBUGLEVEL>1) timer_start(&ti);
2981 224 : rnf = rnfinit(nf, R);
2982 224 : rnfeq = rnf_get_map(rnf);
2983 224 : pol = rnf_zkabs(rnf);
2984 224 : if (DEBUGLEVEL>1) timer_printf(&ti, "absolute basis");
2985 224 : v = polredabs_i(pol, &S, &u, nf_ORIG);
2986 224 : pol = gel(pol,1);
2987 224 : y = gel(v,1); P = findmindisc(y);
2988 224 : a = gel(v,2);
2989 224 : l = lg(y); A = cgetg(l, t_VEC);
2990 1260 : for (i = j = 1; i < l; i++)
2991 1036 : if (ZX_equal(gel(y,i),P))
2992 : {
2993 924 : GEN t = gel(a,i);
2994 924 : if (u) t = RgV_RgC_mul(S.basis, ZM_ZC_mul(u,t));
2995 924 : gel(A,j++) = t;
2996 : }
2997 224 : setlg(A,j); /* mod(A[i], pol) are all roots of P in Q[X]/(pol) */
2998 : }
2999 588 : if (DEBUGLEVEL>1) err_printf("reduced absolute generator: %Ps\n",P);
3000 588 : if (flag & nf_ABSOLUTE)
3001 : {
3002 14 : if (flag & nf_ORIG)
3003 : {
3004 7 : GEN a = gel(rnfeq,2); /* Mod(a,pol) root of T */
3005 7 : GEN k = gel(rnfeq,3); /* Mod(variable(R),R) + k*a root of pol */
3006 7 : if (typ(A) == t_VEC) A = gel(A,1); /* any root will do */
3007 7 : a = RgX_RgXQ_eval(a, lift_shallow(A), P); /* Mod(a, P) root of T */
3008 7 : P = mkvec3(P, mkpolmod(a,P), gsub(A, gmul(k,a)));
3009 : }
3010 14 : return gerepilecopy(av, P);
3011 : }
3012 574 : if (typ(A) != t_VEC)
3013 : {
3014 357 : A = eltabstorel_lift(rnfeq, A);
3015 357 : P = lift_if_rational( RgXQ_charpoly(A, R, varn(R)) );
3016 : }
3017 : else
3018 : { /* canonical factor */
3019 217 : long i, l = lg(A), v = varn(R);
3020 217 : GEN besta = NULL;
3021 1099 : for (i = 1; i < l; i++)
3022 : {
3023 882 : GEN a = eltabstorel_lift(rnfeq, gel(A,i));
3024 882 : GEN p = lift_if_rational( RgXQ_charpoly(a, R, v) );
3025 882 : if (i == 1 || cmp_universal(p, P) < 0) { P = p; besta = a; }
3026 : }
3027 217 : A = besta;
3028 : }
3029 574 : if (flag & nf_ORIG) P = mkvec2(P, mkpolmod(RgXQ_reverse(A,R),P));
3030 574 : return gerepilecopy(av, P);
3031 : }
3032 : GEN
3033 231 : rnfpolredabs(GEN nf, GEN R, long flag)
3034 231 : { return rnfpolred_i(nf,R,flag, 0); }
3035 : GEN
3036 357 : rnfpolredbest(GEN nf, GEN R, long flag)
3037 : {
3038 357 : if (flag < 0 || flag > 3) pari_err_FLAG("rnfpolredbest");
3039 357 : return rnfpolred_i(nf,R,flag, 1);
3040 : }
|
