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