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 281035 : v13checkrnf(GEN v)
30 281035 : { return typ(gel(v,6)) == t_VEC; }
31 : static int
32 50631 : rawcheckbnf(GEN v) { return typ(v)==t_VEC && lg(v)==11; }
33 : static int
34 44401 : 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 37009 : v11checkbnf(GEN v) { return rawchecknf(bnf_get_nf(v)); }
38 : /* v a t_VEC, lg(v) = 12, sanity check for true bnf */
39 : static int
40 4655 : v13checkgchar(GEN v) { return rawcheckbnf(gchar_get_bnf(v)); }
41 : /* v a t_VEC, lg(v) = 10, sanity check for true nf */
42 : static int
43 39004 : 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 290117 : checkrnf_i(GEN rnf)
51 290117 : { return (typ(rnf)==t_VEC && lg(rnf)==13 && v13checkrnf(rnf)); }
52 :
53 : void
54 280461 : checkrnf(GEN rnf)
55 280461 : { if (!checkrnf_i(rnf)) pari_err_TYPE("checkrnf",rnf); }
56 :
57 : GEN
58 2634013 : checkbnf_i(GEN X)
59 : {
60 2634013 : if (typ(X) == t_VEC)
61 2633847 : switch (lg(X))
62 : {
63 2630725 : case 11:
64 2630725 : if (typ(gel(X,6)) != t_INT) return NULL; /* pre-2.2.4 format */
65 2630725 : if (lg(gel(X,10)) != 4) return NULL; /* pre-2.8.1 format */
66 2630725 : return X;
67 2807 : case 7: return checkbnf_i(bnr_get_bnf(X));
68 : }
69 481 : return NULL;
70 : }
71 :
72 : GEN
73 71607668 : checknf_i(GEN X)
74 : {
75 71607668 : if (typ(X)==t_VEC)
76 70958696 : switch(lg(X))
77 : {
78 70478474 : case 10: return X;
79 468595 : case 11: return checknf_i(bnf_get_nf(X));
80 2828 : case 7: return checknf_i(bnr_get_bnf(X));
81 4571 : case 3: if (typ(gel(X,2)) == t_POLMOD) return checknf_i(gel(X,1));
82 : }
83 657764 : return NULL;
84 : }
85 :
86 : GEN
87 1624435 : checkbnf(GEN x)
88 : {
89 1624435 : GEN bnf = checkbnf_i(x);
90 1624431 : if (!bnf) pari_err_TYPE("checkbnf [please apply bnfinit()]",x);
91 1624412 : return bnf;
92 : }
93 :
94 : GEN
95 69742832 : checknf(GEN x)
96 : {
97 69742832 : GEN nf = checknf_i(x);
98 69742799 : if (!nf) pari_err_TYPE("checknf [please apply nfinit()]",x);
99 69742777 : return nf;
100 : }
101 :
102 : GEN
103 1005561 : checkbnr_i(GEN bnr)
104 : {
105 1005561 : if (typ(bnr)!=t_VEC || lg(bnr)!=7 || !checkbnf_i(bnr_get_bnf(bnr)))
106 84 : return NULL;
107 1005476 : return bnr;
108 : }
109 : void
110 822620 : checkbnr(GEN bnr)
111 : {
112 822620 : if (!checkbnr_i(bnr))
113 14 : pari_err_TYPE("checkbnr [please apply bnrinit()]",bnr);
114 822598 : }
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 640918 : checkbid_i(GEN bid)
125 : {
126 : GEN f;
127 640918 : if (typ(bid)!=t_VEC || lg(bid)!=6 || typ(bid_get_U(bid)) != t_VEC)
128 253908 : return NULL;
129 387010 : f = bid_get_mod(bid);
130 387010 : if (typ(f)!=t_VEC || lg(f)!=3) return NULL;
131 387010 : return bid;
132 : }
133 : void
134 384679 : checkbid(GEN bid)
135 : {
136 384679 : if (!checkbid_i(bid)) pari_err_TYPE("checkbid",bid);
137 384672 : }
138 : void
139 19782 : checkabgrp(GEN v)
140 : {
141 19782 : if (typ(v) == t_VEC) switch(lg(v))
142 : {
143 17402 : case 4: if (typ(gel(v,3)) != t_VEC) break;
144 19782 : case 3: if (typ(gel(v,2)) != t_VEC) break;
145 19754 : if (typ(gel(v,1)) != t_INT) break;
146 19754 : return;/*OK*/
147 0 : default: break;
148 : }
149 28 : pari_err_TYPE("checkabgrp",v);
150 : }
151 :
152 : GEN
153 219107 : checknfelt_mod(GEN nf, GEN x, const char *s)
154 : {
155 219107 : GEN T = gel(x,1), a = gel(x,2), Tnf = nf_get_pol(nf);
156 219107 : if (!RgX_equal_var(T, Tnf)) pari_err_MODULUS(s, T, Tnf);
157 219051 : return a;
158 : }
159 :
160 : int
161 9331 : check_ZKmodule_i(GEN M)
162 : {
163 9331 : return (typ(M) ==t_VEC && lg(M) >= 3
164 9331 : && typ(gel(M,1)) == t_MAT
165 9331 : && typ(gel(M,2)) == t_VEC
166 18662 : && lgcols(M) == lg(gel(M,2)));
167 : }
168 : void
169 9275 : check_ZKmodule(GEN M, const char *s)
170 9275 : { 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 21973 : get_bnf(GEN x, long *t)
193 : {
194 21973 : 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 21350 : case t_VEC:
199 21350 : 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 2366 : case 7: *t = typ_BNR;
205 2366 : x = bnr_get_bnf(x);
206 2366 : if (!rawcheckbnf(x)) break;
207 2310 : return x;
208 84 : case 9:
209 84 : if (!v9checkgal(x)) break;
210 84 : *t = typ_GAL; return NULL;
211 420 : case 10:
212 420 : if (!v10checknf(x)) break;
213 420 : *t = typ_NF; return NULL;
214 448 : case 11:
215 448 : if (!v11checkbnf(x)) break;
216 448 : *t = typ_BNF; return x;
217 196 : case 13:
218 196 : 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 115451 : get_nf(GEN x, long *t)
233 : {
234 115451 : 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 112504 : case t_VEC:
239 112504 : 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 1078 : case 10:
256 1078 : if (!v10checknf(x)) break;
257 1078 : *t = typ_NF; return x;
258 224 : case 11:
259 224 : if (!rawchecknf(x = bnf_get_nf(x))) break;
260 224 : *t = typ_BNF; return x;
261 364 : case 13:
262 364 : if (v13checkgchar(x)) { *t = typ_GCHAR; return gchar_get_nf(x); }
263 350 : if (!v13checkrnf(x)) break;
264 350 : *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 83895 : nftyp(GEN x)
278 : {
279 83895 : switch(typ(x))
280 : {
281 21 : case t_POL : return typ_POL;
282 7 : case t_QUAD: return typ_Q;
283 83860 : case t_VEC:
284 83860 : switch(lg(x))
285 : {
286 4095 : case 13:
287 4095 : if (v13checkgchar(x)) return typ_GCHAR;
288 161 : if (!v13checkrnf(x)) break;
289 161 : return typ_RNF;
290 37506 : case 10:
291 37506 : if (!v10checknf(x)) break;
292 37499 : return typ_NF;
293 266 : case 11:
294 266 : if (!v11checkbnf(x)) break;
295 266 : return typ_BNF;
296 36309 : case 7:
297 36309 : x = bnr_get_bnf(x);
298 36309 : if (!rawcheckbnf(x) || !v11checkbnf(x)) break;
299 36295 : 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 7593 : tschirnhaus(GEN x)
319 : {
320 7593 : pari_sp av = avma, av2;
321 7593 : long a, v = varn(x);
322 7593 : GEN u, y = cgetg(5,t_POL);
323 :
324 7593 : if (typ(x)!=t_POL) pari_err_TYPE("tschirnhaus",x);
325 7593 : if (lg(x) < 4) pari_err_CONSTPOL("tschirnhaus");
326 7593 : if (v) { u = leafcopy(x); setvarn(u,0); x=u; }
327 7593 : y[1] = evalsigne(1)|evalvarn(0);
328 : do
329 : {
330 8208 : a = random_bits(2); if (a==0) a = 1; gel(y,4) = stoi(a);
331 8208 : a = random_bits(3); if (a>=4) a -= 8; gel(y,3) = stoi(a);
332 8208 : a = random_bits(3); if (a>=4) a -= 8; gel(y,2) = stoi(a);
333 8208 : u = RgXQ_charpoly(y,x,v); av2 = avma;
334 : }
335 8208 : while (degpol(RgX_gcd(u,RgX_deriv(u)))); /* while u not separable */
336 7593 : if (DEBUGLEVEL>1)
337 0 : err_printf("Tschirnhaus transform. New pol: %Ps",u);
338 7593 : 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 843205 : ZX_Z_normalize(GEN pol, GEN *ptk)
346 : {
347 843205 : long i,j, sk, n = degpol(pol); /* > 0 */
348 : GEN k, fa, P, E, a, POL;
349 :
350 843200 : if (ptk) *ptk = gen_1;
351 843200 : if (!n) return pol;
352 843193 : a = pol + 2; k = gel(a,n-1); /* a[i] = coeff of degree i */
353 1621577 : for (i = n-2; i >= 0; i--)
354 : {
355 1366148 : k = gcdii(k, gel(a,i));
356 1365972 : if (is_pm1(k)) return pol;
357 : }
358 255429 : sk = signe(k);
359 255429 : if (!sk) return pol; /* monomial! */
360 253609 : fa = absZ_factor_limit(k, 0); k = gen_1;
361 253652 : P = gel(fa,1);
362 253652 : E = gel(fa,2);
363 253652 : POL = leafcopy(pol); a = POL+2;
364 544499 : for (i = lg(P)-1; i > 0; i--)
365 : {
366 290843 : GEN p = gel(P,i), pv, pvj;
367 290843 : long vmin = itos(gel(E,i));
368 : /* find v_p(k) = min floor( v_p(a[i]) / (n-i)) */
369 1203699 : for (j=n-1; j>=0; j--)
370 : {
371 : long v;
372 912857 : if (!signe(gel(a,j))) continue;
373 784938 : v = Z_pval(gel(a,j), p) / (n - j);
374 784935 : if (v < vmin) vmin = v;
375 : }
376 290842 : if (!vmin) continue;
377 18073 : pvj = pv = powiu(p,vmin); k = mulii(k, pv);
378 : /* a[j] /= p^(v*(n-j)) */
379 73967 : for (j=n-1; j>=0; j--)
380 : {
381 55893 : if (j < n-1) pvj = mulii(pvj, pv);
382 55893 : gel(a,j) = diviiexact(gel(a,j), pvj);
383 : }
384 : }
385 253656 : if (ptk) *ptk = k;
386 253656 : 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 840564 : ZX_primitive_to_monic(GEN pol, GEN *pL)
394 : {
395 840564 : long i,j, n = degpol(pol);
396 840567 : GEN lc = leading_coeff(pol), L, fa, P, E, a, POL;
397 :
398 840568 : if (is_pm1(lc))
399 : {
400 837697 : if (pL) *pL = gen_1;
401 837697 : return signe(lc) < 0? ZX_neg(pol): pol;
402 : }
403 2870 : if (signe(lc) < 0)
404 35 : POL = ZX_neg(pol);
405 : else
406 2835 : POL = leafcopy(pol);
407 2870 : a = POL+2; lc = gel(a,n);
408 2870 : fa = absZ_factor_limit(lc,0); L = gen_1;
409 2870 : P = gel(fa,1);
410 2870 : E = gel(fa,2);
411 5887 : for (i = lg(P)-1; i > 0; i--)
412 : {
413 3017 : GEN p = gel(P,i), pk, pku;
414 3017 : long v, j0, e = itos(gel(E,i)), k = e/n, d = k*n - e;
415 :
416 3017 : 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 17087 : for (j=n-1; j>0; j--)
419 : {
420 14070 : if (!signe(gel(a,j))) continue;
421 8687 : v = Z_pval(gel(a,j), p);
422 8757 : while (v + d < k * j) { k++; d += n; }
423 : }
424 3017 : pk = powiu(p,k); j0 = d/k;
425 3017 : L = mulii(L, pk);
426 :
427 3017 : pku = powiu(p,d - k*j0);
428 : /* a[j] *= p^(d - kj) */
429 7063 : for (j=j0; j>=0; j--)
430 : {
431 4046 : if (j < j0) pku = mulii(pku, pk);
432 4046 : gel(a,j) = mulii(gel(a,j), pku);
433 : }
434 3017 : j0++;
435 3017 : pku = powiu(p,k*j0 - d);
436 : /* a[j] /= p^(kj - d) */
437 19075 : for (j=j0; j<=n; j++)
438 : {
439 16058 : if (j > j0) pku = mulii(pku, pk);
440 16058 : gel(a,j) = diviiexact(gel(a,j), pku);
441 : }
442 : }
443 2870 : if (pL) *pL = L;
444 2870 : 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 840248 : ZX_Q_normalize(GEN pol, GEN *pL)
452 : {
453 840248 : GEN lc, POL = ZX_primitive_to_monic(pol, &lc);
454 840262 : POL = ZX_Z_normalize(POL, pL);
455 840204 : if (pL) *pL = gdiv(lc, *pL);
456 840276 : return POL;
457 : }
458 :
459 : GEN
460 5182390 : ZX_Q_mul(GEN A, GEN z)
461 : {
462 5182390 : pari_sp av = avma;
463 5182390 : long i, l = lg(A);
464 : GEN d, n, Ad, B, u;
465 5182390 : if (typ(z)==t_INT) return ZX_Z_mul(A,z);
466 4337461 : n = gel(z, 1); d = gel(z, 2);
467 4337461 : Ad = RgX_to_RgC(FpX_red(A, d), l-2);
468 4337461 : u = gcdii(d, FpV_factorback(Ad, NULL, d));
469 4337456 : B = cgetg(l, t_POL);
470 4337461 : B[1] = A[1];
471 4337461 : if (equali1(u))
472 : {
473 4241159 : for(i=2; i<l; i++)
474 2644574 : gel(B, i) = mkfrac(mulii(n, gel(A,i)), d);
475 : } else
476 : {
477 16959945 : for(i=2; i<l; i++)
478 : {
479 14219079 : GEN di = gcdii(gel(Ad, i-1), u);
480 14219043 : GEN ni = mulii(n, diviiexact(gel(A,i), di));
481 14219031 : if (equalii(d, di))
482 3282136 : gel(B, i) = ni;
483 : else
484 10936925 : gel(B, i) = mkfrac(ni, diviiexact(d, di));
485 : }
486 : }
487 4337451 : 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 8778 : ZXX_Q_mul(GEN A, GEN z)
512 : {
513 : long i, l;
514 : GEN B;
515 8778 : if (typ(z)==t_INT) return ZXX_Z_mul(A,z);
516 8330 : B = cgetg_copy(A, &l);
517 8330 : B[1] = A[1];
518 116543 : for (i=2; i<l; i++)
519 : {
520 108213 : GEN Ai = gel(A,i);
521 108213 : gel(B,i) = typ(Ai)==t_POL ? ZX_Q_mul(Ai, z): gmul(Ai, z);
522 : }
523 8330 : 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 94696 : QX_table_nfpoleval(GEN nf, GEN pol, GEN m)
542 : {
543 94696 : pari_sp av = avma;
544 94696 : long i = lg(pol)-1;
545 : GEN res, den;
546 94696 : if (i==1) return gen_0;
547 94696 : pol = Q_remove_denom(pol, &den);
548 94701 : res = scalarcol_shallow(gel(pol,i), nf_get_degree(nf));
549 392780 : for (i-- ; i>=2; i--)
550 298077 : res = ZC_Z_add(ZM_ZC_mul(m, res), gel(pol,i));
551 94703 : if (den) res = RgC_Rg_div(res, den);
552 94702 : 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 94892 : ZC_galoisapply(GEN nf, GEN s, GEN x)
575 : {
576 94892 : x = nf_to_scalar_or_alg(nf, x);
577 94891 : if (typ(x) != t_POL) return scalarcol(x, nf_get_degree(nf));
578 94702 : 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 126152 : galoisapply(GEN nf, GEN aut, GEN x)
673 : {
674 126152 : pari_sp av = avma;
675 : long lx;
676 : GEN y;
677 :
678 126152 : nf = checknf(nf);
679 126154 : 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 1288 : case t_POL:
686 1288 : aut = algtobasis(nf, aut);
687 1288 : y = basistoalg(nf, ZC_galoisapply(nf, aut, x));
688 1288 : 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 88347 : case t_COL:
705 88347 : aut = algtobasis(nf, aut);
706 88348 : 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 94310 : nfgaloismatrix(GEN nf, GEN s)
764 : {
765 94310 : pari_sp av2, av = avma;
766 : GEN zk, D, M, H, m;
767 : long k, n;
768 :
769 94310 : nf = checknf(nf);
770 94310 : zk = nf_get_zkprimpart(nf); n = lg(zk)-1;
771 94310 : M = cgetg(n+1, t_MAT);
772 94310 : gel(M,1) = col_ei(n, 1); /* s(1) = 1 */
773 94310 : if (n == 1) return M;
774 94310 : av2 = avma;
775 94310 : if (typ(s) != t_COL) s = algtobasis(nf, s);
776 94310 : D = nf_get_zkden(nf);
777 94310 : H = RgV_to_RgM(zk, n);
778 94311 : if (n == 2)
779 : {
780 62902 : GEN t = gel(H,2); /* D * s(w_2) */
781 62902 : t = ZC_Z_add(ZC_Z_mul(s, gel(t,2)), gel(t,1));
782 62902 : gel(M,2) = gerepileupto(av2, gdiv(t, D));
783 62902 : return M;
784 : }
785 31409 : m = zk_multable(nf, s);
786 31409 : gel(M,2) = s; /* M[,k] = s(x^(k-1)) */
787 137771 : for (k = 3; k <= n; k++) gel(M,k) = ZM_ZC_mul(m, gel(M,k-1));
788 31406 : M = ZM_mul(M, H);
789 31408 : if (!equali1(D)) M = ZM_Z_divexact(M, D);
790 31407 : 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 586331 : get_nfpol(GEN x, GEN *nf)
1091 : {
1092 586331 : if (typ(x) == t_POL) { *nf = NULL; return x; }
1093 300370 : *nf = checknf(x); return nf_get_pol(*nf);
1094 : }
1095 :
1096 : static GEN
1097 539 : incl_disc(GEN nfa, GEN a, int nolocal)
1098 : {
1099 : GEN d;
1100 539 : if (nfa) return nf_get_disc(nfa);
1101 469 : if (nolocal) return NULL;
1102 462 : d = ZX_disc(a);
1103 462 : if (!signe(d)) pari_err_IRREDPOL("nfisincl",a);
1104 455 : 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 287 : tests_OK(GEN a, GEN nfa, GEN b, GEN nfb, long fliso)
1122 : {
1123 : GEN da2, da, db, fa, P, U;
1124 287 : long i, l, q, m = degpol(a), n = degpol(b);
1125 :
1126 287 : if (m <= 0) pari_err_IRREDPOL("nfisincl",a);
1127 287 : if (n <= 0) pari_err_IRREDPOL("nfisincl",b);
1128 280 : q = n / m; /* relative degree */
1129 280 : if (fliso) { if (n != m) return 0; } else { if (n % m) return 0; }
1130 280 : if (m == 1) return 1;
1131 :
1132 : /*local test expensive if n^2 >> m^4 <=> q = n/m >> m */
1133 273 : db = incl_disc(nfb, b, q > m);
1134 273 : da = db? incl_disc(nfa, a, 0): NULL;
1135 266 : 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 259 : if (!db) return 1;
1142 252 : 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 21 : y = galoisconj(nfb? nfb: b, NULL); settyp(y, t_VEC);
1198 21 : return gerepilecopy(av,y);
1199 : }
1200 56 : if (nfa && !nfb) { swap(a,b); nfb = nfa; nfa = NULL; sw = 1; }
1201 56 : if (!tests_OK(a, nfa, b, nfb, 1)) { set_avma(av); return gen_0; }
1202 :
1203 49 : if (nfb) lb = gen_1; else nfb = b = ZX_Q_normalize(b,&lb);
1204 49 : if (nfa) la = gen_1; else nfa = a = ZX_Q_normalize(a,&la);
1205 49 : va = varn(a); vb = varn(b); newvar = (varncmp(vb,va) <= 0);
1206 49 : if (newvar) { a = leafcopy(a); setvarn(a, fetch_var_higher()); }
1207 49 : y = lift_shallow(nfroots(nfb,a));
1208 49 : if (newvar) (void)delete_var();
1209 49 : lx = lg(y); if (lx==1) { set_avma(av); return gen_0; }
1210 49 : if (sw) { vb = va; b = leafcopy(b); setvarn(b, vb); }
1211 308 : for (i=1; i<lx; i++)
1212 : {
1213 259 : GEN t = gel(y,i);
1214 259 : if (typ(t) == t_POL) setvarn(t, vb); else t = scalarpol(t, vb);
1215 259 : if (lb != gen_1) t = RgX_unscale(t, lb);
1216 259 : if (la != gen_1) t = RgX_Rg_div(t, la);
1217 259 : gel(y,i) = sw? RgXQ_reverse(t, b): t;
1218 : }
1219 49 : return gerepilecopy(av,y);
1220 : }
1221 :
1222 : static GEN
1223 7294 : partmap_reverse(GEN a, GEN b, GEN t, GEN la, GEN lb, long v)
1224 : {
1225 7294 : pari_sp av = avma;
1226 7294 : GEN rnf = rnfequation2(a, t), z;
1227 7294 : if (!RgX_equal(b, gel(rnf,1)))
1228 7 : { setvarn(b,v); pari_err_IRREDPOL("nfisincl", b); }
1229 7287 : z = liftpol_shallow(gel(rnf, 2));
1230 7286 : setvarn(z, v);
1231 7286 : if (!isint1(lb)) z = RgX_unscale(z, lb);
1232 7286 : if (!isint1(la)) z = RgX_Rg_div(z, la);
1233 7286 : 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 8084 : N = gel(L,1); if (!signe(N)) { set_avma(av); return pol_0(v); }
1247 8077 : D = gel(L,2);
1248 8077 : N = RgX_neg(N); setvarn(N, v); setvarn(D, v);
1249 8078 : G = QX_gcd(N,D);
1250 8077 : if (degpol(G)) { N = RgX_div(N,G); D = RgX_div(D,G); }
1251 8077 : if (!isint1(lb)) { N = RgX_unscale(N, lb); D = RgX_unscale(D, lb); }
1252 8077 : if (!isint1(la)) D = RgX_Rg_mul(D, la);
1253 8077 : 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 7175 : nfisincl_from_fact(GEN a, long da, GEN b, GEN la, GEN lb, long vb, GEN y,
1262 : long flag)
1263 : {
1264 7175 : long i, k, l = lg(y), db = degpol(b), d = db / da;
1265 7175 : GEN x = cgetg(l, t_VEC);
1266 7574 : for (i= k = 1; i < l; i++)
1267 : {
1268 7392 : GEN t = gel(y,i);
1269 7392 : if (degpol(t) != d) continue;
1270 7294 : gel(x, k++) = partmap_reverse(a, b, t, la, lb, vb);
1271 7287 : 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 7280 : nfisincl0(GEN fa, GEN fb, long flag)
1298 : {
1299 7280 : pari_sp av = avma;
1300 : long vb;
1301 : GEN a, b, nfa, nfb, x, y, la, lb;
1302 : int newvar;
1303 7280 : if (flag < 0 || flag > 3) pari_err_FLAG("nfisincl");
1304 :
1305 7280 : a = get_nfpol(fa, &nfa);
1306 7280 : b = get_nfpol(fb, &nfb);
1307 7280 : if (!nfa) { a = Q_primpart(a); RgX_check_ZX(a, "nsisincl"); }
1308 7280 : if (!nfb) { b = Q_primpart(b); RgX_check_ZX(b, "nsisincl"); }
1309 7280 : if (ZX_equal(a, b) && flag<=1)
1310 : {
1311 14 : if (flag==1)
1312 : {
1313 7 : x = pol_x(varn(b));
1314 7 : return degpol(b) > 1 ? x: RgX_rem(x,b);
1315 : }
1316 7 : x = galoisconj(fb, NULL); settyp(x, t_VEC);
1317 7 : return gerepilecopy(av,x);
1318 : }
1319 7266 : if (flag==0 && !tests_OK(a, nfa, b, nfb, 0)) { set_avma(av); return gen_0; }
1320 :
1321 7147 : if (nfb) lb = gen_1; else nfb = b = ZX_Q_normalize(b,&lb);
1322 7147 : if (nfa) la = gen_1; else nfa = a = ZX_Q_normalize(a,&la);
1323 7147 : vb = varn(b); newvar = (varncmp(varn(a),vb) <= 0);
1324 7147 : if (newvar) { b = leafcopy(b); setvarn(b, fetch_var_higher()); }
1325 7147 : y = lift_shallow(gel(nffactor(nfa,b),1));
1326 7147 : if (flag==2)
1327 0 : x = nfisincl_from_fact_frac(a, b, la, lb, vb, y);
1328 : else
1329 7147 : x = nfisincl_from_fact(nfa, degpol(a), b, la, lb, vb, y, flag);
1330 7140 : if (newvar) (void)delete_var();
1331 7140 : return gerepilecopy(av,x);
1332 : }
1333 :
1334 : GEN
1335 105 : nfisincl(GEN fa, GEN fb) { return nfisincl0(fa, fb, 0); }
1336 :
1337 : static GEN
1338 14 : lastel(GEN x) { return gel(x, lg(x)-1); }
1339 :
1340 : static GEN
1341 8092 : RgF_to_Flxq(GEN F, GEN T, ulong p)
1342 : {
1343 : GEN N, D, iD;
1344 8092 : if (typ(F)==t_POL) return RgX_to_Flx(F, p);
1345 8085 : N = RgX_to_Flx(gel(F,1), p); D = RgX_to_Flx(gel(F,2), p);
1346 8085 : iD = Flxq_invsafe(D, T, p);
1347 8085 : if (!iD) return NULL;
1348 8078 : return Flxq_mul(N, iD, T, p);
1349 : }
1350 :
1351 : #define pari_APPLY_abort(EXPR)\
1352 : { \
1353 : long i, _l; \
1354 : GEN _y = cgetg_copy(x, &_l);\
1355 : for (i=1; i<_l; i++) \
1356 : { GEN _z = EXPR;\
1357 : if (!_z) return _z;\
1358 : gel(_y,i) = _z;\
1359 : } return _y;\
1360 : }
1361 :
1362 : static GEN
1363 1239 : RgFV_to_FlxqV(GEN x, GEN T, ulong p)
1364 9324 : { pari_APPLY_abort(RgF_to_Flxq(gel(x,i), T, p)) }
1365 :
1366 : static GEN
1367 1231 : nfsplitting_auto(GEN g, GEN R)
1368 : {
1369 1231 : pari_sp av = avma;
1370 : forprime_t T;
1371 1231 : long i, d = degpol(g);
1372 : ulong p;
1373 1231 : GEN P, K, N, G, q, den = Q_denom(R), Rp, Gp;
1374 1232 : u_forprime_init(&T, d*maxss(expu(d)-3, 2), ULONG_MAX);
1375 1232 : av = avma;
1376 : for(;;)
1377 : {
1378 79867 : set_avma(av);
1379 79867 : p = u_forprime_next(&T);
1380 79867 : if (dvdiu(den,p)) continue;
1381 79867 : Gp = ZX_to_Flx(g, p);
1382 79863 : if (!Flx_is_totally_split(Gp, p)) continue;
1383 1239 : P = Flx_roots(Gp, p);
1384 1239 : Rp = RgFV_to_FlxqV(R, Gp, p);
1385 1239 : if (!Rp) { if (DEBUGLEVEL) err_printf("nfsplitting_auto: bad p : %lu\n",p); continue; }
1386 1232 : if (d == 1) return mkvec3(g, mkcol(gel(Rp,1)), utoi(p));
1387 1225 : K = Flm_Flc_invimage(FlxV_to_Flm(Rp, d), vecsmall_ei(d, 2), p);
1388 1225 : if (!K) pari_err_BUG("nfsplitting_auto");
1389 1225 : N = Flm_transpose(FlxV_Flv_multieval(Rp, P, p));
1390 1225 : q = perm_inv(vecsmall_indexsort(gel(N,1)));
1391 1225 : G = cgetg(d+1, t_COL);
1392 35707 : for (i=1; i<=d; i++)
1393 : {
1394 34482 : GEN r = perm_mul(vecsmall_indexsort(gel(N,i)), q);
1395 34482 : gel(G,i) = FlxV_Flc_mul(Rp, vecpermute(K, r), p);
1396 : }
1397 1225 : return mkvec3(g, G, utoi(p));
1398 : }
1399 : }
1400 :
1401 : static GEN
1402 1420 : nfsplitting_composite(GEN P)
1403 : {
1404 1420 : GEN F = gel(ZX_factor(P), 1), Q = NULL;
1405 1421 : long i, n = lg(F)-1;
1406 2849 : for (i = 1; i <= n; i++)
1407 : {
1408 1428 : GEN Fi = gel(F, i);
1409 1428 : if (degpol(Fi) == 1) continue;
1410 1400 : Q = Q ? lastel(compositum(Q, Fi)): Fi;
1411 : }
1412 1421 : return Q ? Q: pol_x(varn(P));
1413 : }
1414 : GEN
1415 1435 : nfsplitting0(GEN T0, GEN D, long flag)
1416 : {
1417 1435 : pari_sp av = avma;
1418 : long d, Ds, v;
1419 1435 : GEN T, F, K, N = NULL, lT = NULL;
1420 1435 : if (flag < 0 || flag > 3) pari_err_FLAG("nfsplitting");
1421 1435 : T = T0 = get_nfpol(T0, &K);
1422 1428 : if (!K)
1423 : {
1424 : GEN c;
1425 1407 : if (typ(T) != t_POL) pari_err_TYPE("nfsplitting",T);
1426 1407 : T = Q_primitive_part(T, &c);
1427 1406 : lT = leading_coeff(T); if (isint1(lT)) lT = NULL;
1428 1406 : if (flag && (c || lT)) pari_err_TYPE("nfsplitting", T0);
1429 1399 : RgX_check_ZX(T,"nfsplitting");
1430 : }
1431 1421 : T = nfsplitting_composite(T);
1432 1421 : if (flag && !ZX_equal(T, T0)) pari_err_IRREDPOL("nfsplitting", T0);
1433 1407 : d = degpol(T); v = varn(T);
1434 1407 : if (d <= 1 && !flag) { set_avma(av); return pol_x(v); }
1435 1379 : if (!K) {
1436 1358 : if (lT) T = polredbest(T,0);
1437 1358 : K = T;
1438 : }
1439 1379 : if (D)
1440 1239 : { if (typ(D) != t_INT || signe(D) < 1) pari_err_TYPE("nfsplitting",D); }
1441 140 : else if (d <= 7 ||
1442 35 : (d <= 11 && pari_is_dir(stack_strcat(pari_datadir, "/galdata"))))
1443 126 : D = gel(polgalois(T,DEFAULTPREC), 1);
1444 : else
1445 14 : D = mpfact(d);
1446 1379 : Ds = itos_or_0(D);
1447 1379 : T = leafcopy(T); setvarn(T, fetch_var_higher());
1448 1379 : for(F = T;;)
1449 1659 : {
1450 3038 : GEN P = gel(nffactor(K, F), 1), Q = gel(P,lg(P)-1);
1451 3038 : if (degpol(gel(P,1)) == degpol(Q))
1452 : {
1453 1288 : if (!flag) break;
1454 1260 : P = liftall_shallow(P);
1455 1260 : if (flag==1)
1456 28 : N = nfisincl_from_fact(K, d, F, gen_1, gen_1, v, P, flag);
1457 : else
1458 1232 : N = nfisincl_from_fact_frac(T0, F, gen_1, gen_1, v, P);
1459 1259 : break;
1460 : }
1461 1750 : F = rnfequation(K,Q);
1462 1750 : if (degpol(F) == Ds && !flag) break;
1463 : }
1464 1378 : if (umodiu(D,degpol(F)))
1465 : {
1466 14 : char *sD = itostr(D);
1467 14 : pari_warn(warner,stack_strcat("ignoring incorrect degree bound ",sD));
1468 : }
1469 1378 : setvarn(F, v); (void)delete_var();
1470 1378 : if (flag) F = flag == 3? nfsplitting_auto(F, N): mkvec2(F, N);
1471 1379 : return gerepilecopy(av, F);
1472 : }
1473 :
1474 : GEN
1475 0 : nfsplitting(GEN T, GEN D) { return nfsplitting0(T, D, 0); }
1476 :
1477 : /*************************************************************************/
1478 : /** **/
1479 : /** INITALG **/
1480 : /** **/
1481 : /*************************************************************************/
1482 : typedef struct {
1483 : GEN T;
1484 : GEN ro; /* roots of T */
1485 : long r1;
1486 : GEN basden;
1487 : long prec;
1488 : long extraprec; /* possibly -1 = irrelevant or not computed */
1489 : GEN M, G; /* possibly NULL = irrelevant or not computed */
1490 : } nffp_t;
1491 :
1492 : static GEN
1493 199800 : get_roots(GEN x, long r1, long prec)
1494 : {
1495 : long i, ru;
1496 : GEN z;
1497 199800 : if (typ(x) != t_POL)
1498 : {
1499 0 : z = leafcopy(x);
1500 0 : ru = (lg(z)-1 + r1) >> 1;
1501 : }
1502 : else
1503 : {
1504 199800 : long n = degpol(x);
1505 199800 : z = (r1 == n)? ZX_realroots_irred(x, prec): QX_complex_roots(x,prec);
1506 199793 : ru = (n+r1)>>1;
1507 : }
1508 436695 : for (i=r1+1; i<=ru; i++) gel(z,i) = gel(z, (i<<1)-r1);
1509 199793 : z[0]=evaltyp(t_VEC)|evallg(ru+1); return z;
1510 : }
1511 :
1512 : GEN
1513 0 : nf_get_allroots(GEN nf)
1514 : {
1515 0 : return embed_roots(nf_get_roots(nf), nf_get_r1(nf));
1516 : }
1517 :
1518 : /* For internal use. compute trace(x mod pol), sym=polsym(pol,deg(pol)-1) */
1519 : static GEN
1520 383792 : quicktrace(GEN x, GEN sym)
1521 : {
1522 383792 : GEN p1 = gen_0;
1523 : long i;
1524 :
1525 383792 : if (typ(x) != t_POL) return gmul(x, gel(sym,1));
1526 383792 : if (signe(x))
1527 : {
1528 383794 : sym--;
1529 2879848 : for (i=lg(x)-1; i>1; i--)
1530 2496130 : p1 = gadd(p1, gmul(gel(x,i),gel(sym,i)));
1531 : }
1532 383716 : return p1;
1533 : }
1534 :
1535 : static GEN
1536 82826 : get_Tr(GEN mul, GEN x, GEN basden)
1537 : {
1538 82826 : GEN t, bas = gel(basden,1), den = gel(basden,2);
1539 82826 : long i, j, n = lg(bas)-1;
1540 82826 : GEN T = cgetg(n+1,t_MAT), TW = cgetg(n+1,t_COL), sym = polsym(x, n-1);
1541 :
1542 82827 : gel(TW,1) = utoipos(n);
1543 260318 : for (i=2; i<=n; i++)
1544 : {
1545 177497 : t = quicktrace(gel(bas,i), sym);
1546 177488 : if (den && gel(den,i)) t = diviiexact(t,gel(den,i));
1547 177490 : gel(TW,i) = t; /* tr(w[i]) */
1548 : }
1549 82821 : gel(T,1) = TW;
1550 260317 : for (i=2; i<=n; i++)
1551 : {
1552 177494 : gel(T,i) = cgetg(n+1,t_COL); gcoeff(T,1,i) = gel(TW,i);
1553 694527 : for (j=2; j<=i; j++) /* Tr(W[i]W[j]) */
1554 517031 : gcoeff(T,i,j) = gcoeff(T,j,i) = ZV_dotproduct(gel(mul,j+(i-1)*n), TW);
1555 : }
1556 82823 : return T;
1557 : }
1558 :
1559 : /* return [bas[i]*denom(bas[i]), denom(bas[i])], denom 1 is given as NULL */
1560 : static GEN
1561 201415 : get_bas_den(GEN bas)
1562 : {
1563 201415 : GEN b,d,den, dbas = leafcopy(bas);
1564 201417 : long i, l = lg(bas);
1565 201417 : int power = 1;
1566 201417 : den = cgetg(l,t_VEC);
1567 934033 : for (i=1; i<l; i++)
1568 : {
1569 732621 : b = Q_remove_denom(gel(bas,i), &d);
1570 732615 : gel(dbas,i) = b;
1571 732615 : gel(den,i) = d; if (d) power = 0;
1572 : }
1573 201412 : if (power) den = NULL; /* power basis */
1574 201412 : return mkvec2(dbas, den);
1575 : }
1576 :
1577 : /* return multiplication table for S->basis */
1578 : static GEN
1579 82827 : nf_multable(nfmaxord_t *S, GEN invbas)
1580 : {
1581 82827 : GEN T = S->T, w = gel(S->basden,1), den = gel(S->basden,2);
1582 82827 : long i,j, n = degpol(T);
1583 82827 : GEN mul = cgetg(n*n+1,t_MAT);
1584 :
1585 : /* i = 1 split for efficiency, assume w[1] = 1 */
1586 343159 : for (j=1; j<=n; j++)
1587 260332 : gel(mul,j) = gel(mul,1+(j-1)*n) = col_ei(n, j);
1588 260323 : for (i=2; i<=n; i++)
1589 694532 : for (j=i; j<=n; j++)
1590 : {
1591 517036 : pari_sp av = avma;
1592 517036 : GEN z = (i == j)? ZXQ_sqr(gel(w,i), T): ZXQ_mul(gel(w,i),gel(w,j), T);
1593 517034 : z = ZM_ZX_mul(invbas, z); /* integral column */
1594 517008 : if (den)
1595 : {
1596 341061 : GEN d = mul_denom(gel(den,i), gel(den,j));
1597 341052 : if (d) z = ZC_Z_divexact(z, d);
1598 : }
1599 517005 : gel(mul,j+(i-1)*n) = gel(mul,i+(j-1)*n) = gerepileupto(av,z);
1600 : }
1601 82819 : return mul;
1602 : }
1603 :
1604 : /* as get_Tr, mul_table not precomputed */
1605 : static GEN
1606 27729 : make_Tr(nfmaxord_t *S)
1607 : {
1608 27729 : GEN T = S->T, w = gel(S->basden,1), den = gel(S->basden,2);
1609 27729 : long i,j, n = degpol(T);
1610 27729 : GEN c, t, d, M = cgetg(n+1,t_MAT), sym = polsym(T, n-1);
1611 :
1612 : /* W[i] = w[i]/den[i]; assume W[1] = 1, case i = 1 split for efficiency */
1613 27729 : c = cgetg(n+1,t_COL); gel(M,1) = c;
1614 27729 : gel(c, 1) = utoipos(n);
1615 83442 : for (j=2; j<=n; j++)
1616 : {
1617 55712 : pari_sp av = avma;
1618 55712 : t = quicktrace(gel(w,j), sym);
1619 55711 : if (den)
1620 : {
1621 35083 : d = gel(den,j);
1622 35083 : if (d) t = diviiexact(t, d);
1623 : }
1624 55710 : gel(c,j) = gerepileuptoint(av, t);
1625 : }
1626 83438 : for (i=2; i<=n; i++)
1627 : {
1628 55712 : c = cgetg(n+1,t_COL); gel(M,i) = c;
1629 206308 : for (j=1; j<i ; j++) gel(c,j) = gcoeff(M,i,j);
1630 206301 : for ( ; j<=n; j++)
1631 : {
1632 150593 : pari_sp av = avma;
1633 150593 : t = (i == j)? ZXQ_sqr(gel(w,i), T): ZXQ_mul(gel(w,i),gel(w,j), T);
1634 150596 : t = quicktrace(t, sym);
1635 150574 : if (den)
1636 : {
1637 117432 : d = mul_denom(gel(den,i),gel(den,j));
1638 117435 : if (d) t = diviiexact(t, d);
1639 : }
1640 150579 : gel(c,j) = gerepileuptoint(av, t); /* Tr (W[i]W[j]) */
1641 : }
1642 : }
1643 27726 : return M;
1644 : }
1645 :
1646 : /* [bas[i]/den[i]]= integer basis. roo = real part of the roots */
1647 : static void
1648 231352 : make_M(nffp_t *F, int trunc)
1649 : {
1650 231352 : GEN bas = gel(F->basden,1), den = gel(F->basden,2), ro = F->ro;
1651 : GEN m, d, M;
1652 231352 : long i, j, l = lg(ro), n = lg(bas);
1653 231352 : M = cgetg(n,t_MAT);
1654 231355 : gel(M,1) = const_col(l-1, gen_1); /* bas[1] = 1 */
1655 824161 : for (j=2; j<n; j++) gel(M,j) = cgetg(l,t_COL);
1656 757692 : for (i=1; i<l; i++)
1657 : {
1658 526350 : GEN r = gel(ro,i), ri;
1659 526350 : ri = (gexpo(r) > 1)? ginv(r): NULL;
1660 2634526 : for (j=2; j<n; j++) gcoeff(M,i,j) = RgX_cxeval(gel(bas,j), r, ri);
1661 : }
1662 231342 : if (den)
1663 515035 : for (j=2; j<n; j++)
1664 : {
1665 399735 : d = gel(den,j); if (!d) continue;
1666 337272 : m = gel(M,j);
1667 1689519 : for (i=1; i<l; i++) gel(m,i) = gdiv(gel(m,i), d);
1668 : }
1669 :
1670 231324 : if (trunc && gprecision(M) > F->prec)
1671 : {
1672 18045 : M = gprec_w(M, F->prec);
1673 18045 : F->ro = gprec_w(ro,F->prec);
1674 : }
1675 231324 : F->M = M;
1676 231324 : }
1677 :
1678 : /* return G real such that G~ * G = T_2 */
1679 : static void
1680 231355 : make_G(nffp_t *F)
1681 : {
1682 231355 : GEN G, M = F->M;
1683 231355 : long i, j, k, r1 = F->r1, l = lg(M);
1684 :
1685 231355 : if (r1 == l-1) { F->G = M; return; }
1686 185340 : G = cgetg(l, t_MAT);
1687 877515 : for (j = 1; j < l; j++)
1688 : {
1689 692237 : GEN g, m = gel(M,j);
1690 692237 : gel(G,j) = g = cgetg(l, t_COL);
1691 1130889 : for (k = i = 1; i <= r1; i++) gel(g,k++) = gel(m,i);
1692 2363044 : for ( ; k < l; i++)
1693 : {
1694 1670869 : GEN r = gel(m,i);
1695 1670869 : if (typ(r) == t_COMPLEX)
1696 : {
1697 1373057 : GEN a = gel(r,1), b = gel(r,2);
1698 1373057 : gel(g,k++) = mpadd(a, b);
1699 1373005 : gel(g,k++) = mpsub(a, b);
1700 : }
1701 : else
1702 : {
1703 297812 : gel(g,k++) = r;
1704 297812 : gel(g,k++) = r;
1705 : }
1706 : }
1707 : }
1708 185278 : F->G = G;
1709 : }
1710 :
1711 : static long
1712 263400 : prec_fix(long prec)
1713 : {
1714 : #ifndef LONG_IS_64BIT
1715 : /* make sure that default accuracy is the same on 32/64bit */
1716 37953 : if (odd(prec)) prec += EXTRAPRECWORD;
1717 : #endif
1718 263400 : return prec;
1719 : }
1720 : static void
1721 231358 : make_M_G(nffp_t *F, int trunc)
1722 : {
1723 : long n, eBD, prec;
1724 231358 : if (F->extraprec < 0)
1725 : { /* not initialized yet; compute roots so that absolute accuracy
1726 : * of M & G >= prec */
1727 : double er;
1728 231338 : n = degpol(F->T);
1729 231338 : eBD = 1 + gexpo(gel(F->basden,1));
1730 231339 : er = F->ro? (1+gexpo(F->ro)): fujiwara_bound(F->T);
1731 231341 : if (er < 0) er = 0;
1732 231341 : F->extraprec = nbits2extraprec(n*er + eBD + log2(n));
1733 : }
1734 231361 : prec = prec_fix(F->prec + F->extraprec);
1735 231361 : if (!F->ro || gprecision(gel(F->ro,1)) < prec)
1736 199800 : F->ro = get_roots(F->T, F->r1, prec);
1737 :
1738 231352 : make_M(F, trunc);
1739 231354 : make_G(F);
1740 231350 : }
1741 :
1742 : static void
1743 173749 : nffp_init(nffp_t *F, nfmaxord_t *S, long prec)
1744 : {
1745 173749 : F->T = S->T;
1746 173749 : F->r1 = S->r1;
1747 173749 : F->basden = S->basden;
1748 173749 : F->ro = NULL;
1749 173749 : F->extraprec = -1;
1750 173749 : F->prec = prec;
1751 173749 : }
1752 :
1753 : /* let bas a t_VEC of QX giving a Z-basis of O_K. Return the index of the
1754 : * basis. Assume bas[1] = 1 and that the leading coefficient of elements
1755 : * of bas are of the form 1/b for a t_INT b */
1756 : static GEN
1757 1498 : get_nfindex(GEN bas)
1758 : {
1759 1498 : pari_sp av = avma;
1760 1498 : long n = lg(bas)-1, i;
1761 : GEN D, d, mat;
1762 :
1763 : /* assume bas[1] = 1 */
1764 1498 : D = gel(bas,1);
1765 1498 : if (! is_pm1(simplify_shallow(D))) pari_err_TYPE("get_nfindex", D);
1766 1498 : D = gen_1;
1767 7917 : for (i = 2; i <= n; i++)
1768 : { /* after nfbasis, basis is upper triangular! */
1769 6426 : GEN B = gel(bas,i), lc;
1770 6426 : if (degpol(B) != i-1) break;
1771 6419 : lc = gel(B, i+1);
1772 6419 : switch (typ(lc))
1773 : {
1774 2471 : case t_INT: continue;
1775 3948 : case t_FRAC: if (is_pm1(gel(lc,1)) ) {D = mulii(D, gel(lc,2)); continue;}
1776 0 : default: pari_err_TYPE("get_nfindex", B);
1777 : }
1778 : }
1779 1498 : if (i <= n)
1780 : { /* not triangular after all */
1781 7 : bas = vecslice(bas,i,n);
1782 7 : bas = Q_remove_denom(bas, &d);
1783 7 : if (!d) return D;
1784 7 : mat = RgV_to_RgM(bas, n);
1785 7 : mat = rowslice(mat, i,n);
1786 7 : D = mulii(D, diviiexact(powiu(d, n-i+1), absi_shallow(ZM_det(mat))));
1787 : }
1788 1498 : return gerepileuptoint(av, D);
1789 : }
1790 : /* make sure all components of S are initialized */
1791 : static void
1792 165653 : nfmaxord_complete(nfmaxord_t *S)
1793 : {
1794 165653 : if (!S->dT) S->dT = ZX_disc(S->T);
1795 165653 : if (!S->index)
1796 : {
1797 1498 : if (S->dK) /* fast */
1798 0 : S->index = sqrti( diviiexact(S->dT, S->dK) );
1799 : else
1800 1498 : S->index = get_nfindex(S->basis);
1801 : }
1802 165653 : if (!S->dK) S->dK = diviiexact(S->dT, sqri(S->index));
1803 165653 : if (S->r1 < 0) S->r1 = ZX_sturm_irred(S->T);
1804 165651 : if (!S->basden) S->basden = get_bas_den(S->basis);
1805 165652 : }
1806 :
1807 : GEN
1808 82826 : nfmaxord_to_nf(nfmaxord_t *S, GEN ro, long prec)
1809 : {
1810 82826 : GEN nf = cgetg(10,t_VEC);
1811 82826 : GEN T = S->T, Tr, D, w, A, dA, MDI, mat = cgetg(9,t_VEC);
1812 82826 : long n = degpol(T);
1813 : nffp_t F;
1814 82826 : nfmaxord_complete(S);
1815 82826 : nffp_init(&F,S,prec);
1816 82826 : F.ro = ro;
1817 82826 : make_M_G(&F, 0);
1818 :
1819 82826 : gel(nf,1) = S->T;
1820 82826 : gel(nf,2) = mkvec2s(S->r1, (n - S->r1)>>1);
1821 82826 : gel(nf,3) = S->dK;
1822 82826 : gel(nf,4) = S->index;
1823 82826 : gel(nf,5) = mat;
1824 82826 : if (gprecision(gel(F.ro,1)) > prec) F.ro = gprec_wtrunc(F.ro, prec);
1825 82827 : gel(nf,6) = F.ro;
1826 82827 : w = S->basis;
1827 82827 : if (!is_pm1(S->index)) w = Q_remove_denom(w, NULL);
1828 82827 : gel(nf,7) = w;
1829 82827 : gel(nf,8) = ZM_inv(RgV_to_RgM(w,n), NULL);
1830 82827 : gel(nf,9) = nf_multable(S, nf_get_invzk(nf));
1831 82826 : gel(mat,1) = F.M;
1832 82826 : gel(mat,2) = F.G;
1833 :
1834 82826 : Tr = get_Tr(gel(nf,9), T, S->basden);
1835 82825 : gel(mat,6) = A = ZM_inv(Tr, &dA); /* dA T^-1, primitive */
1836 82825 : A = ZM_hnfmodid(A, dA);
1837 : /* CAVEAT: nf is not complete yet, but the fields needed for
1838 : * idealtwoelt, zk_scalar_or_multable and idealinv are present ! */
1839 82824 : MDI = idealtwoelt(nf, A);
1840 82825 : gel(MDI,2) = zk_scalar_or_multable(nf, gel(MDI,2));
1841 82823 : gel(mat,7) = MDI;
1842 82823 : if (is_pm1(S->index))
1843 : { /* principal ideal (T'), whose norm is |dK| */
1844 48301 : D = zk_scalar_or_multable(nf, ZX_deriv(T));
1845 48301 : if (typ(D) == t_MAT) D = ZM_hnfmod(D, absi_shallow(S->dK));
1846 : }
1847 : else
1848 : {
1849 34523 : GEN c = diviiexact(dA, gcoeff(A,1,1));
1850 34520 : D = idealHNF_inv_Z(nf, A); /* (A\cap Z) / A */
1851 34524 : if (!is_pm1(c)) D = ZM_Z_mul(D, c);
1852 : }
1853 82826 : gel(mat,3) = RM_round_maxrank(F.G);
1854 82825 : gel(mat,4) = Tr;
1855 82825 : gel(mat,5) = D;
1856 82825 : gel(mat,8) = shallowtrans(S->dKP? S->dKP: gel(absZ_factor(S->dK), 1));
1857 82825 : return nf;
1858 : }
1859 :
1860 : static GEN
1861 3101 : primes_certify(GEN dK, GEN dKP)
1862 : {
1863 3101 : long i, l = lg(dKP);
1864 3101 : GEN v, w, D = dK;
1865 3101 : v = vectrunc_init(l);
1866 3101 : w = vectrunc_init(l);
1867 9695 : for (i = 1; i < l; i++)
1868 : {
1869 6594 : GEN p = gel(dKP,i);
1870 6594 : vectrunc_append(isprime(p)? w: v, p);
1871 6594 : (void)Z_pvalrem(D, p, &D);
1872 : }
1873 3101 : if (!is_pm1(D))
1874 : {
1875 0 : if (signe(D) < 0) D = negi(D);
1876 0 : vectrunc_append(isprime(D)? w: v, D);
1877 : }
1878 3101 : return mkvec2(v,w);
1879 : }
1880 : GEN
1881 3101 : nfcertify(GEN nf)
1882 : {
1883 3101 : pari_sp av = avma;
1884 : GEN vw;
1885 3101 : nf = checknf(nf);
1886 3101 : vw = primes_certify(nf_get_disc(nf), nf_get_ramified_primes(nf));
1887 3101 : return gerepilecopy(av, gel(vw,1));
1888 : }
1889 :
1890 : /* set *pro to roots of S->T */
1891 : static GEN
1892 72863 : get_red_G(nfmaxord_t *S, GEN *pro)
1893 : {
1894 72863 : pari_sp av = avma;
1895 72863 : GEN G, u, u0 = NULL;
1896 72863 : long i, prec, n = degpol(S->T);
1897 : nffp_t F;
1898 :
1899 72863 : prec = nbits2prec(n+32);
1900 72863 : nffp_init(&F, S, prec);
1901 72863 : for (i=1; ; i++)
1902 : {
1903 72863 : F.prec = prec; make_M_G(&F, 0); G = F.G;
1904 72860 : if (u0) G = RgM_mul(G, u0);
1905 72860 : if (DEBUGLEVEL)
1906 0 : err_printf("get_red_G: starting LLL, prec = %ld (%ld + %ld)\n",
1907 0 : prec + F.extraprec, prec, F.extraprec);
1908 72860 : if ((u = lllfp(G, 0.99, LLL_KEEP_FIRST)))
1909 : {
1910 72862 : if (lg(u)-1 == n) break;
1911 : /* singular ==> loss of accuracy */
1912 0 : if (u0) u0 = gerepileupto(av, RgM_mul(u0,u));
1913 0 : else u0 = gerepilecopy(av, u);
1914 : }
1915 0 : prec = precdbl(prec) + nbits2extraprec(gexpo(u0));
1916 0 : F.ro = NULL;
1917 0 : if (DEBUGLEVEL) pari_warn(warnprec,"get_red_G", prec);
1918 : }
1919 72862 : if (u0) u = RgM_mul(u0,u);
1920 72862 : *pro = F.ro; return u;
1921 : }
1922 :
1923 : /* Compute an LLL-reduced basis for the integer basis of nf(T).
1924 : * set *pro = roots of x if computed [NULL if not computed] */
1925 : static void
1926 100620 : set_LLL_basis(nfmaxord_t *S, GEN *pro, long flag, double DELTA)
1927 : {
1928 100620 : GEN B = S->basis;
1929 100620 : long N = degpol(S->T);
1930 100620 : if (S->r1 < 0)
1931 : {
1932 17892 : S->r1 = ZX_sturm_irred(S->T);
1933 17892 : if (odd(N - S->r1)) pari_err_IRREDPOL("set_LLL_basis", S->T);
1934 : }
1935 100613 : if (!S->basden) S->basden = get_bas_den(B);
1936 100613 : *pro = NULL; if (flag & nf_NOLLL) return;
1937 100592 : if (S->r1 == N) {
1938 27729 : pari_sp av = avma;
1939 27729 : GEN u = ZM_lll(make_Tr(S), DELTA, LLL_GRAM|LLL_KEEP_FIRST|LLL_IM);
1940 27730 : B = gerepileupto(av, RgV_RgM_mul(B, u));
1941 : }
1942 : else
1943 72863 : B = RgV_RgM_mul(B, get_red_G(S, pro));
1944 100593 : S->basis = B;
1945 100593 : S->basden = get_bas_den(B);
1946 : }
1947 :
1948 : /* = 1 iff |a| > |b| or equality and a > 0 */
1949 : static int
1950 138775 : cmpii_polred(GEN a, GEN b)
1951 : {
1952 138775 : int fl = abscmpii(a, b);
1953 : long sa, sb;
1954 138775 : if (fl) return fl;
1955 115577 : sa = signe(a);
1956 115577 : sb = signe(b);
1957 115577 : if (sa == sb) return 0;
1958 763 : return sa == 1? 1: -1;
1959 : }
1960 : static int
1961 32599 : ZX_cmp(GEN x, GEN y)
1962 32599 : { return gen_cmp_RgX((void*)cmpii_polred, x, y); }
1963 : /* current best: ZX x of discriminant *dx, is ZX y better than x ?
1964 : * (if so update *dx); both x and y are monic */
1965 : static int
1966 54438 : ZX_is_better(GEN y, GEN x, GEN *dx)
1967 : {
1968 54438 : GEN d = ZX_disc(y);
1969 : int cmp;
1970 54438 : if (!*dx) *dx = ZX_disc(x);
1971 54438 : cmp = abscmpii(d, *dx);
1972 54438 : if (cmp < 0) { *dx = d; return 1; }
1973 49833 : return cmp? 0: (ZX_cmp(y, x) < 0);
1974 : }
1975 :
1976 : static void polredbest_aux(nfmaxord_t *S, GEN *pro, GEN *px, GEN *pdx, GEN *pa);
1977 : /* Seek a simpler, polynomial pol defining the same number field as
1978 : * x (assumed to be monic at this point) */
1979 : static GEN
1980 273 : nfpolred(nfmaxord_t *S, GEN *pro)
1981 : {
1982 273 : GEN x = S->T, dx, b, rev;
1983 273 : long n = degpol(x), v = varn(x);
1984 :
1985 273 : if (n == 1) {
1986 98 : S->T = pol_x(v);
1987 98 : *pro = NULL;
1988 98 : return scalarpol_shallow(negi(gel(x,2)), v);
1989 : }
1990 175 : polredbest_aux(S, pro, &x, &dx, &b);
1991 175 : if (x == S->T) return NULL; /* no improvement */
1992 140 : if (DEBUGLEVEL>1) err_printf("xbest = %Ps\n",x);
1993 :
1994 : /* update T */
1995 140 : rev = QXQ_reverse(b, S->T);
1996 140 : S->basis = QXV_QXQ_eval(S->basis, rev, x);
1997 140 : S->index = sqrti( diviiexact(dx,S->dK) );
1998 140 : S->basden = get_bas_den(S->basis);
1999 140 : S->dT = dx;
2000 140 : S->T = x;
2001 140 : *pro = NULL; /* reset */
2002 140 : return rev;
2003 : }
2004 :
2005 : /* Either nf type or ZX or [monic ZX, data], where data is either an integral
2006 : * basis (deprecated), or listP data (nfbasis input format) to specify
2007 : * a set of primes at with the basis order must be maximal.
2008 : * 1) nf type (or unrecognized): return t_VEC
2009 : * 2) ZX or [ZX, listP]: return t_POL
2010 : * 3) [ZX, order basis]: return 0 (deprecated)
2011 : * incorrect: return -1 */
2012 : static long
2013 84697 : nf_input_type(GEN x)
2014 : {
2015 84697 : GEN T, V, DKP = NULL;
2016 : long i, d, v;
2017 84697 : switch(typ(x))
2018 : {
2019 79286 : case t_POL: return t_POL;
2020 5411 : case t_VEC:
2021 5411 : switch(lg(x))
2022 : {
2023 1645 : case 4: DKP = gel(x,3);
2024 5383 : case 3: break;
2025 28 : default: return t_VEC; /* nf or incorrect */
2026 : }
2027 5383 : T = gel(x,1); V = gel(x,2);
2028 5383 : if (typ(T) != t_POL) return -1;
2029 5383 : switch(typ(V))
2030 : {
2031 224 : case t_INT: case t_MAT: return t_POL;
2032 5159 : case t_VEC: case t_COL:
2033 5159 : if (RgV_is_ZV(V)) return t_POL;
2034 2310 : break;
2035 0 : default: return -1;
2036 : }
2037 2310 : d = degpol(T); v = varn(T);
2038 2310 : if (d<1 || !RgX_is_ZX(T) || !isint1(gel(T,d+2)) || lg(V)-1!=d) return -1;
2039 14847 : for (i = 1; i <= d; i++)
2040 : { /* check integer basis */
2041 12558 : GEN c = gel(V,i);
2042 12558 : switch(typ(c))
2043 : {
2044 28 : case t_INT: break;
2045 12530 : case t_POL: if (varn(c) == v && RgX_is_QX(c) && degpol(c) < d) break;
2046 : /* fall through */
2047 14 : default: return -1;
2048 : }
2049 : }
2050 2289 : if (DKP && (typ(DKP) != t_VEC || !RgV_is_ZV(DKP))) return -1;
2051 2289 : return 0;
2052 : }
2053 0 : return t_VEC; /* nf or incorrect */
2054 : }
2055 :
2056 : /* cater for obsolete nf_PARTIALFACT flag */
2057 : static void
2058 3899 : nfinit_basic_partial(nfmaxord_t *S, GEN T)
2059 : {
2060 3899 : if (typ(T) == t_POL) { nfmaxord(S, mkvec2(T,utoipos(500000)), 0); }
2061 14 : else nfinit_basic(S, T);
2062 3899 : }
2063 : static void
2064 14063 : nfinit_basic_flag(nfmaxord_t *S, GEN x, long flag)
2065 : {
2066 14063 : if (flag & nf_PARTIALFACT)
2067 35 : nfinit_basic_partial(S, x);
2068 : else
2069 14028 : nfinit_basic(S, x);
2070 14056 : }
2071 :
2072 : /* true nf */
2073 : static GEN
2074 57622 : nf_basden(GEN nf)
2075 : {
2076 57622 : GEN zkD = nf_get_zkprimpart(nf), D = nf_get_zkden(nf);
2077 57622 : D = equali1(D)? NULL: const_vec(lg(zkD)-1, D);
2078 57621 : return mkvec2(zkD, D);
2079 : }
2080 : void
2081 84697 : nfinit_basic(nfmaxord_t *S, GEN T)
2082 : {
2083 84697 : switch (nf_input_type(T))
2084 : {
2085 82359 : case t_POL: nfmaxord(S, T, 0); return;
2086 28 : case t_VEC:
2087 : { /* nf, bnf, bnr */
2088 28 : GEN nf = checknf(T);
2089 28 : S->T = S->T0 = nf_get_pol(nf);
2090 28 : S->basis = nf_get_zk(nf); /* probably useless */
2091 28 : S->basden = nf_basden(nf);
2092 28 : S->index = nf_get_index(nf);
2093 28 : S->dK = nf_get_disc(nf);
2094 28 : S->dKP = nf_get_ramified_primes(nf);
2095 28 : S->dT = mulii(S->dK, sqri(S->index));
2096 28 : S->r1 = nf_get_r1(nf); break;
2097 : }
2098 2289 : case 0: /* monic integral polynomial + integer basis (+ ramified primes)*/
2099 2289 : S->T = S->T0 = gel(T,1);
2100 2289 : S->basis = gel(T,2);
2101 2289 : S->basden = NULL;
2102 2289 : S->index = NULL;
2103 2289 : S->dK = NULL;
2104 2289 : S->dKP = lg(T) == 4? gel(T,3): NULL;
2105 2289 : S->dT = NULL;
2106 2289 : S->r1 = -1; break;
2107 20 : default: /* -1 */
2108 20 : pari_err_TYPE("nfinit_basic", T);
2109 : }
2110 2317 : S->dTP = S->dTE = S->dKE = NULL;
2111 2317 : S->unscale = gen_1;
2112 : }
2113 :
2114 : GEN
2115 82824 : nfinit_complete(nfmaxord_t *S, long flag, long prec)
2116 : {
2117 82824 : GEN nf, unscale = S->unscale, rev = NULL;
2118 :
2119 82824 : if (!ZX_is_irred(S->T)) pari_err_IRREDPOL("nfinit",S->T);
2120 82827 : if (!(flag & nf_RED) && !ZX_is_monic(S->T0))
2121 : {
2122 42 : pari_warn(warner,"nonmonic polynomial. Result of the form [nf,c]");
2123 42 : flag |= nf_RED | nf_ORIG;
2124 : }
2125 82827 : if (!(flag & nf_RED) && !isint1(unscale))
2126 : { /* implies lc(x0) = 1 and L := 1/unscale is integral */
2127 4207 : long d = degpol(S->T0);
2128 4207 : GEN L = ginv(unscale); /* x = L^(-deg(x)) x0(L X) */
2129 4207 : GEN f = powiu(L, (d*(d-1)) >> 1);
2130 4207 : S->T = S->T0; /* restore original user-supplied x0, unscale data */
2131 4207 : S->unscale = gen_1;
2132 4207 : S->dT = gmul(S->dT, sqri(f));
2133 4207 : S->basis = RgXV_unscale(S->basis, unscale);
2134 4207 : S->index = gmul(S->index, f);
2135 : }
2136 82827 : nfmaxord_complete(S); /* more expensive after set_LLL_basis */
2137 82826 : if (flag & nf_RED)
2138 : {
2139 : GEN ro;
2140 : /* lie to polred: more efficient to update *after* modreverse, than to
2141 : * unscale in the polred subsystem */
2142 273 : S->unscale = gen_1;
2143 273 : rev = nfpolred(S, &ro);
2144 273 : nf = nfmaxord_to_nf(S, ro, prec);
2145 273 : S->unscale = unscale; /* restore */
2146 : }
2147 : else
2148 : {
2149 82553 : GEN ro; set_LLL_basis(S, &ro, flag, 0.99);
2150 82552 : nf = nfmaxord_to_nf(S, ro, prec);
2151 : }
2152 82825 : if (flag & nf_ORIG)
2153 : {
2154 77 : if (!rev)
2155 : { /* no improvement */
2156 28 : long v = varn(S->T);
2157 28 : rev = degpol(S->T) == 1? pol_0(v): pol_x(v);
2158 : }
2159 77 : if (!isint1(unscale)) rev = RgX_Rg_div(rev, unscale);
2160 77 : nf = mkvec2(nf, mkpolmod(rev, S->T));
2161 : }
2162 82825 : return nf;
2163 : }
2164 : /* Initialize the number field defined by the polynomial x (in variable v)
2165 : * flag & nf_RED: try a polred first.
2166 : * flag & nf_ORIG: return [nfinit(x), Mod(a,red)], where
2167 : * Mod(a,red) = Mod(v,x) (i.e return the base change). */
2168 : GEN
2169 9656 : nfinit0(GEN x, long flag,long prec)
2170 : {
2171 9656 : const pari_sp av = avma;
2172 : nfmaxord_t S;
2173 9656 : if (flag < 0 || flag > 7) pari_err_FLAG("nfinit");
2174 9656 : if (checkrnf_i(x)) return rnf_build_nfabs(x, prec);
2175 9649 : nfinit_basic(&S, x);
2176 9628 : return gerepilecopy(av, nfinit_complete(&S, flag, prec));
2177 : }
2178 : GEN
2179 182 : nfinitred(GEN x, long prec) { return nfinit0(x, nf_RED, prec); }
2180 : GEN
2181 0 : nfinitred2(GEN x, long prec) { return nfinit0(x, nf_RED|nf_ORIG, prec); }
2182 : GEN
2183 6583 : nfinit(GEN x, long prec) { return nfinit0(x, 0, prec); }
2184 :
2185 : /* assume x a bnr/bnf/nf */
2186 : long
2187 901760 : nf_get_prec(GEN x)
2188 : {
2189 901760 : GEN nf = checknf(x), ro = nf_get_roots(nf);
2190 901770 : return (typ(ro)==t_VEC)? precision(gel(ro,1)): DEFAULTPREC;
2191 : }
2192 :
2193 : /* true nf */
2194 : GEN
2195 57594 : nfnewprec_shallow(GEN nf, long prec)
2196 : {
2197 57594 : GEN m, NF = leafcopy(nf);
2198 : nffp_t F;
2199 :
2200 57594 : F.T = nf_get_pol(nf);
2201 57594 : F.ro = NULL;
2202 57594 : F.r1 = nf_get_r1(nf);
2203 57594 : F.basden = nf_basden(nf);
2204 57592 : F.extraprec = -1;
2205 57592 : F.prec = prec; make_M_G(&F, 0);
2206 57594 : gel(NF,5) = m = leafcopy(gel(NF,5));
2207 57594 : gel(m,1) = F.M;
2208 57594 : gel(m,2) = F.G;
2209 57594 : gel(NF,6) = F.ro; return NF;
2210 : }
2211 :
2212 : GEN
2213 154 : nfnewprec(GEN nf, long prec)
2214 : {
2215 : GEN z;
2216 154 : switch(nftyp(nf))
2217 : {
2218 49 : default: pari_err_TYPE("nfnewprec", nf);
2219 14 : case typ_BNF: z = bnfnewprec(nf,prec); break;
2220 7 : case typ_BNR: z = bnrnewprec(nf,prec); break;
2221 84 : case typ_NF: {
2222 84 : pari_sp av = avma;
2223 84 : z = gerepilecopy(av, nfnewprec_shallow(checknf(nf), prec));
2224 84 : break;
2225 : }
2226 : }
2227 105 : return z;
2228 : }
2229 :
2230 : /********************************************************************/
2231 : /** **/
2232 : /** POLRED **/
2233 : /** **/
2234 : /********************************************************************/
2235 : GEN
2236 0 : embednorm_T2(GEN x, long r1)
2237 : {
2238 0 : pari_sp av = avma;
2239 0 : GEN p = RgV_sumpart(x, r1);
2240 0 : GEN q = RgV_sumpart2(x,r1+1, lg(x)-1);
2241 0 : if (q != gen_0) p = gadd(p, gmul2n(q,1));
2242 0 : return avma == av? gcopy(p): gerepileupto(av, p);
2243 : }
2244 :
2245 : /* simplified version of gnorm for scalar, noncomplex inputs, without GC */
2246 : static GEN
2247 164458 : real_norm(GEN x)
2248 : {
2249 164458 : switch(typ(x))
2250 : {
2251 0 : case t_INT: return sqri(x);
2252 164458 : case t_REAL: return sqrr(x);
2253 0 : case t_FRAC: return sqrfrac(x);
2254 : }
2255 0 : pari_err_TYPE("real_norm", x);
2256 : return NULL;/*LCOV_EXCL_LINE*/
2257 : }
2258 : /* simplified version of gnorm, without GC */
2259 : static GEN
2260 23856527 : complex_norm(GEN x)
2261 : {
2262 23856527 : return typ(x) == t_COMPLEX? cxnorm(x): real_norm(x);
2263 : }
2264 : /* return T2(x), argument r1 needed in case x has components whose type
2265 : * is unexpected, e.g. all of them t_INT for embed(gen_1) */
2266 : GEN
2267 71036 : embed_T2(GEN x, long r1)
2268 : {
2269 71036 : pari_sp av = avma;
2270 71036 : long i, l = lg(x);
2271 71036 : GEN c, s = NULL, t = NULL;
2272 71036 : if (typ(gel(x,1)) == t_INT) return muliu(gel(x,1), 2*(l-1)-r1);
2273 235490 : for (i = 1; i <= r1; i++)
2274 : {
2275 164457 : c = real_norm(gel(x,i));
2276 164467 : s = s? gadd(s, c): c;
2277 : }
2278 202755 : for (; i < l; i++)
2279 : {
2280 131724 : c = complex_norm(gel(x,i));
2281 131721 : t = t? gadd(t, c): c;
2282 : }
2283 71031 : if (t) { t = gmul2n(t,1); s = s? gadd(s,t): t; }
2284 71032 : return gerepileupto(av, s);
2285 : }
2286 : /* return N(x) */
2287 : GEN
2288 30164709 : embed_norm(GEN x, long r1)
2289 : {
2290 30164709 : pari_sp av = avma;
2291 30164709 : long i, l = lg(x);
2292 30164709 : GEN c, s = NULL, t = NULL;
2293 30164709 : if (typ(gel(x,1)) == t_INT) return powiu(gel(x,1), 2*(l-1)-r1);
2294 80246012 : for (i = 1; i <= r1; i++)
2295 : {
2296 50207506 : c = gel(x,i);
2297 50207506 : s = s? gmul(s, c): c;
2298 : }
2299 53763053 : for (; i < l; i++)
2300 : {
2301 23724820 : c = complex_norm(gel(x,i));
2302 23723887 : t = t? gmul(t, c): c;
2303 : }
2304 30038233 : if (t) s = s? gmul(s,t): t;
2305 30038200 : return gerepileupto(av, s);
2306 : }
2307 :
2308 : typedef struct {
2309 : long r1, v, prec;
2310 : GEN ZKembed; /* embeddings of fincke-pohst-reduced Zk basis */
2311 : GEN u; /* matrix giving fincke-pohst-reduced Zk basis */
2312 : GEN M; /* embeddings of initial (LLL-reduced) Zk basis */
2313 : GEN bound; /* T2 norm of the polynomial defining nf */
2314 : long expo_best_disc; /* expo(disc(x)), best generator so far */
2315 : } CG_data;
2316 :
2317 : /* characteristic pol of x (given by embeddings) */
2318 : static GEN
2319 259562 : get_pol(CG_data *d, GEN x)
2320 : {
2321 : long e;
2322 259562 : GEN g = grndtoi(roots_to_pol_r1(x, d->v, d->r1), &e);
2323 259560 : return (e > -5)? NULL: g;
2324 : }
2325 :
2326 : /* characteristic pol of x (given as vector on (w_i)) */
2327 : static GEN
2328 100217 : get_polchar(CG_data *d, GEN x)
2329 100217 : { return get_pol(d, RgM_RgC_mul(d->ZKembed,x)); }
2330 :
2331 : /* Choose a canonical polynomial in the pair { Pmin_a, Pmin_{-a} }, i.e.
2332 : * { z(X), (-1)^(deg z) z(-Z) } and keeping the smallest wrt cmpii_polred
2333 : * Either leave z alone (return 1) or set z <- (-1)^n z(-X). In place. */
2334 : int
2335 95836 : ZX_canon_neg(GEN z)
2336 : {
2337 : long i, s;
2338 201395 : for (i = lg(z)-2; i >= 2; i -= 2)
2339 : { /* examine the odd (resp. even) part of z if deg(z) even (resp. odd). */
2340 193842 : s = signe(gel(z,i));
2341 193842 : if (!s) continue;
2342 : /* non trivial */
2343 88283 : if (s < 0) break; /* z(X) < (-1)^n z(-X) */
2344 :
2345 215988 : for (; i>=2; i-=2) gel(z,i) = negi(gel(z,i));
2346 37565 : return 1;
2347 : }
2348 58271 : return 0;
2349 : }
2350 : /* return a defining polynomial for Q(alpha), v = embeddings of alpha.
2351 : * Return NULL on failure: discriminant too large or non primitive */
2352 : static GEN
2353 135953 : try_polmin(CG_data *d, nfmaxord_t *S, GEN v, long flag, GEN *ai)
2354 : {
2355 135953 : const long best = flag & nf_ABSOLUTE;
2356 : long ed;
2357 135953 : pari_sp av = avma;
2358 : GEN g;
2359 135953 : if (best)
2360 : {
2361 135078 : ed = expo(embed_disc(v, d->r1, LOWDEFAULTPREC));
2362 135077 : set_avma(av); if (d->expo_best_disc < ed) return NULL;
2363 : }
2364 : else
2365 875 : ed = 0;
2366 79561 : g = get_pol(d, v);
2367 : /* accuracy too low, compute algebraically */
2368 79563 : if (!g) { set_avma(av); g = ZXQ_charpoly(*ai, S->T, varn(S->T)); }
2369 79563 : g = ZX_radical(g);
2370 79566 : if (best && degpol(g) != degpol(S->T)) return gc_NULL(av);
2371 33549 : g = gerepilecopy(av, g);
2372 33548 : d->expo_best_disc = ed;
2373 33548 : if (flag & nf_ORIG)
2374 : {
2375 29637 : if (ZX_canon_neg(g)) *ai = RgX_neg(*ai);
2376 29637 : if (!isint1(S->unscale)) *ai = RgX_unscale(*ai, S->unscale);
2377 : }
2378 : else
2379 3911 : (void)ZX_canon_neg(g);
2380 33548 : if (DEBUGLEVEL>3) err_printf("polred: generator %Ps\n", g);
2381 33548 : return g;
2382 : }
2383 :
2384 : /* does x generate the correct field ? */
2385 : static GEN
2386 100217 : chk_gen(void *data, GEN x)
2387 : {
2388 100217 : pari_sp av = avma, av1;
2389 100217 : GEN h, g = get_polchar((CG_data*)data,x);
2390 100218 : if (!g) pari_err_PREC("chk_gen");
2391 100218 : av1 = avma;
2392 100218 : h = ZX_gcd(g, ZX_deriv(g));
2393 100218 : if (degpol(h)) return gc_NULL(av);
2394 59243 : if (DEBUGLEVEL>3) err_printf(" generator: %Ps\n",g);
2395 59243 : set_avma(av1); return gerepileupto(av, g);
2396 : }
2397 :
2398 : static long
2399 32039 : chk_gen_prec(long N, long bit)
2400 32039 : { return prec_fix(nbits2prec(10 + (long)log2((double)N) + bit)); }
2401 :
2402 : /* v = [P,A] two vectors (of ZX and ZV resp.) of same length; remove duplicate
2403 : * polynomials in P, updating A, in place. Among elements having the same
2404 : * characteristic pol, choose the smallest according to ZV_abscmp */
2405 : static void
2406 13825 : remove_duplicates(GEN v)
2407 : {
2408 13825 : GEN x, a, P = gel(v,1), A = gel(v,2);
2409 13825 : long k, i, l = lg(P);
2410 13825 : pari_sp av = avma;
2411 :
2412 13825 : if (l < 2) return;
2413 13825 : (void)sort_factor_pol(mkvec2(P, A), cmpii);
2414 13825 : x = gel(P,1); a = gel(A,1);
2415 58011 : for (k=1,i=2; i<l; i++)
2416 44186 : if (ZX_equal(gel(P,i), x))
2417 : {
2418 21666 : if (ZV_abscmp(gel(A,i), a) < 0) a = gel(A,i);
2419 : }
2420 : else
2421 : {
2422 22520 : gel(A,k) = a;
2423 22520 : gel(P,k) = x;
2424 22520 : k++;
2425 22520 : x = gel(P,i); a = gel(A,i);
2426 : }
2427 13825 : l = k+1;
2428 13825 : gel(A,k) = a; setlg(A,l);
2429 13825 : gel(P,k) = x; setlg(P,l); set_avma(av);
2430 : }
2431 :
2432 : static void
2433 18067 : polred_init(nfmaxord_t *S, nffp_t *F, CG_data *d)
2434 : {
2435 18067 : long e, prec, n = degpol(S->T);
2436 : double log2rho;
2437 : GEN ro;
2438 18067 : set_LLL_basis(S, &ro, 0, 0.9999);
2439 : /* || polchar ||_oo < 2^e ~ 2 (n * rho)^n, rho = max modulus of root */
2440 18060 : log2rho = ro ? (double)gexpo(ro): fujiwara_bound(S->T);
2441 18060 : e = n * (long)(log2rho + log2((double)n)) + 1;
2442 18060 : if (e < 0) e = 0; /* can occur if n = 1 */
2443 18060 : prec = chk_gen_prec(n, e);
2444 18060 : nffp_init(F,S,prec);
2445 18060 : F->ro = ro;
2446 18060 : make_M_G(F, 1);
2447 :
2448 18060 : d->v = varn(S->T);
2449 18060 : d->expo_best_disc = -1;
2450 18060 : d->ZKembed = NULL;
2451 18060 : d->M = NULL;
2452 18060 : d->u = NULL;
2453 18060 : d->r1= S->r1;
2454 18060 : }
2455 : static GEN
2456 13986 : findmindisc(GEN y)
2457 : {
2458 13986 : GEN x = gel(y,1), dx = NULL;
2459 13986 : long i, l = lg(y);
2460 35043 : for (i = 2; i < l; i++)
2461 : {
2462 21057 : GEN yi = gel(y,i);
2463 21057 : if (ZX_is_better(yi,x,&dx)) x = yi;
2464 : }
2465 13986 : return x;
2466 : }
2467 : /* filter [y,b] from polred_aux: keep a single polynomial of degree n in y
2468 : * [ the best wrt discriminant ordering ], but keep all imprimitive
2469 : * polynomials */
2470 : static void
2471 4081 : filter(GEN y, GEN b, long n)
2472 : {
2473 : GEN x, a, dx;
2474 4081 : long i, k = 1, l = lg(y);
2475 4081 : a = x = dx = NULL;
2476 37693 : for (i = 1; i < l; i++)
2477 : {
2478 33612 : GEN yi = gel(y,i), ai = gel(b,i);
2479 33612 : if (degpol(yi) == n)
2480 : {
2481 33444 : pari_sp av = avma;
2482 33444 : if (!dx) dx = ZX_disc(yi);
2483 29538 : else if (!ZX_is_better(yi,x,&dx)) { set_avma(av); continue; }
2484 5850 : x = yi; a = ai; continue;
2485 : }
2486 168 : gel(y,k) = yi;
2487 168 : gel(b,k) = ai; k++;
2488 : }
2489 4081 : if (dx)
2490 : {
2491 3906 : gel(y,k) = x;
2492 3906 : gel(b,k) = a; k++;
2493 : }
2494 4081 : setlg(y, k);
2495 4081 : setlg(b, k);
2496 4081 : }
2497 :
2498 : static GEN
2499 4116 : polred_aux(nfmaxord_t *S, GEN *pro, long flag)
2500 : { /* only keep polynomials of max degree and best discriminant */
2501 4116 : const long best = flag & nf_ABSOLUTE;
2502 4116 : const long orig = flag & nf_ORIG;
2503 4116 : GEN M, b, y, x = S->T;
2504 4116 : long maxi, i, j, k, v = varn(x), n = lg(S->basis)-1;
2505 : nffp_t F;
2506 : CG_data d;
2507 :
2508 4116 : if (n == 1)
2509 : {
2510 28 : if (!best)
2511 : {
2512 14 : GEN X = pol_x(v);
2513 14 : return orig? mkmat2(mkcol(X),mkcol(gen_1)): mkvec(X);
2514 : }
2515 : else
2516 14 : return orig? trivial_fact(): cgetg(1,t_VEC);
2517 : }
2518 :
2519 4088 : polred_init(S, &F, &d);
2520 4081 : if (pro) *pro = F.ro;
2521 4081 : M = F.M;
2522 4081 : if (best)
2523 : {
2524 4018 : if (!S->dT) S->dT = ZX_disc(S->T);
2525 4018 : d.expo_best_disc = expi(S->dT);
2526 : }
2527 :
2528 : /* n + 2 sum_{1 <= i <= n} n-i = n + n(n-1) = n*n */
2529 4081 : y = cgetg(n*n + 1, t_VEC);
2530 4081 : b = cgetg(n*n + 1, t_COL);
2531 4081 : k = 1;
2532 4081 : if (!best) { gel(y,1) = pol_x(v); gel(b,1) = gen_0; k++; }
2533 26894 : for (i = 2; i <= n; i++)
2534 : {
2535 : GEN ch, ai;
2536 22813 : ai = gel(S->basis,i);
2537 22813 : ch = try_polmin(&d, S, gel(M,i), flag, &ai);
2538 22813 : if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
2539 : }
2540 4081 : maxi = minss(n, 3);
2541 15946 : for (i = 1; i <= maxi; i++)
2542 68437 : for (j = i+1; j <= n; j++)
2543 : {
2544 : GEN ch, ai, v;
2545 56572 : ai = gadd(gel(S->basis,i), gel(S->basis,j));
2546 56572 : v = RgV_add(gel(M,i), gel(M,j));
2547 : /* defining polynomial for Q(w_i+w_j) */
2548 56570 : ch = try_polmin(&d, S, v, flag, &ai);
2549 56572 : if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
2550 :
2551 56572 : ai = gsub(gel(S->basis,i), gel(S->basis,j));
2552 56574 : v = RgV_sub(gel(M,i), gel(M,j));
2553 : /* defining polynomial for Q(w_i-w_j) */
2554 56574 : ch = try_polmin(&d, S, v, flag, &ai);
2555 56572 : if (ch) { gel(y,k) = ch; gel(b,k) = ai; k++; }
2556 : }
2557 4081 : setlg(y, k);
2558 4081 : setlg(b, k); filter(y, b, n);
2559 4081 : if (!orig) return gen_sort_uniq(y, (void*)cmpii, &gen_cmp_RgX);
2560 3402 : settyp(y, t_COL);
2561 3402 : (void)sort_factor_pol(mkmat2(y, b), cmpii);
2562 3402 : return mkmat2(b, y);
2563 : }
2564 :
2565 : /* FIXME: obsolete */
2566 : static GEN
2567 84 : Polred(GEN x, long flag, GEN fa)
2568 : {
2569 84 : pari_sp av = avma;
2570 : nfmaxord_t S;
2571 84 : if (fa)
2572 14 : nfinit_basic(&S, mkvec2(x,fa));
2573 : else
2574 70 : nfinit_basic_flag(&S, x, flag);
2575 77 : return gerepilecopy(av, polred_aux(&S, NULL, flag));
2576 : }
2577 :
2578 : /* finds "best" polynomial in polred_aux list, defaulting to S->T if none of
2579 : * them is primitive. *px is the ZX, characteristic polynomial of Mod(*pb,S->T),
2580 : * *pdx its discriminant if pdx != NULL. Set *pro = polroots(S->T) */
2581 : static void
2582 4039 : polredbest_aux(nfmaxord_t *S, GEN *pro, GEN *px, GEN *pdx, GEN *pb)
2583 : {
2584 4039 : GEN y, dx, x = S->T; /* default value */
2585 : long i, l;
2586 4039 : y = polred_aux(S, pro, pb? nf_ORIG|nf_ABSOLUTE: nf_ABSOLUTE);
2587 4032 : dx = S->dT;
2588 4032 : if (pb)
2589 : {
2590 3388 : GEN a, b = deg1pol_shallow(S->unscale, gen_0, varn(x));
2591 3388 : a = gel(y,1); l = lg(a);
2592 3388 : y = gel(y,2);
2593 6594 : for (i=1; i<l; i++)
2594 : {
2595 3206 : GEN yi = gel(y,i);
2596 3206 : pari_sp av = avma;
2597 3206 : if (ZX_is_better(yi,x,&dx)) { x = yi; b = gel(a,i); } else set_avma(av);
2598 : }
2599 3388 : *pb = b;
2600 : }
2601 : else
2602 : {
2603 644 : l = lg(y);
2604 1281 : for (i=1; i<l; i++)
2605 : {
2606 637 : GEN yi = gel(y,i);
2607 637 : pari_sp av = avma;
2608 637 : if (ZX_is_better(yi,x,&dx)) x = yi; else set_avma(av);
2609 : }
2610 : }
2611 4032 : if (pdx) { if (!dx) dx = ZX_disc(x); *pdx = dx; }
2612 4032 : *px = x;
2613 4032 : }
2614 : static GEN
2615 3864 : polredbest_i(GEN T0, long flag)
2616 : {
2617 3864 : GEN T = T0, a;
2618 : nfmaxord_t S;
2619 3864 : nfinit_basic_partial(&S, T);
2620 3864 : polredbest_aux(&S, NULL, &T, NULL, flag? &a: NULL);
2621 3857 : if (flag == 2)
2622 350 : T = mkvec2(T, a);
2623 3507 : else if (flag == 1)
2624 : {
2625 2863 : GEN b = (T0 == T)? pol_x(varn(T)): QXQ_reverse(a, T0);
2626 : /* charpoly(Mod(a,T0)) = T; charpoly(Mod(b,T)) = S.x */
2627 2863 : if (degpol(T) == 1) b = grem(b,T);
2628 2863 : T = mkvec2(T, mkpolmod(b,T));
2629 : }
2630 3857 : return T;
2631 : }
2632 : GEN
2633 3507 : polredbest(GEN T, long flag)
2634 : {
2635 3507 : pari_sp av = avma;
2636 3507 : if (flag < 0 || flag > 1) pari_err_FLAG("polredbest");
2637 3507 : 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 22175 : set_mulid(GEN V, GEN M, GEN Mi, long r1, long r2, long N, long k)
2703 : {
2704 22175 : GEN v, Mk = cgetg(N+1, t_MAT);
2705 : long i, e;
2706 35097 : for (i = 1; i < k; i++) gel(Mk,i) = gmael(V, i, k);
2707 146179 : for ( ; i <=N; i++)
2708 : {
2709 124004 : v = vecmul(gel(M,k), gel(M,i));
2710 124008 : v = RgM_RgC_mul(Mi, split_realimag(v, r1, r2));
2711 124005 : gel(Mk,i) = grndtoi(v, &e);
2712 124004 : if (e > -5) return NULL;
2713 : }
2714 22175 : gel(V,k) = Mk; return Mk;
2715 : }
2716 :
2717 : static GEN
2718 12446 : ZM_image_shallow(GEN M, long *pr)
2719 : {
2720 : long j, k, r;
2721 12446 : GEN y, d = ZM_pivots(M, &k);
2722 12446 : r = lg(M)-1 - k;
2723 12446 : y = cgetg(r+1,t_MAT);
2724 55709 : for (j=k=1; j<=r; k++)
2725 43263 : if (d[k]) gel(y,j++) = gel(M,k);
2726 12446 : *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 13987 : chk_gen_init(FP_chk_fun *chk, GEN R, GEN U)
2733 : {
2734 13987 : CG_data *d = (CG_data*)chk->data;
2735 : GEN P, V, D, inv, bound, S, M;
2736 13987 : long N = lg(U)-1, r1 = d->r1, r2 = (N-r1)>>1;
2737 13987 : long i, j, prec, firstprim = 0, skipfirst = 0;
2738 : pari_sp av;
2739 :
2740 13987 : d->u = U;
2741 13987 : d->ZKembed = M = RgM_mul(d->M, U);
2742 :
2743 13987 : av = avma; bound = d->bound;
2744 13987 : D = cgetg(N+1, t_VECSMALL);
2745 93726 : for (i = 1; i <= N; i++)
2746 : {
2747 79747 : pari_sp av2 = avma;
2748 79747 : P = get_pol(d, gel(M,i));
2749 79746 : if (!P) pari_err_PREC("chk_gen_init");
2750 79738 : P = gerepilecopy(av2, ZX_radical(P));
2751 79739 : D[i] = degpol(P);
2752 79739 : if (D[i] == N)
2753 : { /* primitive element */
2754 56777 : GEN B = embed_T2(gel(M,i), r1);
2755 56777 : if (!firstprim) firstprim = i; /* index of first primitive element */
2756 56777 : if (DEBUGLEVEL>2) err_printf("chk_gen_init: generator %Ps\n",P);
2757 56777 : if (gcmp(B,bound) < 0) bound = gerepileuptoleaf(av2, B);
2758 : }
2759 : else
2760 : {
2761 22962 : if (DEBUGLEVEL>2) err_printf("chk_gen_init: subfield %Ps\n",P);
2762 22962 : if (firstprim)
2763 : { /* cycle basis vectors so that primitive elements come last */
2764 2812 : GEN u = d->u, e = M;
2765 2812 : GEN te = gel(e,i), tu = gel(u,i), tR = gel(R,i);
2766 2812 : long tS = D[i];
2767 8511 : for (j = i; j > firstprim; j--)
2768 : {
2769 5699 : u[j] = u[j-1];
2770 5699 : e[j] = e[j-1];
2771 5699 : R[j] = R[j-1];
2772 5699 : D[j] = D[j-1];
2773 : }
2774 2812 : gel(u,firstprim) = tu;
2775 2812 : gel(e,firstprim) = te;
2776 2812 : gel(R,firstprim) = tR;
2777 2812 : D[firstprim] = tS; firstprim++;
2778 : }
2779 : }
2780 : }
2781 13979 : 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 35 : P = get_pol(d, e); if (!P) pari_err_PREC( "chk_gen_init");
2794 35 : if (!ZX_is_squarefree(P)) continue;
2795 35 : if (DEBUGLEVEL>2) err_printf("chk_gen_init: generator %Ps\n",P);
2796 35 : bound = B ;
2797 : }
2798 : }
2799 :
2800 13979 : if (firstprim != 1)
2801 : {
2802 13748 : inv = ginv( split_realimag(M, r1, r2) ); /*TODO: use QR?*/
2803 13748 : V = gel(inv,1);
2804 53487 : for (i = 2; i <= r1+r2; i++) V = gadd(V, gel(inv,i));
2805 : /* V corresponds to 1_Z */
2806 13748 : V = grndtoi(V, &j);
2807 13748 : if (j > -5) pari_err_BUG("precision too low in chk_gen_init");
2808 13748 : S = mkmat(V); /* 1 */
2809 :
2810 13748 : V = cgetg(N+1, t_VEC);
2811 35314 : 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 35314 : long j, k, h, rkM, dP = D[i];
2815 :
2816 35314 : if (dP == N) break; /* primitive */
2817 22175 : Mx = set_mulid(V, M, inv, r1, r2, N, i);
2818 22175 : if (!Mx) break; /* prec. problem. Stop */
2819 22175 : if (dP == 1) continue;
2820 9128 : rkM = lg(S)-1;
2821 9128 : M2 = cgetg(N+1, t_MAT); /* we will add to S the elts of M2 */
2822 9128 : gel(M2,1) = col_ei(N, i); /* nf.zk[i] */
2823 9128 : k = 2;
2824 18252 : for (h = 1; h < dP; h++)
2825 : {
2826 : long r; /* add to M2 the elts of S * nf.zk[i] */
2827 38745 : for (j = 1; j <= rkM; j++) gel(M2,k++) = ZM_ZC_mul(Mx, gel(S,j));
2828 12446 : setlg(M2, k); k = 1;
2829 12446 : S = ZM_image_shallow(shallowconcat(S,M2), &r);
2830 13055 : if (r == rkM) break;
2831 9733 : if (r > rkM)
2832 : {
2833 9733 : rkM = r;
2834 9733 : if (rkM == N) break;
2835 : }
2836 : }
2837 9128 : 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 13979 : chk->skipfirst = skipfirst;
2843 13979 : 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 13979 : bound = gerepileuptoleaf(av, bound);
2847 13979 : prec = chk_gen_prec(N, (gexpo(bound)*N)/2);
2848 13979 : if (DEBUGLEVEL)
2849 0 : err_printf("chk_gen_init: new prec = %ld (initially %ld)\n", prec, d->prec);
2850 13979 : if (prec > d->prec) pari_err_BUG("polredabs (precision problem)");
2851 13979 : if (prec < d->prec) d->ZKembed = gprec_w(M, prec);
2852 13979 : return bound;
2853 : }
2854 :
2855 : static GEN
2856 13993 : polredabs_i(GEN x, nfmaxord_t *S, GEN *u, long flag)
2857 : {
2858 13993 : 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 13993 : nfinit_basic_flag(S, x, flag);
2865 13993 : x = S->T0;
2866 13993 : 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 13979 : chk.data = (void*)&d;
2874 13979 : polred_init(S, &F, &d);
2875 13979 : d.bound = embed_T2(F.ro, d.r1);
2876 13979 : if (realprec(d.bound) > F.prec) d.bound = rtor(d.bound, F.prec);
2877 : for (;;)
2878 20 : {
2879 13999 : GEN R = R_from_QR(F.G, F.prec);
2880 13999 : if (R)
2881 : {
2882 13987 : d.prec = F.prec;
2883 13987 : d.M = F.M;
2884 13987 : v = fincke_pohst(mkvec(R),NULL,-1, 0, &chk);
2885 13987 : 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 13979 : y = gel(v,1);
2893 13979 : a = gel(v,2); l = lg(a);
2894 72802 : for (i = 1; i < l; i++) /* normalize wrt z -> -z */
2895 58823 : if (ZX_canon_neg(gel(y,i)) && (flag & (nf_ORIG|nf_RAW)))
2896 490 : gel(a,i) = ZC_neg(gel(a,i));
2897 13979 : *u = d.u; return v;
2898 : }
2899 :
2900 : GEN
2901 13825 : polredabs0(GEN x, long flag)
2902 : {
2903 13825 : pari_sp av = avma;
2904 : GEN Y, A, u, v;
2905 : nfmaxord_t S;
2906 : long i, l;
2907 :
2908 13825 : v = polredabs_i(x, &S, &u, flag);
2909 13825 : remove_duplicates(v);
2910 13825 : Y = gel(v,1);
2911 13825 : A = gel(v,2);
2912 13825 : l = lg(A); if (l == 1) pari_err_BUG("polredabs (missing vector)");
2913 13825 : if (DEBUGLEVEL) err_printf("Found %ld minimal polynomials.\n",l-1);
2914 13825 : if (!(flag & nf_ALL))
2915 : {
2916 13818 : GEN y = findmindisc(Y);
2917 28099 : for (i = 1; i < l; i++)
2918 28099 : if (ZX_equal(gel(Y,i), y)) break;
2919 13818 : Y = mkvec(gel(Y,i));
2920 13818 : A = mkvec(gel(A,i)); l = 2;
2921 : }
2922 14091 : 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 13825 : 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 532 : rnfpolred_i(GEN nf, GEN R, long flag, long best)
2949 : {
2950 532 : const char *f = best? "rnfpolredbest": "rnfpolredabs";
2951 532 : const long abs = ((flag & nf_ORIG) && (flag & nf_ABSOLUTE));
2952 532 : GEN listP = NULL, red, pol, A, P, T, rnfeq;
2953 532 : pari_sp av = avma;
2954 :
2955 532 : 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 532 : if (typ(R) != t_POL) pari_err_TYPE(f,R);
2961 532 : nf = checknf(nf);
2962 532 : T = nf_get_pol(nf);
2963 532 : R = RgX_nffix(f, T, R, 0);
2964 532 : 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 168 : if (DEBUGLEVEL>1) timer_start(&ti);
2981 168 : rnf = rnfinit(nf, R);
2982 168 : rnfeq = rnf_get_map(rnf);
2983 168 : pol = rnf_zkabs(rnf);
2984 168 : if (DEBUGLEVEL>1) timer_printf(&ti, "absolute basis");
2985 168 : v = polredabs_i(pol, &S, &u, nf_ORIG);
2986 168 : pol = gel(pol,1);
2987 168 : y = gel(v,1); P = findmindisc(y);
2988 168 : a = gel(v,2);
2989 168 : l = lg(y); A = cgetg(l, t_VEC);
2990 994 : for (i = j = 1; i < l; i++)
2991 826 : if (ZX_equal(gel(y,i),P))
2992 : {
2993 714 : GEN t = gel(a,i);
2994 714 : if (u) t = RgV_RgC_mul(S.basis, ZM_ZC_mul(u,t));
2995 714 : gel(A,j++) = t;
2996 : }
2997 168 : setlg(A,j); /* mod(A[i], pol) are all roots of P in Q[X]/(pol) */
2998 : }
2999 532 : if (DEBUGLEVEL>1) err_printf("reduced absolute generator: %Ps\n",P);
3000 532 : 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 518 : 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 161 : long i, l = lg(A), v = varn(R);
3020 161 : GEN besta = NULL;
3021 833 : for (i = 1; i < l; i++)
3022 : {
3023 672 : GEN a = eltabstorel_lift(rnfeq, gel(A,i));
3024 672 : GEN p = lift_if_rational( RgXQ_charpoly(a, R, v) );
3025 672 : if (i == 1 || cmp_universal(p, P) < 0) { P = p; besta = a; }
3026 : }
3027 161 : A = besta;
3028 : }
3029 518 : if (flag & nf_ORIG) P = mkvec2(P, mkpolmod(RgXQ_reverse(A,R),P));
3030 518 : return gerepilecopy(av, P);
3031 : }
3032 : GEN
3033 175 : rnfpolredabs(GEN nf, GEN R, long flag)
3034 175 : { 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 : }
|
