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 : Check the License for details. You should have received a copy of it, along
11 : with the package; see the file 'COPYING'. If not, write to the Free Software
12 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
13 :
14 : /*******************************************************************/
15 : /* */
16 : /* MAXIMAL ORDERS */
17 : /* */
18 : /*******************************************************************/
19 : #include "pari.h"
20 : #include "paripriv.h"
21 :
22 : #define DEBUGLEVEL DEBUGLEVEL_nf
23 :
24 : /* allow p = -1 from factorizations, avoid oo loop on p = 1 */
25 : static long
26 15029 : safe_Z_pvalrem(GEN x, GEN p, GEN *z)
27 : {
28 15029 : if (is_pm1(p))
29 : {
30 28 : if (signe(p) > 0) return gvaluation(x,p); /*error*/
31 21 : *z = absi(x); return 1;
32 : }
33 15001 : return Z_pvalrem(x, p, z);
34 : }
35 : /* D an integer, P a ZV, return a factorization matrix for D over P, removing
36 : * entries with 0 exponent. */
37 : static GEN
38 4235 : fact_from_factors(GEN D, GEN P, long flag)
39 : {
40 4235 : long i, l = lg(P), iq = 1;
41 4235 : GEN Q = cgetg(l+1,t_COL);
42 4235 : GEN E = cgetg(l+1,t_COL);
43 19257 : for (i=1; i<l; i++)
44 : {
45 15029 : GEN p = gel(P,i);
46 : long k;
47 15029 : if (flag && !equalim1(p))
48 : {
49 14 : p = gcdii(p, D);
50 14 : if (is_pm1(p)) continue;
51 : }
52 15029 : k = safe_Z_pvalrem(D, p, &D);
53 15022 : if (k) { gel(Q,iq) = p; gel(E,iq) = utoipos(k); iq++; }
54 : }
55 4228 : D = absi_shallow(D);
56 4228 : if (!equali1(D))
57 : {
58 819 : long k = Z_isanypower(D, &D);
59 819 : if (!k) k = 1;
60 819 : gel(Q,iq) = D; gel(E,iq) = utoipos(k); iq++;
61 : }
62 4228 : setlg(Q,iq);
63 4228 : setlg(E,iq); return mkmat2(Q,E);
64 : }
65 :
66 : /* d a t_INT; f a t_MAT factorisation of some t_INT sharing some divisors
67 : * with d, or a prime (t_INT). Return a factorization F of d: "primes"
68 : * entries in f _may_ be composite, and are included as is in d. */
69 : static GEN
70 2355 : update_fact(GEN d, GEN f)
71 : {
72 : GEN P;
73 2355 : switch (typ(f))
74 : {
75 2341 : case t_INT: case t_VEC: case t_COL: return f;
76 14 : case t_MAT:
77 14 : if (lg(f) == 3) { P = gel(f,1); break; }
78 : /*fall through*/
79 : default:
80 7 : pari_err_TYPE("nfbasis [factorization expected]",f);
81 : return NULL;/*LCOV_EXCL_LINE*/
82 : }
83 7 : return fact_from_factors(d, P, 1);
84 : }
85 :
86 : /* T = C T0(X/L); C = L^d / lt(T0), d = deg(T)
87 : * disc T = C^2(d - 1) L^-(d(d-1)) disc T0 = (L^d / lt(T0)^2)^(d-1) disc T0 */
88 : static GEN
89 816222 : set_disc(nfmaxord_t *S)
90 : {
91 : GEN L, dT;
92 : long d;
93 816222 : if (S->T0 == S->T) return ZX_disc(S->T);
94 249992 : d = degpol(S->T0);
95 249992 : L = S->unscale;
96 249992 : if (typ(L) == t_FRAC && abscmpii(gel(L,1), gel(L,2)) < 0)
97 11778 : dT = ZX_disc(S->T); /* more efficient */
98 : else
99 : {
100 238214 : GEN l0 = leading_coeff(S->T0);
101 238214 : GEN a = gpowgs(gdiv(gpowgs(L, d), sqri(l0)), d-1);
102 238214 : dT = gmul(a, ZX_disc(S->T0)); /* more efficient */
103 : }
104 249992 : return S->dT = dT;
105 : }
106 :
107 : /* dT != 0 */
108 : static GEN
109 791395 : poldiscfactors_i(GEN T, GEN dT, long flag)
110 : {
111 : GEN U, fa, Z, E, P, Tp;
112 : long i, l;
113 :
114 791395 : fa = absZ_factor_limit_strict(dT, minuu(tridiv_bound(dT), maxprime()), &U);
115 791395 : if (!U) return fa;
116 788 : Z = mkcol(gel(U,1)); P = gel(fa,1); Tp = NULL;
117 1702 : while (lg(Z) != 1)
118 : { /* pop and handle last element of Z */
119 914 : GEN p = veclast(Z), r;
120 914 : setlg(Z, lg(Z)-1);
121 914 : if (!Tp) /* first time: p is composite and not a power */
122 788 : Tp = ZX_deriv(T);
123 : else
124 : {
125 126 : (void)Z_isanypower(p, &p);
126 126 : if ((flag || lgefint(p)==3) && BPSW_psp(p))
127 96 : { P = vec_append(P, p); continue; }
128 : }
129 818 : r = FpX_gcd_check(T, Tp, p);
130 818 : if (r)
131 63 : Z = shallowconcat(Z, Z_cba(r, diviiexact(p,r)));
132 755 : else if (flag)
133 7 : P = shallowconcat(P, gel(Z_factor(p),1));
134 : else
135 748 : P = vec_append(P, p);
136 : }
137 788 : ZV_sort_inplace(P); l = lg(P); E = cgetg(l, t_COL);
138 7104 : for (i = 1; i < l; i++) gel(E,i) = utoipos(Z_pvalrem(dT, gel(P,i), &dT));
139 788 : return mkmat2(P,E);
140 : }
141 :
142 : GEN
143 42 : poldiscfactors(GEN T, long flag)
144 : {
145 42 : pari_sp av = avma;
146 : GEN dT;
147 42 : if (typ(T) != t_POL || !RgX_is_ZX(T)) pari_err_TYPE("poldiscfactors",T);
148 42 : if (flag < 0 || flag > 1) pari_err_FLAG("poldiscfactors");
149 42 : dT = ZX_disc(T);
150 42 : if (!signe(dT)) retmkvec2(gen_0, Z_factor(gen_0));
151 35 : return gc_GEN(av, mkvec2(dT, poldiscfactors_i(T, dT, flag)));
152 : }
153 :
154 : static void
155 816222 : nfmaxord_check_args(nfmaxord_t *S, GEN T, long flag)
156 : {
157 816222 : GEN dT, L, E, P, fa = NULL;
158 : pari_timer t;
159 816222 : long l, ty = typ(T);
160 :
161 816222 : if (DEBUGLEVEL) timer_start(&t);
162 816222 : if (ty == t_VEC) {
163 24862 : if (lg(T) != 3) pari_err_TYPE("nfmaxord",T);
164 24862 : fa = gel(T,2); T = gel(T,1); ty = typ(T);
165 : }
166 816222 : if (ty != t_POL) pari_err_TYPE("nfmaxord",T);
167 816222 : T = Q_primpart(T);
168 816222 : if (degpol(T) <= 0) pari_err_CONSTPOL("nfmaxord");
169 816222 : RgX_check_ZX(T, "nfmaxord");
170 816222 : S->T0 = T;
171 816222 : S->T = T = ZX_Q_normalize(T, &L);
172 816222 : S->unscale = L;
173 816222 : S->dT = dT = set_disc(S);
174 816222 : S->certify = 1;
175 816222 : if (!signe(dT)) pari_err_IRREDPOL("nfmaxord",T);
176 816222 : if (fa)
177 : {
178 24862 : const long MIN = 100; /* include at least all p < 101 */
179 24862 : GEN P0 = NULL, U;
180 24862 : S->certify = 0;
181 24862 : if (!isint1(L)) fa = update_fact(dT, fa);
182 24855 : switch(typ(fa))
183 : {
184 224 : case t_MAT:
185 224 : if (!is_Z_factornon0(fa)) pari_err_TYPE("nfmaxord",fa);
186 217 : P0 = gel(fa,1); /* fall through */
187 4228 : case t_VEC: case t_COL:
188 4228 : if (!P0)
189 : {
190 4011 : if (!RgV_is_ZV(fa)) pari_err_TYPE("nfmaxord",fa);
191 4011 : P0 = fa;
192 : }
193 4228 : P = gel(absZ_factor_limit_strict(dT, MIN, &U), 1);
194 4228 : if (lg(P) != 0) { settyp(P, typ(P0)); P0 = shallowconcat(P0,P); }
195 4228 : P0 = ZV_sort_uniq_shallow(P0);
196 4228 : fa = fact_from_factors(dT, P0, 0);
197 4221 : break;
198 20613 : case t_INT:
199 20613 : fa = absZ_factor_limit(dT, (signe(fa) <= 0)? 1: maxuu(itou(fa), MIN));
200 20613 : break;
201 7 : default:
202 7 : pari_err_TYPE("nfmaxord",fa);
203 : }
204 : }
205 : else
206 : {
207 791360 : S->certify = !(flag & nf_PARTIALFACT);
208 791360 : fa = poldiscfactors_i(T, dT, 0);
209 : }
210 816194 : P = gel(fa,1); l = lg(P);
211 816194 : E = gel(fa,2);
212 816194 : if (l > 1 && is_pm1(gel(P,1)))
213 : {
214 21 : l--;
215 21 : P = vecslice(P, 2, l);
216 21 : E = vecslice(E, 2, l);
217 : }
218 816194 : S->dTP = P;
219 816194 : S->dTE = vec_to_vecsmall(E);
220 816194 : if (DEBUGLEVEL>2) timer_printf(&t, "disc. factorisation");
221 816194 : }
222 :
223 : static int
224 224693 : fnz(GEN x,long j)
225 : {
226 : long i;
227 729245 : for (i=1; i<j; i++)
228 555768 : if (signe(gel(x,i))) return 0;
229 173477 : return 1;
230 : }
231 : /* return list u[i], 2 by 2 coprime with the same prime divisors as ab */
232 : static GEN
233 294 : get_coprimes(GEN a, GEN b)
234 : {
235 294 : long i, k = 1;
236 294 : GEN u = cgetg(3, t_COL);
237 294 : gel(u,1) = a;
238 294 : gel(u,2) = b;
239 : /* u1,..., uk 2 by 2 coprime */
240 1071 : while (k+1 < lg(u))
241 : {
242 777 : GEN d, c = gel(u,k+1);
243 777 : if (is_pm1(c)) { k++; continue; }
244 1309 : for (i=1; i<=k; i++)
245 : {
246 840 : GEN ui = gel(u,i);
247 840 : if (is_pm1(ui)) continue;
248 483 : d = gcdii(c, ui);
249 483 : if (d == gen_1) continue;
250 483 : c = diviiexact(c, d);
251 483 : gel(u,i) = diviiexact(ui, d);
252 483 : u = vec_append(u, d);
253 : }
254 469 : gel(u,++k) = c;
255 : }
256 1365 : for (i = k = 1; i < lg(u); i++)
257 1071 : if (!is_pm1(gel(u,i))) gel(u,k++) = gel(u,i);
258 294 : setlg(u, k); return u;
259 : }
260 :
261 : /*******************************************************************/
262 : /* */
263 : /* ROUND 4 */
264 : /* */
265 : /*******************************************************************/
266 : typedef struct {
267 : /* constants */
268 : long pisprime; /* -1: unknown, 1: prime, 0: composite */
269 : GEN p, f; /* goal: factor f p-adically */
270 : long df;
271 : GEN pdf; /* p^df = reduced discriminant of f */
272 : long mf; /* */
273 : GEN psf, pmf; /* stability precision for f, wanted precision for f */
274 : long vpsf; /* v_p(p_f) */
275 : /* these are updated along the way */
276 : GEN phi; /* a p-integer, in Q[X] */
277 : GEN phi0; /* a p-integer, in Q[X] from testb2 / testc2, to be composed with
278 : * phi when correct precision is known */
279 : GEN chi; /* characteristic polynomial of phi (mod psc) in Z[X] */
280 : GEN nu; /* irreducible divisor of chi mod p, in Z[X] */
281 : GEN invnu; /* numerator ( 1/ Mod(nu, chi) mod pmr ) */
282 : GEN Dinvnu;/* denominator ( ... ) */
283 : long vDinvnu; /* v_p(Dinvnu) */
284 : GEN prc, psc; /* reduced discriminant of chi, stability precision for chi */
285 : long vpsc; /* v_p(p_c) */
286 : GEN ns, nsf, precns; /* cached Newton sums for nsf and their precision */
287 : } decomp_t;
288 : static GEN maxord_i(decomp_t *S, GEN p, GEN f, long mf, GEN w, long flag);
289 : static GEN dbasis(GEN p, GEN f, long mf, GEN alpha, GEN U);
290 : static GEN maxord(GEN p,GEN f,long mf);
291 : static GEN ZX_Dedekind(GEN F, GEN *pg, GEN p);
292 :
293 : static void
294 512 : fix_PE(GEN *pP, GEN *pE, long i, GEN u, GEN N)
295 : {
296 : GEN P, E;
297 512 : long k, l = lg(u), lP = lg(*pP);
298 : pari_sp av;
299 :
300 512 : *pP = P = shallowconcat(*pP, vecslice(u, 2, l-1));
301 512 : *pE = E = vecsmall_lengthen(*pE, lP + l-2);
302 512 : gel(P,i) = gel(u,1); av = avma;
303 512 : E[i] = Z_pvalrem(N, gel(P,i), &N);
304 1031 : for (k=lP, lP=lg(P); k < lP; k++) E[k] = Z_pvalrem(N, gel(P,k), &N);
305 512 : set_avma(av);
306 512 : }
307 : static long
308 686746 : diag_denomval(GEN M, GEN p)
309 : {
310 : long j, v, l;
311 686746 : if (typ(M) != t_MAT) return 0;
312 463431 : v = 0; l = lg(M);
313 2072780 : for (j=1; j<l; j++)
314 : {
315 1609349 : GEN t = gcoeff(M,j,j);
316 1609349 : if (typ(t) == t_FRAC) v += Z_pval(gel(t,2), p);
317 : }
318 463431 : return v;
319 : }
320 :
321 : /* n > 1 is composite, not a pure power, and has no prime divisor < 2^14;
322 : * return a BPSW divisor of n and smallest k-th root of largest coprime cofactor */
323 : static GEN
324 197 : Z_fac(GEN n)
325 : {
326 197 : GEN p = icopy(n), part = ifac_start(p, 0);
327 : long e;
328 197 : ifac_next(&part, &p, &e); n = diviiexact(n, powiu(p, e));
329 197 : (void)Z_isanypower(n, &n); return mkvec2(p, n);
330 : }
331 :
332 : /* Warning: data computed for T = ZX_Q_normalize(T0). If S.unscale !=
333 : * gen_1, caller must take steps to correct the components if it wishes
334 : * to stick to the original T0. Return a vector of p-maximal orders, for
335 : * those p s.t p^2 | disc(T) [ = S->dTP ]*/
336 : static GEN
337 816222 : get_maxord(nfmaxord_t *S, GEN T0, long flag)
338 : {
339 : GEN P, E;
340 : VOLATILE GEN O;
341 : VOLATILE long lP, i;
342 :
343 816222 : nfmaxord_check_args(S, T0, flag);
344 816194 : P = S->dTP; lP = lg(P);
345 816194 : E = S->dTE;
346 816194 : O = cgetg(1, t_VEC);
347 3154858 : for (i=1; i<lP; i=i+1)
348 : {
349 : VOLATILE pari_sp av;
350 : /* includes the silly case where P[i] = -1 */
351 2338664 : if (E[i] <= 1)
352 : {
353 1307182 : if (S->certify)
354 : {
355 1298560 : GEN p = gel(P,i);
356 1298560 : if (signe(p) > 0 && !BPSW_psp(p))
357 : {
358 197 : fix_PE(&P, &E, i, Z_fac(p), S->dT);
359 197 : lP = lg(P); i=i-1; continue;
360 : }
361 : }
362 1306985 : O = vec_append(O, gen_1); continue;
363 : }
364 1031482 : av = avma;
365 1031482 : pari_CATCH(CATCH_ALL) {
366 294 : GEN u, err = pari_err_last();
367 : long l;
368 294 : switch(err_get_num(err))
369 : {
370 294 : case e_INV:
371 : {
372 294 : GEN p, x = err_get_compo(err, 2);
373 : long k;
374 294 : if (typ(x) == t_INTMOD)
375 : { /* caught false prime, update factorization */
376 294 : p = gcdii(gel(x,1), gel(x,2));
377 294 : u = diviiexact(gel(x,1),p);
378 294 : if (DEBUGLEVEL) pari_warn(warner,"impossible inverse: %Ps", x);
379 294 : (void)gc_all(av, 2, &p, &u);
380 :
381 294 : u = get_coprimes(p, u); l = lg(u);
382 : /* no small factors, but often a prime power */
383 882 : for (k = 1; k < l; k++) (void)Z_isanypower(gel(u,k), &gel(u,k));
384 294 : break;
385 : }
386 : /* fall through */
387 : }
388 : case e_PRIME: case e_IRREDPOL:
389 : { /* we're here because we failed BPSW_isprime(), no point in
390 : * reporting a possible counter-example to the BPSW test */
391 0 : GEN p = gel(P,i);
392 0 : set_avma(av);
393 0 : if (DEBUGLEVEL)
394 0 : pari_warn(warner,"large composite in nfmaxord:loop(), %Ps", p);
395 0 : if (expi(p) < 100)
396 0 : u = gel(Z_factor(p), 1); /* p < 2^100 should take ~20ms */
397 0 : else if (S->certify)
398 0 : u = Z_fac(p);
399 : else
400 0 : { /* give up, probably not maximal */
401 0 : GEN B, g, k = ZX_Dedekind(S->T, &g, p);
402 0 : k = FpX_normalize(k, p);
403 0 : B = dbasis(p, S->T, E[i], NULL, FpX_div(S->T,k,p));
404 0 : O = vec_append(O, B);
405 0 : pari_CATCH_reset(); continue;
406 : }
407 0 : break;
408 : }
409 0 : default: pari_err(0, err);
410 : return NULL;/*LCOV_EXCL_LINE*/
411 : }
412 294 : fix_PE(&P, &E, i, u, S->dT);
413 294 : lP = lg(P); av = avma;
414 1031776 : } pari_RETRY {
415 1031776 : GEN p = gel(P,i), O2;
416 1031776 : if (DEBUGLEVEL>2) err_printf("Treating p^k = %Ps^%ld\n",p,E[i]);
417 1031776 : O2 = maxord(p,S->T,E[i]);
418 1031482 : if (S->certify && (odd(E[i]) || E[i] != 2*diag_denomval(O2, p))
419 611202 : && !BPSW_psp(p))
420 : {
421 21 : fix_PE(&P, &E, i, gel(Z_factor(p), 1), S->dT);
422 21 : lP = lg(P); i=i-1;
423 : }
424 : else
425 1031461 : O = vec_append(O, O2);
426 1031482 : } pari_ENDCATCH;
427 : }
428 816194 : S->dTP = P; S->dTE = E; return O;
429 : }
430 :
431 : /* M a QM, return denominator of diagonal. All denominators are powers of
432 : * a given integer */
433 : static GEN
434 100813 : diag_denom(GEN M)
435 : {
436 100813 : GEN d = gen_1;
437 100813 : long j, l = lg(M);
438 703978 : for (j=1; j<l; j++)
439 : {
440 603165 : GEN t = gcoeff(M,j,j);
441 603165 : if (typ(t) == t_INT) continue;
442 214709 : t = gel(t,2);
443 214709 : if (abscmpii(t,d) > 0) d = t;
444 : }
445 100813 : return d;
446 : }
447 : static void
448 747321 : setPE(GEN D, GEN P, GEN *pP, GEN *pE)
449 : {
450 747321 : long k, j, l = lg(P);
451 : GEN P2, E2;
452 747321 : *pP = P2 = cgetg(l, t_VEC);
453 747321 : *pE = E2 = cgetg(l, t_VECSMALL);
454 2874067 : for (k = j = 1; j < l; j++)
455 : {
456 2126746 : long v = Z_pvalrem(D, gel(P,j), &D);
457 2126746 : if (v) { gel(P2,k) = gel(P,j); E2[k] = v; k++; }
458 : }
459 747321 : setlg(P2, k);
460 747321 : setlg(E2, k);
461 747321 : }
462 : void
463 102530 : nfmaxord(nfmaxord_t *S, GEN T0, long flag)
464 : {
465 102530 : GEN O = get_maxord(S, T0, flag);
466 102523 : GEN f = S->T, P = S->dTP, a = NULL, da = NULL;
467 102523 : long n = degpol(f), lP = lg(P), i, j, k;
468 102523 : int centered = 0;
469 102523 : pari_sp av = avma;
470 : /* r1 & basden not initialized here */
471 102523 : S->r1 = -1;
472 102523 : S->basden = NULL;
473 360052 : for (i=1; i<lP; i++)
474 : {
475 257529 : GEN M, db, b = gel(O,i);
476 257529 : if (b == gen_1) continue;
477 100813 : db = diag_denom(b);
478 100813 : if (db == gen_1) continue;
479 :
480 : /* db = denom(b), (da,db) = 1. Compute da Im(b) + db Im(a) */
481 100813 : b = Q_muli_to_int(b,db);
482 100813 : if (!da) { da = db; a = b; }
483 : else
484 : { /* optimization: easy as long as both matrix are diagonal */
485 136449 : j=2; while (j<=n && fnz(gel(a,j),j) && fnz(gel(b,j),j)) j++;
486 51230 : k = j-1; M = cgetg(2*n-k+1,t_MAT);
487 187679 : for (j=1; j<=k; j++)
488 : {
489 136449 : gel(M,j) = gel(a,j);
490 136449 : gcoeff(M,j,j) = mulii(gcoeff(a,j,j),gcoeff(b,j,j));
491 : }
492 : /* could reduce mod M(j,j) but not worth it: usually close to da*db */
493 282999 : for ( ; j<=n; j++) gel(M,j) = ZC_Z_mul(gel(a,j), db);
494 282999 : for ( ; j<=2*n-k; j++) gel(M,j) = ZC_Z_mul(gel(b,j+k-n), da);
495 51230 : da = mulii(da,db);
496 51230 : a = ZM_hnfmodall_i(M, da, hnf_MODID|hnf_CENTER);
497 51230 : (void)gc_all(av, 2, &a, &da);
498 51230 : centered = 1;
499 : }
500 : }
501 102523 : if (da)
502 : {
503 49583 : GEN index = diviiexact(da, gcoeff(a,1,1));
504 234947 : for (j=2; j<=n; j++) index = mulii(index, diviiexact(da, gcoeff(a,j,j)));
505 49583 : if (!centered) a = ZM_hnfcenter(a);
506 49583 : a = RgM_Rg_div(a, da);
507 49583 : S->index = index;
508 49583 : S->dK = diviiexact(S->dT, sqri(index));
509 : }
510 : else
511 : {
512 52940 : S->index = gen_1;
513 52940 : S->dK = S->dT;
514 52940 : a = matid(n);
515 : }
516 102523 : setPE(S->dK, P, &S->dKP, &S->dKE);
517 102523 : S->basis = RgM_to_RgXV(a, varn(f));
518 102523 : }
519 : GEN
520 938 : nfbasis(GEN x, GEN *pdK)
521 : {
522 938 : pari_sp av = avma;
523 : nfmaxord_t S;
524 : GEN B;
525 938 : nfmaxord(&S, x, 0);
526 938 : B = RgXV_unscale(S.basis, S.unscale);
527 938 : if (pdK) *pdK = S.dK;
528 938 : return gc_all(av, pdK? 2: 1, &B, pdK);
529 : }
530 : /* field discriminant: faster than nfmaxord, use local data only */
531 : static GEN
532 713692 : maxord_disc(nfmaxord_t *S, GEN x)
533 : {
534 713692 : GEN O = get_maxord(S, x, 0), I = gen_1;
535 713671 : long n = degpol(S->T), lP = lg(O), i, j;
536 2794588 : for (i = 1; i < lP; i++)
537 : {
538 2080917 : GEN b = gel(O,i);
539 2080917 : if (b == gen_1) continue;
540 2708615 : for (j = 1; j <= n; j++)
541 : {
542 2116877 : GEN c = gcoeff(b,j,j);
543 2116877 : if (typ(c) == t_FRAC) I = mulii(I, gel(c,2)) ;
544 : }
545 : }
546 713671 : return diviiexact(S->dT, sqri(I));
547 : }
548 : GEN
549 68901 : nfdisc(GEN x)
550 : {
551 68901 : pari_sp av = avma;
552 : nfmaxord_t S;
553 68901 : return gc_INT(av, maxord_disc(&S, x));
554 : }
555 : GEN
556 644798 : nfdiscfactors(GEN x)
557 : {
558 644798 : pari_sp av = avma;
559 644798 : GEN E, P, D, nf = checknf_i(x);
560 644798 : if (nf)
561 : {
562 7 : D = nf_get_disc(nf);
563 7 : P = nf_get_ramified_primes(nf);
564 : }
565 : else
566 : {
567 : nfmaxord_t S;
568 644791 : D = maxord_disc(&S, x);
569 644791 : P = S.dTP;
570 : }
571 644798 : setPE(D, P, &P, &E); settyp(P, t_COL);
572 644798 : return gc_GEN(av, mkvec2(D, mkmat2(P, zc_to_ZC(E))));
573 : }
574 :
575 : static ulong
576 1604583 : Flx_checkdeflate(GEN x)
577 : {
578 1604583 : ulong d = 0, i, lx = (ulong)lg(x);
579 2573035 : for (i=3; i<lx; i++)
580 1722490 : if (x[i]) { d = ugcd(d,i-2); if (d == 1) break; }
581 1604583 : return d;
582 : }
583 :
584 : /* product of (monic) irreducible factors of f over Fp[X]
585 : * Assume f reduced mod p, otherwise valuation at x may be wrong */
586 : static GEN
587 1604583 : Flx_radical(GEN f, ulong p)
588 : {
589 1604583 : long v0 = Flx_valrem(f, &f);
590 : ulong du, d, e;
591 : GEN u;
592 :
593 1604583 : d = Flx_checkdeflate(f);
594 1604583 : if (!d) return v0? polx_Flx(f[1]): pol1_Flx(f[1]);
595 1006695 : if (u_lvalrem(d,p, &e)) f = Flx_deflate(f, d/e); /* f(x^p^i) -> f(x) */
596 1006695 : u = Flx_gcd(f, Flx_deriv(f, p), p); /* (f,f') */
597 1006689 : du = degpol(u);
598 1006689 : if (du)
599 : {
600 318507 : if (du == (ulong)degpol(f))
601 0 : f = Flx_radical(Flx_deflate(f,p), p);
602 : else
603 : {
604 318507 : u = Flx_normalize(u, p);
605 318507 : f = Flx_div(f, u, p);
606 318507 : if (p <= du)
607 : {
608 66844 : GEN w = (degpol(f) >= degpol(u))? Flx_rem(f, u, p): f;
609 66844 : w = Flxq_powu(w, du, u, p);
610 66844 : w = Flx_div(u, Flx_gcd(w,u,p), p); /* u / gcd(u, v^(deg u-1)) */
611 66844 : f = Flx_mul(f, Flx_radical(Flx_deflate(w,p), p), p);
612 : }
613 : }
614 : }
615 1006689 : if (v0) f = Flx_shift(f, 1);
616 1006689 : return f;
617 : }
618 : /* Assume f reduced mod p, otherwise valuation at x may be wrong */
619 : static GEN
620 5697 : FpX_radical(GEN f, GEN p)
621 : {
622 : GEN u;
623 : long v0;
624 5697 : if (lgefint(p) == 3)
625 : {
626 1754 : ulong q = p[2];
627 1754 : return Flx_to_ZX( Flx_radical(ZX_to_Flx(f, q), q) );
628 : }
629 3943 : v0 = ZX_valrem(f, &f);
630 3943 : u = FpX_gcd(f,FpX_deriv(f, p), p);
631 3655 : if (degpol(u)) f = FpX_div(f, u, p);
632 3655 : if (v0) f = RgX_shift(f, 1);
633 3655 : return f;
634 : }
635 : /* f / a */
636 : static GEN
637 1535985 : zx_z_div(GEN f, ulong a)
638 : {
639 1535985 : long i, l = lg(f);
640 1535985 : GEN g = cgetg(l, t_VECSMALL);
641 1535985 : g[1] = f[1];
642 5219702 : for (i = 2; i < l; i++) g[i] = f[i] / a;
643 1535985 : return g;
644 : }
645 : /* Dedekind criterion; return k = gcd(g,h, (f-gh)/p), where
646 : * f = \prod f_i^e_i, g = \prod f_i, h = \prod f_i^{e_i-1}
647 : * k = 1 iff Z[X]/(f) is p-maximal */
648 : static GEN
649 1541682 : ZX_Dedekind(GEN F, GEN *pg, GEN p)
650 : {
651 : GEN k, h, g, f, f2;
652 1541682 : ulong q = p[2];
653 1541682 : if (lgefint(p) == 3 && q < (1UL << BITS_IN_HALFULONG))
654 1535985 : {
655 1535985 : ulong q2 = q*q;
656 1535985 : f2 = ZX_to_Flx(F, q2);
657 1535985 : f = Flx_red(f2, q);
658 1535985 : g = Flx_radical(f, q);
659 1535985 : h = Flx_div(f, g, q);
660 1535985 : k = zx_z_div(Flx_sub(f2, Flx_mul(g,h,q2), q2), q);
661 1535985 : k = Flx_gcd(k, Flx_gcd(g,h,q), q);
662 1535985 : k = Flx_to_ZX(k);
663 1535985 : g = Flx_to_ZX(g);
664 : }
665 : else
666 : {
667 5697 : f2 = FpX_red(F, sqri(p));
668 5697 : f = FpX_red(f2, p);
669 5697 : g = FpX_radical(f, p);
670 5403 : h = FpX_div(f, g, p);
671 5403 : k = ZX_Z_divexact(ZX_sub(f2, ZX_mul(g,h)), p);
672 5403 : k = FpX_gcd(FpX_red(k, p), FpX_gcd(g,h,p), p);
673 : }
674 1541388 : *pg = g; return k;
675 : }
676 :
677 : /* p-maximal order of Z[x]/f; mf = v_p(Disc(f)) or < 0 [unknown].
678 : * Return gen_1 if p-maximal */
679 : static GEN
680 1541682 : maxord(GEN p, GEN f, long mf)
681 : {
682 1541682 : const pari_sp av = avma;
683 1541682 : GEN res, g, k = ZX_Dedekind(f, &g, p);
684 1541388 : long dk = degpol(k);
685 1541388 : if (DEBUGLEVEL>2) err_printf(" ZX_Dedekind: gcd has degree %ld\n", dk);
686 1541388 : if (!dk) { set_avma(av); return gen_1; }
687 876450 : if (mf < 0) mf = ZpX_disc_val(f, p);
688 876450 : k = FpX_normalize(k, p);
689 876450 : if (2*dk >= mf-1)
690 421400 : res = dbasis(p, f, mf, NULL, FpX_div(f,k,p));
691 : else
692 : {
693 : GEN w, F1, F2;
694 : decomp_t S;
695 455050 : F1 = FpX_factor(k,p);
696 455050 : F2 = FpX_factor(FpX_div(g,k,p),p);
697 455050 : w = merge_sort_uniq(gel(F1,1),gel(F2,1),(void*)cmpii,&gen_cmp_RgX);
698 455050 : res = maxord_i(&S, p, f, mf, w, 0);
699 : }
700 876450 : return gc_GEN(av,res);
701 : }
702 : /* T monic separable ZX, p prime */
703 : GEN
704 0 : ZpX_primedec(GEN T, GEN p)
705 : {
706 0 : const pari_sp av = avma;
707 0 : GEN w, F1, F2, res, g, k = ZX_Dedekind(T, &g, p);
708 : decomp_t S;
709 0 : if (!degpol(k)) return zm_to_ZM(FpX_degfact(T, p));
710 0 : k = FpX_normalize(k, p);
711 0 : F1 = FpX_factor(k,p);
712 0 : F2 = FpX_factor(FpX_div(g,k,p),p);
713 0 : w = merge_sort_uniq(gel(F1,1),gel(F2,1),(void*)cmpii,&gen_cmp_RgX);
714 0 : res = maxord_i(&S, p, T, ZpX_disc_val(T, p), w, -1);
715 0 : if (!res)
716 : {
717 0 : long f = degpol(S.nu), e = degpol(T) / f;
718 0 : set_avma(av); retmkmat2(mkcols(f), mkcols(e));
719 : }
720 0 : return gc_GEN(av,res);
721 : }
722 :
723 : static GEN
724 4808736 : Zlx_sylvester_echelon(GEN f1, GEN f2, long early_abort, ulong p, ulong pm)
725 : {
726 4808736 : long j, n = degpol(f1);
727 4808736 : GEN h, a = cgetg(n+1,t_MAT);
728 4808736 : f1 = Flx_get_red(f1, pm);
729 4808736 : h = Flx_rem(f2,f1,pm);
730 4808736 : for (j=1;; j++)
731 : {
732 16831170 : gel(a,j) = Flx_to_Flv(h, n);
733 16831170 : if (j == n) break;
734 12022434 : h = Flx_rem(Flx_shift(h, 1), f1, pm);
735 : }
736 4808736 : return zlm_echelon(a, early_abort, p, pm);
737 : }
738 : /* Sylvester's matrix, mod p^m (assumes f1 monic). If early_abort
739 : * is set, return NULL if one pivot is 0 mod p^m */
740 : static GEN
741 74449 : ZpX_sylvester_echelon(GEN f1, GEN f2, long early_abort, GEN p, GEN pm)
742 : {
743 74449 : long j, n = degpol(f1);
744 74449 : GEN h, a = cgetg(n+1,t_MAT);
745 74449 : h = FpXQ_red(f2,f1,pm);
746 74449 : for (j=1;; j++)
747 : {
748 432401 : gel(a,j) = RgX_to_RgC(h, n);
749 432401 : if (j == n) break;
750 357952 : h = FpX_rem(RgX_shift_shallow(h, 1), f1, pm);
751 : }
752 74449 : return ZpM_echelon(a, early_abort, p, pm);
753 : }
754 :
755 : /* polynomial gcd mod p^m (assumes f1 monic). Return a QpX ! */
756 : static GEN
757 246519 : Zlx_gcd(GEN f1, GEN f2, ulong p, ulong pm)
758 : {
759 246519 : pari_sp av = avma;
760 246519 : GEN a = Zlx_sylvester_echelon(f1,f2,0,p,pm);
761 246519 : long c, l = lg(a), sv = f1[1];
762 762541 : for (c = 1; c < l; c++)
763 : {
764 762541 : ulong t = ucoeff(a,c,c);
765 762541 : if (t)
766 : {
767 246519 : a = Flx_to_ZX(Flv_to_Flx(gel(a,c), sv));
768 246519 : if (t == 1) return gc_GEN(av, a);
769 74803 : return gc_upto(av, RgX_Rg_div(a, utoipos(t)));
770 : }
771 : }
772 0 : set_avma(av);
773 0 : a = cgetg(2,t_POL); a[1] = sv; return a;
774 : }
775 : GEN
776 255254 : ZpX_gcd(GEN f1, GEN f2, GEN p, GEN pm)
777 : {
778 255254 : pari_sp av = avma;
779 : GEN a;
780 : long c, l, v;
781 255254 : if (lgefint(pm) == 3)
782 : {
783 246519 : ulong q = pm[2];
784 246519 : return Zlx_gcd(ZX_to_Flx(f1, q), ZX_to_Flx(f2,q), p[2], q);
785 : }
786 8735 : a = ZpX_sylvester_echelon(f1,f2,0,p,pm);
787 8735 : l = lg(a); v = varn(f1);
788 60280 : for (c = 1; c < l; c++)
789 : {
790 60280 : GEN t = gcoeff(a,c,c);
791 60280 : if (signe(t))
792 : {
793 8735 : a = RgV_to_RgX(gel(a,c), v);
794 8735 : if (equali1(t)) return gc_GEN(av, a);
795 2469 : return gc_upto(av, RgX_Rg_div(a, t));
796 : }
797 : }
798 0 : set_avma(av); return pol_0(v);
799 : }
800 :
801 : /* Return m > 0, such that p^m ~ 2^16 for initial value of m; assume p prime */
802 : static long
803 4519710 : init_m(GEN p)
804 : {
805 : ulong pp;
806 4519710 : if (lgefint(p) > 3) return 1;
807 4518555 : pp = p[2]; /* m ~ 16 / log2(pp) */
808 4518555 : if (pp < 41) switch(pp)
809 : {
810 1225600 : case 2: return 16;
811 352745 : case 3: return 10;
812 262389 : case 5: return 6;
813 153344 : case 7: return 5;
814 209491 : case 11: case 13: return 4;
815 331904 : default: return 3;
816 : }
817 1983082 : return pp < 257? 2: 1;
818 : }
819 :
820 : /* reduced resultant mod p^m (assumes x monic) */
821 : GEN
822 995793 : ZpX_reduced_resultant(GEN x, GEN y, GEN p, GEN pm)
823 : {
824 995793 : pari_sp av = avma;
825 : GEN z;
826 995793 : if (lgefint(pm) == 3)
827 : {
828 983944 : ulong q = pm[2];
829 983944 : z = Zlx_sylvester_echelon(ZX_to_Flx(x,q), ZX_to_Flx(y,q),0,p[2],q);
830 983944 : if (lg(z) > 1)
831 : {
832 983944 : ulong c = ucoeff(z,1,1);
833 983944 : if (c) return gc_utoipos(av, c);
834 : }
835 : }
836 : else
837 : {
838 11849 : z = ZpX_sylvester_echelon(x,y,0,p,pm);
839 11849 : if (lg(z) > 1)
840 : {
841 11849 : GEN c = gcoeff(z,1,1);
842 11849 : if (signe(c)) return gc_INT(av, c);
843 : }
844 : }
845 129582 : set_avma(av); return gen_0;
846 : }
847 : /* Assume Res(f,g) divides p^M. Return Res(f, g), using dynamic p-adic
848 : * precision (until result is nonzero or p^M). */
849 : GEN
850 933998 : ZpX_reduced_resultant_fast(GEN f, GEN g, GEN p, long M)
851 : {
852 933998 : GEN R, q = NULL;
853 : long m;
854 933998 : m = init_m(p); if (m < 1) m = 1;
855 61795 : for(;; m <<= 1) {
856 995793 : if (M < 2*m) break;
857 95159 : q = q? sqri(q): powiu(p, m); /* p^m */
858 95159 : R = ZpX_reduced_resultant(f,g, p, q); if (signe(R)) return R;
859 : }
860 900634 : q = powiu(p, M);
861 900634 : R = ZpX_reduced_resultant(f,g, p, q); return signe(R)? R: q;
862 : }
863 :
864 : /* v_p(Res(x,y) mod p^m), assumes (lc(x),p) = 1 */
865 : static long
866 3632138 : ZpX_resultant_val_i(GEN x, GEN y, GEN p, GEN pm)
867 : {
868 3632138 : pari_sp av = avma;
869 : GEN z;
870 : long i, l, v;
871 3632138 : if (lgefint(pm) == 3)
872 : {
873 3578273 : ulong q = pm[2], pp = p[2];
874 3578273 : z = Zlx_sylvester_echelon(ZX_to_Flx(x,q), ZX_to_Flx(y,q), 1, pp, q);
875 3578273 : if (!z) return gc_long(av,-1); /* failure */
876 3375681 : v = 0; l = lg(z);
877 14023249 : for (i = 1; i < l; i++) v += u_lval(ucoeff(z,i,i), pp);
878 : }
879 : else
880 : {
881 53865 : z = ZpX_sylvester_echelon(x, y, 1, p, pm);
882 53865 : if (!z) return gc_long(av,-1); /* failure */
883 53052 : v = 0; l = lg(z);
884 196356 : for (i = 1; i < l; i++) v += Z_pval(gcoeff(z,i,i), p);
885 : }
886 3428733 : return v;
887 : }
888 :
889 : /* assume (lc(f),p) = 1; no assumption on g */
890 : long
891 3585712 : ZpX_resultant_val(GEN f, GEN g, GEN p, long M)
892 : {
893 3585712 : pari_sp av = avma;
894 3585712 : GEN q = NULL;
895 : long v, m;
896 3585712 : m = init_m(p); if (m < 2) m = 2;
897 46426 : for(;; m <<= 1) {
898 3632138 : if (m > M) m = M;
899 3632138 : q = q? sqri(q): powiu(p, m); /* p^m */
900 3632138 : v = ZpX_resultant_val_i(f,g, p, q); if (v >= 0) return gc_long(av,v);
901 203405 : if (m == M) return gc_long(av,M);
902 : }
903 : }
904 :
905 : /* assume f separable and (lc(f),p) = 1 */
906 : long
907 185054 : ZpX_disc_val(GEN f, GEN p)
908 : {
909 185054 : pari_sp av = avma;
910 : long v;
911 185054 : if (degpol(f) == 1) return 0;
912 185054 : v = ZpX_resultant_val(f, ZX_deriv(f), p, LONG_MAX);
913 185054 : return gc_long(av,v);
914 : }
915 :
916 : /* *e a ZX, *d, *z in Z, *d = p^(*vd). Simplify e / d by cancelling a
917 : * common factor p^v; if z!=NULL, update it by cancelling the same power of p */
918 : static void
919 3583684 : update_den(GEN p, GEN *e, GEN *d, long *vd, GEN *z)
920 : {
921 : GEN newe;
922 3583684 : long ve = ZX_pvalrem(*e, p, &newe);
923 3583684 : if (ve) {
924 : GEN newd;
925 1764035 : long v = minss(*vd, ve);
926 1764035 : if (v) {
927 1764035 : if (v == *vd)
928 : { /* rare, denominator cancelled */
929 384756 : if (ve != v) newe = ZX_Z_mul(newe, powiu(p, ve - v));
930 384756 : newd = gen_1;
931 384756 : *vd = 0;
932 384756 : if (z) *z =diviiexact(*z, powiu(p, v));
933 : }
934 : else
935 : { /* v = ve < vd, generic case */
936 1379279 : GEN q = powiu(p, v);
937 1379279 : newd = diviiexact(*d, q);
938 1379279 : *vd -= v;
939 1379279 : if (z) *z = diviiexact(*z, q);
940 : }
941 1764035 : *e = newe;
942 1764035 : *d = newd;
943 : }
944 : }
945 3583684 : }
946 :
947 : /* return denominator, a power of p */
948 : static GEN
949 2770381 : QpX_denom(GEN x)
950 : {
951 2770381 : long i, l = lg(x);
952 2770381 : GEN maxd = gen_1;
953 9831563 : for (i=2; i<l; i++)
954 : {
955 7061182 : GEN d = gel(x,i);
956 7061182 : if (typ(d) == t_FRAC && cmpii(gel(d,2), maxd) > 0) maxd = gel(d,2);
957 : }
958 2770381 : return maxd;
959 : }
960 : static GEN
961 509906 : QpXV_denom(GEN x)
962 : {
963 509906 : long l = lg(x), i;
964 509906 : GEN maxd = gen_1;
965 1533059 : for (i = 1; i < l; i++)
966 : {
967 1023153 : GEN d = QpX_denom(gel(x,i));
968 1023153 : if (cmpii(d, maxd) > 0) maxd = d;
969 : }
970 509906 : return maxd;
971 : }
972 :
973 : static GEN
974 1747228 : QpX_remove_denom(GEN x, GEN p, GEN *pdx, long *pv)
975 : {
976 1747228 : *pdx = QpX_denom(x);
977 1747228 : if (*pdx == gen_1) { *pv = 0; *pdx = NULL; }
978 : else {
979 1269802 : x = Q_muli_to_int(x,*pdx);
980 1269802 : *pv = Z_pval(*pdx, p);
981 : }
982 1747228 : return x;
983 : }
984 :
985 : /* p^v * f o g mod (T,q). q = p^vq */
986 : static GEN
987 287745 : compmod(GEN p, GEN f, GEN g, GEN T, GEN q, long v)
988 : {
989 287745 : GEN D = NULL, z, df, dg, qD;
990 287745 : long vD = 0, vdf, vdg;
991 :
992 287745 : f = QpX_remove_denom(f, p, &df, &vdf);
993 287745 : if (typ(g) == t_VEC) /* [num,den,v_p(den)] */
994 0 : { vdg = itos(gel(g,3)); dg = gel(g,2); g = gel(g,1); }
995 : else
996 287745 : g = QpX_remove_denom(g, p, &dg, &vdg);
997 287745 : if (df) { D = df; vD = vdf; }
998 287745 : if (dg) {
999 56197 : long degf = degpol(f);
1000 56197 : D = mul_content(D, powiu(dg, degf));
1001 56197 : vD += degf * vdg;
1002 : }
1003 287745 : qD = D ? mulii(q, D): q;
1004 287745 : if (dg) f = FpX_rescale(f, dg, qD);
1005 287745 : z = FpX_FpXQ_eval(f, g, T, qD);
1006 287745 : if (!D) {
1007 0 : if (v) {
1008 0 : if (v > 0)
1009 0 : z = ZX_Z_mul(z, powiu(p, v));
1010 : else
1011 0 : z = RgX_Rg_div(z, powiu(p, -v));
1012 : }
1013 0 : return z;
1014 : }
1015 287745 : update_den(p, &z, &D, &vD, NULL);
1016 287745 : qD = mulii(D,q);
1017 287745 : if (v) vD -= v;
1018 287745 : z = FpX_center_i(z, qD, shifti(qD,-1));
1019 287745 : if (vD > 0)
1020 287745 : z = RgX_Rg_div(z, powiu(p, vD));
1021 0 : else if (vD < 0)
1022 0 : z = ZX_Z_mul(z, powiu(p, -vD));
1023 287745 : return z;
1024 : }
1025 :
1026 : /* fast implementation of ZM_hnfmodid(M, D) / D, D = p^k */
1027 : static GEN
1028 455050 : ZpM_hnfmodid(GEN M, GEN p, GEN D)
1029 : {
1030 455050 : long i, l = lg(M);
1031 455050 : M = RgM_Rg_div(ZpM_echelon(M,0,p,D), D);
1032 2048425 : for (i = 1; i < l; i++)
1033 1593375 : if (gequal0(gcoeff(M,i,i))) gcoeff(M,i,i) = gen_1;
1034 455050 : return M;
1035 : }
1036 :
1037 : /* Return Z-basis for Z[a] + U(a)/p Z[a] in Z[t]/(f), mf = v_p(disc f), U
1038 : * a ZX. Special cases: a = t is coded as NULL, U = 0 is coded as NULL */
1039 : static GEN
1040 621497 : dbasis(GEN p, GEN f, long mf, GEN a, GEN U)
1041 : {
1042 621497 : long n = degpol(f), i, dU;
1043 : GEN b, h;
1044 :
1045 621497 : if (n == 1) return matid(1);
1046 621497 : if (a && gequalX(a)) a = NULL;
1047 621497 : if (DEBUGLEVEL>5)
1048 : {
1049 0 : err_printf(" entering Dedekind Basis with parameters p=%Ps\n",p);
1050 0 : err_printf(" f = %Ps,\n a = %Ps\n",f, a? a: pol_x(varn(f)));
1051 : }
1052 621497 : if (a)
1053 : {
1054 200097 : GEN pd = powiu(p, mf >> 1);
1055 200097 : GEN da, pdp = mulii(pd,p), D = pdp;
1056 : long vda;
1057 200097 : dU = U ? degpol(U): 0;
1058 200097 : b = cgetg(n+1, t_MAT);
1059 200097 : h = scalarpol(pd, varn(f));
1060 200097 : a = QpX_remove_denom(a, p, &da, &vda);
1061 200097 : if (da) D = mulii(D, da);
1062 200097 : gel(b,1) = scalarcol_shallow(pd, n);
1063 570222 : for (i=2; i<=n; i++)
1064 : {
1065 370125 : if (i == dU+1)
1066 0 : h = compmod(p, U, mkvec3(a,da,stoi(vda)), f, pdp, (mf>>1) - 1);
1067 : else
1068 : {
1069 370125 : h = FpXQ_mul(h, a, f, D);
1070 370125 : if (da) h = ZX_Z_divexact(h, da);
1071 : }
1072 370125 : gel(b,i) = RgX_to_RgC(h,n);
1073 : }
1074 200097 : return ZpM_hnfmodid(b, p, pd);
1075 : }
1076 : else
1077 : {
1078 421400 : if (!U) return matid(n);
1079 421400 : dU = degpol(U);
1080 421400 : if (dU == n) return matid(n);
1081 421400 : U = FpX_normalize(U, p);
1082 421400 : b = cgetg(n+1, t_MAT);
1083 1634461 : for (i = 1; i <= dU; i++) gel(b,i) = vec_ei(n, i);
1084 421400 : h = RgX_Rg_div(U, p);
1085 474119 : for ( ; i <= n; i++)
1086 : {
1087 474119 : gel(b, i) = RgX_to_RgC(h,n);
1088 474119 : if (i == n) break;
1089 52719 : h = RgX_shift_shallow(h,1);
1090 : }
1091 421400 : return b;
1092 : }
1093 : }
1094 :
1095 : static GEN
1096 509906 : get_partial_order_as_pols(GEN p, GEN f)
1097 : {
1098 509906 : GEN O = maxord(p, f, -1);
1099 509906 : long v = varn(f);
1100 509906 : return O == gen_1? pol_x_powers(degpol(f), v): RgM_to_RgXV(O, v);
1101 : }
1102 :
1103 : static long
1104 2292 : p_is_prime(decomp_t *S)
1105 : {
1106 2292 : if (S->pisprime < 0) S->pisprime = BPSW_psp(S->p);
1107 2292 : return S->pisprime;
1108 : }
1109 : static GEN ZpX_monic_factor_squarefree(GEN f, GEN p, long prec);
1110 :
1111 : /* if flag = 0, maximal order, else factorization to precision r = flag */
1112 : static GEN
1113 255254 : Decomp(decomp_t *S, long flag)
1114 : {
1115 255254 : pari_sp av = avma;
1116 : GEN fred, pr2, pr, pk, ph2, ph, b1, b2, a, e, de, f1, f2, dt, th, chip;
1117 255254 : GEN p = S->p;
1118 255254 : long vde, vdt, k, r = maxss(flag, 2*S->df + 1);
1119 :
1120 255254 : if (DEBUGLEVEL>5) err_printf(" entering Decomp: %Ps^%ld\n f = %Ps\n",
1121 : p, r, S->f);
1122 255254 : else if (DEBUGLEVEL>2) err_printf(" entering Decomp\n");
1123 255254 : chip = FpX_red(S->chi, p);
1124 255254 : if (!FpX_valrem(chip, S->nu, p, &b1))
1125 : {
1126 0 : if (!p_is_prime(S)) pari_err_PRIME("Decomp",p);
1127 0 : pari_err_BUG("Decomp (not a factor)");
1128 : }
1129 255254 : b2 = FpX_div(chip, b1, p);
1130 255254 : a = FpX_mul(FpXQ_inv(b2, b1, p), b2, p);
1131 : /* E = e / de, e in Z[X], de in Z, E = a(phi) mod (f, p) */
1132 255254 : th = QpX_remove_denom(S->phi, p, &dt, &vdt);
1133 255254 : if (dt)
1134 : {
1135 122925 : long dega = degpol(a);
1136 122925 : vde = dega * vdt;
1137 122925 : de = powiu(dt, dega);
1138 122925 : pr = mulii(p, de);
1139 122925 : a = FpX_rescale(a, dt, pr);
1140 : }
1141 : else
1142 : {
1143 132329 : vde = 0;
1144 132329 : de = gen_1;
1145 132329 : pr = p;
1146 : }
1147 255254 : e = FpX_FpXQ_eval(a, th, S->f, pr);
1148 255254 : update_den(p, &e, &de, &vde, NULL);
1149 :
1150 255254 : pk = p; k = 1;
1151 : /* E, (1 - E) tend to orthogonal idempotents in Zp[X]/(f) */
1152 1181519 : while (k < r + vde)
1153 : { /* E <-- E^2(3-2E) mod p^2k, with E = e/de */
1154 : GEN D;
1155 926265 : pk = sqri(pk); k <<= 1;
1156 926265 : e = ZX_mul(ZX_sqr(e), Z_ZX_sub(mului(3,de), gmul2n(e,1)));
1157 926265 : de= mulii(de, sqri(de));
1158 926265 : vde *= 3;
1159 926265 : D = mulii(pk, de);
1160 926265 : e = FpX_rem(e, centermod(S->f, D), D); /* e/de defined mod pk */
1161 926265 : update_den(p, &e, &de, &vde, NULL);
1162 : }
1163 : /* required precision of the factors */
1164 255254 : pr = powiu(p, r); pr2 = shifti(pr, -1);
1165 255254 : ph = mulii(de,pr);ph2 = shifti(ph, -1);
1166 255254 : e = FpX_center_i(FpX_red(e, ph), ph, ph2);
1167 255254 : fred = FpX_red(S->f, ph);
1168 :
1169 255254 : f1 = ZpX_gcd(fred, Z_ZX_sub(de, e), p, ph); /* p-adic gcd(f, 1-e) */
1170 255254 : if (!is_pm1(de))
1171 : {
1172 122925 : fred = FpX_red(fred, pr);
1173 122925 : f1 = FpX_red(f1, pr);
1174 : }
1175 255254 : f2 = FpX_div(fred,f1, pr);
1176 255254 : f1 = FpX_center_i(f1, pr, pr2);
1177 255254 : f2 = FpX_center_i(f2, pr, pr2);
1178 :
1179 255254 : if (DEBUGLEVEL>5)
1180 0 : err_printf(" leaving Decomp: f1 = %Ps\nf2 = %Ps\ne = %Ps\nde= %Ps\n", f1,f2,e,de);
1181 :
1182 255254 : if (flag < 0)
1183 : {
1184 0 : GEN m = vconcat(ZpX_primedec(f1, p), ZpX_primedec(f2, p));
1185 0 : return sort_factor(m, (void*)&cmpii, &cmp_nodata);
1186 : }
1187 255254 : else if (flag)
1188 : {
1189 301 : (void)gc_all(av, 2, &f1, &f2);
1190 301 : return shallowconcat(ZpX_monic_factor_squarefree(f1, p, flag),
1191 : ZpX_monic_factor_squarefree(f2, p, flag));
1192 : } else {
1193 : GEN D, d1, d2, B1, B2, M;
1194 : long n, n1, n2, i;
1195 254953 : (void)gc_all(av, 4, &f1, &f2, &e, &de);
1196 254953 : D = de;
1197 254953 : B1 = get_partial_order_as_pols(p,f1); n1 = lg(B1)-1;
1198 254953 : B2 = get_partial_order_as_pols(p,f2); n2 = lg(B2)-1; n = n1+n2;
1199 254953 : d1 = QpXV_denom(B1);
1200 254953 : d2 = QpXV_denom(B2); if (cmpii(d1, d2) < 0) d1 = d2;
1201 254953 : if (d1 != gen_1) {
1202 157511 : B1 = Q_muli_to_int(B1, d1);
1203 157511 : B2 = Q_muli_to_int(B2, d1);
1204 157511 : D = mulii(d1, D);
1205 : }
1206 254953 : fred = centermod_i(S->f, D, shifti(D,-1));
1207 254953 : M = cgetg(n+1, t_MAT);
1208 821813 : for (i=1; i<=n1; i++)
1209 566860 : gel(M,i) = RgX_to_RgC(FpX_rem(FpX_mul(gel(B1,i),e,D), fred, D), n);
1210 254953 : e = Z_ZX_sub(de, e); B2 -= n1;
1211 711246 : for ( ; i<=n; i++)
1212 456293 : gel(M,i) = RgX_to_RgC(FpX_rem(FpX_mul(gel(B2,i),e,D), fred, D), n);
1213 254953 : return ZpM_hnfmodid(M, p, D);
1214 : }
1215 : }
1216 :
1217 : /* minimum extension valuation: L/E */
1218 : static void
1219 625431 : vstar(GEN p,GEN h, long *L, long *E)
1220 : {
1221 625431 : long first, j, k, v, w, m = degpol(h);
1222 :
1223 625431 : first = 1; k = 1; v = 0;
1224 2587714 : for (j=1; j<=m; j++)
1225 : {
1226 1962283 : GEN c = gel(h, m-j+2);
1227 1962283 : if (signe(c))
1228 : {
1229 1887254 : w = Z_pval(c,p);
1230 1887254 : if (first || w*k < v*j) { v = w; k = j; }
1231 1887254 : first = 0;
1232 : }
1233 : }
1234 : /* v/k = min_j ( v_p(h_{m-j}) / j ) */
1235 625431 : w = (long)ugcd(v,k);
1236 625431 : *L = v/w;
1237 625431 : *E = k/w;
1238 625431 : }
1239 :
1240 : static GEN
1241 64666 : redelt_i(GEN a, GEN N, GEN p, GEN *pda, long *pvda)
1242 : {
1243 : GEN z;
1244 64666 : a = Q_remove_denom(a, pda);
1245 64666 : *pvda = 0;
1246 64666 : if (*pda)
1247 : {
1248 64666 : long v = Z_pvalrem(*pda, p, &z);
1249 64666 : if (v) {
1250 64666 : *pda = powiu(p, v);
1251 64666 : *pvda = v;
1252 64666 : N = mulii(*pda, N);
1253 : }
1254 : else
1255 0 : *pda = NULL;
1256 64666 : if (!is_pm1(z)) a = ZX_Z_mul(a, Fp_inv(z, N));
1257 : }
1258 64666 : return centermod(a, N);
1259 : }
1260 : /* reduce the element a modulo N [ a power of p ], taking first care of the
1261 : * denominators */
1262 : static GEN
1263 48702 : redelt(GEN a, GEN N, GEN p)
1264 : {
1265 : GEN da;
1266 : long vda;
1267 48702 : a = redelt_i(a, N, p, &da, &vda);
1268 48702 : if (da) a = RgX_Rg_div(a, da);
1269 48702 : return a;
1270 : }
1271 :
1272 : /* compute the c first Newton sums modulo pp of the
1273 : characteristic polynomial of a/d mod chi, d > 0 power of p (NULL = gen_1),
1274 : a, chi in Zp[X], vda = v_p(da)
1275 : ns = Newton sums of chi */
1276 : static GEN
1277 708014 : newtonsums(GEN p, GEN a, GEN da, long vda, GEN chi, long c, GEN pp, GEN ns)
1278 : {
1279 : GEN va, pa, dpa, s;
1280 708014 : long j, k, vdpa, lns = lg(ns);
1281 : pari_sp av;
1282 :
1283 708014 : a = centermod(a, pp); av = avma;
1284 708014 : dpa = pa = NULL; /* -Wall */
1285 708014 : vdpa = 0;
1286 708014 : va = zerovec(c);
1287 2927872 : for (j = 1; j <= c; j++)
1288 : { /* pa/dpa = (a/d)^(j-1) mod (chi, pp), dpa = p^vdpa */
1289 : long l;
1290 2226643 : pa = j == 1? a: FpXQ_mul(pa, a, chi, pp);
1291 2226643 : l = lg(pa); if (l == 2) break;
1292 2226643 : if (lns < l) l = lns;
1293 :
1294 2226643 : if (da) {
1295 2091285 : dpa = j == 1? da: mulii(dpa, da);
1296 2091285 : vdpa += vda;
1297 2091285 : update_den(p, &pa, &dpa, &vdpa, &pp);
1298 : }
1299 2226643 : s = mulii(gel(pa,2), gel(ns,2)); /* k = 2 */
1300 11075313 : for (k = 3; k < l; k++) s = addii(s, mulii(gel(pa,k), gel(ns,k)));
1301 2226643 : if (da) {
1302 : GEN r;
1303 2091285 : s = dvmdii(s, dpa, &r);
1304 2091285 : if (r != gen_0) return NULL;
1305 : }
1306 2219858 : gel(va,j) = centermodii(s, pp, shifti(pp,-1));
1307 :
1308 2219858 : if (gc_needed(av, 1))
1309 : {
1310 7 : if(DEBUGMEM>1) pari_warn(warnmem, "newtonsums");
1311 7 : (void)gc_all(av, dpa?4:3, &pa, &va, &pp, &dpa);
1312 : }
1313 : }
1314 701229 : for (; j <= c; j++) gel(va,j) = gen_0;
1315 701229 : return va;
1316 : }
1317 :
1318 : /* compute the characteristic polynomial of a/da mod chi (a in Z[X]), given
1319 : * by its Newton sums to a precision of pp using Newton sums */
1320 : static GEN
1321 701229 : newtoncharpoly(GEN pp, GEN p, GEN NS)
1322 : {
1323 701229 : long n = lg(NS)-1, j, k;
1324 701229 : GEN c = cgetg(n + 2, t_VEC), pp2 = shifti(pp,-1);
1325 :
1326 701229 : gel(c,1) = (n & 1 ? gen_m1: gen_1);
1327 2910581 : for (k = 2; k <= n+1; k++)
1328 : {
1329 2209450 : pari_sp av2 = avma;
1330 2209450 : GEN s = gen_0;
1331 : ulong z;
1332 2209450 : long v = u_pvalrem(k - 1, p, &z);
1333 9361761 : for (j = 1; j < k; j++)
1334 : {
1335 7152311 : GEN t = mulii(gel(NS,j), gel(c,k-j));
1336 7152311 : if (!odd(j)) t = negi(t);
1337 7152311 : s = addii(s, t);
1338 : }
1339 2209450 : if (v) {
1340 846243 : s = gdiv(s, powiu(p, v));
1341 846243 : if (typ(s) != t_INT) return NULL;
1342 : }
1343 2209352 : s = mulii(s, Fp_inv(utoipos(z), pp));
1344 2209352 : gel(c,k) = gc_INT(av2, Fp_center_i(s, pp, pp2));
1345 : }
1346 1869696 : for (k = odd(n)? 1: 2; k <= n+1; k += 2) gel(c,k) = negi(gel(c,k));
1347 701131 : return gtopoly(c, 0);
1348 : }
1349 :
1350 : static void
1351 708014 : manage_cache(decomp_t *S, GEN f, GEN pp)
1352 : {
1353 708014 : GEN t = S->precns;
1354 :
1355 708014 : if (!t) t = mulii(S->pmf, powiu(S->p, S->df));
1356 708014 : if (cmpii(t, pp) < 0) t = pp;
1357 :
1358 708014 : if (!S->precns || !RgX_equal(f, S->nsf) || cmpii(S->precns, t) < 0)
1359 : {
1360 522574 : if (DEBUGLEVEL>4)
1361 0 : err_printf(" Precision for cached Newton sums for %Ps: %Ps -> %Ps\n",
1362 0 : f, S->precns? S->precns: gen_0, t);
1363 522574 : S->nsf = f;
1364 522574 : S->ns = FpX_Newton(f, degpol(f), t);
1365 522574 : S->precns = t;
1366 : }
1367 708014 : }
1368 :
1369 : /* return NULL if a mod f is not an integer
1370 : * The denominator of any integer in Zp[X]/(f) divides pdr */
1371 : static GEN
1372 708014 : mycaract(decomp_t *S, GEN f, GEN a, GEN pp, GEN pdr)
1373 : {
1374 : pari_sp av;
1375 : GEN d, chi, prec1, prec2, prec3, ns;
1376 708014 : long vd, n = degpol(f);
1377 :
1378 708014 : if (gequal0(a)) return pol_0(varn(f));
1379 :
1380 708014 : a = QpX_remove_denom(a, S->p, &d, &vd);
1381 708014 : prec1 = pp;
1382 708014 : if (lgefint(S->p) == 3)
1383 707962 : prec1 = mulii(prec1, powiu(S->p, factorial_lval(n, itou(S->p))));
1384 708014 : if (d)
1385 : {
1386 643167 : GEN p1 = powiu(d, n);
1387 643167 : prec2 = mulii(prec1, p1);
1388 643167 : prec3 = mulii(prec1, gmin_shallow(mulii(p1, d), pdr));
1389 : }
1390 : else
1391 64847 : prec2 = prec3 = prec1;
1392 708014 : manage_cache(S, f, prec3);
1393 :
1394 708014 : av = avma;
1395 708014 : ns = newtonsums(S->p, a, d, vd, f, n, prec2, S->ns);
1396 708014 : if (!ns) return NULL;
1397 701229 : chi = newtoncharpoly(prec1, S->p, ns);
1398 701229 : if (!chi) return NULL;
1399 701131 : setvarn(chi, varn(f));
1400 701131 : return gc_upto(av, centermod(chi, pp));
1401 : }
1402 :
1403 : static GEN
1404 642684 : get_nu(GEN chi, GEN p, long *ptl)
1405 : { /* split off powers of x first for efficiency */
1406 642684 : long v = ZX_valrem(FpX_red(chi,p), &chi), n;
1407 : GEN P;
1408 642684 : if (!degpol(chi)) { *ptl = 1; return pol_x(varn(chi)); }
1409 476911 : P = gel(FpX_factor(chi,p), 1); n = lg(P)-1;
1410 476911 : *ptl = v? n+1: n; return gel(P,n);
1411 : }
1412 :
1413 : /* Factor characteristic polynomial chi of phi mod p. If it splits, update
1414 : * S->{phi, chi, nu} and return 1. In any case, set *nu to an irreducible
1415 : * factor mod p of chi */
1416 : static int
1417 481769 : split_char(decomp_t *S, GEN chi, GEN phi, GEN phi0, GEN *nu)
1418 : {
1419 : long l;
1420 481769 : *nu = get_nu(chi, S->p, &l);
1421 481769 : if (l == 1) return 0; /* single irreducible factor: doesn't split */
1422 : /* phi o phi0 mod (p, f) */
1423 122925 : S->phi = compmod(S->p, phi, phi0, S->f, S->p, 0);
1424 122925 : S->chi = chi;
1425 122925 : S->nu = *nu; return 1;
1426 : }
1427 :
1428 : /* Return the prime element in Zp[phi], a t_INT (iff *Ep = 1) or QX;
1429 : * nup, chip are ZX. phi = NULL codes X
1430 : * If *Ep < oE or Ep divides Ediv (!=0) return NULL (uninteresting) */
1431 : static GEN
1432 564749 : getprime(decomp_t *S, GEN phi, GEN chip, GEN nup, long *Lp, long *Ep,
1433 : long oE, long Ediv)
1434 : {
1435 : GEN z, chin, q, qp;
1436 : long r, s;
1437 :
1438 564749 : if (phi && dvdii(constant_coeff(chip), S->psc))
1439 : {
1440 1700 : chip = mycaract(S, S->chi, phi, S->pmf, S->prc);
1441 1700 : if (dvdii(constant_coeff(chip), S->pmf))
1442 1246 : chip = ZXQ_charpoly(phi, S->chi, varn(chip));
1443 : }
1444 564749 : if (degpol(nup) == 1)
1445 : {
1446 525542 : GEN c = gel(nup,2); /* nup = X + c */
1447 525542 : chin = signe(c)? RgX_Rg_translate(chip, negi(c)): chip;
1448 : }
1449 : else
1450 39207 : chin = ZXQ_charpoly(nup, chip, varn(chip));
1451 :
1452 564749 : vstar(S->p, chin, Lp, Ep);
1453 564749 : if (*Ep < oE || (Ediv && Ediv % *Ep == 0)) return NULL;
1454 :
1455 443690 : if (*Ep == 1) return S->p;
1456 305914 : (void)cbezout(*Lp, -*Ep, &r, &s); /* = 1 */
1457 305914 : if (r <= 0)
1458 : {
1459 60094 : long t = 1 + ((-r) / *Ep);
1460 60094 : r += t * *Ep;
1461 60094 : s += t * *Lp;
1462 : }
1463 : /* r > 0 minimal such that r L/E - s = 1/E
1464 : * pi = nu^r / p^s is an element of valuation 1/E,
1465 : * so is pi + O(p) since 1/E < 1. May compute nu^r mod p^(s+1) */
1466 305914 : q = powiu(S->p, s); qp = mulii(q, S->p);
1467 305914 : nup = FpXQ_powu(nup, r, S->chi, qp);
1468 305914 : if (!phi) return RgX_Rg_div(nup, q); /* phi = X : no composition */
1469 48702 : z = compmod(S->p, nup, phi, S->chi, qp, -s);
1470 48702 : return signe(z)? z: NULL;
1471 : }
1472 :
1473 : static int
1474 277033 : update_phi(decomp_t *S)
1475 : {
1476 277033 : GEN PHI = NULL, prc, psc, X = pol_x(varn(S->f));
1477 : long k, vpsc;
1478 277033 : for (k = 1;; k++)
1479 : {
1480 279399 : prc = ZpX_reduced_resultant_fast(S->chi, ZX_deriv(S->chi), S->p, S->vpsc);
1481 : /* if prc == S->psc then either chi is not separable or
1482 : the reduced discriminant of chi is too large */
1483 279399 : if (!equalii(prc, S->psc)) break;
1484 :
1485 : /* increase precision */
1486 2366 : S->vpsc = maxss(S->vpsf, S->vpsc + 1);
1487 2366 : S->psc = (S->vpsc == S->vpsf)? S->psf: mulii(S->psc, S->p);
1488 :
1489 2366 : PHI = S->phi;
1490 2366 : if (S->phi0) PHI = compmod(S->p, PHI, S->phi0, S->f, S->psc, 0);
1491 : /* change phi (in case not separable) */
1492 2366 : PHI = gadd(PHI, ZX_Z_mul(X, mului(k, S->p)));
1493 2366 : S->chi = mycaract(S, S->f, PHI, S->psc, S->pdf);
1494 : }
1495 277033 : psc = mulii(sqri(prc), S->p);
1496 277033 : vpsc = 2*Z_pval(prc, S->p) + 1;
1497 :
1498 277033 : if (!PHI) /* break out of above loop immediately (k = 1) */
1499 : {
1500 274667 : PHI = S->phi;
1501 274667 : if (S->phi0) PHI = compmod(S->p, PHI, S->phi0, S->f, psc, 0);
1502 274667 : if (S->phi0 || cmpii(psc,S->psc) > 0)
1503 : {
1504 : for(;;)
1505 : {
1506 114333 : S->chi = mycaract(S, S->f, PHI, psc, S->pdf);
1507 114333 : prc = ZpX_reduced_resultant_fast(S->chi, ZX_deriv(S->chi), S->p, vpsc);
1508 114333 : if (!equalii(prc, psc)) break;
1509 497 : psc = mulii(psc, S->p); vpsc++;
1510 : /* it can happen that S->chi is never squarefree: then change PHI */
1511 497 : if (vpsc > 2*S->mf) PHI = gadd(PHI, ZX_Z_mul(X, S->p));
1512 : }
1513 113836 : psc = mulii(sqri(prc), S->p);
1514 113836 : vpsc = 2*Z_pval(prc, S->p) + 1;
1515 : }
1516 : }
1517 277033 : S->phi = PHI;
1518 277033 : S->chi = FpX_red(S->chi, psc);
1519 :
1520 : /* may happen if p is unramified */
1521 277033 : if (is_pm1(prc)) return 0;
1522 232286 : S->prc = prc;
1523 232286 : S->psc = psc;
1524 232286 : S->vpsc = vpsc; return 1;
1525 : }
1526 :
1527 : /* return 1 if at least 2 factors mod p ==> chi splits
1528 : * Replace S->phi such that F increases (to D) */
1529 : static int
1530 67416 : testb2(decomp_t *S, long D, GEN theta)
1531 : {
1532 67416 : long v = varn(S->chi), dlim = degpol(S->chi)-1;
1533 67416 : GEN T0 = S->phi, chi, phi, nu;
1534 67416 : if (DEBUGLEVEL>4) err_printf(" Increasing Fa\n");
1535 : for (;;)
1536 : {
1537 67476 : phi = gadd(theta, random_FpX(dlim, v, S->p));
1538 67476 : chi = mycaract(S, S->chi, phi, S->psf, S->prc);
1539 : /* phi nonprimary ? */
1540 67476 : if (split_char(S, chi, phi, T0, &nu)) return 1;
1541 67476 : if (degpol(nu) == D) break;
1542 : }
1543 : /* F_phi=lcm(F_alpha, F_theta)=D and E_phi=E_alpha */
1544 67416 : S->phi0 = T0;
1545 67416 : S->chi = chi;
1546 67416 : S->phi = phi;
1547 67416 : S->nu = nu; return 0;
1548 : }
1549 :
1550 : /* return 1 if at least 2 factors mod p ==> chi can be split.
1551 : * compute a new S->phi such that E = lcm(Ea, Et);
1552 : * A a ZX, T a t_INT (iff Et = 1, probably impossible ?) or QX */
1553 : static int
1554 48702 : testc2(decomp_t *S, GEN A, long Ea, GEN T, long Et)
1555 : {
1556 48702 : GEN c, chi, phi, nu, T0 = S->phi;
1557 :
1558 48702 : if (DEBUGLEVEL>4) err_printf(" Increasing Ea\n");
1559 48702 : if (Et == 1) /* same as other branch, split for efficiency */
1560 0 : c = A; /* Et = 1 => s = 1, r = 0, t = 0 */
1561 : else
1562 : {
1563 : long r, s, t;
1564 48702 : (void)cbezout(Ea, Et, &r, &s); t = 0;
1565 48800 : while (r < 0) { r = r + Et; t++; }
1566 48912 : while (s < 0) { s = s + Ea; t++; }
1567 :
1568 : /* A^s T^r / p^t */
1569 48702 : c = RgXQ_mul(RgXQ_powu(A, s, S->chi), RgXQ_powu(T, r, S->chi), S->chi);
1570 48702 : c = RgX_Rg_div(c, powiu(S->p, t));
1571 48702 : c = redelt(c, S->psc, S->p);
1572 : }
1573 48702 : phi = RgX_add(c, pol_x(varn(S->chi)));
1574 48702 : chi = mycaract(S, S->chi, phi, S->psf, S->prc);
1575 48702 : if (split_char(S, chi, phi, T0, &nu)) return 1;
1576 : /* E_phi = lcm(E_alpha,E_theta) */
1577 48702 : S->phi0 = T0;
1578 48702 : S->chi = chi;
1579 48702 : S->phi = phi;
1580 48702 : S->nu = nu; return 0;
1581 : }
1582 :
1583 : /* Return h^(-degpol(P)) P(x * h) if result is integral, NULL otherwise */
1584 : static GEN
1585 59952 : ZX_rescale_inv(GEN P, GEN h)
1586 : {
1587 59952 : long i, l = lg(P);
1588 59952 : GEN Q = cgetg(l,t_POL), hi = h;
1589 59952 : gel(Q,l-1) = gel(P,l-1);
1590 174287 : for (i=l-2; i>=2; i--)
1591 : {
1592 : GEN r;
1593 174287 : gel(Q,i) = dvmdii(gel(P,i), hi, &r);
1594 174287 : if (signe(r)) return NULL;
1595 174287 : if (i == 2) break;
1596 114335 : hi = mulii(hi,h);
1597 : }
1598 59952 : Q[1] = P[1]; return Q;
1599 : }
1600 :
1601 : /* x p^-eq nu^-er mod p */
1602 : static GEN
1603 304149 : get_gamma(decomp_t *S, GEN x, long eq, long er)
1604 : {
1605 304149 : GEN q, g = x, Dg = powiu(S->p, eq);
1606 304149 : long vDg = eq;
1607 304149 : if (er)
1608 : {
1609 23135 : if (!S->invnu)
1610 : {
1611 15964 : while (gdvd(S->chi, S->nu)) S->nu = RgX_Rg_add(S->nu, S->p);
1612 15964 : S->invnu = QXQ_inv(S->nu, S->chi);
1613 15964 : S->invnu = redelt_i(S->invnu, S->psc, S->p, &S->Dinvnu, &S->vDinvnu);
1614 : }
1615 23135 : if (S->Dinvnu) {
1616 23135 : Dg = mulii(Dg, powiu(S->Dinvnu, er));
1617 23135 : vDg += er * S->vDinvnu;
1618 : }
1619 23135 : q = mulii(S->p, Dg);
1620 23135 : g = ZX_mul(g, FpXQ_powu(S->invnu, er, S->chi, q));
1621 23135 : g = FpX_rem(g, S->chi, q);
1622 23135 : update_den(S->p, &g, &Dg, &vDg, NULL);
1623 23135 : g = centermod(g, mulii(S->p, Dg));
1624 : }
1625 304149 : if (!is_pm1(Dg)) g = RgX_Rg_div(g, Dg);
1626 304149 : return g;
1627 : }
1628 : static GEN
1629 357218 : get_g(decomp_t *S, long Ea, long L, long E, GEN beta, GEN *pchig,
1630 : long *peq, long *per)
1631 : {
1632 : long eq, er;
1633 357218 : GEN g, chig, chib = NULL;
1634 : for(;;) /* at most twice */
1635 : {
1636 364101 : if (L < 0)
1637 : {
1638 60682 : chib = ZXQ_charpoly(beta, S->chi, varn(S->chi));
1639 60682 : vstar(S->p, chib, &L, &E);
1640 : }
1641 364101 : eq = L / E; er = L*Ea / E - eq*Ea;
1642 : /* floor(L Ea/E) = eq Ea + er */
1643 364101 : if (er || !chib)
1644 : { /* g might not be an integer ==> chig = NULL */
1645 304149 : g = get_gamma(S, beta, eq, er);
1646 304149 : chig = mycaract(S, S->chi, g, S->psc, S->prc);
1647 : }
1648 : else
1649 : { /* g = beta/p^eq, special case of the above */
1650 59952 : GEN h = powiu(S->p, eq);
1651 59952 : g = RgX_Rg_div(beta, h);
1652 59952 : chig = ZX_rescale_inv(chib, h); /* chib(x h) / h^N */
1653 59952 : if (chig) chig = FpX_red(chig, S->pmf);
1654 : }
1655 : /* either success or second consecutive failure */
1656 364101 : if (chig || chib) break;
1657 : /* if g fails the v*-test, v(beta) was wrong. Retry once */
1658 6883 : L = -1;
1659 : }
1660 357218 : *pchig = chig; *peq = eq; *per = er; return g;
1661 : }
1662 :
1663 : /* return 1 if at least 2 factors mod p ==> chi can be split */
1664 : static int
1665 239043 : loop(decomp_t *S, long Ea)
1666 : {
1667 239043 : pari_sp av = avma;
1668 239043 : GEN beta = FpXQ_powu(S->nu, Ea, S->chi, S->p);
1669 239043 : long N = degpol(S->f), v = varn(S->f);
1670 239043 : S->invnu = NULL;
1671 : for (;;)
1672 118175 : { /* beta tends to a factor of chi */
1673 : long L, i, Fg, eq, er;
1674 357218 : GEN chig = NULL, d, g, nug;
1675 :
1676 357218 : if (DEBUGLEVEL>4) err_printf(" beta = %Ps\n", beta);
1677 357218 : L = ZpX_resultant_val(S->chi, beta, S->p, S->mf+1);
1678 357218 : if (L > S->mf) L = -1; /* from scratch */
1679 357218 : g = get_g(S, Ea, L, N, beta, &chig, &eq, &er);
1680 357218 : if (DEBUGLEVEL>4) err_printf(" (eq,er) = (%ld,%ld)\n", eq,er);
1681 : /* g = beta p^-eq nu^-er (a unit), chig = charpoly(g) */
1682 475406 : if (split_char(S, chig, g,S->phi, &nug)) return 1;
1683 :
1684 236363 : Fg = degpol(nug);
1685 236363 : if (Fg == 1)
1686 : { /* frequent special case nug = x - d */
1687 : long Le, Ee;
1688 : GEN chie, nue, e, pie;
1689 160574 : d = negi(gel(nug,2));
1690 160574 : chie = RgX_Rg_translate(chig, d);
1691 160574 : nue = pol_x(v);
1692 160574 : e = RgX_Rg_sub(g, d);
1693 160574 : pie = getprime(S, e, chie, nue, &Le, &Ee, 0,Ea);
1694 160574 : if (pie) return testc2(S, S->nu, Ea, pie, Ee);
1695 : }
1696 : else
1697 : {
1698 75789 : long Fa = degpol(S->nu), vdeng;
1699 : GEN deng, numg, nume;
1700 79368 : if (Fa % Fg) return testb2(S, ulcm(Fa,Fg), g);
1701 : /* nu & nug irreducible mod p, deg nug | deg nu. To improve beta, look
1702 : * for a root d of nug in Fp[phi] such that v_p(g - d) > 0 */
1703 8373 : if (ZX_equal(nug, S->nu))
1704 6081 : d = pol_x(v);
1705 : else
1706 : {
1707 2292 : if (!p_is_prime(S)) pari_err_PRIME("FpX_ffisom",S->p);
1708 2292 : d = FpX_ffisom(nug, S->nu, S->p);
1709 : }
1710 : /* write g = numg / deng, e = nume / deng */
1711 8373 : numg = QpX_remove_denom(g, S->p, &deng, &vdeng);
1712 13218 : for (i = 1; i <= Fg; i++)
1713 : {
1714 : GEN chie, nue, e;
1715 13218 : if (i != 1) d = FpXQ_pow(d, S->p, S->nu, S->p); /* next root */
1716 13218 : nume = ZX_sub(numg, ZX_Z_mul(d, deng));
1717 : /* test e = nume / deng */
1718 13218 : if (ZpX_resultant_val(S->chi, nume, S->p, vdeng*N+1) <= vdeng*N)
1719 4845 : continue;
1720 8373 : e = RgX_Rg_div(nume, deng);
1721 8373 : chie = mycaract(S, S->chi, e, S->psc, S->prc);
1722 9882 : if (split_char(S, chie, e,S->phi, &nue)) return 1;
1723 6303 : if (RgX_is_monomial(nue))
1724 : { /* v_p(e) = v_p(g - d) > 0 */
1725 : long Le, Ee;
1726 : GEN pie;
1727 6303 : pie = getprime(S, e, chie, nue, &Le, &Ee, 0,Ea);
1728 6303 : if (pie) return testc2(S, S->nu, Ea, pie, Ee);
1729 4794 : break;
1730 : }
1731 : }
1732 4794 : if (i > Fg)
1733 : {
1734 0 : if (!p_is_prime(S)) pari_err_PRIME("nilord",S->p);
1735 0 : pari_err_BUG("nilord (no root)");
1736 : }
1737 : }
1738 118175 : if (eq) d = gmul(d, powiu(S->p, eq));
1739 118175 : if (er) d = gmul(d, gpowgs(S->nu, er));
1740 118175 : beta = gsub(beta, d);
1741 :
1742 118175 : if (gc_needed(av,1))
1743 : {
1744 0 : if (DEBUGMEM > 1) pari_warn(warnmem, "nilord");
1745 0 : (void)gc_all(av, S->invnu? 6: 4, &beta, &(S->precns), &(S->ns), &(S->nsf), &(S->invnu), &(S->Dinvnu));
1746 : }
1747 : }
1748 : }
1749 :
1750 : /* E and F cannot decrease; return 1 if O = Zp[phi], 2 if we can get a
1751 : * decomposition and 0 otherwise */
1752 : static long
1753 394988 : progress(decomp_t *S, GEN *ppa, long *pE)
1754 : {
1755 394988 : long E = *pE, F;
1756 394988 : GEN pa = *ppa;
1757 394988 : S->phi0 = NULL; /* no delayed composition */
1758 : for(;;)
1759 2884 : {
1760 : long l, La, Ea; /* N.B If E = 0, getprime cannot return NULL */
1761 397872 : GEN pia = getprime(S, NULL, S->chi, S->nu, &La, &Ea, E,0);
1762 397872 : if (pia) { /* success, we break out in THIS loop */
1763 394988 : pa = (typ(pia) == t_POL)? RgX_RgXQ_eval(pia, S->phi, S->f): pia;
1764 394988 : E = Ea;
1765 394988 : if (La == 1) break; /* no need to change phi so that nu = pia */
1766 : }
1767 : /* phi += prime elt */
1768 65095 : S->phi = typ(pa) == t_INT? RgX_Rg_add_shallow(S->phi, pa)
1769 160915 : : RgX_add(S->phi, pa);
1770 : /* recompute char. poly. chi from scratch */
1771 160915 : S->chi = mycaract(S, S->f, S->phi, S->psf, S->pdf);
1772 160915 : S->nu = get_nu(S->chi, S->p, &l);
1773 160915 : if (l > 1) return 2;
1774 160915 : if (!update_phi(S)) return 1; /* unramified */
1775 160915 : if (pia) break;
1776 : }
1777 394988 : *pE = E; *ppa = pa; F = degpol(S->nu);
1778 394988 : if (DEBUGLEVEL>4) err_printf(" (E, F) = (%ld,%ld)\n", E, F);
1779 394988 : if (E * F == degpol(S->f)) return 1;
1780 239043 : if (loop(S, E)) return 2;
1781 116118 : if (!update_phi(S)) return 1;
1782 71371 : return 0;
1783 : }
1784 :
1785 : /* flag != 0 iff we're looking for the p-adic factorization,
1786 : in which case it is the p-adic precision we want */
1787 : static GEN
1788 455946 : maxord_i(decomp_t *S, GEN p, GEN f, long mf, GEN w, long flag)
1789 : {
1790 455946 : long oE, n = lg(w)-1; /* factor of largest degree */
1791 455946 : GEN opa, D = ZpX_reduced_resultant_fast(f, ZX_deriv(f), p, mf);
1792 455946 : S->pisprime = -1;
1793 455946 : S->p = p;
1794 455946 : S->mf = mf;
1795 455946 : S->nu = gel(w,n);
1796 455946 : S->df = Z_pval(D, p);
1797 455946 : S->pdf = powiu(p, S->df);
1798 455946 : S->phi = pol_x(varn(f));
1799 455946 : S->chi = S->f = f;
1800 455946 : if (n > 1) return Decomp(S, flag); /* FIXME: use bezout_lift_fact */
1801 :
1802 323617 : if (DEBUGLEVEL>4)
1803 0 : err_printf(" entering Nilord: %Ps^%ld\n f = %Ps, nu = %Ps\n",
1804 : p, S->df, S->f, S->nu);
1805 323617 : else if (DEBUGLEVEL>2) err_printf(" entering Nilord\n");
1806 323617 : S->psf = S->psc = mulii(sqri(D), p);
1807 323617 : S->vpsf = S->vpsc = 2*S->df + 1;
1808 323617 : S->prc = mulii(D, p);
1809 323617 : S->chi = FpX_red(S->f, S->psc);
1810 323617 : S->pmf = powiu(p, S->mf+1);
1811 323617 : S->precns = NULL;
1812 323617 : for(opa = NULL, oE = 0;;)
1813 71371 : {
1814 394988 : long n = progress(S, &opa, &oE);
1815 394988 : if (n == 1) return flag? NULL: dbasis(p, S->f, S->mf, S->phi, S->chi);
1816 194296 : if (n == 2) return Decomp(S, flag);
1817 : }
1818 : }
1819 :
1820 : static int
1821 889 : expo_is_squarefree(GEN e)
1822 : {
1823 889 : long i, l = lg(e);
1824 1260 : for (i=1; i<l; i++)
1825 1022 : if (e[i] != 1) return 0;
1826 238 : return 1;
1827 : }
1828 : /* pure round 4 */
1829 : static GEN
1830 896 : ZpX_round4(GEN f, GEN p, GEN w, long prec)
1831 : {
1832 : decomp_t S;
1833 896 : GEN L = maxord_i(&S, p, f, ZpX_disc_val(f,p), w, prec);
1834 896 : return L? L: mkvec(f);
1835 : }
1836 : /* f a squarefree ZX with leading_coeff 1, degree > 0. Return list of
1837 : * irreducible factors in Zp[X] (computed mod p^prec) */
1838 : static GEN
1839 1155 : ZpX_monic_factor_squarefree(GEN f, GEN p, long prec)
1840 : {
1841 1155 : pari_sp av = avma;
1842 : GEN L, fa, w, e;
1843 : long i, l;
1844 1155 : if (degpol(f) == 1) return mkvec(f);
1845 889 : fa = FpX_factor(f,p); w = gel(fa,1); e = gel(fa,2);
1846 : /* no repeated factors: Hensel lift */
1847 889 : if (expo_is_squarefree(e)) return ZpX_liftfact(f, w, powiu(p,prec), p, prec);
1848 651 : l = lg(w);
1849 651 : if (l == 2)
1850 : {
1851 392 : L = ZpX_round4(f,p,w,prec);
1852 392 : if (lg(L) == 2) { set_avma(av); return mkvec(f); }
1853 : }
1854 : else
1855 : { /* >= 2 factors mod p: partial Hensel lift */
1856 259 : GEN D = ZpX_reduced_resultant_fast(f, ZX_deriv(f), p, ZpX_disc_val(f,p));
1857 259 : long r = maxss(2*Z_pval(D,p)+1, prec);
1858 259 : GEN W = cgetg(l, t_VEC);
1859 833 : for (i = 1; i < l; i++)
1860 574 : gel(W,i) = e[i] == 1? gel(w,i): FpX_powu(gel(w,i), e[i], p);
1861 259 : L = ZpX_liftfact(f, W, powiu(p,r), p, r);
1862 833 : for (i = 1; i < l; i++)
1863 574 : gel(L,i) = e[i] == 1? mkvec(gel(L,i))
1864 574 : : ZpX_round4(gel(L,i), p, mkvec(gel(w,i)), prec);
1865 259 : L = shallowconcat1(L);
1866 : }
1867 399 : return gc_GEN(av, L);
1868 : }
1869 :
1870 : /* assume T a ZX with leading_coeff 1, degree > 0 */
1871 : GEN
1872 546 : ZpX_monic_factor(GEN T, GEN p, long prec)
1873 : {
1874 : GEN Q, P, E, F;
1875 : long L, l, i, v;
1876 :
1877 546 : if (degpol(T) == 1) return mkmat2(mkcol(T), mkcol(gen_1));
1878 546 : v = ZX_valrem(T, &T);
1879 546 : Q = ZX_squff(T, &F); l = lg(Q); L = v? l + 1: l;
1880 546 : P = cgetg(L, t_VEC);
1881 546 : E = cgetg(L, t_VEC);
1882 1099 : for (i = 1; i < l; i++)
1883 : {
1884 553 : GEN w = ZpX_monic_factor_squarefree(gel(Q,i), p, prec);
1885 553 : gel(P,i) = w; settyp(w, t_COL);
1886 553 : gel(E,i) = const_col(lg(w)-1, utoipos(F[i]));
1887 : }
1888 546 : if (v) { gel(P,i) = pol_x(varn(T)); gel(E,i) = utoipos(v); }
1889 546 : return mkmat2(shallowconcat1(P), shallowconcat1(E));
1890 : }
1891 :
1892 : /* DT = multiple of disc(T) or NULL
1893 : * Return a multiple of the denominator of an algebraic integer (in Q[X]/(T))
1894 : * when expressed in terms of the power basis */
1895 : GEN
1896 44126 : indexpartial(GEN T, GEN DT)
1897 : {
1898 44126 : pari_sp av = avma;
1899 : long i, nb;
1900 44126 : GEN fa, E, P, U, res = gen_1, dT = ZX_deriv(T);
1901 :
1902 44126 : if (!DT) DT = ZX_disc(T);
1903 44126 : fa = absZ_factor_limit_strict(DT, 0, &U);
1904 44126 : P = gel(fa,1);
1905 44126 : E = gel(fa,2); nb = lg(P)-1;
1906 212162 : for (i = 1; i <= nb; i++)
1907 : {
1908 168036 : long e = itou(gel(E,i)), e2 = e >> 1;
1909 168036 : GEN p = gel(P,i), q = p;
1910 168036 : if (e2 >= 2) q = ZpX_reduced_resultant_fast(T, dT, p, e2);
1911 168036 : res = mulii(res, q);
1912 : }
1913 44126 : if (U)
1914 : {
1915 1923 : long e = itou(gel(U,2)), e2 = e >> 1;
1916 1923 : GEN p = gel(U,1), q = powiu(p, odd(e)? e2+1: e2);
1917 1923 : res = mulii(res, q);
1918 : }
1919 44126 : return gc_INT(av,res);
1920 : }
1921 :
1922 : /*******************************************************************/
1923 : /* */
1924 : /* 2-ELT REPRESENTATION FOR PRIME IDEALS (dividing index) */
1925 : /* */
1926 : /*******************************************************************/
1927 : /* to compute norm of elt in basis form */
1928 : typedef struct {
1929 : long r1;
1930 : GEN M; /* via embed_norm */
1931 :
1932 : GEN D, w, T; /* via resultant if M = NULL */
1933 : } norm_S;
1934 :
1935 : static GEN
1936 505223 : get_norm(norm_S *S, GEN a)
1937 : {
1938 505223 : if (S->M)
1939 : {
1940 : long e;
1941 503817 : GEN N = grndtoi( embed_norm(RgM_RgC_mul(S->M, a), S->r1), &e );
1942 503817 : if (e > -5) pari_err_PREC( "get_norm");
1943 503817 : return N;
1944 : }
1945 1406 : if (S->w) a = RgV_RgC_mul(S->w, a);
1946 1406 : return ZX_resultant_all(S->T, a, S->D, 0);
1947 : }
1948 : static void
1949 216656 : init_norm(norm_S *S, GEN nf, GEN p)
1950 : {
1951 216656 : GEN T = nf_get_pol(nf), M = nf_get_M(nf);
1952 216656 : long N = degpol(T), ex = gexpo(M) + gexpo(mului(8 * N, p));
1953 :
1954 216656 : S->r1 = nf_get_r1(nf);
1955 216656 : if (N * ex <= gprecision(M) - 20)
1956 : { /* enough prec to use embed_norm */
1957 216465 : S->M = M;
1958 216465 : S->D = NULL;
1959 216465 : S->w = NULL;
1960 216465 : S->T = NULL;
1961 : }
1962 : else
1963 : {
1964 191 : GEN w = leafcopy(nf_get_zkprimpart(nf)), D = nf_get_zkden(nf), Dp = sqri(p);
1965 : long i;
1966 191 : if (!equali1(D))
1967 : {
1968 191 : GEN w1 = D;
1969 191 : long v = Z_pval(D, p);
1970 191 : D = powiu(p, v);
1971 191 : Dp = mulii(D, Dp);
1972 191 : gel(w, 1) = remii(w1, Dp);
1973 : }
1974 3969 : for (i=2; i<=N; i++) gel(w,i) = FpX_red(gel(w,i), Dp);
1975 191 : S->M = NULL;
1976 191 : S->D = D;
1977 191 : S->w = w;
1978 191 : S->T = T;
1979 : }
1980 216656 : }
1981 : /* f = f(pr/p), q = p^(f+1), a in pr.
1982 : * Return 1 if v_pr(a) = 1, and 0 otherwise */
1983 : static int
1984 505223 : is_uniformizer(GEN a, GEN q, norm_S *S) { return !dvdii(get_norm(S,a), q); }
1985 :
1986 : /* Return x * y, x, y are t_MAT (Fp-basis of in O_K/p), assume (x,y)=1.
1987 : * Either x or y may be NULL (= O_K), not both */
1988 : static GEN
1989 711843 : mul_intersect(GEN x, GEN y, GEN p)
1990 : {
1991 711843 : if (!x) return y;
1992 375864 : if (!y) return x;
1993 263871 : return FpM_intersect_i(x, y, p);
1994 : }
1995 : /* Fp-basis of (ZK/pr): applied to the primes found in primedec_aux()
1996 : * true nf */
1997 : static GEN
1998 311943 : Fp_basis(GEN nf, GEN pr)
1999 : {
2000 : long i, j, l;
2001 : GEN x, y;
2002 : /* already in basis form (from Buchman-Lenstra) ? */
2003 311943 : if (typ(pr) == t_MAT) return pr;
2004 : /* ordinary prid (from Kummer) */
2005 74083 : x = pr_hnf(nf, pr);
2006 74083 : l = lg(x);
2007 74083 : y = cgetg(l, t_MAT);
2008 620401 : for (i=j=1; i<l; i++)
2009 546318 : if (gequal1(gcoeff(x,i,i))) gel(y,j++) = gel(x,i);
2010 74083 : setlg(y, j); return y;
2011 : }
2012 : /* Let Ip = prod_{ P | p } P be the p-radical. The list L contains the
2013 : * P (mod Ip) seen as sub-Fp-vector spaces of ZK/Ip.
2014 : * Return the list of (Ip / P) (mod Ip).
2015 : * N.B: All ideal multiplications are computed as intersections of Fp-vector
2016 : * spaces. true nf */
2017 : static GEN
2018 216656 : get_LV(GEN nf, GEN L, GEN p, long N)
2019 : {
2020 216656 : long i, l = lg(L)-1;
2021 : GEN LV, LW, A, B;
2022 :
2023 216656 : LV = cgetg(l+1, t_VEC);
2024 216656 : if (l == 1) { gel(LV,1) = matid(N); return LV; }
2025 111993 : LW = cgetg(l+1, t_VEC);
2026 423936 : for (i=1; i<=l; i++) gel(LW,i) = Fp_basis(nf, gel(L,i));
2027 :
2028 : /* A[i] = L[1]...L[i-1], i = 2..l */
2029 111993 : A = cgetg(l+1, t_VEC); gel(A,1) = NULL;
2030 311943 : for (i=1; i < l; i++) gel(A,i+1) = mul_intersect(gel(A,i), gel(LW,i), p);
2031 : /* B[i] = L[i+1]...L[l], i = 1..(l-1) */
2032 111993 : B = cgetg(l+1, t_VEC); gel(B,l) = NULL;
2033 311943 : for (i=l; i>=2; i--) gel(B,i-1) = mul_intersect(gel(B,i), gel(LW,i), p);
2034 423936 : for (i=1; i<=l; i++) gel(LV,i) = mul_intersect(gel(A,i), gel(B,i), p);
2035 111993 : return LV;
2036 : }
2037 :
2038 : static void
2039 0 : errprime(GEN p) { pari_err_PRIME("idealprimedec",p); }
2040 :
2041 : /* P = Fp-basis (over O_K/p) for pr.
2042 : * V = Z-basis for I_p/pr. ramif != 0 iff some pr|p is ramified.
2043 : * Return a p-uniformizer for pr. Assume pr not inert, i.e. m > 0 */
2044 : static GEN
2045 299589 : uniformizer(GEN nf, norm_S *S, GEN P, GEN V, GEN p, int ramif)
2046 : {
2047 299589 : long i, l, f, m = lg(P)-1, N = nf_get_degree(nf);
2048 : GEN u, Mv, x, q;
2049 :
2050 299589 : f = N - m; /* we want v_p(Norm(x)) = p^f */
2051 299589 : q = powiu(p,f+1);
2052 :
2053 299589 : u = FpM_FpC_invimage(shallowconcat(P, V), col_ei(N,1), p);
2054 299589 : setlg(u, lg(P));
2055 299589 : u = centermod(ZM_ZC_mul(P, u), p);
2056 299589 : if (is_uniformizer(u, q, S)) return u;
2057 159021 : if (signe(gel(u,1)) <= 0) /* make sure u[1] in ]-p,p] */
2058 135014 : gel(u,1) = addii(gel(u,1), p); /* try u + p */
2059 : else
2060 24007 : gel(u,1) = subii(gel(u,1), p); /* try u - p */
2061 159021 : if (!ramif || is_uniformizer(u, q, S)) return u;
2062 :
2063 : /* P/p ramified, u in P^2, not in Q for all other Q|p */
2064 85791 : Mv = zk_multable(nf, Z_ZC_sub(gen_1,u));
2065 85791 : l = lg(P);
2066 115125 : for (i=1; i<l; i++)
2067 : {
2068 115125 : x = centermod(ZC_add(u, ZM_ZC_mul(Mv, gel(P,i))), p);
2069 115125 : if (is_uniformizer(x, q, S)) return x;
2070 : }
2071 0 : errprime(p);
2072 : return NULL; /* LCOV_EXCL_LINE */
2073 : }
2074 :
2075 : /*******************************************************************/
2076 : /* */
2077 : /* BUCHMANN-LENSTRA ALGORITHM */
2078 : /* */
2079 : /*******************************************************************/
2080 : static GEN
2081 4398853 : mk_pr(GEN p, GEN u, long e, long f, GEN t)
2082 4398853 : { return mkvec5(p, u, utoipos(e), utoipos(f), t); }
2083 :
2084 : /* nf a true nf, u in Z[X]/(T); pr = p Z_K + u Z_K of ramification index e */
2085 : GEN
2086 3938728 : idealprimedec_kummer(GEN nf,GEN u,long e,GEN p)
2087 : {
2088 3938728 : GEN t, T = nf_get_pol(nf);
2089 3938728 : long f = degpol(u), N = degpol(T);
2090 :
2091 3938728 : if (f == N)
2092 : { /* inert */
2093 625344 : u = scalarcol_shallow(p,N);
2094 625344 : t = gen_1;
2095 : }
2096 : else
2097 : {
2098 3313384 : t = centermod(poltobasis(nf, FpX_div(T, u, p)), p);
2099 3313384 : u = centermod(poltobasis(nf, u), p);
2100 3313384 : if (e == 1)
2101 : { /* make sure v_pr(u) = 1 (automatic if e>1) */
2102 3030222 : GEN cw, w = Q_primitive_part(nf_to_scalar_or_alg(nf, u), &cw);
2103 3030222 : long v = cw? f - Q_pval(cw, p) * N: f;
2104 3030222 : if (ZpX_resultant_val(T, w, p, v + 1) > v)
2105 : {
2106 124157 : GEN c = gel(u,1);
2107 124157 : gel(u,1) = signe(c) > 0? subii(c, p): addii(c, p);
2108 : }
2109 : }
2110 3313384 : t = zk_multable(nf, t);
2111 : }
2112 3938728 : return mk_pr(p,u,e,f,t);
2113 : }
2114 :
2115 : typedef struct {
2116 : GEN nf, p;
2117 : long I;
2118 : } eltmod_muldata;
2119 :
2120 : static GEN
2121 830756 : sqr_mod(void *data, GEN x)
2122 : {
2123 830756 : eltmod_muldata *D = (eltmod_muldata*)data;
2124 830756 : return FpC_red(nfsqri(D->nf, x), D->p);
2125 : }
2126 : static GEN
2127 352490 : ei_msqr_mod(void *data, GEN x)
2128 : {
2129 352490 : GEN x2 = sqr_mod(data, x);
2130 352490 : eltmod_muldata *D = (eltmod_muldata*)data;
2131 352490 : return FpC_red(zk_ei_mul(D->nf, x2, D->I), D->p);
2132 : }
2133 : /* nf a true nf; compute lift(nf.zk[I]^p mod p) */
2134 : static GEN
2135 748050 : pow_ei_mod_p(GEN nf, long I, GEN p)
2136 : {
2137 748050 : pari_sp av = avma;
2138 : eltmod_muldata D;
2139 748050 : long N = nf_get_degree(nf);
2140 748050 : GEN y = col_ei(N,I);
2141 748050 : if (I == 1) return y;
2142 529063 : D.nf = nf;
2143 529063 : D.p = p;
2144 529063 : D.I = I;
2145 529063 : y = gen_pow_fold(y, p, (void*)&D, &sqr_mod, &ei_msqr_mod);
2146 529063 : return gc_upto(av,y);
2147 : }
2148 :
2149 : /* nf a true nf; return a Z basis of Z_K's p-radical, phi = x--> x^p-x */
2150 : static GEN
2151 216656 : pradical(GEN nf, GEN p, GEN *phi)
2152 : {
2153 216656 : long i, N = nf_get_degree(nf);
2154 : GEN q,m,frob,rad;
2155 :
2156 : /* matrix of Frob: x->x^p over Z_K/p */
2157 216656 : frob = cgetg(N+1,t_MAT);
2158 955074 : for (i=1; i<=N; i++) gel(frob,i) = pow_ei_mod_p(nf,i,p);
2159 :
2160 216656 : m = frob; q = p;
2161 308516 : while (abscmpiu(q,N) < 0) { q = mulii(q,p); m = FpM_mul(m, frob, p); }
2162 216656 : rad = FpM_ker(m, p); /* m = Frob^k, s.t p^k >= N */
2163 955074 : for (i=1; i<=N; i++) gcoeff(frob,i,i) = subiu(gcoeff(frob,i,i), 1);
2164 216656 : *phi = frob; return rad;
2165 : }
2166 :
2167 : /* return powers of a: a^0, ... , a^d, d = dim A */
2168 : static GEN
2169 161193 : get_powers(GEN mul, GEN p)
2170 : {
2171 161193 : long i, d = lgcols(mul);
2172 161193 : GEN z, pow = cgetg(d+2,t_MAT), P = pow+1;
2173 :
2174 161193 : gel(P,0) = scalarcol_shallow(gen_1, d-1);
2175 161193 : z = gel(mul,1);
2176 765756 : for (i=1; i<=d; i++)
2177 : {
2178 604563 : gel(P,i) = z; /* a^i */
2179 604563 : if (i!=d) z = FpM_FpC_mul(mul, z, p);
2180 : }
2181 161193 : return pow;
2182 : }
2183 :
2184 : /* minimal polynomial of a in A (dim A = d).
2185 : * mul = multiplication table by a in A */
2186 : static GEN
2187 118741 : pol_min(GEN mul, GEN p)
2188 : {
2189 118741 : pari_sp av = avma;
2190 118741 : GEN z = FpM_deplin(get_powers(mul, p), p);
2191 118741 : return gc_GEN(av, RgV_to_RgX(z,0));
2192 : }
2193 :
2194 : static GEN
2195 416284 : get_pr(GEN nf, norm_S *S, GEN p, GEN P, GEN V, int ramif, long N, long flim)
2196 : {
2197 : GEN u, t;
2198 : long e, f;
2199 :
2200 416284 : if (typ(P) == t_VEC)
2201 : { /* already done (Kummer) */
2202 74083 : f = pr_get_f(P);
2203 74083 : if (flim > 0 && f > flim) return NULL;
2204 72757 : if (flim == -2) return (GEN)f;
2205 72750 : return P;
2206 : }
2207 342201 : f = N - (lg(P)-1);
2208 342201 : if (flim > 0 && f > flim) return NULL;
2209 340168 : if (flim == -2) return (GEN)f;
2210 : /* P = (p,u) prime. t is an anti-uniformizer: Z_K + t/p Z_K = P^(-1),
2211 : * so that v_P(t) = e(P/p)-1 */
2212 339790 : if (f == N) {
2213 40201 : u = scalarcol_shallow(p,N);
2214 40201 : t = gen_1;
2215 40201 : e = 1;
2216 : } else {
2217 : GEN mt;
2218 299589 : u = uniformizer(nf, S, P, V, p, ramif);
2219 299589 : t = FpM_deplin(zk_multable(nf,u), p);
2220 299589 : mt = zk_multable(nf, t);
2221 299589 : e = ramif? 1 + ZC_nfval(t,mk_pr(p,u,0,0,mt)): 1;
2222 299589 : t = mt;
2223 : }
2224 339790 : return mk_pr(p,u,e,f,t);
2225 : }
2226 :
2227 : /* true nf */
2228 : static GEN
2229 216656 : primedec_end(GEN nf, GEN L, GEN p, long flim)
2230 : {
2231 216656 : long i, j, l = lg(L), N = nf_get_degree(nf);
2232 216656 : GEN LV = get_LV(nf, L,p,N);
2233 216656 : int ramif = dvdii(nf_get_disc(nf), p);
2234 216656 : norm_S S; init_norm(&S, nf, p);
2235 632541 : for (i = j = 1; i < l; i++)
2236 : {
2237 416284 : GEN P = get_pr(nf, &S, p, gel(L,i), gel(LV,i), ramif, N, flim);
2238 416284 : if (!P) continue;
2239 412925 : gel(L,j++) = P;
2240 412925 : if (flim == -1) return P;
2241 : }
2242 216257 : setlg(L, j); return L;
2243 : }
2244 :
2245 : /* prime ideal decomposition of p; if flim>0, restrict to f(P,p) <= flim
2246 : * if flim = -1 return only the first P
2247 : * if flim = -2 return only the f(P/p) in a t_VECSMALL; true nf */
2248 : static GEN
2249 2831772 : primedec_aux(GEN nf, GEN p, long flim)
2250 : {
2251 2831772 : const long TYP = (flim == -2)? t_VECSMALL: t_VEC;
2252 2831772 : GEN E, F, L, Ip, phi, f, g, h, UN, T = nf_get_pol(nf);
2253 : long i, k, c, iL, N;
2254 : int kummer;
2255 :
2256 2831772 : F = FpX_factor(T, p);
2257 2831772 : E = gel(F,2);
2258 2831772 : F = gel(F,1);
2259 :
2260 2831772 : k = lg(F); if (k == 1) errprime(p);
2261 2831772 : if ( !dvdii(nf_get_index(nf),p) ) /* p doesn't divide index */
2262 : {
2263 2613723 : L = cgetg(k, TYP);
2264 6334847 : for (i=1; i<k; i++)
2265 : {
2266 4335159 : GEN t = gel(F,i);
2267 4335159 : long f = degpol(t);
2268 4335159 : if (flim > 0 && f > flim) { setlg(L, i); break; }
2269 3725933 : if (flim == -2)
2270 0 : L[i] = f;
2271 : else
2272 3725933 : gel(L,i) = idealprimedec_kummer(nf, t, E[i],p);
2273 3725933 : if (flim == -1) return gel(L,1);
2274 : }
2275 2608914 : return L;
2276 : }
2277 :
2278 218049 : kummer = 0;
2279 218049 : g = FpXV_prod(F, p);
2280 218049 : h = FpX_div(T,g,p);
2281 218049 : f = FpX_red(ZX_Z_divexact(ZX_sub(ZX_mul(g,h), T), p), p);
2282 :
2283 218049 : N = degpol(T);
2284 218049 : L = cgetg(N+1,TYP);
2285 218049 : iL = 1;
2286 531921 : for (i=1; i<k; i++)
2287 315265 : if (E[i] == 1 || signe(FpX_rem(f,gel(F,i),p)))
2288 74083 : {
2289 75476 : GEN t = gel(F,i);
2290 75476 : kummer = 1;
2291 75476 : gel(L,iL++) = idealprimedec_kummer(nf, t, E[i],p);
2292 75476 : if (flim == -1) return gel(L,1);
2293 : }
2294 : else /* F[i] | (f,g,h), happens at least once by Dedekind criterion */
2295 239789 : E[i] = 0;
2296 :
2297 : /* phi matrix of x -> x^p - x in algebra Z_K/p */
2298 216656 : Ip = pradical(nf,p,&phi);
2299 :
2300 : /* split etale algebra Z_K / (p,Ip) */
2301 216656 : h = cgetg(N+1,t_VEC);
2302 216656 : if (kummer)
2303 : { /* split off Kummer factors */
2304 47391 : GEN mb, b = NULL;
2305 176374 : for (i=1; i<k; i++)
2306 128983 : if (!E[i]) b = b? FpX_mul(b, gel(F,i), p): gel(F,i);
2307 47391 : if (!b) errprime(p);
2308 47391 : b = FpC_red(poltobasis(nf,b), p);
2309 47391 : mb = FpM_red(zk_multable(nf,b), p);
2310 : /* Fp-base of ideal (Ip, b) in ZK/p */
2311 47391 : gel(h,1) = FpM_image(shallowconcat(mb,Ip), p);
2312 : }
2313 : else
2314 169265 : gel(h,1) = Ip;
2315 :
2316 216656 : UN = col_ei(N, 1);
2317 472685 : for (c=1; c; c--)
2318 : { /* Let A:= (Z_K/p) / Ip etale; split A2 := A / Im H ~ Im M2
2319 : H * ? + M2 * Mi2 = Id_N ==> M2 * Mi2 projector A --> A2 */
2320 256029 : GEN M, Mi, M2, Mi2, phi2, mat1, H = gel(h,c); /* maximal rank */
2321 256029 : long dim, r = lg(H)-1;
2322 :
2323 256029 : M = FpM_suppl(shallowconcat(H,UN), p);
2324 256029 : Mi = FpM_inv(M, p);
2325 256029 : M2 = vecslice(M, r+1,N); /* M = (H|M2) invertible */
2326 256029 : Mi2 = rowslice(Mi,r+1,N);
2327 : /* FIXME: FpM_mul(,M2) could be done with vecpermute */
2328 256029 : phi2 = FpM_mul(Mi2, FpM_mul(phi,M2, p), p);
2329 256029 : mat1 = FpM_ker(phi2, p);
2330 256029 : dim = lg(mat1)-1; /* A2 product of 'dim' fields */
2331 256029 : if (dim > 1)
2332 : { /* phi2 v = 0 => a = M2 v in Ker phi, a not in Fp.1 + H */
2333 118741 : GEN R, a, mula, mul2, v = gel(mat1,2);
2334 : long n;
2335 :
2336 118741 : a = FpM_FpC_mul(M2,v, p); /* not a scalar */
2337 118741 : mula = FpM_red(zk_multable(nf,a), p);
2338 118741 : mul2 = FpM_mul(Mi2, FpM_mul(mula,M2, p), p);
2339 118741 : R = FpX_roots(pol_min(mul2,p), p); /* totally split mod p */
2340 118741 : n = lg(R)-1;
2341 363349 : for (i=1; i<=n; i++)
2342 : {
2343 244608 : GEN I = RgM_Rg_sub_shallow(mula, gel(R,i));
2344 244608 : gel(h,c++) = FpM_image(shallowconcat(H, I), p);
2345 : }
2346 118741 : if (n == dim)
2347 306339 : for (i=1; i<=n; i++) gel(L,iL++) = gel(h,--c);
2348 : }
2349 : else /* A2 field ==> H maximal, f = N-r = dim(A2) */
2350 137288 : gel(L,iL++) = H;
2351 : }
2352 216656 : setlg(L, iL);
2353 216656 : return primedec_end(nf, L, p, flim);
2354 : }
2355 :
2356 : GEN
2357 2824800 : idealprimedec_limit_f(GEN nf, GEN p, long f)
2358 : {
2359 2824800 : pari_sp av = avma;
2360 : GEN v;
2361 2824800 : if (typ(p) != t_INT) pari_err_TYPE("idealprimedec",p);
2362 2824800 : if (f < 0) pari_err_DOMAIN("idealprimedec", "f", "<", gen_0, stoi(f));
2363 2824800 : v = primedec_aux(checknf(nf), p, f);
2364 2824800 : v = gen_sort(v, (void*)&cmp_prime_over_p, &cmp_nodata);
2365 2824800 : return gc_upto(av,v);
2366 : }
2367 : /* true nf */
2368 : GEN
2369 6601 : idealprimedec_galois(GEN nf, GEN p)
2370 : {
2371 6601 : pari_sp av = avma;
2372 6601 : GEN v = primedec_aux(nf, p, -1);
2373 6601 : return gc_GEN(av,v);
2374 : }
2375 : /* true nf */
2376 : GEN
2377 371 : idealprimedec_degrees(GEN nf, GEN p)
2378 : {
2379 371 : pari_sp av = avma;
2380 371 : GEN v = primedec_aux(nf, p, -2);
2381 371 : vecsmall_sort(v); return gc_leaf(av, v);
2382 : }
2383 : GEN
2384 509474 : idealprimedec_limit_norm(GEN nf, GEN p, GEN B)
2385 509474 : { return idealprimedec_limit_f(nf, p, logint(B,p)); }
2386 : GEN
2387 1300131 : idealprimedec(GEN nf, GEN p)
2388 1300131 : { return idealprimedec_limit_f(nf, p, 0); }
2389 : static GEN
2390 26614 : nf_pV_to_prVV(GEN nf, GEN x)
2391 89649 : { pari_APPLY_same(idealprimedec(nf, gel(x,i))); }
2392 : GEN
2393 38605 : nf_pV_to_prV(GEN nf, GEN x)
2394 : {
2395 38605 : if (lg(x) == 1) return leafcopy(x);
2396 26614 : return shallowconcat1(nf_pV_to_prVV(nf, x));
2397 : }
2398 :
2399 : /* return [Fp[x]: Fp] */
2400 : static long
2401 4109 : ffdegree(GEN x, GEN frob, GEN p)
2402 : {
2403 4109 : pari_sp av = avma;
2404 4109 : long d, f = lg(frob)-1;
2405 4109 : GEN y = x;
2406 :
2407 13209 : for (d=1; d < f; d++)
2408 : {
2409 10878 : y = FpM_FpC_mul(frob, y, p);
2410 10878 : if (ZV_equal(y, x)) break;
2411 : }
2412 4109 : return gc_long(av,d);
2413 : }
2414 :
2415 : static GEN
2416 93325 : lift_to_zk(GEN v, GEN c, long N)
2417 : {
2418 93325 : GEN w = zerocol(N);
2419 93325 : long i, l = lg(c);
2420 310943 : for (i=1; i<l; i++) gel(w,c[i]) = gel(v,i);
2421 93325 : return w;
2422 : }
2423 :
2424 : /* return t = 1 mod pr, t = 0 mod p / pr^e(pr/p) */
2425 : GEN
2426 1005886 : pr_anti_uniformizer(GEN nf, GEN pr)
2427 : {
2428 1005886 : long N = nf_get_degree(nf), e = pr_get_e(pr);
2429 : GEN p, b, z;
2430 :
2431 1005886 : if (e * pr_get_f(pr) == N) return gen_1;
2432 498379 : p = pr_get_p(pr);
2433 498379 : b = pr_get_tau(pr); /* ZM */
2434 498379 : if (e != 1)
2435 : {
2436 23597 : GEN q = powiu(pr_get_p(pr), e-1);
2437 23597 : b = ZM_Z_divexact(ZM_powu(b,e), q);
2438 : }
2439 : /* b = tau^e / p^(e-1), v_pr(b) = 0, v_Q(b) >= e(Q/p) for other Q | p */
2440 498379 : z = ZM_hnfmodid(FpM_red(b,p), p); /* ideal (p) / pr^e, coprime to pr */
2441 498379 : z = idealaddtoone_raw(nf, pr, z);
2442 498379 : return Z_ZC_sub(gen_1, FpC_center(FpC_red(z,p), p, shifti(p,-1)));
2443 : }
2444 :
2445 : #define mpr_TAU 1
2446 : #define mpr_FFP 2
2447 : #define mpr_NFP 5
2448 : #define SMALLMODPR 4
2449 : #define LARGEMODPR 6
2450 : static GEN
2451 4050650 : modpr_TAU(GEN modpr)
2452 : {
2453 4050650 : GEN tau = gel(modpr,mpr_TAU);
2454 4050650 : return isintzero(tau)? NULL: tau;
2455 : }
2456 :
2457 : /* H = HNF matrix, which is identity but for the first line. Return a
2458 : * projector to Z^n / H ~ Z/qZ, with q = H[1,1] */
2459 : GEN
2460 1084728 : hnf_Znproj(GEN H)
2461 : {
2462 1084728 : long i, l = lg(H);
2463 1084728 : GEN p = cgetg(l, t_VEC), q = gcoeff(H,1,1);
2464 1084728 : gel(p,1) = gen_1;
2465 2488384 : for (i = 2; i < l; i++) gel(p,i) = Fp_neg(gcoeff(H,1,i), q);
2466 1084728 : return p;
2467 : }
2468 :
2469 : /* p not necessarily prime, but coprime to denom(basis) */
2470 : GEN
2471 210 : QXQV_to_FpM(GEN basis, GEN T, GEN p)
2472 : {
2473 210 : long i, l = lg(basis), f = degpol(T);
2474 210 : GEN z = cgetg(l, t_MAT);
2475 4942 : for (i = 1; i < l; i++)
2476 : {
2477 4732 : GEN w = gel(basis,i);
2478 4732 : if (typ(w) == t_INT)
2479 0 : w = scalarcol_shallow(w, f);
2480 : else
2481 : {
2482 : GEN dx;
2483 4732 : w = Q_remove_denom(w, &dx);
2484 4732 : w = FpXQ_red(w, T, p);
2485 4732 : if (dx)
2486 : {
2487 0 : dx = Fp_inv(dx, p);
2488 0 : if (!equali1(dx)) w = FpX_Fp_mul(w, dx, p);
2489 : }
2490 4732 : w = RgX_to_RgC(w, f);
2491 : }
2492 4732 : gel(z,i) = w; /* w_i mod (T,p) */
2493 : }
2494 210 : return z;
2495 : }
2496 :
2497 : /* initialize reduction mod pr; if zk = 1, will only init data required to
2498 : * reduce *integral* element. Realize (O_K/pr) as Fp[X] / (T), for a
2499 : * *monic* T; use variable vT for varn(T) */
2500 : static GEN
2501 1212622 : modprinit(GEN nf, GEN pr, int zk, long vT)
2502 : {
2503 1212622 : pari_sp av = avma;
2504 : GEN res, tau, mul, x, p, T, pow, ffproj, nfproj, prh, c;
2505 : long N, i, k, f;
2506 :
2507 1212622 : nf = checknf(nf); checkprid(pr);
2508 1212608 : if (vT < 0) vT = nf_get_varn(nf);
2509 1212608 : f = pr_get_f(pr);
2510 1212608 : N = nf_get_degree(nf);
2511 1212608 : prh = pr_hnf(nf, pr);
2512 1212608 : tau = zk? gen_0: pr_anti_uniformizer(nf, pr);
2513 1212608 : p = pr_get_p(pr);
2514 :
2515 1212608 : if (f == 1)
2516 : {
2517 1065881 : res = cgetg(SMALLMODPR, t_COL);
2518 1065881 : gel(res,mpr_TAU) = tau;
2519 1065881 : gel(res,mpr_FFP) = hnf_Znproj(prh);
2520 1065881 : gel(res,3) = pr; return gc_GEN(av, res);
2521 : }
2522 :
2523 146727 : c = cgetg(f+1, t_VECSMALL);
2524 146727 : ffproj = cgetg(N+1, t_MAT);
2525 604219 : for (k=i=1; i<=N; i++)
2526 : {
2527 457492 : x = gcoeff(prh, i,i);
2528 457492 : if (!is_pm1(x)) { c[k] = i; gel(ffproj,i) = col_ei(N, i); k++; }
2529 : else
2530 130144 : gel(ffproj,i) = ZC_neg(gel(prh,i));
2531 : }
2532 146727 : ffproj = rowpermute(ffproj, c);
2533 146727 : if (! dvdii(nf_get_index(nf), p))
2534 : {
2535 104275 : GEN basis = nf_get_zkprimpart(nf), D = nf_get_zkden(nf);
2536 104275 : if (N == f)
2537 : { /* pr inert */
2538 45227 : T = nf_get_pol(nf);
2539 45227 : T = FpX_red(T,p);
2540 45227 : ffproj = RgV_to_RgM(basis, lg(basis)-1);
2541 : }
2542 : else
2543 : {
2544 59048 : T = RgV_RgC_mul(basis, pr_get_gen(pr));
2545 59048 : T = FpX_normalize(FpX_red(T,p),p);
2546 59048 : basis = FqV_red(vecpermute(basis,c), T, p);
2547 59048 : basis = RgV_to_RgM(basis, lg(basis)-1);
2548 59048 : ffproj = ZM_mul(basis, ffproj);
2549 : }
2550 104275 : setvarn(T, vT);
2551 104275 : ffproj = FpM_red(ffproj, p);
2552 104275 : if (!equali1(D))
2553 : {
2554 33201 : D = modii(D,p);
2555 33201 : if (!equali1(D)) ffproj = FpM_Fp_mul(ffproj, Fp_inv(D,p), p);
2556 : }
2557 :
2558 104275 : res = cgetg(SMALLMODPR+1, t_COL);
2559 104275 : gel(res,mpr_TAU) = tau;
2560 104275 : gel(res,mpr_FFP) = ffproj;
2561 104275 : gel(res,3) = pr;
2562 104275 : gel(res,4) = T; return gc_GEN(av, res);
2563 : }
2564 :
2565 42452 : if (uisprime(f))
2566 : {
2567 40121 : mul = ei_multable(nf, c[2]);
2568 40121 : mul = vecpermute(mul, c);
2569 : }
2570 : else
2571 : {
2572 : GEN v, u, u2, frob;
2573 : long deg,deg1,deg2;
2574 :
2575 : /* matrix of Frob: x->x^p over Z_K/pr = < w[c1], ..., w[cf] > over Fp */
2576 2331 : frob = cgetg(f+1, t_MAT);
2577 11963 : for (i=1; i<=f; i++)
2578 : {
2579 9632 : x = pow_ei_mod_p(nf,c[i],p);
2580 9632 : gel(frob,i) = FpM_FpC_mul(ffproj, x, p);
2581 : }
2582 2331 : u = col_ei(f,2); k = 2;
2583 2331 : deg1 = ffdegree(u, frob, p);
2584 4102 : while (deg1 < f)
2585 : {
2586 1771 : k++; u2 = col_ei(f, k);
2587 1771 : deg2 = ffdegree(u2, frob, p);
2588 1771 : deg = ulcm(deg1,deg2);
2589 1771 : if (deg == deg1) continue;
2590 1764 : if (deg == deg2) { deg1 = deg2; u = u2; continue; }
2591 7 : u = ZC_add(u, u2);
2592 7 : while (ffdegree(u, frob, p) < deg) u = ZC_add(u, u2);
2593 7 : deg1 = deg;
2594 : }
2595 2331 : v = lift_to_zk(u,c,N);
2596 :
2597 2331 : mul = cgetg(f+1,t_MAT);
2598 2331 : gel(mul,1) = v; /* assume w_1 = 1 */
2599 9632 : for (i=2; i<=f; i++) gel(mul,i) = zk_ei_mul(nf,v,c[i]);
2600 : }
2601 :
2602 : /* Z_K/pr = Fp(v), mul = mul by v */
2603 42452 : mul = FpM_red(mul, p);
2604 42452 : mul = FpM_mul(ffproj, mul, p);
2605 :
2606 42452 : pow = get_powers(mul, p);
2607 42452 : T = RgV_to_RgX(FpM_deplin(pow, p), vT);
2608 42452 : nfproj = cgetg(f+1, t_MAT);
2609 133446 : for (i=1; i<=f; i++) gel(nfproj,i) = lift_to_zk(gel(pow,i), c, N);
2610 :
2611 42452 : setlg(pow, f+1);
2612 42452 : ffproj = FpM_mul(FpM_inv(pow, p), ffproj, p);
2613 :
2614 42452 : res = cgetg(LARGEMODPR, t_COL);
2615 42452 : gel(res,mpr_TAU) = tau;
2616 42452 : gel(res,mpr_FFP) = ffproj;
2617 42452 : gel(res,3) = pr;
2618 42452 : gel(res,4) = T;
2619 42452 : gel(res,mpr_NFP) = nfproj; return gc_GEN(av, res);
2620 : }
2621 :
2622 : GEN
2623 7 : nfmodprinit(GEN nf, GEN pr) { return modprinit(nf, pr, 0, -1); }
2624 : GEN
2625 175369 : zkmodprinit(GEN nf, GEN pr) { return modprinit(nf, pr, 1, -1); }
2626 : GEN
2627 77 : nfmodprinit0(GEN nf, GEN pr, long v) { return modprinit(nf, pr, 0, v); }
2628 :
2629 : /* x may be a modpr */
2630 : static int
2631 4916921 : ok_modpr(GEN x)
2632 4916921 : { return typ(x) == t_COL && lg(x) >= SMALLMODPR && lg(x) <= LARGEMODPR; }
2633 : void
2634 210 : checkmodpr(GEN x)
2635 : {
2636 210 : if (!ok_modpr(x)) pari_err_TYPE("checkmodpr [use nfmodprinit]", x);
2637 210 : checkprid(modpr_get_pr(x));
2638 210 : }
2639 : GEN
2640 137753 : get_modpr(GEN x)
2641 137753 : { return ok_modpr(x)? x: NULL; }
2642 :
2643 : int
2644 23629249 : checkprid_i(GEN x)
2645 : {
2646 22809656 : return (typ(x) == t_VEC && lg(x) == 6
2647 22736611 : && typ(gel(x,2)) == t_COL && typ(gel(x,3)) == t_INT
2648 46438905 : && typ(gel(x,5)) != t_COL); /* tau changed to t_MAT/t_INT in 2.6 */
2649 : }
2650 : void
2651 18388261 : checkprid(GEN x)
2652 18388261 : { if (!checkprid_i(x)) pari_err_TYPE("checkprid",x); }
2653 : GEN
2654 940002 : get_prid(GEN x)
2655 : {
2656 940002 : long lx = lg(x);
2657 940002 : if (lx == 3 && typ(x) == t_VEC) x = gel(x,1);
2658 940002 : if (checkprid_i(x)) return x;
2659 695380 : if (ok_modpr(x)) {
2660 108465 : x = modpr_get_pr(x);
2661 108465 : if (checkprid_i(x)) return x;
2662 : }
2663 586915 : return NULL;
2664 : }
2665 :
2666 : static GEN
2667 4083578 : to_ff_init(GEN nf, GEN *pr, GEN *T, GEN *p, int zk)
2668 : {
2669 4083578 : GEN modpr = ok_modpr(*pr)? *pr: modprinit(nf, *pr, zk, -1);
2670 4083564 : *T = modpr_get_T(modpr);
2671 4083564 : *pr = modpr_get_pr(modpr);
2672 4083564 : *p = pr_get_p(*pr); return modpr;
2673 : }
2674 :
2675 : /* Return an element of O_K which is set to x Mod T */
2676 : GEN
2677 4508 : modpr_genFq(GEN modpr)
2678 : {
2679 4508 : switch(lg(modpr))
2680 : {
2681 917 : case SMALLMODPR: /* Fp */
2682 917 : return gen_1;
2683 1568 : case LARGEMODPR: /* painful case, p \mid index */
2684 1568 : return gmael(modpr,mpr_NFP, 2);
2685 2023 : default: /* trivial case : p \nmid index */
2686 : {
2687 2023 : long v = varn( modpr_get_T(modpr) );
2688 2023 : return pol_x(v);
2689 : }
2690 : }
2691 : }
2692 :
2693 : GEN
2694 4049467 : nf_to_Fq_init(GEN nf, GEN *pr, GEN *T, GEN *p) {
2695 4049467 : GEN modpr = to_ff_init(nf,pr,T,p,0);
2696 4049453 : GEN tau = modpr_TAU(modpr);
2697 4049453 : if (!tau) gel(modpr,mpr_TAU) = pr_anti_uniformizer(nf, *pr);
2698 4049453 : return modpr;
2699 : }
2700 : GEN
2701 34111 : zk_to_Fq_init(GEN nf, GEN *pr, GEN *T, GEN *p) {
2702 34111 : return to_ff_init(nf,pr,T,p,1);
2703 : }
2704 :
2705 : /* assume x in 'basis' form (t_COL) */
2706 : GEN
2707 6521292 : zk_to_Fq(GEN x, GEN modpr)
2708 : {
2709 6521292 : GEN pr = modpr_get_pr(modpr), p = pr_get_p(pr);
2710 6521292 : GEN ffproj = gel(modpr,mpr_FFP);
2711 6521292 : GEN T = modpr_get_T(modpr);
2712 6521292 : return T? FpM_FpC_mul_FpX(ffproj,x, p, varn(T)): FpV_dotproduct(ffproj,x, p);
2713 : }
2714 :
2715 : /* REDUCTION Modulo a prime ideal */
2716 :
2717 : /* nf a true nf, not GC-clean, OK for gc_upto */
2718 : static GEN
2719 14189338 : nf_to_Fq_i(GEN nf, GEN x0, GEN modpr)
2720 : {
2721 14189338 : GEN x = x0, den, pr = modpr_get_pr(modpr), p = pr_get_p(pr);
2722 14189338 : long tx = typ(x);
2723 :
2724 14189338 : if (tx == t_POLMOD) { x = gel(x,2); tx = typ(x); }
2725 14189338 : switch(tx)
2726 : {
2727 7756972 : case t_INT: return modii(x, p);
2728 5574 : case t_FRAC: return Rg_to_Fp(x, p);
2729 204077 : case t_POL:
2730 204077 : switch(lg(x))
2731 : {
2732 21 : case 2: return gen_0;
2733 25760 : case 3: return Rg_to_Fp(gel(x,2), p);
2734 : }
2735 178296 : x = Q_remove_denom(x, &den);
2736 178296 : x = poltobasis(nf, x);
2737 : /* content(x) and den may not be coprime */
2738 178184 : break;
2739 6222715 : case t_COL:
2740 6222715 : x = Q_remove_denom(x, &den);
2741 : /* content(x) and den are coprime */
2742 6222715 : if (lg(x)-1 == nf_get_degree(nf)) break;
2743 56 : default: pari_err_TYPE("Rg_to_ff",x);
2744 : return NULL;/*LCOV_EXCL_LINE*/
2745 : }
2746 6400843 : if (den)
2747 : {
2748 49013 : long v = Z_pvalrem(den, p, &den);
2749 49013 : if (v)
2750 : {
2751 1799 : if (tx == t_POL) v -= ZV_pvalrem(x, p, &x);
2752 : /* now v = valuation(true denominator of x) */
2753 1799 : if (v > 0)
2754 : {
2755 1197 : GEN tau = modpr_TAU(modpr);
2756 1197 : if (!tau) pari_err_TYPE("zk_to_ff", x0);
2757 1197 : x = nfmuli(nf,x, nfpow_u(nf, tau, v));
2758 1197 : v -= ZV_pvalrem(x, p, &x);
2759 : }
2760 1799 : if (v > 0) pari_err_INV("Rg_to_ff", mkintmod(gen_0,p));
2761 1771 : if (v) return gen_0;
2762 1232 : if (is_pm1(den)) den = NULL;
2763 : }
2764 48446 : x = FpC_red(x, p);
2765 : }
2766 6400276 : x = zk_to_Fq(x, modpr);
2767 6400276 : if (den)
2768 : {
2769 47921 : GEN c = Fp_inv(den, p);
2770 47921 : x = typ(x) == t_INT? Fp_mul(x,c,p): FpX_Fp_mul(x,c,p);
2771 : }
2772 6400276 : return x;
2773 : }
2774 :
2775 : GEN
2776 210 : nfreducemodpr(GEN nf, GEN x, GEN modpr)
2777 : {
2778 210 : pari_sp av = avma;
2779 210 : nf = checknf(nf); checkmodpr(modpr);
2780 210 : return gc_upto(av, algtobasis(nf, Fq_to_nf(nf_to_Fq_i(nf,x,modpr),modpr)));
2781 : }
2782 :
2783 : GEN
2784 350 : nfmodpr(GEN nf, GEN x, GEN pr)
2785 : {
2786 350 : pari_sp av = avma;
2787 : GEN T, p, modpr;
2788 350 : nf = checknf(nf);
2789 350 : modpr = nf_to_Fq_init(nf, &pr, &T, &p);
2790 343 : if (typ(x) == t_MAT && lg(x) == 3)
2791 42 : {
2792 49 : GEN y, v = famat_nfvalrem(nf, x, pr, &y);
2793 49 : long s = signe(v);
2794 49 : if (s < 0) pari_err_INV("nfmodpr", mkintmod(gen_0,p));
2795 42 : if (s > 0)
2796 28 : x = gen_0;
2797 : else
2798 14 : x = FqV_factorback(nfV_to_FqV(gel(y,1), nf, modpr), gel(y,2), T, p);
2799 : }
2800 : else
2801 294 : x = nf_to_Fq_i(nf, x, modpr);
2802 224 : if (!T) T = pol_x(nf_get_varn(nf));
2803 224 : x = Fq_to_FF(x, Tp_to_FF(T,p));
2804 224 : return gc_GEN(av, x);
2805 : }
2806 : GEN
2807 77 : nfmodprlift(GEN nf, GEN x, GEN pr)
2808 : {
2809 77 : pari_sp av = avma;
2810 : GEN T, p, modpr;
2811 : long d;
2812 77 : nf = checknf(nf);
2813 77 : switch(typ(x))
2814 : {
2815 7 : case t_INT: return icopy(x);
2816 0 : case t_INTMOD: return icopy(gel(x,2));
2817 42 : case t_FFELT: break;
2818 28 : case t_VEC: case t_COL: case t_MAT:
2819 63 : pari_APPLY_same(nfmodprlift(nf,gel(x,i),pr));
2820 0 : default: pari_err_TYPE("nfmodprlit",x);
2821 : }
2822 42 : x = FF_to_FpXQ(x);
2823 42 : setvarn(x, nf_get_varn(nf));
2824 42 : d = degpol(x);
2825 42 : if (d <= 0) { set_avma(av); return d? gen_0: icopy(gel(x,2)); }
2826 14 : modpr = nf_to_Fq_init(nf, &pr, &T, &p);
2827 14 : return gc_GEN(av, Fq_to_nf(x, modpr));
2828 : }
2829 :
2830 : /* lift A from residue field to nf */
2831 : GEN
2832 3087513 : Fq_to_nf(GEN A, GEN modpr)
2833 : {
2834 : long dA;
2835 3087513 : if (typ(A) == t_INT || lg(modpr) < LARGEMODPR) return A;
2836 45342 : dA = degpol(A);
2837 45342 : if (dA <= 0) return dA ? gen_0: gel(A,2);
2838 41149 : return ZM_ZX_mul(gel(modpr,mpr_NFP), A);
2839 : }
2840 : GEN
2841 0 : FqV_to_nfV(GEN x, GEN modpr)
2842 0 : { pari_APPLY_same(Fq_to_nf(gel(x,i), modpr)) }
2843 : GEN
2844 2934 : FqM_to_nfM(GEN A, GEN modpr)
2845 : {
2846 2934 : long i,j,h,l = lg(A);
2847 2934 : GEN B = cgetg(l, t_MAT);
2848 :
2849 2934 : if (l == 1) return B;
2850 2633 : h = lgcols(A);
2851 11010 : for (j=1; j<l; j++)
2852 : {
2853 8377 : GEN Aj = gel(A,j), Bj = cgetg(h,t_COL); gel(B,j) = Bj;
2854 52732 : for (i=1; i<h; i++) gel(Bj,i) = Fq_to_nf(gel(Aj,i), modpr);
2855 : }
2856 2633 : return B;
2857 : }
2858 : GEN
2859 10226 : FqX_to_nfX(GEN x, GEN modpr)
2860 : {
2861 10226 : if (typ(x) != t_POL) return icopy(x); /* scalar */
2862 43271 : pari_APPLY_pol(Fq_to_nf(gel(x,i), modpr));
2863 : }
2864 :
2865 : /* true nf */
2866 : static GEN
2867 14188834 : gc_nf_to_Fq(GEN nf, GEN A, GEN modpr)
2868 : {
2869 14188834 : pari_sp av = avma;
2870 14188834 : return gc_upto(av, nf_to_Fq_i(nf, A, modpr));
2871 : }
2872 : GEN
2873 12922213 : nf_to_Fq(GEN nf, GEN A, GEN modpr)
2874 12922213 : { return gc_nf_to_Fq(checknf(nf), A, modpr); }
2875 : /* A t_VEC/t_COL */
2876 : GEN
2877 167410 : nfV_to_FqV(GEN x, GEN nf, GEN modpr)
2878 : {
2879 167410 : nf = checknf(nf);
2880 916830 : pari_APPLY_same(gc_nf_to_Fq(nf, gel(x,i), modpr));
2881 : }
2882 : /* A t_MAT */
2883 : GEN
2884 1872 : nfM_to_FqM(GEN A, GEN nf, GEN modpr)
2885 : {
2886 1872 : long i,j,h,l = lg(A);
2887 1872 : GEN B = cgetg(l,t_MAT);
2888 :
2889 1872 : if (l == 1) return B;
2890 1872 : h = lgcols(A); nf = checknf(nf);
2891 42624 : for (j=1; j<l; j++)
2892 : {
2893 40752 : GEN Aj = gel(A,j), Bj = cgetg(h,t_COL); gel(B,j) = Bj;
2894 305810 : for (i=1; i<h; i++) gel(Bj,i) = gc_nf_to_Fq(nf, gel(Aj,i), modpr);
2895 : }
2896 1872 : return B;
2897 : }
2898 : /* A t_POL */
2899 : GEN
2900 9415 : nfX_to_FqX(GEN x, GEN nf, GEN modpr)
2901 : {
2902 9415 : nf = checknf(nf);
2903 51988 : pari_APPLY_pol(gc_nf_to_Fq(nf, gel(x,i), modpr));
2904 : }
2905 :
2906 : /*******************************************************************/
2907 : /* */
2908 : /* RELATIVE ROUND 2 */
2909 : /* */
2910 : /*******************************************************************/
2911 : /* Shallow functions */
2912 : /* FIXME: use a bb_field and export the nfX_* routines */
2913 : static GEN
2914 4144 : nfX_sub(GEN nf, GEN x, GEN y)
2915 : {
2916 4144 : long i, lx = lg(x), ly = lg(y);
2917 : GEN z;
2918 4144 : if (ly <= lx) {
2919 4144 : z = cgetg(lx,t_POL); z[1] = x[1];
2920 24871 : for (i=2; i < ly; i++) gel(z,i) = nfsub(nf,gel(x,i),gel(y,i));
2921 4144 : for ( ; i < lx; i++) gel(z,i) = gel(x,i);
2922 4144 : z = normalizepol_lg(z, lx);
2923 : } else {
2924 0 : z = cgetg(ly,t_POL); z[1] = y[1];
2925 0 : for (i=2; i < lx; i++) gel(z,i) = nfsub(nf,gel(x,i),gel(y,i));
2926 0 : for ( ; i < ly; i++) gel(z,i) = gneg(gel(y,i));
2927 0 : z = normalizepol_lg(z, ly);
2928 : }
2929 4144 : return z;
2930 : }
2931 : /* FIXME: quadratic multiplication */
2932 : static GEN
2933 20667 : nfX_mul(GEN nf, GEN a, GEN b)
2934 : {
2935 20667 : long da = degpol(a), db = degpol(b), dc, lc, k;
2936 : GEN c;
2937 20667 : if (da < 0 || db < 0) return gen_0;
2938 20667 : dc = da + db;
2939 20667 : if (dc == 0) return nfmul(nf, gel(a,2),gel(b,2));
2940 20667 : lc = dc+3;
2941 20667 : c = cgetg(lc, t_POL); c[1] = a[1];
2942 170070 : for (k = 0; k <= dc; k++)
2943 : {
2944 149403 : long i, I = minss(k, da);
2945 149403 : GEN d = NULL;
2946 541730 : for (i = maxss(k-db, 0); i <= I; i++)
2947 : {
2948 392327 : GEN e = nfmul(nf, gel(a, i+2), gel(b, k-i+2));
2949 392327 : d = d? nfadd(nf, d, e): e;
2950 : }
2951 149403 : gel(c, k+2) = d;
2952 : }
2953 20667 : return normalizepol_lg(c, lc);
2954 : }
2955 : /* assume b monic */
2956 : static GEN
2957 16523 : nfX_rem(GEN nf, GEN a, GEN b)
2958 : {
2959 16523 : long da = degpol(a), db = degpol(b);
2960 16523 : if (da < 0) return gen_0;
2961 16523 : a = leafcopy(a);
2962 39759 : while (da >= db)
2963 : {
2964 23236 : long i, k = da;
2965 23236 : GEN A = gel(a, k+2);
2966 190325 : for (i = db-1, k--; i >= 0; i--, k--)
2967 167089 : gel(a,k+2) = nfsub(nf, gel(a,k+2), nfmul(nf, A, gel(b,i+2)));
2968 23236 : a = normalizepol_lg(a, lg(a)-1);
2969 23236 : da = degpol(a);
2970 : }
2971 16523 : return a;
2972 : }
2973 : static GEN
2974 16523 : nfXQ_mul(GEN nf, GEN a, GEN b, GEN T)
2975 : {
2976 16523 : GEN c = nfX_mul(nf, a, b);
2977 16523 : if (typ(c) != t_POL) return c;
2978 16523 : return nfX_rem(nf, c, T);
2979 : }
2980 :
2981 : static void
2982 4222 : fill(long l, GEN H, GEN Hx, GEN I, GEN Ix)
2983 : {
2984 : long i;
2985 4222 : if (typ(Ix) == t_VEC) /* standard */
2986 15335 : for (i=1; i<l; i++) { gel(H,i) = gel(Hx,i); gel(I,i) = gel(Ix,i); }
2987 : else /* constant ideal */
2988 3830 : for (i=1; i<l; i++) { gel(H,i) = gel(Hx,i); gel(I,i) = Ix; }
2989 4222 : }
2990 :
2991 : /* given MODULES x and y by their pseudo-bases, returns a pseudo-basis of the
2992 : * module generated by x and y. */
2993 : static GEN
2994 2111 : rnfjoinmodules_i(GEN nf, GEN Hx, GEN Ix, GEN Hy, GEN Iy)
2995 : {
2996 2111 : long lx = lg(Hx), ly = lg(Hy), l = lx+ly-1;
2997 2111 : GEN H = cgetg(l, t_MAT), I = cgetg(l, t_VEC);
2998 2111 : fill(lx, H , Hx, I , Ix);
2999 2111 : fill(ly, H+lx-1, Hy, I+lx-1, Iy); return nfhnf(nf, mkvec2(H, I));
3000 : }
3001 : static GEN
3002 1329 : rnfjoinmodules(GEN nf, GEN x, GEN y)
3003 : {
3004 1329 : if (!x) return y;
3005 644 : if (!y) return x;
3006 644 : return rnfjoinmodules_i(nf, gel(x,1), gel(x,2), gel(y,1), gel(y,2));
3007 : }
3008 :
3009 : typedef struct {
3010 : GEN multab, T,p;
3011 : long h;
3012 : } rnfeltmod_muldata;
3013 :
3014 : static GEN
3015 16190 : _sqr(void *data, GEN x)
3016 : {
3017 16190 : rnfeltmod_muldata *D = (rnfeltmod_muldata *) data;
3018 10356 : GEN z = x? tablesqr(D->multab,x)
3019 16190 : : tablemul_ei_ej(D->multab,D->h,D->h);
3020 16190 : return FqC_red(z,D->T,D->p);
3021 : }
3022 : static GEN
3023 4158 : _msqr(void *data, GEN x)
3024 : {
3025 4158 : GEN x2 = _sqr(data, x), z;
3026 4158 : rnfeltmod_muldata *D = (rnfeltmod_muldata *) data;
3027 4158 : z = tablemul_ei(D->multab, x2, D->h);
3028 4158 : return FqC_red(z,D->T,D->p);
3029 : }
3030 :
3031 : /* Compute W[h]^n mod (T,p) in the extension, assume n >= 0. T a ZX */
3032 : static GEN
3033 5834 : rnfeltid_powmod(GEN multab, long h, GEN n, GEN T, GEN p)
3034 : {
3035 5834 : pari_sp av = avma;
3036 : GEN y;
3037 : rnfeltmod_muldata D;
3038 :
3039 5834 : if (!signe(n)) return gen_1;
3040 :
3041 5834 : D.multab = multab;
3042 5834 : D.h = h;
3043 5834 : D.T = T;
3044 5834 : D.p = p;
3045 5834 : y = gen_pow_fold(NULL, n, (void*)&D, &_sqr, &_msqr);
3046 5834 : return gc_GEN(av, y);
3047 : }
3048 :
3049 : /* P != 0 has at most degpol(P) roots. Look for an element in Fq which is not
3050 : * a root, cf repres() */
3051 : static GEN
3052 21 : FqX_non_root(GEN P, GEN T, GEN p)
3053 : {
3054 21 : long dP = degpol(P), f, vT;
3055 : long i, j, k, pi, pp;
3056 : GEN v;
3057 :
3058 21 : if (dP == 0) return gen_1;
3059 21 : pp = is_bigint(p) ? dP+1: itos(p);
3060 21 : v = cgetg(dP + 2, t_VEC);
3061 21 : gel(v,1) = gen_0;
3062 21 : if (T)
3063 0 : { f = degpol(T); vT = varn(T); }
3064 : else
3065 21 : { f = 1; vT = 0; }
3066 42 : for (i=pi=1; i<=f; i++,pi*=pp)
3067 : {
3068 21 : GEN gi = i == 1? gen_1: pol_xn(i-1, vT), jgi = gi;
3069 42 : for (j=1; j<pp; j++)
3070 : {
3071 42 : for (k=1; k<=pi; k++)
3072 : {
3073 21 : GEN z = Fq_add(gel(v,k), jgi, T,p);
3074 21 : if (!gequal0(FqX_eval(P, z, T,p))) return z;
3075 21 : gel(v, j*pi+k) = z;
3076 : }
3077 21 : if (j < pp-1) jgi = Fq_add(jgi, gi, T,p); /* j*g[i] */
3078 : }
3079 : }
3080 21 : return NULL;
3081 : }
3082 :
3083 : /* true nf, x t_POL */
3084 : static int
3085 8253 : nfpolisintegral_i(GEN nf, GEN x)
3086 : {
3087 8253 : GEN d, T = nf_get_pol(nf);
3088 8253 : long l = lg(x);
3089 8253 : if (varn(x) != varn(T)) pari_err_VAR("nfisintegral", x,T);
3090 8253 : if (l >= lg(T)) { x = RgX_rem(x, T); l = lg(x); }
3091 8253 : if (l == 2) return 1;
3092 8253 : if (l == 3)
3093 : {
3094 0 : switch(typ(gel(x,2)))
3095 : {
3096 0 : case t_INT: return 1;
3097 0 : case t_FRAC: return 0;
3098 0 : default: pari_err_TYPE("nfisintegral",x);
3099 : }
3100 : }
3101 8253 : x = Q_remove_denom(x, &d);
3102 8253 : if (!RgX_is_ZX(x)) pari_err_TYPE("nfisintegral",x);
3103 8253 : if (!d) return 1;
3104 672 : x = ZM_ZX_mul(nf_get_invzk(nf), x);
3105 672 : return ZV_Z_dvd(x, d);
3106 : }
3107 : static int
3108 8253 : nfpolisintegral(GEN nf, GEN x)
3109 8253 : { pari_sp av = avma; return gc_int(av, nfpolisintegral_i(nf, x)); }
3110 :
3111 : /* true nf */
3112 : static int
3113 32235 : nfisintegral(GEN nf, GEN x)
3114 : {
3115 32235 : switch(typ(x))
3116 : {
3117 23975 : case t_INT: return 1;
3118 7 : case t_FRAC: return 0;
3119 0 : case t_POLMOD:
3120 0 : x = checknfelt_mod(nf,x,"nfisintegral");
3121 0 : switch(typ(x))
3122 : {
3123 0 : case t_INT: return 1;
3124 0 : case t_FRAC: return 0;
3125 0 : case t_POL: return nfpolisintegral(nf,x);
3126 : }
3127 0 : break;
3128 8253 : case t_POL: return nfpolisintegral(nf,x);
3129 0 : case t_COL:
3130 0 : if (lg(x)-1 != nf_get_degree(nf)) break;
3131 0 : return RgV_is_ZV(x);
3132 : }
3133 0 : pari_err_TYPE("nfisintegral",x);
3134 : return 0; /* LCOV_EXCL_LINE */
3135 : }
3136 : /* true nf */
3137 : static int
3138 7329 : nfXisintegral(GEN nf, GEN x)
3139 : {
3140 7329 : long i, l = lg(x);
3141 39550 : for (i = 2; i < l; i++)
3142 32235 : if (!nfisintegral(nf, gel(x,i))) return 0;
3143 7315 : return 1;
3144 : }
3145 :
3146 : /* Relative Dedekind criterion over (true) nf, applied to the order defined by a
3147 : * root of monic irreducible polynomial P, modulo the prime ideal pr. Assume
3148 : * vdisc = v_pr( disc(P) ).
3149 : * Return NULL if nf[X]/P is pr-maximal. Otherwise, return [flag, O, v]:
3150 : * O = enlarged order, given by a pseudo-basis
3151 : * flag = 1 if O is proven pr-maximal (may be 0 and O nevertheless pr-maximal)
3152 : * v = v_pr(disc(O)). */
3153 : static GEN
3154 4172 : rnfdedekind_i(GEN nf, GEN P, GEN pr, long vdisc, long only_maximal)
3155 : {
3156 : GEN Ppr, A, I, p, tau, g, h, k, base, T, gzk, hzk, prinvp, pal, nfT, modpr;
3157 : long m, vt, r, d, i, j, mpr;
3158 :
3159 4172 : if (vdisc < 0) pari_err_TYPE("rnfdedekind [non integral pol]", P);
3160 4165 : if (vdisc == 1) return NULL; /* pr-maximal */
3161 4165 : if (!only_maximal && !gequal1(leading_coeff(P)))
3162 0 : pari_err_IMPL( "the full Dedekind criterion in the nonmonic case");
3163 4165 : if (!nfXisintegral(nf, P))
3164 0 : pari_err_IMPL("non integral polynomial in rnfdedekind");
3165 : /* either monic OR only_maximal = 1 */
3166 4165 : m = degpol(P);
3167 4165 : nfT = nf_get_pol(nf);
3168 4165 : modpr = nf_to_Fq_init(nf,&pr, &T, &p);
3169 4165 : Ppr = nfX_to_FqX(P, nf, modpr);
3170 4165 : mpr = degpol(Ppr);
3171 4165 : if (mpr < m) /* nonmonic => only_maximal = 1 */
3172 : {
3173 21 : if (mpr < 0) return NULL;
3174 21 : if (! RgX_valrem(Ppr, &Ppr))
3175 : { /* nonzero constant coefficient */
3176 0 : Ppr = RgX_shift_shallow(RgX_recip_i(Ppr), m - mpr);
3177 0 : P = RgX_recip_i(P);
3178 : }
3179 : else
3180 : {
3181 21 : GEN z = FqX_non_root(Ppr, T, p);
3182 21 : if (!z) pari_err_IMPL( "Dedekind in the difficult case");
3183 0 : z = Fq_to_nf(z, modpr);
3184 0 : if (typ(z) == t_INT)
3185 0 : P = RgX_Rg_translate(P, z);
3186 : else
3187 0 : P = RgXQX_RgXQ_translate(P, z, T);
3188 0 : P = RgX_recip_i(P);
3189 0 : Ppr = nfX_to_FqX(P, nf, modpr); /* degpol(P) = degpol(Ppr) = m */
3190 : }
3191 : }
3192 4144 : A = gel(FqX_factor(Ppr,T,p),1);
3193 4144 : r = lg(A); /* > 1 */
3194 4144 : g = gel(A,1);
3195 7266 : for (i=2; i<r; i++) g = FqX_mul(g, gel(A,i), T, p);
3196 4144 : h = FqX_div(Ppr,g, T, p);
3197 4144 : gzk = FqX_to_nfX(g, modpr);
3198 4144 : hzk = FqX_to_nfX(h, modpr);
3199 4144 : k = nfX_sub(nf, P, nfX_mul(nf, gzk,hzk));
3200 4144 : tau = pr_get_tau(pr);
3201 4144 : switch(typ(tau))
3202 : {
3203 2086 : case t_INT: k = gdiv(k, p); break;
3204 2058 : case t_MAT: k = RgX_Rg_div(tablemulvec(NULL,tau, k), p); break;
3205 : }
3206 4144 : k = nfX_to_FqX(k, nf, modpr);
3207 4144 : k = FqX_normalize(FqX_gcd(FqX_gcd(g,h, T,p), k, T,p), T,p);
3208 4144 : d = degpol(k); /* <= m */
3209 4144 : if (!d) return NULL; /* pr-maximal */
3210 1952 : if (only_maximal) return gen_0; /* not maximal */
3211 :
3212 1931 : A = cgetg(m+d+1,t_MAT);
3213 1931 : I = cgetg(m+d+1,t_VEC); base = mkvec2(A, I);
3214 : /* base[2] temporarily multiplied by p, for the final nfhnfmod,
3215 : * which requires integral ideals */
3216 1931 : prinvp = pr_inv_p(pr); /* again multiplied by p */
3217 10427 : for (j=1; j<=m; j++)
3218 : {
3219 8496 : gel(A,j) = col_ei(m, j);
3220 8496 : gel(I,j) = p;
3221 : }
3222 1931 : pal = FqX_to_nfX(FqX_div(Ppr,k, T,p), modpr);
3223 4198 : for ( ; j<=m+d; j++)
3224 : {
3225 2267 : gel(A,j) = RgX_to_RgC(pal,m);
3226 2267 : gel(I,j) = prinvp;
3227 2267 : if (j < m+d) pal = RgXQX_rem(RgX_shift_shallow(pal,1),P,nfT);
3228 : }
3229 : /* the modulus is integral */
3230 1931 : base = nfhnfmod(nf,base, idealmulpowprime(nf, powiu(p,m), pr, utoineg(d)));
3231 1931 : gel(base,2) = gdiv(gel(base,2), p); /* cancel the factor p */
3232 1931 : vt = vdisc - 2*d;
3233 1931 : return mkvec3(vt < 2? gen_1: gen_0, base, stoi(vt));
3234 : }
3235 :
3236 : /* [L:K] = n */
3237 : static GEN
3238 1513 : triv_order(long n)
3239 : {
3240 1513 : GEN z = cgetg(3, t_VEC);
3241 1513 : gel(z,1) = matid(n);
3242 1513 : gel(z,2) = const_vec(n, gen_1); return z;
3243 : }
3244 :
3245 : /* if flag is set, return gen_1 (resp. gen_0) if the order K[X]/(P)
3246 : * is pr-maximal (resp. not pr-maximal). */
3247 : GEN
3248 91 : rnfdedekind(GEN nf, GEN P, GEN pr, long flag)
3249 : {
3250 91 : pari_sp av = avma;
3251 : GEN z, dP;
3252 : long v;
3253 :
3254 91 : nf = checknf(nf);
3255 91 : P = RgX_nffix("rnfdedekind", nf_get_pol(nf), P, 1);
3256 91 : dP = nfX_disc(nf, P);
3257 91 : if (gequal0(dP))
3258 7 : pari_err_DOMAIN("rnfdedekind","issquarefree(pol)","=",gen_0,P);
3259 84 : if (!pr)
3260 : {
3261 21 : GEN fa = idealfactor(nf, dP);
3262 21 : GEN Q = gel(fa,1), E = gel(fa,2);
3263 21 : pari_sp av2 = avma;
3264 21 : long i, l = lg(Q);
3265 21 : for (i = 1; i < l; i++, set_avma(av2))
3266 : {
3267 21 : v = itos(gel(E,i));
3268 21 : if (rnfdedekind_i(nf,P,gel(Q,i),v,1)) { set_avma(av); return gen_0; }
3269 0 : set_avma(av2);
3270 : }
3271 0 : set_avma(av); return gen_1;
3272 : }
3273 63 : else if (typ(pr) == t_VEC)
3274 : { /* flag = 1 is implicit */
3275 63 : if (lg(pr) == 1) { set_avma(av); return gen_1; }
3276 63 : if (typ(gel(pr,1)) == t_VEC)
3277 : { /* list of primes */
3278 14 : GEN Q = pr;
3279 14 : pari_sp av2 = avma;
3280 14 : long i, l = lg(Q);
3281 14 : for (i = 1; i < l; i++, set_avma(av2))
3282 : {
3283 14 : v = nfval(nf, dP, gel(Q,i));
3284 14 : if (rnfdedekind_i(nf,P,gel(Q,i),v,1)) { set_avma(av); return gen_0; }
3285 : }
3286 0 : set_avma(av); return gen_1;
3287 : }
3288 : }
3289 : /* single prime */
3290 49 : v = nfval(nf, dP, pr);
3291 49 : z = rnfdedekind_i(nf, P, pr, v, flag);
3292 42 : if (z)
3293 : {
3294 21 : if (flag) { set_avma(av); return gen_0; }
3295 14 : z = gc_GEN(av, z);
3296 : }
3297 : else
3298 : {
3299 21 : set_avma(av); if (flag) return gen_1;
3300 7 : z = cgetg(4, t_VEC);
3301 7 : gel(z,1) = gen_1;
3302 7 : gel(z,2) = triv_order(degpol(P));
3303 7 : gel(z,3) = stoi(v);
3304 : }
3305 21 : return z;
3306 : }
3307 :
3308 : static int
3309 8049 : ideal_is1(GEN x) {
3310 8049 : switch(typ(x))
3311 : {
3312 4119 : case t_INT: return is_pm1(x);
3313 3391 : case t_MAT: return RgM_isidentity(x);
3314 : }
3315 539 : return 0;
3316 : }
3317 :
3318 : /* return a in ideal A such that v_pr(a) = v_pr(A) */
3319 : static GEN
3320 3776 : minval(GEN nf, GEN A, GEN pr)
3321 : {
3322 3776 : GEN ab = idealtwoelt(nf,A), a = gel(ab,1), b = gel(ab,2);
3323 3776 : if (nfval(nf,a,pr) > nfval(nf,b,pr)) a = b;
3324 3776 : return a;
3325 : }
3326 :
3327 : /* nf a true nf. Return NULL if power order is pr-maximal */
3328 : static GEN
3329 4088 : rnfmaxord(GEN nf, GEN pol, GEN pr, long vdisc)
3330 : {
3331 4088 : pari_sp av = avma, av1;
3332 : long i, j, k, n, nn, vpol, cnt, sep;
3333 : GEN q, q1, p, T, modpr, W, I, p1;
3334 : GEN prhinv, mpi, Id;
3335 :
3336 4088 : if (DEBUGLEVEL>1) err_printf(" treating %Ps^%ld\n", pr, vdisc);
3337 4088 : modpr = nf_to_Fq_init(nf,&pr,&T,&p);
3338 4088 : av1 = avma;
3339 4088 : p1 = rnfdedekind_i(nf, pol, modpr, vdisc, 0);
3340 4088 : if (!p1) return gc_NULL(av);
3341 1917 : if (is_pm1(gel(p1,1))) return gc_GEN(av,gel(p1,2));
3342 832 : sep = itos(gel(p1,3));
3343 832 : W = gmael(p1,2,1);
3344 832 : I = gmael(p1,2,2);
3345 832 : (void)gc_all(av1, 2, &W, &I);
3346 :
3347 832 : mpi = zk_multable(nf, pr_get_gen(pr));
3348 832 : n = degpol(pol); nn = n*n;
3349 832 : vpol = varn(pol);
3350 832 : q1 = q = pr_norm(pr);
3351 1021 : while (abscmpiu(q1,n) < 0) q1 = mulii(q1,q);
3352 832 : Id = matid(n);
3353 832 : prhinv = pr_inv(pr);
3354 832 : av1 = avma;
3355 832 : for(cnt=1;; cnt++)
3356 1040 : {
3357 1872 : GEN I0 = leafcopy(I), W0 = leafcopy(W);
3358 : GEN Wa, Winv, Ip, A, MW, MWmod, F, pseudo, C, G;
3359 1872 : GEN Tauinv = cgetg(n+1, t_VEC), Tau = cgetg(n+1, t_VEC);
3360 :
3361 1872 : if (DEBUGLEVEL>1) err_printf(" pass no %ld\n",cnt);
3362 9578 : for (j=1; j<=n; j++)
3363 : {
3364 : GEN tau, tauinv;
3365 7706 : if (ideal_is1(gel(I,j)))
3366 : {
3367 3930 : gel(I,j) = gel(Tau,j) = gel(Tauinv,j) = gen_1;
3368 3930 : continue;
3369 : }
3370 3776 : gel(Tau,j) = tau = minval(nf, gel(I,j), pr);
3371 3776 : gel(Tauinv,j) = tauinv = nfinv(nf, tau);
3372 3776 : gel(W,j) = nfC_nf_mul(nf, gel(W,j), tau);
3373 3776 : gel(I,j) = idealmul(nf, tauinv, gel(I,j)); /* v_pr(I[j]) = 0 */
3374 : }
3375 : /* W = (Z_K/pr)-basis of O/pr. O = (W0,I0) ~ (W, I) */
3376 :
3377 : /* compute MW: W_i*W_j = sum MW_k,(i,j) W_k */
3378 1872 : Wa = RgM_to_RgXV(W,vpol);
3379 1872 : Winv = nfM_inv(nf, W);
3380 1872 : MW = cgetg(nn+1, t_MAT);
3381 : /* W_1 = 1 */
3382 9578 : for (j=1; j<=n; j++) gel(MW, j) = gel(MW, (j-1)*n+1) = gel(Id,j);
3383 7706 : for (i=2; i<=n; i++)
3384 22357 : for (j=i; j<=n; j++)
3385 : {
3386 16523 : GEN z = nfXQ_mul(nf, gel(Wa,i), gel(Wa,j), pol);
3387 16523 : if (typ(z) != t_POL)
3388 0 : z = nfC_nf_mul(nf, gel(Winv,1), z);
3389 : else
3390 : {
3391 16523 : z = RgX_to_RgC(z, lg(Winv)-1);
3392 16523 : z = nfM_nfC_mul(nf, Winv, z);
3393 : }
3394 16523 : gel(MW, (i-1)*n+j) = gel(MW, (j-1)*n+i) = z;
3395 : }
3396 :
3397 : /* compute Ip = pr-radical [ could use Ker(trace) if q large ] */
3398 1872 : MWmod = nfM_to_FqM(MW,nf,modpr);
3399 1872 : F = cgetg(n+1, t_MAT); gel(F,1) = gel(Id,1);
3400 7706 : for (j=2; j<=n; j++) gel(F,j) = rnfeltid_powmod(MWmod, j, q1, T,p);
3401 1872 : Ip = FqM_ker(F,T,p);
3402 1872 : if (lg(Ip) == 1) { W = W0; I = I0; break; }
3403 :
3404 : /* Fill C: W_k A_j = sum_i C_(i,j),k A_i */
3405 1467 : A = FqM_to_nfM(FqM_suppl(Ip,T,p), modpr);
3406 3795 : for (j = lg(Ip); j<=n; j++) gel(A,j) = nfC_multable_mul(gel(A,j), mpi);
3407 1467 : MW = nfM_mul(nf, nfM_inv(nf,A), MW);
3408 1467 : C = cgetg(n+1, t_MAT);
3409 7481 : for (k=1; k<=n; k++)
3410 : {
3411 6014 : GEN mek = vecslice(MW, (k-1)*n+1, k*n), Ck;
3412 6014 : gel(C,k) = Ck = cgetg(nn+1, t_COL);
3413 38020 : for (j = 1; j <= n; j++)
3414 : {
3415 32006 : GEN z = nfM_nfC_mul(nf, mek, gel(A,j));
3416 241576 : for (i = 1; i <= n; i++)
3417 209570 : gel(Ck, (j-1)*n+i) = gc_nf_to_Fq(nf,gel(z,i),modpr);
3418 : }
3419 : }
3420 1467 : G = FqM_to_nfM(FqM_ker(C,T,p), modpr);
3421 :
3422 1467 : pseudo = rnfjoinmodules_i(nf, G,prhinv, Id,I);
3423 : /* express W in terms of the power basis */
3424 1467 : W = nfM_mul(nf, W, gel(pseudo,1));
3425 1467 : I = gel(pseudo,2);
3426 : /* restore the HNF property W[i,i] = 1. NB: W upper triangular, with
3427 : * W[i,i] = Tau[i] */
3428 7481 : for (j=1; j<=n; j++)
3429 6014 : if (gel(Tau,j) != gen_1)
3430 : {
3431 2699 : gel(W,j) = nfC_nf_mul(nf, gel(W,j), gel(Tauinv,j));
3432 2699 : gel(I,j) = idealmul(nf, gel(Tau,j), gel(I,j));
3433 : }
3434 1467 : if (DEBUGLEVEL>3) err_printf(" new order:\n%Ps\n%Ps\n", W, I);
3435 1467 : if (sep <= 3 || gequal(I,I0)) break;
3436 :
3437 1040 : if (gc_needed(av1,2))
3438 : {
3439 0 : if(DEBUGMEM>1) pari_warn(warnmem,"rnfmaxord");
3440 0 : (void)gc_all(av1,2, &W,&I);
3441 : }
3442 : }
3443 832 : return gc_GEN(av, mkvec2(W, I));
3444 : }
3445 :
3446 : GEN
3447 969544 : Rg_nffix(const char *f, GEN T, GEN c, int lift)
3448 : {
3449 969544 : switch(typ(c))
3450 : {
3451 537520 : case t_INT: case t_FRAC: return c;
3452 71775 : case t_POL:
3453 71775 : if (lg(c) >= lg(T)) c = RgX_rem(c,T);
3454 71775 : break;
3455 360242 : case t_POLMOD:
3456 360242 : if (!RgX_equal_var(gel(c,1), T)) pari_err_MODULUS(f, gel(c,1),T);
3457 359661 : c = gel(c,2);
3458 359661 : switch(typ(c))
3459 : {
3460 313830 : case t_POL: break;
3461 45831 : case t_INT: case t_FRAC: return c;
3462 0 : default: pari_err_TYPE(f, c);
3463 : }
3464 313830 : break;
3465 7 : default: pari_err_TYPE(f,c);
3466 : }
3467 : /* typ(c) = t_POL */
3468 385605 : if (varn(c) != varn(T)) pari_err_VAR(f, c,T);
3469 385591 : switch(lg(c))
3470 : {
3471 13346 : case 2: return gen_0;
3472 30099 : case 3:
3473 30099 : c = gel(c,2); if (is_rational_t(typ(c))) return c;
3474 0 : pari_err_TYPE(f,c);
3475 : }
3476 342146 : RgX_check_QX(c, f);
3477 342125 : return lift? c: mkpolmod(c, T);
3478 : }
3479 : /* check whether x is a polynomials with coeffs in number field Q[y]/(T)
3480 : * and returned a normalized copy. If 'lift' is set return lifted coefs
3481 : * (t_POL/t_FRAC/t_INT) else t_POLMOD/t_FRAC/t_INT */
3482 : GEN
3483 329976 : RgX_nffix(const char *f, GEN T, GEN x, int lift)
3484 : {
3485 329976 : long vT = varn(T);
3486 329976 : if (typ(x) != t_POL) pari_err_TYPE(stack_strcat(f," [t_POL expected]"), x);
3487 329976 : if (varncmp(varn(x), vT) >= 0) pari_err_PRIORITY(f, x, ">=", vT);
3488 1233441 : pari_APPLY_pol_normalized(Rg_nffix(f, T, gel(x,i), lift));
3489 : }
3490 : GEN
3491 49 : RgV_nffix(const char *f, GEN T, GEN x, int lift)
3492 119 : { pari_APPLY_same(Rg_nffix(f, T, gel(x,i), lift)); }
3493 :
3494 : static GEN
3495 3220 : get_d(GEN nf, GEN d)
3496 : {
3497 3220 : GEN b = idealredmodpower(nf, d, 2, 100000);
3498 3220 : return nfmul(nf, d, nfsqr(nf,b));
3499 : }
3500 :
3501 : /* true nf */
3502 : static GEN
3503 4277 : pr_factorback(GEN nf, GEN fa)
3504 : {
3505 4277 : GEN P = gel(fa,1), E = gel(fa,2), z = gen_1;
3506 4277 : long i, l = lg(P);
3507 8357 : for (i = 1; i < l; i++) z = idealmulpowprime(nf, z, gel(P,i), gel(E,i));
3508 4277 : return z;
3509 : }
3510 : /* true nf */
3511 : static GEN
3512 4277 : pr_factorback_scal(GEN nf, GEN fa)
3513 : {
3514 4277 : GEN D = pr_factorback(nf,fa);
3515 4277 : if (typ(D) == t_MAT && RgM_isscalar(D,NULL)) D = gcoeff(D,1,1);
3516 4277 : return D;
3517 : }
3518 :
3519 : /* nf = base field K
3520 : * pol= monic polynomial in Z_K[X] defining a relative extension L = K[X]/(pol).
3521 : * Returns a pseudo-basis [A,I] of Z_L, set *pD to [D,d] and *pf to the
3522 : * index-ideal; rnf is used when lim != 0 and may be NULL */
3523 : GEN
3524 3171 : rnfallbase(GEN nf, GEN pol, GEN lim, GEN rnf, GEN *pD, GEN *pf, GEN *pDKP)
3525 : {
3526 : long i, j, jf, l;
3527 : GEN fa, E, P, Ef, Pf, z, disc;
3528 :
3529 3171 : nf = checknf(nf); pol = liftpol_shallow(pol);
3530 3171 : if (!gequal1(leading_coeff(pol)))
3531 7 : pari_err_IMPL("nonmonic relative polynomials in rnfallbase");
3532 3164 : if (!nfXisintegral(nf, pol))
3533 14 : pari_err_IMPL("non integral polynomial in rnfallbase");
3534 3150 : disc = nf_to_scalar_or_basis(nf, nfX_disc(nf, pol));
3535 3150 : if (gequal0(disc))
3536 7 : pari_err_DOMAIN("rnfpseudobasis","issquarefree(pol)","=",gen_0, pol);
3537 3143 : if (lim)
3538 : {
3539 : GEN rnfeq, zknf, dzknf, U, vU, dA, A, MB, dB, BdB, vj, B, Tabs;
3540 1015 : GEN D = idealhnf_shallow(nf, disc), extendP = NULL;
3541 1015 : long rU, m = nf_get_degree(nf), n = degpol(pol), N = n*m;
3542 : nfmaxord_t S;
3543 :
3544 1015 : if (typ(lim) == t_INT)
3545 133 : P = ZV_union_shallow(nf_get_ramified_primes(nf),
3546 133 : gel(Z_factor_limit(gcoeff(D,1,1), itou(lim)), 1));
3547 : else
3548 : {
3549 882 : P = cgetg_copy(lim, &l);
3550 3101 : for (i = 1; i < l; i++)
3551 : {
3552 2219 : GEN p = gel(lim,i);
3553 2219 : if (typ(p) != t_INT) p = pr_get_p(p);
3554 2219 : gel(P,i) = p;
3555 : }
3556 882 : P = ZV_sort_uniq_shallow(P);
3557 : }
3558 1015 : if (rnf)
3559 : {
3560 966 : rnfeq = rnf_get_map(rnf);
3561 966 : zknf = rnf_get_nfzk(rnf);
3562 : }
3563 : else
3564 : {
3565 49 : rnfeq = nf_rnfeq(nf, pol);
3566 49 : zknf = nf_nfzk(nf, rnfeq);
3567 : }
3568 1015 : dzknf = gel(zknf,1);
3569 1015 : if (gequal1(dzknf)) dzknf = NULL;
3570 882 : RESTART:
3571 1036 : if (extendP)
3572 : {
3573 21 : GEN oldP = P;
3574 21 : if (typ(extendP)==t_POL)
3575 : {
3576 0 : long l = lg(extendP);
3577 0 : for (i = 2; i < l; i++)
3578 : {
3579 0 : GEN q = gel(extendP,i);
3580 0 : if (typ(q) == t_FRAC) P = ZV_cba_extend(P, gel(q,2));
3581 : }
3582 : } else /*t_FRAC*/
3583 21 : P = ZV_cba_extend(P, gel(extendP,2));
3584 21 : if (ZV_equal(P, oldP))
3585 0 : pari_err(e_MISC, "rnfpseudobasis fails, try increasing B");
3586 21 : extendP = NULL;
3587 : }
3588 1036 : Tabs = gel(rnfeq,1);
3589 1036 : nfmaxord(&S, mkvec2(Tabs,P), 0);
3590 1036 : B = RgXV_unscale(S.basis, S.unscale);
3591 1036 : BdB = Q_remove_denom(B, &dB);
3592 1036 : MB = RgXV_to_RgM(BdB, N); /* HNF */
3593 :
3594 1036 : vU = cgetg(N+1, t_VEC);
3595 1036 : vj = cgetg(N+1, t_VECSMALL);
3596 1036 : gel(vU,1) = U = cgetg(m+1, t_MAT);
3597 1036 : gel(U,1) = col_ei(N, 1);
3598 1036 : A = dB? (dzknf? gdiv(dB,dzknf): dB): NULL;
3599 1036 : if (A)
3600 : {
3601 987 : if (typ(A) != t_INT) { extendP = A; goto RESTART; }
3602 966 : if (equali1(A)) A = NULL;
3603 : }
3604 2065 : for (j = 2; j <= m; j++)
3605 : {
3606 1050 : GEN t = gel(zknf,j);
3607 1050 : if (!RgX_is_ZX(t)) { extendP = t; goto RESTART; }
3608 1050 : if (A) t = ZX_Z_mul(t, A);
3609 1050 : gel(U,j) = hnf_solve(MB, RgX_to_RgC(t, N));
3610 : }
3611 6321 : for (i = 2; i <= N; i++)
3612 : {
3613 5306 : GEN b = gel(BdB,i);
3614 5306 : gel(vU,i) = U = cgetg(m+1, t_MAT);
3615 5306 : gel(U,1) = hnf_solve(MB, RgX_to_RgC(b, N));
3616 11508 : for (j = 2; j <= m; j++)
3617 : {
3618 6202 : GEN t = ZX_rem(ZX_mul(b, gel(zknf,j)), Tabs);
3619 6202 : if (dzknf)
3620 : {
3621 5586 : t = RgX_Rg_div(t, dzknf);
3622 5586 : if (!RgX_is_ZX(t)) { extendP = t; goto RESTART; }
3623 : }
3624 6202 : gel(U,j) = hnf_solve(MB, RgX_to_RgC(t, N));
3625 : }
3626 : }
3627 1015 : vj[1] = 1; U = gel(vU,1); rU = m;
3628 2156 : for (i = j = 2; i <= N; i++)
3629 : {
3630 2149 : GEN V = shallowconcat(U, gel(vU,i));
3631 2149 : if (ZM_rank(V) != rU)
3632 : {
3633 2149 : U = V; rU += m; vj[j++] = i;
3634 2149 : if (rU == N) break;
3635 : }
3636 : }
3637 1015 : if (dB) for(;;)
3638 1316 : {
3639 2282 : GEN c = gen_1, H = ZM_hnfmodid(U, dB);
3640 2282 : long ic = 0;
3641 19957 : for (i = 1; i <= N; i++)
3642 17675 : if (cmpii(gcoeff(H,i,i), c) > 0) { c = gcoeff(H,i,i); ic = i; }
3643 2282 : if (!ic) break;
3644 1316 : vj[j++] = ic;
3645 1316 : U = shallowconcat(H, gel(vU, ic));
3646 : }
3647 1015 : setlg(vj, j);
3648 1015 : B = vecpermute(B, vj);
3649 :
3650 1015 : l = lg(B);
3651 1015 : A = cgetg(l,t_MAT);
3652 5495 : for (j = 1; j < l; j++)
3653 : {
3654 4480 : GEN t = eltabstorel_lift(rnfeq, gel(B,j));
3655 4480 : gel(A,j) = Rg_to_RgC(t, n);
3656 : }
3657 1015 : A = RgM_to_nfM(nf, A);
3658 1015 : A = Q_remove_denom(A, &dA);
3659 1015 : if (!dA)
3660 : { /* order is maximal */
3661 63 : z = triv_order(n);
3662 63 : if (pf) *pf = gen_1;
3663 : }
3664 : else
3665 : {
3666 : GEN fi;
3667 : /* the first n columns of A are probably in HNF already */
3668 952 : A = shallowconcat(vecslice(A,n+1,lg(A)-1), vecslice(A,1,n));
3669 952 : A = mkvec2(A, const_vec(l-1,gen_1));
3670 952 : if (DEBUGLEVEL > 2) err_printf("rnfallbase: nfhnf in dim %ld\n", l-1);
3671 952 : z = nfhnfmod(nf, A, nfdetint(nf,A));
3672 952 : gel(z,2) = gdiv(gel(z,2), dA);
3673 952 : fi = idealprod(nf,gel(z,2));
3674 952 : D = idealmul(nf, D, idealsqr(nf, fi));
3675 952 : if (pf) *pf = idealinv(nf, fi);
3676 : }
3677 1015 : if (RgM_isscalar(D,NULL)) D = gcoeff(D,1,1);
3678 1015 : if (pDKP) *pDKP = S.dKP;
3679 1015 : *pD = mkvec2(D, get_d(nf, disc)); return z;
3680 : }
3681 2128 : fa = idealfactor(nf, disc);
3682 2128 : P = gel(fa,1); l = lg(P); z = NULL;
3683 2128 : E = gel(fa,2);
3684 2128 : Pf = cgetg(l, t_COL);
3685 2128 : Ef = cgetg(l, t_COL);
3686 5998 : for (i = j = jf = 1; i < l; i++)
3687 : {
3688 3870 : GEN pr = gel(P,i);
3689 3870 : long e = itos(gel(E,i));
3690 3870 : if (e > 1)
3691 : {
3692 2877 : GEN vD = rnfmaxord(nf, pol, pr, e);
3693 2877 : if (vD)
3694 : {
3695 1329 : long ef = idealprodval(nf, gel(vD,2), pr);
3696 1329 : z = rnfjoinmodules(nf, z, vD);
3697 1329 : if (ef) { gel(Pf, jf) = pr; gel(Ef, jf++) = stoi(-ef); }
3698 1329 : e += 2 * ef;
3699 : }
3700 : }
3701 3870 : if (e) { gel(P, j) = pr; gel(E, j++) = stoi(e); }
3702 : }
3703 2128 : setlg(P,j);
3704 2128 : setlg(E,j);
3705 2128 : if (pDKP) *pDKP = prV_primes(P);
3706 2128 : if (pf)
3707 : {
3708 2072 : setlg(Pf, jf);
3709 2072 : setlg(Ef, jf); *pf = pr_factorback_scal(nf, mkmat2(Pf,Ef));
3710 : }
3711 2128 : *pD = mkvec2(pr_factorback_scal(nf,fa), get_d(nf, disc));
3712 2128 : return z? z: triv_order(degpol(pol));
3713 : }
3714 :
3715 : static GEN
3716 1666 : RgX_to_algX(GEN nf, GEN x)
3717 10143 : { pari_APPLY_pol_normalized(nf_to_scalar_or_alg(nf, gel(x,i))); }
3718 :
3719 : GEN
3720 1680 : nfX_to_monic(GEN nf, GEN T, GEN *pL)
3721 : {
3722 : GEN lT, g, a;
3723 1680 : long i, l = lg(T);
3724 1680 : if (l == 2) return pol_0(varn(T));
3725 1680 : if (l == 3) return pol_1(varn(T));
3726 1680 : nf = checknf(nf);
3727 1680 : T = Q_primpart(RgX_to_nfX(nf, T));
3728 1680 : lT = leading_coeff(T); if (pL) *pL = lT;
3729 1680 : if (isint1(T)) return T;
3730 1680 : g = cgetg_copy(T, &l); g[1] = T[1]; a = lT;
3731 1680 : gel(g, l-1) = gen_1;
3732 1680 : gel(g, l-2) = gel(T,l-2);
3733 1680 : if (l == 4) { gel(g,l-2) = nf_to_scalar_or_alg(nf, gel(g,l-2)); return g; }
3734 1666 : if (typ(lT) == t_INT)
3735 : {
3736 1652 : gel(g, l-3) = gmul(a, gel(T,l-3));
3737 5124 : for (i = l-4; i > 1; i--) { a = mulii(a,lT); gel(g,i) = gmul(a, gel(T,i)); }
3738 : }
3739 : else
3740 : {
3741 14 : gel(g, l-3) = nfmul(nf, a, gel(T,l-3));
3742 35 : for (i = l-3; i > 1; i--)
3743 : {
3744 21 : a = nfmul(nf,a,lT);
3745 21 : gel(g,i) = nfmul(nf, a, gel(T,i));
3746 : }
3747 : }
3748 1666 : return RgX_to_algX(nf, g);
3749 : }
3750 :
3751 : GEN
3752 868 : rnfdisc_factored(GEN nf, GEN pol, GEN *pd)
3753 : {
3754 : long i, j, l;
3755 : GEN fa, E, P, disc, lim;
3756 :
3757 868 : pol = rnfdisc_get_T(nf, pol, &lim);
3758 868 : disc = nf_to_scalar_or_basis(nf, nfX_disc(nf, pol));
3759 868 : if (gequal0(disc))
3760 0 : pari_err_DOMAIN("rnfdisc","issquarefree(pol)","=",gen_0, pol);
3761 868 : pol = nfX_to_monic(nf, pol, NULL);
3762 868 : fa = idealfactor_partial(nf, disc, lim);
3763 868 : P = gel(fa,1); l = lg(P);
3764 868 : E = gel(fa,2);
3765 2352 : for (i = j = 1; i < l; i++)
3766 : {
3767 1484 : long e = itos(gel(E,i));
3768 1484 : GEN pr = gel(P,i);
3769 1484 : if (e > 1)
3770 : {
3771 1211 : GEN vD = rnfmaxord(nf, pol, pr, e);
3772 1211 : if (vD) e += 2*idealprodval(nf, gel(vD,2), pr);
3773 : }
3774 1484 : if (e) { gel(P, j) = pr; gel(E, j++) = stoi(e); }
3775 : }
3776 868 : if (pd) *pd = get_d(nf, disc);
3777 868 : setlg(P, j);
3778 868 : setlg(E, j); return fa;
3779 : }
3780 : GEN
3781 77 : rnfdiscf(GEN nf, GEN pol)
3782 : {
3783 77 : pari_sp av = avma;
3784 : GEN d, fa;
3785 77 : nf = checknf(nf); fa = rnfdisc_factored(nf, pol, &d);
3786 77 : return gc_GEN(av, mkvec2(pr_factorback_scal(nf,fa), d));
3787 : }
3788 :
3789 : GEN
3790 35 : gen_if_principal(GEN bnf, GEN x)
3791 : {
3792 35 : pari_sp av = avma;
3793 35 : GEN z = bnfisprincipal0(bnf,x, nf_GEN_IF_PRINCIPAL | nf_FORCE);
3794 35 : return isintzero(z)? gc_NULL(av): z;
3795 : }
3796 :
3797 : /* given bnf and a HNF pseudo-basis of a proj. module, simplify the HNF as
3798 : * much as possible. The resulting matrix will be upper triangular but the
3799 : * diagonal coefficients will not be equal to 1. The ideals are integral and
3800 : * primitive. */
3801 : GEN
3802 0 : rnfsimplifybasis(GEN bnf, GEN M)
3803 : {
3804 0 : pari_sp av = avma;
3805 : long i, l;
3806 : GEN y, Az, Iz, nf, A, I;
3807 :
3808 0 : bnf = checkbnf(bnf); nf = bnf_get_nf(bnf);
3809 0 : if (!check_ZKmodule_i(M)) pari_err_TYPE("rnfsimplifybasis",M);
3810 0 : A = gel(M,1);
3811 0 : I = gel(M,2); l = lg(I);
3812 0 : Az = cgetg(l, t_MAT);
3813 0 : Iz = cgetg(l, t_VEC); y = mkvec2(Az, Iz);
3814 0 : for (i = 1; i < l; i++)
3815 : {
3816 : GEN c, d;
3817 0 : if (ideal_is1(gel(I,i)))
3818 : {
3819 0 : gel(Iz,i) = gen_1;
3820 0 : gel(Az,i) = gel(A,i); continue;
3821 : }
3822 :
3823 0 : gel(Iz,i) = Q_primitive_part(gel(I,i), &c);
3824 0 : gel(Az,i) = c? RgC_Rg_mul(gel(A,i),c): gel(A,i);
3825 0 : if (c && ideal_is1(gel(Iz,i))) continue;
3826 :
3827 0 : d = gen_if_principal(bnf, gel(Iz,i));
3828 0 : if (d)
3829 : {
3830 0 : gel(Iz,i) = gen_1;
3831 0 : gel(Az,i) = nfC_nf_mul(nf, gel(Az,i), d);
3832 : }
3833 : }
3834 0 : return gc_GEN(av, y);
3835 : }
3836 :
3837 : static GEN
3838 63 : get_module(GEN nf, GEN O, const char *s)
3839 : {
3840 63 : if (typ(O) == t_POL) return rnfpseudobasis(nf, O);
3841 56 : if (!check_ZKmodule_i(O)) pari_err_TYPE(s, O);
3842 56 : return shallowcopy(O);
3843 : }
3844 :
3845 : GEN
3846 14 : rnfdet(GEN nf, GEN M)
3847 : {
3848 14 : pari_sp av = avma;
3849 : GEN D;
3850 14 : nf = checknf(nf);
3851 14 : M = get_module(nf, M, "rnfdet");
3852 14 : D = idealmul(nf, nfM_det(nf, gel(M,1)), idealprod(nf, gel(M,2)));
3853 14 : return gc_upto(av, D);
3854 : }
3855 :
3856 : /* Given two fractional ideals a and b, gives x in a, y in b, z in b^-1,
3857 : t in a^-1 such that xt-yz=1. In the present version, z is in Z. */
3858 : static void
3859 63 : nfidealdet1(GEN nf, GEN a, GEN b, GEN *px, GEN *py, GEN *pz, GEN *pt)
3860 : {
3861 : GEN x, uv, y, da, db;
3862 :
3863 63 : a = idealinv(nf,a);
3864 63 : a = Q_remove_denom(a, &da);
3865 63 : b = Q_remove_denom(b, &db);
3866 63 : x = idealcoprime(nf,a,b);
3867 63 : uv = idealaddtoone(nf, idealmul(nf,x,a), b);
3868 63 : y = gel(uv,2);
3869 63 : if (da) x = gmul(x,da);
3870 63 : if (db) y = gdiv(y,db);
3871 63 : *px = x;
3872 63 : *py = y;
3873 63 : *pz = db ? negi(db): gen_m1;
3874 63 : *pt = nfdiv(nf, gel(uv,1), x);
3875 63 : }
3876 :
3877 : /* given a pseudo-basis of a proj. module in HNF [A,I] (or [A,I,D,d]), gives
3878 : * an n x n matrix (not HNF) of a pseudo-basis and an ideal vector
3879 : * [1,...,1,I] such that M ~ Z_K^(n-1) x I. Uses the approximation theorem.*/
3880 : GEN
3881 28 : rnfsteinitz(GEN nf, GEN M)
3882 : {
3883 28 : pari_sp av = avma;
3884 : long i, n;
3885 : GEN A, I;
3886 :
3887 28 : nf = checknf(nf);
3888 28 : M = get_module(nf, M, "rnfsteinitz");
3889 28 : A = RgM_to_nfM(nf, gel(M,1));
3890 28 : I = leafcopy(gel(M,2)); n = lg(A)-1;
3891 189 : for (i = 1; i < n; i++)
3892 : {
3893 161 : GEN c1, c2, b, a = gel(I,i);
3894 161 : gel(I,i) = gen_1;
3895 161 : if (ideal_is1(a)) continue;
3896 :
3897 63 : c1 = gel(A,i);
3898 63 : c2 = gel(A,i+1);
3899 63 : b = gel(I,i+1);
3900 63 : if (ideal_is1(b))
3901 : {
3902 0 : gel(A,i) = c2;
3903 0 : gel(A,i+1) = gneg(c1);
3904 0 : gel(I,i+1) = a;
3905 : }
3906 : else
3907 : {
3908 63 : pari_sp av2 = avma;
3909 : GEN x, y, z, t, c;
3910 63 : nfidealdet1(nf,a,b, &x,&y,&z,&t);
3911 63 : x = RgC_add(nfC_nf_mul(nf, c1, x), nfC_nf_mul(nf, c2, y));
3912 63 : y = RgC_add(nfC_nf_mul(nf, c1, z), nfC_nf_mul(nf, c2, t));
3913 63 : (void)gc_all(av2, 2, &x,&y);
3914 63 : gel(A,i) = x;
3915 63 : gel(A,i+1) = y;
3916 63 : gel(I,i+1) = Q_primitive_part(idealmul(nf,a,b), &c);
3917 63 : if (c) gel(A,i+1) = nfC_nf_mul(nf, gel(A,i+1), c);
3918 : }
3919 : }
3920 28 : gel(M,1) = A;
3921 28 : gel(M,2) = I; return gc_GEN(av, M);
3922 : }
3923 :
3924 : /* Given bnf and a proj. module (or a t_POL -> rnfpseudobasis), and outputs a
3925 : * basis if it is free, an n+1-generating set if it is not */
3926 : GEN
3927 21 : rnfbasis(GEN bnf, GEN M)
3928 : {
3929 21 : pari_sp av = avma;
3930 : long j, n;
3931 : GEN nf, A, I, cl, col, a;
3932 :
3933 21 : bnf = checkbnf(bnf); nf = bnf_get_nf(bnf);
3934 21 : M = get_module(nf, M, "rnfbasis");
3935 21 : I = gel(M,2); n = lg(I)-1;
3936 98 : j = 1; while (j < n && ideal_is1(gel(I,j))) j++;
3937 21 : if (j < n) { M = rnfsteinitz(nf,M); I = gel(M,2); }
3938 21 : A = gel(M,1);
3939 21 : col= gel(A,n); A = vecslice(A, 1, n-1);
3940 21 : cl = gel(I,n);
3941 21 : a = gen_if_principal(bnf, cl);
3942 21 : if (!a)
3943 : {
3944 7 : GEN v = idealtwoelt(nf, cl);
3945 7 : A = vec_append(A, gmul(gel(v,1), col));
3946 7 : a = gel(v,2);
3947 : }
3948 21 : A = vec_append(A, nfC_nf_mul(nf, col, a));
3949 21 : return gc_GEN(av, A);
3950 : }
3951 :
3952 : /* Given a Z_K-module M (or a polynomial => rnfpseudobasis) outputs a
3953 : * Z_K-basis in HNF if it exists, zero if not */
3954 : GEN
3955 7 : rnfhnfbasis(GEN bnf, GEN M)
3956 : {
3957 7 : pari_sp av = avma;
3958 : long j, l;
3959 : GEN nf, A, I, a;
3960 :
3961 7 : bnf = checkbnf(bnf); nf = bnf_get_nf(bnf);
3962 7 : if (typ(M) == t_POL) M = rnfpseudobasis(nf, M);
3963 : else
3964 : {
3965 7 : if (typ(M) != t_VEC) pari_err_TYPE("rnfhnfbasis", M);
3966 7 : if (lg(M) == 5) M = mkvec2(gel(M,1), gel(M,2));
3967 7 : M = nfhnf(nf, M); /* in case M is not in HNF */
3968 : }
3969 7 : A = shallowcopy(gel(M,1));
3970 7 : I = gel(M,2); l = lg(A);
3971 42 : for (j = 1; j < l; j++)
3972 : {
3973 35 : if (ideal_is1(gel(I,j))) continue;
3974 14 : a = gen_if_principal(bnf, gel(I,j));
3975 14 : if (!a) return gc_const(av, gen_0);
3976 14 : gel(A,j) = nfC_nf_mul(nf, gel(A,j), a);
3977 : }
3978 7 : return gc_GEN(av,A);
3979 : }
3980 :
3981 : long
3982 7 : rnfisfree(GEN bnf, GEN M)
3983 : {
3984 7 : pari_sp av = avma;
3985 : GEN nf, P, I;
3986 : long l, j;
3987 :
3988 7 : bnf = checkbnf(bnf);
3989 7 : if (is_pm1( bnf_get_no(bnf) )) return 1;
3990 0 : nf = bnf_get_nf(bnf);
3991 0 : M = get_module(nf, M, "rnfisfree");
3992 0 : I = gel(M,2); l = lg(I); P = NULL;
3993 0 : for (j = 1; j < l; j++)
3994 0 : if (!ideal_is1(gel(I,j))) P = P? idealmul(nf, P, gel(I,j)): gel(I,j);
3995 0 : return gc_long(av, P? gequal0( isprincipal(bnf,P) ): 1);
3996 : }
3997 :
3998 : /**********************************************************************/
3999 : /** **/
4000 : /** COMPOSITUM OF TWO NUMBER FIELDS **/
4001 : /** **/
4002 : /**********************************************************************/
4003 : static GEN
4004 26698 : compositum_fix(GEN nf, GEN A)
4005 : {
4006 : int ok;
4007 26698 : if (nf)
4008 : {
4009 980 : A = RgXQX_red(A, nf_get_pol(nf));
4010 980 : A = Q_primpart(liftpol_shallow(A)); RgX_check_ZXX(A,"polcompositum");
4011 980 : ok = nfissquarefree(nf,A);
4012 : }
4013 : else
4014 : {
4015 25718 : A = Q_primpart(A); RgX_check_ZX(A,"polcompositum");
4016 25718 : ok = ZX_is_squarefree(A);
4017 : }
4018 26698 : if (!ok) pari_err_DOMAIN("polcompositum","issquarefree(arg)","=",gen_0,A);
4019 26691 : return A;
4020 : }
4021 : #define next_lambda(a) (a>0 ? -a : 1-a)
4022 :
4023 : static long
4024 511 : nfcompositum_lambda(GEN nf, GEN A, GEN B, long lambda)
4025 : {
4026 511 : pari_sp av = avma;
4027 : forprime_t S;
4028 511 : GEN T = nf_get_pol(nf);
4029 511 : long vT = varn(T);
4030 : ulong p;
4031 511 : init_modular_big(&S);
4032 511 : p = u_forprime_next(&S);
4033 : while (1)
4034 42 : {
4035 : GEN Hp, Tp, a;
4036 553 : if (DEBUGLEVEL>4) err_printf("Trying lambda = %ld\n", lambda);
4037 553 : a = ZXX_to_FlxX(RgX_rescale(A, stoi(-lambda)), p, vT);
4038 553 : Tp = ZX_to_Flx(T, p);
4039 553 : Hp = FlxqX_composedsum(a, ZXX_to_FlxX(B, p, vT), Tp, p);
4040 553 : if (!FlxqX_is_squarefree(Hp, Tp, p))
4041 42 : { lambda = next_lambda(lambda); continue; }
4042 511 : if (DEBUGLEVEL>4) err_printf("Final lambda = %ld\n", lambda);
4043 511 : return gc_long(av, lambda);
4044 : }
4045 : }
4046 :
4047 : /* modular version */
4048 : GEN
4049 13461 : nfcompositum(GEN nf, GEN A, GEN B, long flag)
4050 : {
4051 13461 : pari_sp av = avma;
4052 : int same;
4053 : long v, k;
4054 : GEN C, D, LPRS;
4055 :
4056 13461 : if (typ(A)!=t_POL) pari_err_TYPE("polcompositum",A);
4057 13461 : if (typ(B)!=t_POL) pari_err_TYPE("polcompositum",B);
4058 13461 : if (degpol(A)<=0 || degpol(B)<=0) pari_err_CONSTPOL("polcompositum");
4059 13461 : v = varn(A);
4060 13461 : if (varn(B) != v) pari_err_VAR("polcompositum", A,B);
4061 13461 : if (nf)
4062 : {
4063 553 : nf = checknf(nf);
4064 546 : if (varncmp(v,nf_get_varn(nf))>=0) pari_err_PRIORITY("polcompositum", nf, ">=", v);
4065 : }
4066 13419 : same = (A == B || RgX_equal(A,B));
4067 13419 : A = compositum_fix(nf,A);
4068 13412 : B = same ? A: compositum_fix(nf,B);
4069 :
4070 13412 : D = LPRS = NULL; /* -Wall */
4071 13412 : k = same? -1: 1;
4072 13412 : if (nf)
4073 : {
4074 511 : long v0 = fetch_var();
4075 511 : GEN q, T = nf_get_pol(nf);
4076 511 : A = liftpol_shallow(A);
4077 511 : B = liftpol_shallow(B);
4078 511 : k = nfcompositum_lambda(nf, A, B, k);
4079 511 : if (flag&1)
4080 : {
4081 : GEN H0, H1;
4082 196 : GEN chgvar = deg1pol_shallow(stoi(k),pol_x(v0),v);
4083 196 : GEN B1 = poleval(QXQX_to_mod_shallow(B, T), chgvar);
4084 196 : C = RgX_resultant_all(QXQX_to_mod_shallow(A, T), B1, &q);
4085 196 : C = gsubst(C,v0,pol_x(v));
4086 196 : C = lift_if_rational(C);
4087 196 : H0 = gsubst(gel(q,2),v0,pol_x(v));
4088 196 : H1 = gsubst(gel(q,3),v0,pol_x(v));
4089 196 : if (typ(H0) != t_POL) H0 = scalarpol_shallow(H0,v);
4090 196 : if (typ(H1) != t_POL) H1 = scalarpol_shallow(H1,v);
4091 196 : H0 = lift_if_rational(H0);
4092 196 : H1 = lift_if_rational(H1);
4093 196 : LPRS = mkvec2(H0,H1);
4094 : }
4095 : else
4096 : {
4097 315 : C = nf_direct_compositum(nf, RgX_rescale(A,stoi(-k)), B);
4098 315 : setvarn(C, v); C = QXQX_to_mod_shallow(C, T);
4099 : }
4100 511 : C = RgX_normalize(C);
4101 : }
4102 : else
4103 : {
4104 12901 : B = leafcopy(B); setvarn(B,fetch_var_higher());
4105 3185 : C = (flag&1)? ZX_ZXY_resultant_all(A, B, &k, &LPRS)
4106 12901 : : ZX_compositum(A, B, &k);
4107 12901 : setvarn(C, v);
4108 : }
4109 : /* C = Res_Y (A(Y), B(X + kY)) guaranteed squarefree */
4110 13412 : if (flag & 2)
4111 10311 : C = mkvec(C);
4112 : else
4113 : {
4114 3101 : if (same)
4115 : {
4116 105 : D = RgX_rescale(A, stoi(1 - k));
4117 105 : if (nf) D = RgX_normalize(QXQX_to_mod_shallow(D, nf_get_pol(nf)));
4118 105 : C = RgX_div(C, D);
4119 105 : if (degpol(C) <= 0)
4120 0 : C = mkvec(D);
4121 : else
4122 105 : C = shallowconcat(nf? gel(nffactor(nf,C),1): ZX_DDF(C), D);
4123 : }
4124 : else
4125 2996 : C = nf? gel(nffactor(nf,C),1): ZX_DDF(C);
4126 : }
4127 13412 : gen_sort_inplace(C, (void*)(nf?&cmp_RgX: &cmpii), &gen_cmp_RgX, NULL);
4128 13412 : if (flag&1)
4129 : { /* a,b,c root of A,B,C = compositum, c = b - k a */
4130 3381 : long i, l = lg(C);
4131 3381 : GEN a, b, mH0 = RgX_neg(gel(LPRS,1)), H1 = gel(LPRS,2);
4132 3381 : setvarn(mH0,v);
4133 3381 : setvarn(H1,v);
4134 6839 : for (i=1; i<l; i++)
4135 : {
4136 3458 : GEN D = gel(C,i);
4137 3458 : a = RgXQ_mul(mH0, nf? RgXQ_inv(H1,D): QXQ_inv(H1,D), D);
4138 3458 : b = gadd(pol_x(v), gmulsg(k,a));
4139 3458 : if (degpol(D) == 1) b = RgX_rem(b,D);
4140 3458 : gel(C,i) = mkvec4(D, mkpolmod(a,D), mkpolmod(b,D), stoi(-k));
4141 : }
4142 : }
4143 13412 : (void)delete_var();
4144 13412 : settyp(C, t_VEC);
4145 13412 : if (flag&2) C = gel(C,1);
4146 13412 : return gc_GEN(av, C);
4147 : }
4148 : GEN
4149 12908 : polcompositum0(GEN A, GEN B, long flag)
4150 12908 : { return nfcompositum(NULL,A,B,flag); }
4151 :
4152 : GEN
4153 91 : compositum(GEN pol1,GEN pol2) { return polcompositum0(pol1,pol2,0); }
4154 : GEN
4155 2870 : compositum2(GEN pol1,GEN pol2) { return polcompositum0(pol1,pol2,1); }
|
