Line data Source code
1 : /* Copyright (C) 2000-2004 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /*******************************************************************/
16 : /* */
17 : /* POLYNOMIAL FACTORIZATION IN A NUMBER FIELD */
18 : /* */
19 : /*******************************************************************/
20 : #include "pari.h"
21 : #include "paripriv.h"
22 :
23 : #define DEBUGLEVEL DEBUGLEVEL_nffactor
24 :
25 : static GEN nfsqff(GEN nf,GEN pol,long fl,GEN den);
26 : static int nfsqff_use_Trager(long n, long dpol);
27 :
28 : enum { FACTORS = 0, ROOTS, ROOTS_SPLIT };
29 :
30 : /* for nf_bestlift: reconstruction of algebraic integers known mod P^k,
31 : * P maximal ideal above p */
32 : typedef struct {
33 : long k; /* input known mod P^k */
34 : GEN p, pk; /* p^k = denom(prk^-1) [ assume pr unramified ]*/
35 : GEN prk; /* |.|^2 LLL-reduced basis (b_i) of P^k (NOT T2-reduced) */
36 : GEN iprk; /* den * prk^-1 */
37 : GEN GSmin; /* min |b_i^*|^2 */
38 :
39 : GEN Tp; /* Tpk mod p */
40 : GEN Tpk;
41 : GEN ZqProj;/* projector to Zp / P^k = Z/p^k[X] / Tpk */
42 :
43 : GEN tozk;
44 : GEN topow;
45 : GEN topowden; /* topow x / topowden = basistoalg(x) */
46 : GEN dn; /* NULL (we trust nf.zk) or a t_INT > 1 (an alg. integer has
47 : denominator dividing dn, when expressed on nf.zk */
48 : } nflift_t;
49 :
50 : typedef struct
51 : {
52 : nflift_t *L;
53 : GEN nf;
54 : GEN pol, polbase; /* leading coeff is a t_INT */
55 : GEN fact;
56 : GEN Br, bound, ZC, BS_2;
57 : } nfcmbf_t;
58 :
59 : /*******************************************************************/
60 : /* RATIONAL RECONSTRUCTION (use ratlift) */
61 : /*******************************************************************/
62 : static GEN
63 27603034 : lift_to_frac(GEN t, GEN N, GEN amax, GEN bmax, GEN den, GEN tden)
64 : {
65 27603034 : pari_sp av = avma;
66 : GEN a, b;
67 27603034 : if (signe(t) < 0) /* in case t is a centerlift */
68 : {
69 3069535 : if (abscmpii(t, amax) < 0) return icopy(t);
70 2805895 : t = addii(t, N);
71 : }
72 27339394 : if (tden)
73 : {
74 7662947 : a = Fp_center_i(Fp_mul(t, tden, N), N, shifti(N,-1));
75 7662947 : if (abscmpii(a, amax) < 0) return gc_upto(av, Qdivii(a, tden));
76 : }
77 20545312 : if (!Fp_ratlift(t, N, amax,bmax, &a,&b)) return gc_NULL(av);
78 20224285 : if (is_pm1(b)) return gc_GEN(av, a);
79 3245725 : if ((den && !dvdii(den,b)) || !is_pm1(gcdii(a,b))) return gc_NULL(av);
80 3239144 : return gc_GEN(av, mkfrac(a, b));
81 : }
82 :
83 : /* Compute rational lifting for all the components of P modulo N. Assume
84 : * Fp_ratlift preconditions are met; we allow centerlifts. If one component
85 : * fails, return NULL. If den != NULL, check that the deninators divide den;
86 : * assume (N, den) = 1. */
87 : GEN
88 5347515 : FpC_ratlift(GEN P, GEN N, GEN amax, GEN bmax, GEN den)
89 : {
90 5347515 : pari_sp av = avma;
91 : long j, l;
92 5347515 : GEN tden = NULL, Q = cgetg_copy(P, &l);
93 5347515 : if (l==1) return Q;
94 5347515 : if (den && cmpii(bmax, den) > 0) bmax = den;
95 31768473 : for (j = 1; j < l; ++j)
96 : {
97 26638017 : GEN a = lift_to_frac(gel(P,j), N, amax, bmax, den, tden);
98 26638017 : if (!a) return gc_NULL(av);
99 26420958 : if (typ(a) == t_FRAC)
100 : {
101 6802354 : GEN d = gel(a,2);
102 6802354 : tden = tden? (cmpii(tden, d) < 0? d: tden): d;
103 : }
104 26420958 : gel(Q,j) = a;
105 : }
106 5130456 : return Q;
107 : }
108 : GEN
109 364212 : FpX_ratlift(GEN P, GEN N, GEN amax, GEN bmax, GEN den)
110 : {
111 364212 : pari_sp av = avma;
112 : long j, l;
113 364212 : GEN tden = NULL, Q = cgetg_copy(P, &l);
114 364212 : Q[1] = P[1];
115 364212 : if (den && cmpii(bmax, den) > 0) bmax = den;
116 1218680 : for (j = 2; j < l; ++j)
117 : {
118 965017 : GEN a = lift_to_frac(gel(P,j), N, amax, bmax, den, tden);
119 965017 : if (!a) return gc_NULL(av);
120 854468 : if (typ(a) == t_FRAC)
121 : {
122 584380 : GEN d = gel(a,2);
123 584380 : tden = tden? (cmpii(tden, d) < 0? d: tden): d;
124 : }
125 854468 : gel(Q,j) = a;
126 : }
127 253663 : return Q;
128 : }
129 :
130 : GEN
131 2785305 : FpM_ratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN den)
132 : {
133 2785305 : pari_sp av = avma;
134 2785305 : long j, l = lg(M);
135 2785305 : GEN N = cgetg_copy(M, &l);
136 2785305 : if (l == 1) return N;
137 7831788 : for (j = 1; j < l; ++j)
138 : {
139 5253942 : GEN a = FpC_ratlift(gel(M, j), mod, amax, bmax, den);
140 5253942 : if (!a) return gc_NULL(av);
141 5046483 : gel(N,j) = a;
142 : }
143 2577846 : return N;
144 : }
145 :
146 : /*******************************************************************/
147 : /* GCD in K[X], K NUMBER FIELD */
148 : /*******************************************************************/
149 : /* P a nonzero ZXQX */
150 : static GEN
151 121395 : lead_simplify(GEN P)
152 : {
153 121395 : GEN x = leading_coeff(P); /* x a nonzero ZX or t_INT */
154 121395 : if (typ(x) == t_POL)
155 : {
156 3864 : if (degpol(x)) return x;
157 3577 : x = gel(x,2);
158 : }
159 121108 : return is_pm1(x)? NULL: x;
160 : }
161 : /* P,Q in Z[X,Y], T in Z[Y] irreducible. compute GCD in Q[Y]/(T)[X].
162 : *
163 : * M. Encarnacion "On a modular Algorithm for computing GCDs of polynomials
164 : * over number fields" (ISSAC'94).
165 : *
166 : * We procede as follows
167 : * 1:compute the gcd modulo primes discarding bad primes as they are detected.
168 : * 2:reconstruct the result via FpM_ratlift, stoping as soon as we get weird
169 : * denominators.
170 : * 3:if FpM_ratlift succeeds, try the full division.
171 : * Suppose accuracy is insufficient to get the result right: FpM_ratlift will
172 : * rarely succeed, and even if it does the polynomial we get has sensible
173 : * coefficients, so the full division will not be too costly.
174 : *
175 : * If not NULL, den must be a multiple of the denominator of the gcd,
176 : * for example the discriminant of T.
177 : *
178 : * NOTE: if T is not irreducible, nfgcd may loop forever, esp. if gcd | T */
179 : GEN
180 90330 : nfgcd_all(GEN P, GEN Q, GEN T, GEN den, GEN *Pnew)
181 : {
182 90330 : pari_sp btop, ltop = avma;
183 90330 : GEN lP, lQ, M, dsol, R, bo, sol, mod = NULL, lden = NULL;
184 90330 : long vP = varn(P), vT = varn(T), dT = degpol(T), dM = 0, dR;
185 : forprime_t S;
186 :
187 90330 : if (!signe(P)) { if (Pnew) *Pnew = pol_0(vT); return gcopy(Q); }
188 90295 : if (!signe(Q)) { if (Pnew) *Pnew = pol_1(vT); return gcopy(P); }
189 : /* Compute denominators */
190 90148 : if ((lP = lead_simplify(P)) && (lQ = lead_simplify(Q)))
191 : {
192 14839 : if (typ(lP) == t_INT && typ(lQ) == t_INT)
193 14692 : lden = powiu(gcdii(lP, lQ), dT);
194 147 : else if (typ(lP) == t_INT)
195 7 : lden = gcdii(powiu(lP, dT), ZX_resultant(lQ, T));
196 140 : else if (typ(lQ) == t_INT)
197 0 : lden = gcdii(powiu(lQ, dT), ZX_resultant(lP, T));
198 : else
199 140 : lden = gcdii(ZX_resultant(lP, T), ZX_resultant(lQ, T));
200 14839 : if (is_pm1(lden)) lden = NULL;
201 14839 : if (den && lden) den = mulii(den, lden);
202 : }
203 90148 : init_modular_small(&S);
204 90148 : btop = avma;
205 : for(;;)
206 7775 : {
207 97923 : ulong p = u_forprime_next(&S);
208 : GEN Tp;
209 97923 : if (!p) pari_err_OVERFLOW("nfgcd [ran out of primes]");
210 : /*Discard primes dividing disc(T) or lc(PQ) */
211 97923 : if (lden && !umodiu(lden, p)) continue;
212 97923 : Tp = ZX_to_Flx(T,p);
213 97923 : if (!Flx_is_squarefree(Tp, p)) continue;
214 : /*Discard primes when modular gcd does not exist*/
215 97923 : if ((R = FlxqX_safegcd(ZXX_to_FlxX(P,p,vT),
216 : ZXX_to_FlxX(Q,p,vT),
217 0 : Tp, p)) == NULL) continue;
218 97923 : dR = degpol(R);
219 97923 : if (dR == 0) { set_avma(ltop); if (Pnew) *Pnew = P; return pol_1(vP); }
220 10946 : if (mod && dR > dM) continue; /* p divides Res(P/gcd, Q/gcd). Discard. */
221 :
222 10946 : R = FlxX_to_Flm(R, dT);
223 : /* previous primes divided Res(P/gcd, Q/gcd)? Discard them. */
224 10946 : if (!mod || dR < dM) { M = ZM_init_CRT(R, p); mod = utoipos(p); dM = dR; continue; }
225 7775 : (void)ZM_incremental_CRT(&M,R, &mod,p);
226 7775 : if (gc_needed(btop, 1))
227 : {
228 0 : if (DEBUGMEM>1) pari_warn(warnmem,"nfgcd");
229 0 : (void)gc_all(btop, 2, &M, &mod);
230 : }
231 : /* I suspect it must be better to take amax > bmax*/
232 7775 : bo = sqrti(shifti(mod, -1));
233 7775 : if ((sol = FpM_ratlift(M, mod, bo, bo, den)) == NULL) continue;
234 3211 : sol = RgM_to_RgXX(sol,vP,vT);
235 3211 : dsol = Q_primpart(sol);
236 :
237 3211 : if (!ZXQX_dvd(Q, dsol, T)) continue;
238 3171 : if (Pnew)
239 : {
240 364 : *Pnew = RgXQX_pseudodivrem(P, dsol, T, &R);
241 364 : if (signe(R)) continue;
242 : }
243 2807 : else if (!ZXQX_dvd(P, dsol, T)) continue;
244 3171 : return gc_all(ltop, Pnew? 2: 1, &dsol, Pnew); /* both remainders are 0 */
245 : }
246 : }
247 : GEN
248 48048 : nfgcd(GEN P, GEN Q, GEN T, GEN den)
249 48048 : { return nfgcd_all(P, Q, T, den, NULL); }
250 :
251 : GEN
252 5544 : ZXQX_gcd(GEN P, GEN Q, GEN T)
253 5544 : { return nfgcd_all(P, Q, T, NULL, NULL); }
254 :
255 : GEN
256 3045 : QXQX_gcd(GEN P, GEN Q, GEN T)
257 : {
258 3045 : pari_sp av = avma;
259 3045 : GEN P1 = Q_remove_denom(P, NULL);
260 3045 : GEN Q1 = Q_remove_denom(Q, NULL);
261 3045 : return gc_upto(av, ZXQX_gcd(P1, Q1, T));
262 : }
263 :
264 : int
265 49826 : nfissquarefree(GEN nf, GEN x)
266 : {
267 49826 : pari_sp av = avma;
268 49826 : GEN g, y = RgX_deriv(x);
269 49826 : if (RgX_is_rational(x)) g = QX_gcd(x, y);
270 : else
271 : {
272 46774 : GEN T = get_nfpol(nf,&nf);
273 46774 : x = Q_primpart( liftpol_shallow(x) );
274 46774 : y = Q_primpart( liftpol_shallow(y) );
275 46774 : g = nfgcd(x, y, T, nf? nf_get_index(nf): NULL);
276 : }
277 49826 : return gc_bool(av, degpol(g) == 0);
278 : }
279 :
280 : /*******************************************************************/
281 : /* FACTOR OVER (Z_K/pr)[X] --> FqX_factor */
282 : /*******************************************************************/
283 : GEN
284 7 : nffactormod(GEN nf, GEN x, GEN pr)
285 : {
286 7 : long j, l, vx = varn(x), vn;
287 7 : pari_sp av = avma;
288 : GEN F, E, rep, xrd, modpr, T, p;
289 :
290 7 : nf = checknf(nf);
291 7 : vn = nf_get_varn(nf);
292 7 : if (typ(x)!=t_POL) pari_err_TYPE("nffactormod",x);
293 7 : if (varncmp(vx,vn) >= 0) pari_err_PRIORITY("nffactormod", x, ">=", vn);
294 :
295 7 : modpr = nf_to_Fq_init(nf, &pr, &T, &p);
296 7 : xrd = nfX_to_FqX(x, nf, modpr);
297 7 : rep = FqX_factor(xrd,T,p);
298 7 : settyp(rep, t_MAT);
299 7 : F = gel(rep,1); l = lg(F);
300 7 : E = gel(rep,2); settyp(E, t_COL);
301 14 : for (j = 1; j < l; j++) {
302 7 : gel(F,j) = FqX_to_nfX(gel(F,j), modpr);
303 7 : gel(E,j) = stoi(E[j]);
304 : }
305 7 : return gc_GEN(av, rep);
306 : }
307 :
308 : /*******************************************************************/
309 : /* MAIN ROUTINES nfroots / nffactor */
310 : /*******************************************************************/
311 : static GEN
312 38369 : QXQX_normalize(GEN P, GEN T)
313 : {
314 38369 : GEN P0 = leading_coeff(P);
315 38369 : long t = typ(P0);
316 38369 : if (t == t_POL)
317 : {
318 623 : if (degpol(P0)) return RgXQX_RgXQ_mul(P, QXQ_inv(P0,T), T);
319 518 : P0 = gel(P0,2); t = typ(P0);
320 : }
321 : /* t = t_INT/t_FRAC */
322 38264 : if (t == t_INT && is_pm1(P0) && signe(P0) > 0) return P; /* monic */
323 7476 : return RgX_Rg_div(P, P0);
324 : }
325 : /* assume leading term of P is an integer */
326 : static GEN
327 36645 : RgX_int_normalize(GEN P)
328 : {
329 36645 : GEN P0 = leading_coeff(P);
330 : /* cater for t_POL */
331 36645 : if (typ(P0) == t_POL)
332 : {
333 0 : P0 = gel(P0,2); /* nonzero constant */
334 0 : P = shallowcopy(P);
335 0 : gel(P,lg(P)-1) = P0; /* now leading term is a t_INT */
336 : }
337 36645 : if (typ(P0) != t_INT) pari_err_BUG("RgX_int_normalize");
338 36645 : if (is_pm1(P0)) return signe(P0) > 0? P: RgX_neg(P);
339 16464 : return RgX_Rg_div(P, P0);
340 : }
341 :
342 : /* discard change of variable if nf is of the form [nf,c] as return by nfinit
343 : * for nonmonic polynomials */
344 : static GEN
345 12752 : proper_nf(GEN nf)
346 12752 : { return (lg(nf) == 3)? gel(nf,1): nf; }
347 :
348 : /* if *pnf = NULL replace if by a "quick" K = nfinit(T), ensuring maximality
349 : * by small primes only. Return a multiplicative bound for the denominator of
350 : * algebraic integers in Z_K in terms of K.zk */
351 : static GEN
352 36451 : fix_nf(GEN *pnf, GEN *pT, GEN *pA)
353 : {
354 36451 : GEN nf, NF, P, q, D, T = *pT;
355 : nfmaxord_t S;
356 : long i, l, lim;
357 :
358 36451 : if (*pnf) return gen_1;
359 12752 : lim = GP_DATA->factorlimit + 1;
360 12752 : nfmaxord(&S, mkvec2(T, utoipos(lim)), nf_PARTIALFACT);
361 12752 : NF = nfinit_complete(&S, 0, DEFAULTPREC);
362 12752 : *pnf = nf = proper_nf(NF);
363 12752 : if (nf != NF) { /* t_POL defining base field changed (not monic) */
364 14 : GEN A = *pA, a = cgetg_copy(A, &l);
365 14 : GEN rev = gel(NF,2), pow, dpow;
366 :
367 14 : *pT = T = nf_get_pol(nf); /* need to update T */
368 14 : pow = QXQ_powers(lift_shallow(rev), degpol(T)-1, T);
369 14 : pow = Q_remove_denom(pow, &dpow);
370 14 : a[1] = A[1];
371 56 : for (i=2; i<l; i++) {
372 42 : GEN c = gel(A,i);
373 42 : if (typ(c) == t_POL) c = QX_ZXQV_eval(c, pow, dpow);
374 42 : gel(a,i) = c;
375 : }
376 14 : *pA = Q_primpart(a); /* need to update A */
377 : }
378 :
379 12752 : D = nf_get_disc(nf); if (is_pm1(D)) return gen_1;
380 : /* D may be incorrect */
381 12752 : P = nf_get_ramified_primes(nf); l = lg(P);
382 32845 : for (i = 1, q = gen_1; i < l; i++)
383 : {
384 20093 : GEN p = gel(P,i);
385 20093 : if (cmpiu(p, lim) >= 0 && !BPSW_psp(p)) q = mulii(q, p);
386 : }
387 12752 : return q;
388 : }
389 :
390 : /* lt(A) is an integer; ensure it is not a constant t_POL. In place */
391 : static void
392 36619 : ensure_lt_INT(GEN A)
393 : {
394 36619 : long n = lg(A)-1;
395 36619 : GEN lt = gel(A,n);
396 36808 : while (typ(lt) != t_INT) gel(A,n) = lt = gel(lt,2);
397 36619 : }
398 :
399 : /* set B = A/gcd(A,A'), squarefree */
400 : static GEN
401 36605 : get_nfsqff_data(GEN *pnf, GEN *pT, GEN *pA, GEN *pB, GEN *ptbad)
402 : {
403 36605 : GEN den, bad, D, B, A = *pA, T = *pT;
404 36605 : long n = degpol(T);
405 :
406 36605 : A = Q_primpart( QXQX_normalize(A, T) );
407 36605 : if (nfsqff_use_Trager(n, degpol(A)))
408 : {
409 175 : *pnf = T;
410 175 : bad = den = absi_shallow(ZX_disc(T));
411 175 : if (is_pm1(leading_coeff(T))) den = indexpartial(T, den);
412 : }
413 : else
414 : {
415 36430 : den = fix_nf(pnf, &T, &A);
416 36430 : bad = nf_get_index(*pnf);
417 36430 : if (den != gen_1) bad = mulii(bad, den);
418 : }
419 36605 : D = nfgcd_all(A, RgX_deriv(A), T, bad, &B);
420 36605 : if (degpol(D)) B = Q_primpart( QXQX_normalize(B, T) );
421 36605 : if (ptbad) *ptbad = bad;
422 36605 : *pA = A;
423 36605 : *pB = B; ensure_lt_INT(B);
424 36605 : *pT = T; return den;
425 : }
426 :
427 : /* return the roots of pol in nf */
428 : GEN
429 43213 : nfroots(GEN nf,GEN pol)
430 : {
431 43213 : pari_sp av = avma;
432 : GEN z, A, B, T, den;
433 : long d, dT;
434 :
435 43213 : if (!nf) return nfrootsQ(pol);
436 29997 : T = get_nfpol(nf, &nf);
437 29997 : RgX_check_ZX(T,"nfroots");
438 29997 : A = RgX_nffix("nfroots", T,pol,1);
439 29997 : d = degpol(A);
440 29997 : if (d < 0) pari_err_ROOTS0("nfroots");
441 29997 : if (d == 0) return cgetg(1,t_COL);
442 29997 : if (d == 1)
443 : {
444 315 : A = QXQX_normalize(A,T);
445 315 : A = mkpolmod(gneg_i(gel(A,2)), T);
446 315 : return gc_GEN(av, mkcol(A));
447 : }
448 29682 : dT = degpol(T);
449 29682 : if (dT == 1) return gc_upto(av, nfrootsQ(simplify_shallow(A)));
450 :
451 24047 : den = get_nfsqff_data(&nf, &T, &A, &B, NULL);
452 24047 : if (RgX_is_ZX(B))
453 : {
454 8594 : GEN v = gel(ZX_factor(B), 1);
455 8594 : long i, l = lg(v), p = mael(factoru(dT),1,1); /* smallest prime divisor */
456 8594 : z = cgetg(1, t_VEC);
457 22935 : for (i = 1; i < l; i++)
458 : {
459 14341 : GEN b = gel(v,i); /* irreducible / Q */
460 14341 : long db = degpol(b);
461 14341 : if (db != 1 && degpol(b) < p) continue;
462 14229 : z = shallowconcat(z, nfsqff(nf, b, ROOTS, den));
463 : }
464 : }
465 : else
466 15453 : z = nfsqff(nf,B, ROOTS, den);
467 24047 : z = gc_upto(av, QXQV_to_mod(z, T));
468 24047 : gen_sort_inplace(z, (void*)&cmp_RgX, &cmp_nodata, NULL);
469 24047 : settyp(z, t_COL); return z;
470 : }
471 :
472 : static GEN
473 477366 : _norml2(GEN x) { return RgC_fpnorml2(x, DEFAULTPREC); }
474 :
475 : /* return a minimal lift of elt modulo id, as a ZC */
476 : static GEN
477 263399 : nf_bestlift(GEN elt, GEN bound, nflift_t *L)
478 : {
479 : GEN u;
480 263399 : long i,l = lg(L->prk), t = typ(elt);
481 263399 : if (t != t_INT)
482 : {
483 16069 : if (t == t_POL) elt = ZM_ZX_mul(L->tozk, elt);
484 16069 : u = ZM_ZC_mul(L->iprk,elt);
485 335637 : for (i=1; i<l; i++) gel(u,i) = diviiround(gel(u,i), L->pk);
486 : }
487 : else
488 : {
489 247330 : u = ZC_Z_mul(gel(L->iprk,1), elt);
490 1925056 : for (i=1; i<l; i++) gel(u,i) = diviiround(gel(u,i), L->pk);
491 247330 : elt = scalarcol(elt, l-1);
492 : }
493 263399 : u = ZC_sub(elt, ZM_ZC_mul(L->prk, u));
494 263399 : if (bound && gcmp(_norml2(u), bound) > 0) u = NULL;
495 263399 : return u;
496 : }
497 :
498 : /* Warning: return L->topowden * (best lift). */
499 : static GEN
500 158555 : nf_bestlift_to_pol(GEN elt, GEN bound, nflift_t *L)
501 : {
502 158555 : pari_sp av = avma;
503 158555 : GEN u,v = nf_bestlift(elt,bound,L);
504 158555 : if (!v) return NULL;
505 149026 : if (ZV_isscalar(v))
506 : {
507 65269 : if (L->topowden)
508 65269 : u = mulii(L->topowden, gel(v,1));
509 : else
510 0 : u = icopy(gel(v,1));
511 65269 : u = gc_INT(av, u);
512 : }
513 : else
514 : {
515 83757 : v = gclone(v); set_avma(av);
516 83757 : u = RgV_dotproduct(L->topow, v);
517 83757 : gunclone(v);
518 : }
519 149026 : return u;
520 : }
521 :
522 : /* return the T->powden * (lift of pol with coefficients of T2-norm <= C)
523 : * if it exists. */
524 : static GEN
525 36624 : nf_pol_lift(GEN pol, GEN bound, nflift_t *L)
526 : {
527 36624 : long i, l = lg(pol);
528 36624 : GEN x = cgetg(l,t_POL);
529 :
530 36624 : x[1] = pol[1];
531 171700 : for (i=l-1; i>1; i--)
532 : {
533 144605 : GEN t = nf_bestlift_to_pol(gel(pol,i), bound, L);
534 144605 : if (!t) return NULL;
535 135076 : gel(x,i) = t;
536 : }
537 27095 : return x;
538 : }
539 :
540 : static GEN
541 0 : zerofact(long v)
542 : {
543 0 : GEN z = cgetg(3, t_MAT);
544 0 : gel(z,1) = mkcol(pol_0(v));
545 0 : gel(z,2) = mkcol(gen_1); return z;
546 : }
547 :
548 : /* Return the factorization of A in Q[X]/(T) in rep [pre-allocated with
549 : * cgetg(3,t_MAT)], reclaiming all memory between avma and rep.
550 : * y is the vector of irreducible factors of B = Q_primpart( A/gcd(A,A') ).
551 : * Bad primes divide 'bad' */
552 : static void
553 12572 : fact_from_sqff(GEN rep, GEN A, GEN B, GEN y, GEN T, GEN bad)
554 : {
555 12572 : pari_sp av = (pari_sp)rep;
556 12572 : long n = lg(y)-1;
557 : GEN ex;
558 :
559 12572 : if (A != B)
560 : { /* not squarefree */
561 140 : if (n == 1)
562 : { /* perfect power, simple ! */
563 28 : long e = degpol(A) / degpol(gel(y,1));
564 28 : y = gc_upto(av, QXQXV_to_mod(y, T));
565 28 : ex = mkcol(utoipos(e));
566 : }
567 : else
568 : { /* compute valuations mod a prime of degree 1 (avoid coeff explosion) */
569 112 : GEN quo, p, r, Bp, lb = leading_coeff(B), E = cgetalloc(n+1,t_VECSMALL);
570 112 : pari_sp av1 = avma;
571 : ulong pp;
572 : long j;
573 : forprime_t S;
574 112 : u_forprime_init(&S, degpol(T), ULONG_MAX);
575 280 : for (; ; set_avma(av1))
576 : {
577 392 : pp = u_forprime_next(&S);
578 392 : if (! umodiu(bad,pp) || !umodiu(lb, pp)) continue;
579 371 : p = utoipos(pp);
580 371 : r = FpX_oneroot(T, p);
581 371 : if (!r) continue;
582 203 : Bp = FpXY_evalx(B, r, p);
583 203 : if (FpX_is_squarefree(Bp, p)) break;
584 : }
585 :
586 112 : quo = FpXY_evalx(Q_primpart(A), r, p);
587 259 : for (j=n; j>=2; j--)
588 : {
589 147 : GEN junk, fact = Q_remove_denom(gel(y,j), &junk);
590 147 : long e = 0;
591 147 : fact = FpXY_evalx(fact, r, p);
592 308 : for(;; e++)
593 308 : {
594 455 : GEN q = FpX_divrem(quo,fact,p,ONLY_DIVIDES);
595 455 : if (!q) break;
596 308 : quo = q;
597 : }
598 147 : E[j] = e;
599 : }
600 112 : E[1] = degpol(quo) / degpol(gel(y,1));
601 112 : y = gc_upto(av, QXQXV_to_mod(y, T));
602 112 : ex = zc_to_ZC(E); pari_free((void*)E);
603 : }
604 : }
605 : else
606 : {
607 12432 : y = gc_upto(av, QXQXV_to_mod(y, T));
608 12432 : ex = const_col(n, gen_1);
609 : }
610 12572 : gel(rep,1) = y; settyp(y, t_COL);
611 12572 : gel(rep,2) = ex;
612 12572 : }
613 :
614 : /* return the factorization of polynomial pol in nf */
615 : static GEN
616 12901 : nffactor_i(GEN nf,GEN T,GEN pol)
617 : {
618 12901 : GEN bad, A, B, y, den, rep = cgetg(3, t_MAT);
619 12901 : pari_sp av = avma;
620 : long dA;
621 : pari_timer ti;
622 :
623 12901 : if (DEBUGLEVEL>2) { timer_start(&ti); err_printf("\nEntering nffactor:\n"); }
624 12901 : A = RgX_nffix("nffactor",T,pol,1);
625 12901 : dA = degpol(A);
626 12901 : if (dA <= 0) {
627 7 : set_avma((pari_sp)(rep + 3));
628 7 : return (dA == 0)? trivial_fact(): zerofact(varn(pol));
629 : }
630 12894 : if (dA == 1) {
631 : GEN c;
632 210 : A = Q_primpart( QXQX_normalize(A, T) );
633 210 : A = gc_GEN(av, A); c = gel(A,2);
634 210 : if (typ(c) == t_POL && degpol(c) > 0) gel(A,2) = mkpolmod(c, ZX_copy(T));
635 210 : gel(rep,1) = mkcol(A);
636 210 : gel(rep,2) = mkcol(gen_1); return rep;
637 : }
638 12684 : if (degpol(T) == 1) return gc_upto(av, QX_factor(simplify_shallow(A)));
639 :
640 12558 : den = get_nfsqff_data(&nf, &T, &A, &B, &bad);
641 12558 : if (DEBUGLEVEL>2) timer_printf(&ti, "squarefree test");
642 12558 : if (RgX_is_ZX(B))
643 : {
644 11970 : GEN v = gel(ZX_factor(B), 1);
645 11970 : long i, l = lg(v);
646 11970 : y = cgetg(1, t_VEC);
647 23996 : for (i = 1; i < l; i++)
648 : {
649 12026 : GEN b = gel(v,i); /* irreducible / Q */
650 12026 : y = shallowconcat(y, nfsqff(nf, b, 0, den));
651 : }
652 : }
653 : else
654 588 : y = nfsqff(nf,B, 0, den);
655 12558 : if (DEBUGLEVEL>3) err_printf("number of factor(s) found: %ld\n", lg(y)-1);
656 :
657 12558 : fact_from_sqff(rep, A, B, y, T, bad);
658 12558 : return rep;
659 : }
660 :
661 : /* return the factorization of P in nf */
662 : GEN
663 12894 : nffactor(GEN nf, GEN P)
664 : {
665 12894 : GEN y, T = get_nfpol(nf, &nf);
666 12894 : if (!nf) RgX_check_ZX(T,"nffactor");
667 12894 : if (typ(P) == t_RFRAC)
668 : {
669 14 : pari_sp av = avma;
670 14 : GEN a = gel(P, 1), b = gel(P, 2);
671 14 : y = famat_inv_shallow(nffactor_i(nf, T, b));
672 14 : if (typ(a) == t_POL && varn(a) == varn(b))
673 7 : y = famat_mul_shallow(nffactor_i(nf, T, a), y);
674 14 : y = gc_GEN(av, y);
675 : }
676 : else
677 12880 : y = nffactor_i(nf, T, P);
678 12894 : return sort_factor_pol(y, cmp_RgX);
679 : }
680 :
681 : /* assume x scalar or t_COL, G t_MAT */
682 : static GEN
683 202270 : arch_for_T2(GEN G, GEN x)
684 : {
685 2478 : return (typ(x) == t_COL)? RgM_RgC_mul(G,x)
686 204748 : : RgC_Rg_mul(gel(G,1),x);
687 : }
688 :
689 : /* polbase a zkX with t_INT leading coeff; return a bound for T_2(P),
690 : * P | polbase in C[X]. NB: Mignotte bound: A | S ==>
691 : * |a_i| <= binom(d-1, i-1) || S ||_2 + binom(d-1, i) lc(S)
692 : *
693 : * Apply to sigma(S) for all embeddings sigma, then take the L_2 norm over
694 : * sigma, then take the sup over i */
695 : static GEN
696 12075 : nf_Mignotte_bound(GEN nf, GEN polbase)
697 12075 : { GEN lS = leading_coeff(polbase); /* t_INT */
698 : GEN p1, C, N2, binlS, bin;
699 12075 : long prec = nf_get_prec(nf), n = nf_get_degree(nf), r1 = nf_get_r1(nf);
700 12075 : long i, j, d = degpol(polbase);
701 :
702 12075 : binlS = bin = vecbinomial(d-1);
703 12075 : if (!isint1(lS)) binlS = ZC_Z_mul(bin,lS);
704 :
705 12075 : N2 = cgetg(n+1, t_VEC);
706 : for (;;)
707 0 : {
708 12075 : GEN G = nf_get_G(nf), matGS = cgetg(d+2, t_MAT);
709 :
710 133721 : for (j=0; j<=d; j++) gel(matGS,j+1) = arch_for_T2(G, gel(polbase,j+2));
711 12075 : matGS = shallowtrans(matGS);
712 19733 : for (j=1; j <= r1; j++) /* N2[j] = || sigma_j(S) ||_2 */
713 : {
714 7658 : GEN c = sqrtr( _norml2(gel(matGS,j)) );
715 7658 : gel(N2,j) = c; if (!signe(c)) goto PRECPB;
716 : }
717 31248 : for ( ; j <= n; j+=2)
718 : {
719 19173 : GEN q1 = _norml2(gel(matGS, j));
720 19173 : GEN q2 = _norml2(gel(matGS, j+1));
721 19173 : GEN c = sqrtr( gmul2n(addrr(q1, q2), -1) );
722 19173 : gel(N2,j) = gel(N2,j+1) = c; if (!signe(c)) goto PRECPB;
723 : }
724 12075 : break; /* done */
725 0 : PRECPB:
726 0 : prec = precdbl(prec);
727 0 : nf = nfnewprec_shallow(nf, prec);
728 0 : if (DEBUGLEVEL>1) pari_warn(warnprec, "nf_factor_bound", prec);
729 : }
730 :
731 : /* Take sup over 0 <= i <= d of
732 : * sum_j | binom(d-1, i-1) ||sigma_j(S)||_2 + binom(d-1,i) lc(S) |^2 */
733 :
734 : /* i = 0: n lc(S)^2 */
735 12075 : C = mului(n, sqri(lS));
736 : /* i = d: sum_sigma ||sigma(S)||_2^2 */
737 12075 : p1 = gnorml2(N2); if (gcmp(C, p1) < 0) C = p1;
738 109571 : for (i = 1; i < d; i++)
739 : {
740 97496 : GEN B = gel(bin,i), L = gel(binlS,i+1);
741 97496 : GEN s = sqrr(addri(mulir(B, gel(N2,1)), L)); /* j=1 */
742 627193 : for (j = 2; j <= n; j++) s = addrr(s, sqrr(addri(mulir(B, gel(N2,j)), L)));
743 97496 : if (mpcmp(C, s) < 0) C = s;
744 : }
745 12075 : return C;
746 : }
747 :
748 : /* return a bound for T_2(P), P | polbase
749 : * max |b_i|^2 <= 3^{3/2 + d} / (4 \pi d) [P]_2,
750 : * where [P]_2 is Bombieri's 2-norm
751 : * Sum over conjugates */
752 : static GEN
753 12075 : nf_Beauzamy_bound(GEN nf, GEN polbase)
754 : {
755 : GEN lt, C, s, POL, bin;
756 12075 : long d = degpol(polbase), n = nf_get_degree(nf), prec = nf_get_prec(nf);
757 12075 : bin = vecbinomial(d);
758 12075 : POL = polbase + 2;
759 : /* compute [POL]_2 */
760 : for (;;)
761 0 : {
762 12075 : GEN G = nf_get_G(nf);
763 : long i;
764 :
765 12075 : s = real_0(prec);
766 133721 : for (i=0; i<=d; i++)
767 : {
768 121646 : GEN c = gel(POL,i);
769 121646 : if (gequal0(c)) continue;
770 80624 : c = _norml2(arch_for_T2(G,c));
771 80624 : if (!signe(c)) goto PRECPB;
772 : /* s += T2(POL[i]) / binomial(d,i) */
773 80624 : s = addrr(s, divri(c, gel(bin,i+1)));
774 : }
775 12075 : break;
776 0 : PRECPB:
777 0 : prec = precdbl(prec);
778 0 : nf = nfnewprec_shallow(nf, prec);
779 0 : if (DEBUGLEVEL>1) pari_warn(warnprec, "nf_factor_bound", prec);
780 : }
781 12075 : lt = leading_coeff(polbase);
782 12075 : s = mulri(s, muliu(sqri(lt), n));
783 12075 : C = powruhalf(utor(3,DEFAULTPREC), 3 + 2*d); /* 3^{3/2 + d} */
784 12075 : return divrr(mulrr(C, s), mulur(d, mppi(DEFAULTPREC)));
785 : }
786 :
787 : static GEN
788 12075 : nf_factor_bound(GEN nf, GEN polbase)
789 : {
790 12075 : pari_sp av = avma;
791 12075 : GEN a = nf_Mignotte_bound(nf, polbase);
792 12075 : GEN b = nf_Beauzamy_bound(nf, polbase);
793 12075 : if (DEBUGLEVEL>2)
794 : {
795 0 : err_printf("Mignotte bound: %Ps\n",a);
796 0 : err_printf("Beauzamy bound: %Ps\n",b);
797 : }
798 12075 : return gc_upto(av, gmin(a, b));
799 : }
800 :
801 : /* True nf; return Bs: if r a root of sigma_i(P), |r| < Bs[i] */
802 : static GEN
803 19015 : nf_root_bounds(GEN nf, GEN P)
804 : {
805 : long lR, i, j, l, prec, r1;
806 : GEN Ps, R, V;
807 :
808 19015 : if (RgX_is_rational(P)) return polrootsbound(P, NULL);
809 5598 : r1 = nf_get_r1(nf);
810 5598 : P = Q_primpart(P);
811 5598 : prec = ZXX_max_lg(P) + 1;
812 5598 : l = lg(P);
813 5598 : if (nf_get_prec(nf) >= prec)
814 5590 : R = nf_get_roots(nf);
815 : else
816 8 : R = QX_complex_roots(nf_get_pol(nf), prec);
817 5598 : lR = lg(R);
818 5598 : V = cgetg(lR, t_VEC);
819 5598 : Ps = cgetg(l, t_POL); /* sigma (P) */
820 5598 : Ps[1] = P[1];
821 19727 : for (j=1; j<lg(R); j++)
822 : {
823 14129 : GEN r = gel(R,j);
824 101981 : for (i=2; i<l; i++) gel(Ps,i) = poleval(gel(P,i), r);
825 14129 : gel(V,j) = polrootsbound(Ps, NULL);
826 : }
827 5598 : return mkvec2(vecslice(V,1,r1), vecslice(V,r1+1,lg(V)-1));
828 : }
829 :
830 : /* return B such that, if x = sum x_i K.zk[i] in O_K, then ||x||_2^2 <= B T_2(x)
831 : * den = multiplicative bound for denom(x) [usually NULL, for 1, but when we
832 : * use nf_PARTIALFACT K.zk may not generate O_K] */
833 : GEN
834 23337 : nf_L2_bound(GEN nf, GEN den, GEN *pL)
835 : {
836 23337 : GEN M, L, prep, T = nf_get_pol(nf), tozk = nf_get_invzk(nf);
837 23337 : long prec = ZM_max_lg(tozk) + ZX_max_lg(T) + nbits2prec(degpol(T));
838 23337 : (void)initgaloisborne(nf, den? den: gen_1, prec, &L, &prep, NULL);
839 23337 : M = vandermondeinverse(L, RgX_gtofp(T,prec), den, prep);
840 23337 : if (pL) *pL = L;
841 23337 : return RgM_fpnorml2(RgM_mul(tozk,M), DEFAULTPREC);
842 : }
843 :
844 : /* sum_i L[i]^p */
845 : static GEN
846 21734 : normlp(GEN L, long p)
847 : {
848 21734 : long i, l = lg(L);
849 : GEN z;
850 21734 : if (l == 1) return gen_0;
851 11336 : z = gpowgs(gel(L,1), p);
852 27670 : for (i=2; i<l; i++) z = gadd(z, gpowgs(gel(L,i), p));
853 11336 : return z;
854 : }
855 : /* \sum_i deg(sigma_i) L[i]^p in dimension n (L may be a scalar
856 : * or [L1,L2], where Ld corresponds to the archimedean places of degree d) */
857 : static GEN
858 26368 : normTp(GEN L, long p, long n)
859 : {
860 26368 : if (typ(L) != t_VEC) return gmulsg(n, gpowgs(L, p));
861 10867 : return gadd(normlp(gel(L,1),p), gmul2n(normlp(gel(L,2),p), 1));
862 : }
863 :
864 : /* S = S0 + tS1, P = P0 + tP1 (Euclidean div. by t integer). For a true
865 : * factor (vS, vP), we have:
866 : * | S vS + P vP |^2 < Btra
867 : * This implies | S1 vS + P1 vP |^2 < Bhigh, assuming t > sqrt(Btra).
868 : * d = dimension of low part (= [nf:Q])
869 : * n0 = bound for |vS|^2
870 : * */
871 : static double
872 12166 : get_Bhigh(long n0, long d)
873 : {
874 12166 : double sqrtd = sqrt((double)d);
875 12166 : double z = n0*sqrtd + sqrtd/2 * (d * (n0+1));
876 12166 : z = 1. + 0.5 * z; return z * z;
877 : }
878 :
879 : typedef struct {
880 : GEN d;
881 : GEN dPinvS; /* d P^(-1) S [ integral ] */
882 : double **PinvSdbl; /* P^(-1) S as double */
883 : GEN S1, P1; /* S = S0 + S1 q, idem P */
884 : } trace_data;
885 :
886 : /* S1 * u - P1 * round(P^-1 S u). K nonzero coords in u given by ind */
887 : static GEN
888 206133 : get_trace(GEN ind, trace_data *T)
889 : {
890 206133 : long i, j, l, K = lg(ind)-1;
891 : GEN z, s, v;
892 :
893 206133 : s = gel(T->S1, ind[1]);
894 206133 : if (K == 1) return s;
895 :
896 : /* compute s = S1 u */
897 503043 : for (j=2; j<=K; j++) s = ZC_add(s, gel(T->S1, ind[j]));
898 :
899 : /* compute v := - round(P^1 S u) */
900 179124 : l = lg(s);
901 179124 : v = cgetg(l, t_VECSMALL);
902 2621083 : for (i=1; i<l; i++)
903 : {
904 2441959 : double r, t = 0.;
905 : /* quick approximate computation */
906 9301921 : for (j=1; j<=K; j++) t += T->PinvSdbl[ ind[j] ][i];
907 2441959 : r = floor(t + 0.5);
908 2441959 : if (fabs(t + 0.5 - r) < 0.0001)
909 : { /* dubious, compute exactly */
910 362 : z = gen_0;
911 1391 : for (j=1; j<=K; j++) z = addii(z, ((GEN**)T->dPinvS)[ ind[j] ][i]);
912 362 : v[i] = - itos( diviiround(z, T->d) );
913 : }
914 : else
915 2441597 : v[i] = - (long)r;
916 : }
917 179124 : return ZC_add(s, ZM_zc_mul(T->P1, v));
918 : }
919 :
920 : static trace_data *
921 24150 : init_trace(trace_data *T, GEN S, nflift_t *L, GEN q)
922 : {
923 24150 : long e = gexpo(S), i,j, l,h;
924 : GEN qgood, S1, invd;
925 :
926 24150 : if (e < 0) return NULL; /* S = 0 */
927 :
928 20383 : qgood = int2n(e - 32); /* single precision check */
929 20383 : if (cmpii(qgood, q) > 0) q = qgood;
930 :
931 20383 : S1 = gdivround(S, q);
932 20383 : if (gequal0(S1)) return NULL;
933 :
934 4052 : invd = invr(itor(L->pk, DEFAULTPREC));
935 :
936 4052 : T->dPinvS = ZM_mul(L->iprk, S);
937 4052 : l = lg(S);
938 4052 : h = lgcols(T->dPinvS);
939 4052 : T->PinvSdbl = (double**)cgetg(l, t_MAT);
940 34723 : for (j = 1; j < l; j++)
941 : {
942 30671 : double *t = (double *) stack_malloc_align(h * sizeof(double), sizeof(double));
943 30671 : GEN c = gel(T->dPinvS,j);
944 30671 : pari_sp av = avma;
945 30671 : T->PinvSdbl[j] = t;
946 322530 : for (i=1; i < h; i++) t[i] = rtodbl(mulri(invd, gel(c,i)));
947 30671 : set_avma(av);
948 : }
949 :
950 4052 : T->d = L->pk;
951 4052 : T->P1 = gdivround(L->prk, q);
952 4052 : T->S1 = S1; return T;
953 : }
954 :
955 : static void
956 208864 : update_trace(trace_data *T, long k, long i)
957 : {
958 208864 : if (!T) return;
959 92230 : gel(T->S1,k) = gel(T->S1,i);
960 92230 : gel(T->dPinvS,k) = gel(T->dPinvS,i);
961 92230 : T->PinvSdbl[k] = T->PinvSdbl[i];
962 : }
963 :
964 : /* reduce coeffs mod (T,pk), then center mod pk */
965 : static GEN
966 52955 : FqX_centermod(GEN z, GEN T, GEN pk, GEN pks2)
967 : {
968 : long i, l;
969 : GEN y;
970 52955 : if (!T) return centermod_i(z, pk, pks2);
971 15435 : y = FpXQX_red(z, T, pk); l = lg(y);
972 138537 : for (i = 2; i < l; i++)
973 : {
974 123102 : GEN c = gel(y,i);
975 123102 : if (typ(c) == t_INT)
976 81417 : c = Fp_center_i(c, pk, pks2);
977 : else
978 41685 : c = FpX_center_i(c, pk, pks2);
979 123102 : gel(y,i) = c;
980 : }
981 15435 : return y;
982 : }
983 :
984 : typedef struct {
985 : GEN lt, C, Clt, C2lt, C2ltpol;
986 : } div_data;
987 :
988 : static void
989 19092 : init_div_data(div_data *D, GEN pol, nflift_t *L)
990 : {
991 19092 : GEN C2lt, Clt, C = mul_content(L->topowden, L->dn);
992 19092 : GEN lc = leading_coeff(pol), lt = is_pm1(lc)? NULL: absi_shallow(lc);
993 19092 : if (C)
994 : {
995 19092 : GEN C2 = sqri(C);
996 19092 : if (lt) {
997 3588 : C2lt = mulii(C2, lt);
998 3588 : Clt = mulii(C,lt);
999 : } else {
1000 15504 : C2lt = C2;
1001 15504 : Clt = C;
1002 : }
1003 : }
1004 : else
1005 0 : C2lt = Clt = lt;
1006 19092 : D->lt = lt;
1007 19092 : D->C = C;
1008 19092 : D->Clt = Clt;
1009 19092 : D->C2lt = C2lt;
1010 19092 : D->C2ltpol = C2lt? RgX_Rg_mul(pol, C2lt): pol;
1011 19092 : }
1012 : static void
1013 24409 : update_target(div_data *D, GEN pol)
1014 24409 : { D->C2ltpol = D->Clt? RgX_Rg_mul(pol, D->Clt): pol; }
1015 :
1016 : /* nb = number of modular factors; return a "good" K such that naive
1017 : * recombination of up to maxK modular factors is not too costly */
1018 : long
1019 114772 : cmbf_maxK(long nb)
1020 : {
1021 114772 : if (nb > 10) return 3;
1022 110712 : return nb-1;
1023 : }
1024 : /* Naive recombination of modular factors: combine up to maxK modular
1025 : * factors, degree <= klim
1026 : *
1027 : * target = polynomial we want to factor
1028 : * famod = array of modular factors. Product should be congruent to
1029 : * target/lc(target) modulo p^a
1030 : * For true factors: S1,S2 <= p^b, with b <= a and p^(b-a) < 2^31 */
1031 : /* set *done = 1 if factorisation is known to be complete */
1032 : static GEN
1033 12075 : nfcmbf(nfcmbf_t *T, long klim, long *pmaxK, int *done)
1034 : {
1035 12075 : GEN nf = T->nf, famod = T->fact, bound = T->bound;
1036 12075 : GEN ltdn, nfpol = nf_get_pol(nf);
1037 12075 : long K = 1, cnt = 1, i,j,k, curdeg, lfamod = lg(famod)-1, dnf = degpol(nfpol);
1038 12075 : pari_sp av0 = avma;
1039 12075 : GEN Tpk = T->L->Tpk, pk = T->L->pk, pks2 = shifti(pk,-1);
1040 12075 : GEN ind = cgetg(lfamod+1, t_VECSMALL);
1041 12075 : GEN deg = cgetg(lfamod+1, t_VECSMALL);
1042 12075 : GEN degsofar = cgetg(lfamod+1, t_VECSMALL);
1043 12075 : GEN fa = cgetg(lfamod+1, t_VEC);
1044 12075 : const double Bhigh = get_Bhigh(lfamod, dnf);
1045 : trace_data _T1, _T2, *T1, *T2;
1046 : div_data D;
1047 : pari_timer ti;
1048 :
1049 12075 : timer_start(&ti);
1050 :
1051 12075 : *pmaxK = cmbf_maxK(lfamod);
1052 12075 : init_div_data(&D, T->pol, T->L);
1053 12075 : ltdn = mul_content(D.lt, T->L->dn);
1054 : {
1055 12075 : GEN q = ceil_safe(sqrtr(T->BS_2));
1056 12075 : GEN t1,t2, lt2dn = mul_content(ltdn, D.lt);
1057 12075 : GEN trace1 = cgetg(lfamod+1, t_MAT);
1058 12075 : GEN trace2 = cgetg(lfamod+1, t_MAT);
1059 61270 : for (i=1; i <= lfamod; i++)
1060 : {
1061 49195 : pari_sp av = avma;
1062 49195 : GEN P = gel(famod,i);
1063 49195 : long d = degpol(P);
1064 :
1065 49195 : deg[i] = d; P += 2;
1066 49195 : t1 = gel(P,d-1);/* = - S_1 */
1067 49195 : t2 = Fq_sqr(t1, Tpk, pk);
1068 49195 : if (d > 1) t2 = Fq_sub(t2, gmul2n(gel(P,d-2), 1), Tpk, pk);
1069 : /* t2 = S_2 Newton sum */
1070 49195 : if (ltdn)
1071 : {
1072 518 : t1 = Fq_Fp_mul(t1, ltdn, Tpk, pk);
1073 518 : t2 = Fq_Fp_mul(t2, lt2dn, Tpk, pk);
1074 : }
1075 49195 : gel(trace1,i) = gclone( nf_bestlift(t1, NULL, T->L) );
1076 49195 : gel(trace2,i) = gclone( nf_bestlift(t2, NULL, T->L) ); set_avma(av);
1077 : }
1078 12075 : T1 = init_trace(&_T1, trace1, T->L, q);
1079 12075 : T2 = init_trace(&_T2, trace2, T->L, q);
1080 61270 : for (i=1; i <= lfamod; i++) {
1081 49195 : gunclone(gel(trace1,i));
1082 49195 : gunclone(gel(trace2,i));
1083 : }
1084 : }
1085 12075 : degsofar[0] = 0; /* sentinel */
1086 :
1087 : /* ind runs through strictly increasing sequences of length K,
1088 : * 1 <= ind[i] <= lfamod */
1089 16344 : nextK:
1090 16344 : if (K > *pmaxK || 2*K > lfamod) goto END;
1091 14021 : if (DEBUGLEVEL > 3)
1092 0 : err_printf("\n### K = %d, %Ps combinations\n", K,binomial(utoipos(lfamod), K));
1093 14021 : setlg(ind, K+1); ind[1] = 1;
1094 14021 : i = 1; curdeg = deg[ind[1]];
1095 : for(;;)
1096 : { /* try all combinations of K factors */
1097 253625 : for (j = i; j < K; j++)
1098 : {
1099 29715 : degsofar[j] = curdeg;
1100 29715 : ind[j+1] = ind[j]+1; curdeg += deg[ind[j+1]];
1101 : }
1102 223910 : if (curdeg <= klim) /* trial divide */
1103 14545 : {
1104 : GEN t, y, q;
1105 : pari_sp av;
1106 :
1107 223910 : av = avma;
1108 223910 : if (T1)
1109 : { /* d-1 test */
1110 100723 : t = get_trace(ind, T1);
1111 100723 : if (rtodbl(_norml2(t)) > Bhigh)
1112 : {
1113 92348 : if (DEBUGLEVEL>6) err_printf(".");
1114 92348 : set_avma(av); goto NEXT;
1115 : }
1116 : }
1117 131562 : if (T2)
1118 : { /* d-2 test */
1119 105410 : t = get_trace(ind, T2);
1120 105410 : if (rtodbl(_norml2(t)) > Bhigh)
1121 : {
1122 95162 : if (DEBUGLEVEL>3) err_printf("|");
1123 95162 : set_avma(av); goto NEXT;
1124 : }
1125 : }
1126 36400 : set_avma(av);
1127 36400 : y = ltdn; /* full computation */
1128 89355 : for (i=1; i<=K; i++)
1129 : {
1130 52955 : GEN q = gel(famod, ind[i]);
1131 52955 : if (y) q = gmul(y, q);
1132 52955 : y = FqX_centermod(q, Tpk, pk, pks2);
1133 : }
1134 36400 : y = nf_pol_lift(y, bound, T->L);
1135 36400 : if (!y)
1136 : {
1137 9515 : if (DEBUGLEVEL>3) err_printf("@");
1138 9515 : set_avma(av); goto NEXT;
1139 : }
1140 : /* y = topowden*dn*lt*\prod_{i in ind} famod[i] is apparently in O_K[X],
1141 : * in fact in (Z[Y]/nf.pol)[X] due to multiplication by C = topowden*dn.
1142 : * Try out this candidate factor */
1143 26885 : q = RgXQX_divrem(D.C2ltpol, y, nfpol, ONLY_DIVIDES);
1144 26885 : if (!q)
1145 : {
1146 2588 : if (DEBUGLEVEL>3) err_printf("*");
1147 2588 : set_avma(av); goto NEXT;
1148 : }
1149 : /* Original T->pol in O_K[X] with leading coeff lt in Z,
1150 : * y = C*lt \prod famod[i] is in O_K[X] with leading coeff in Z
1151 : * q = C^2*lt*pol / y = C * (lt*pol) / (lt*\prod famod[i]) is a
1152 : * K-rational factor, in fact in Z[Y]/nf.pol)[X] as above, with
1153 : * leading term C*lt. */
1154 24297 : update_target(&D, q);
1155 24297 : gel(fa,cnt++) = D.C2lt? RgX_int_normalize(y): y; /* make monic */
1156 160124 : for (i=j=k=1; i <= lfamod; i++)
1157 : { /* remove used factors */
1158 135827 : if (j <= K && i == ind[j]) j++;
1159 : else
1160 : {
1161 104432 : gel(famod,k) = gel(famod,i);
1162 104432 : update_trace(T1, k, i);
1163 104432 : update_trace(T2, k, i);
1164 104432 : deg[k] = deg[i]; k++;
1165 : }
1166 : }
1167 24297 : lfamod -= K;
1168 24297 : *pmaxK = cmbf_maxK(lfamod);
1169 24297 : if (lfamod < 2*K) goto END;
1170 14545 : i = 1; curdeg = deg[ind[1]];
1171 14545 : if (DEBUGLEVEL > 2)
1172 : {
1173 0 : err_printf("\n"); timer_printf(&ti, "to find factor %Ps",gel(fa,cnt-1));
1174 0 : err_printf("remaining modular factor(s): %ld\n", lfamod);
1175 : }
1176 14545 : continue;
1177 : }
1178 :
1179 0 : NEXT:
1180 199613 : for (i = K+1;;)
1181 : {
1182 226499 : if (--i == 0) { K++; goto nextK; }
1183 222230 : if (++ind[i] <= lfamod - K + i)
1184 : {
1185 195344 : curdeg = degsofar[i-1] + deg[ind[i]];
1186 195344 : if (curdeg <= klim) break;
1187 : }
1188 : }
1189 : }
1190 12075 : END:
1191 12075 : *done = 1;
1192 12075 : if (degpol(D.C2ltpol) > 0)
1193 : { /* leftover factor */
1194 12075 : GEN q = D.C2ltpol;
1195 12075 : if (D.C2lt) q = RgX_int_normalize(q);
1196 12075 : if (lfamod >= 2*K)
1197 : { /* restore leading coefficient [#930] */
1198 91 : if (D.lt) q = RgX_Rg_mul(q, D.lt);
1199 91 : *done = 0; /* ... may still be reducible */
1200 : }
1201 12075 : setlg(famod, lfamod+1);
1202 12075 : gel(fa,cnt++) = q;
1203 : }
1204 12075 : if (DEBUGLEVEL>6) err_printf("\n");
1205 12075 : setlg(fa, cnt);
1206 12075 : return gc_GEN(av0, fa);
1207 : }
1208 :
1209 : static GEN
1210 98 : nf_chk_factors(nfcmbf_t *T, GEN P, GEN M_L, GEN famod, GEN pk)
1211 : {
1212 98 : GEN nf = T->nf, bound = T->bound;
1213 98 : GEN nfT = nf_get_pol(nf);
1214 : long i, r;
1215 98 : GEN pol = P, list, piv, y;
1216 98 : GEN Tpk = T->L->Tpk;
1217 : div_data D;
1218 :
1219 98 : piv = ZM_hnf_knapsack(M_L);
1220 98 : if (!piv) return NULL;
1221 77 : if (DEBUGLEVEL>3) err_printf("ZM_hnf_knapsack output:\n%Ps\n",piv);
1222 :
1223 77 : r = lg(piv)-1;
1224 77 : list = cgetg(r+1, t_VEC);
1225 77 : init_div_data(&D, pol, T->L);
1226 77 : for (i = 1;;)
1227 112 : {
1228 189 : pari_sp av = avma;
1229 189 : if (DEBUGLEVEL) err_printf("nf_LLL_cmbf: checking factor %ld\n", i);
1230 189 : y = chk_factors_get(D.lt, famod, gel(piv,i), Tpk, pk);
1231 :
1232 189 : if (! (y = nf_pol_lift(y, bound, T->L)) ) return NULL;
1233 182 : y = gc_GEN(av, y);
1234 : /* y is the candidate factor */
1235 182 : pol = RgXQX_divrem(D.C2ltpol, y, nfT, ONLY_DIVIDES);
1236 182 : if (!pol) return NULL;
1237 :
1238 182 : if (D.C2lt) y = RgX_int_normalize(y);
1239 182 : gel(list,i) = y;
1240 182 : if (++i >= r) break;
1241 :
1242 112 : update_target(&D, pol);
1243 : }
1244 70 : gel(list,i) = RgX_int_normalize(pol); return list;
1245 : }
1246 :
1247 : static GEN
1248 165205 : nf_to_Zq(GEN x, GEN T, GEN pk, GEN pks2, GEN proj)
1249 : {
1250 : GEN y;
1251 165205 : if (typ(x) != t_COL) return centermodii(x, pk, pks2);
1252 23224 : if (!T)
1253 : {
1254 22951 : y = ZV_dotproduct(proj, x);
1255 22951 : return centermodii(y, pk, pks2);
1256 : }
1257 273 : y = ZM_ZC_mul(proj, x);
1258 273 : y = RgV_to_RgX(y, varn(T));
1259 273 : return FpX_center_i(FpX_rem(y, T, pk), pk, pks2);
1260 : }
1261 :
1262 : /* Assume P in nfX form, lc(P) != 0 mod p. Reduce P to Zp[X]/(T) mod p^a */
1263 : static GEN
1264 19057 : ZqX(GEN P, GEN pk, GEN T, GEN proj)
1265 : {
1266 19057 : long i, l = lg(P);
1267 19057 : GEN z, pks2 = shifti(pk,-1);
1268 :
1269 19057 : z = cgetg(l,t_POL); z[1] = P[1];
1270 184227 : for (i=2; i<l; i++) gel(z,i) = nf_to_Zq(gel(P,i),T,pk,pks2,proj);
1271 19057 : return normalizepol_lg(z, l);
1272 : }
1273 :
1274 : static GEN
1275 19022 : ZqX_normalize(GEN P, GEN lt, nflift_t *L)
1276 : {
1277 19022 : GEN R = lt? RgX_Rg_mul(P, Fp_inv(lt, L->pk)): P;
1278 19022 : return ZqX(R, L->pk, L->Tpk, L->ZqProj);
1279 : }
1280 :
1281 : /* k allowing to reconstruct x, T_2(x)^2 < C, from x mod pr^k */
1282 : /* return log [ 2sqrt(C/d) * ( (3/2)sqrt(gamma) )^(d-1) ] ^d / log N(pr)
1283 : * cf. Belabas relative van Hoeij algorithm, lemma 3.12 */
1284 : static double
1285 19057 : bestlift_bound(GEN C, long d, double alpha, GEN p, long f)
1286 : {
1287 19057 : const double g = 1 / (alpha - 0.25); /* = 2 if alpha = 3/4 */
1288 19057 : GEN C4 = shiftr(gtofp(C,DEFAULTPREC), 2);
1289 19057 : double t, logp = log(gtodouble(p));
1290 19057 : if (f == d)
1291 : { /* p inert, no LLL fudge factor: p^(2k) / 4 > C */
1292 35 : t = 0.5 * rtodbl(mplog(C4));
1293 35 : return ceil(t / logp);
1294 : }
1295 : /* (1/2)log (4C/d) + (d-1)(log 3/2 sqrt(gamma)) */
1296 19022 : t = 0.5 * rtodbl(mplog(divru(C4,d))) + (d-1) * log(1.5 * sqrt(g));
1297 19022 : return ceil((t * d) / (logp * f));
1298 : }
1299 :
1300 : /* M invertible ZM (would also work with a mix of t_INT/t_REAL entries).
1301 : * Return upper triangular matrix R from QR decomposition of
1302 : * M and set B = squared lengths of orthogonalized vectors.
1303 : * R and B have t_INT/t_REAL coefs */
1304 : static GEN
1305 19275 : get_R(GEN M, GEN *pB)
1306 : {
1307 : GEN B, R;
1308 19275 : long i, l, prec = nbits2prec( gexpo(M) + 64 );
1309 :
1310 : for(;;)
1311 : {
1312 19275 : R = gaussred_from_QR(M, prec);
1313 19275 : if (R) break;
1314 0 : prec = precdbl(prec);
1315 : }
1316 19275 : l = lg(R); B = cgetg(l, t_VEC);
1317 92121 : for (i=1; i<l; i++) { gel(B,i) = gcoeff(R,i,i); gcoeff(R,i,i) = gen_1; }
1318 19275 : *pB = B; return R;
1319 : }
1320 :
1321 : static void
1322 19057 : init_proj(nflift_t *L, GEN prkHNF, GEN nfT)
1323 : {
1324 19057 : if (degpol(L->Tp)>1)
1325 : {
1326 210 : GEN coTp = FpX_div(FpX_red(nfT, L->p), L->Tp, L->p); /* Tp's cofactor */
1327 : GEN z, proj;
1328 210 : z = ZpX_liftfact(nfT, mkvec2(L->Tp, coTp), L->pk, L->p, L->k);
1329 210 : L->Tpk = gel(z,1);
1330 210 : proj = QXQV_to_FpM(L->topow, L->Tpk, L->pk);
1331 210 : if (L->topowden)
1332 210 : proj = FpM_red(ZM_Z_mul(proj, Fp_inv(L->topowden, L->pk)), L->pk);
1333 210 : L->ZqProj = proj;
1334 : }
1335 : else
1336 : {
1337 18847 : L->Tpk = NULL;
1338 18847 : L->ZqProj = hnf_Znproj(prkHNF);
1339 : }
1340 19057 : }
1341 :
1342 : /* Square of the radius of largest ball inscript in PRK's fundamental domain,
1343 : * Rmax ^2 = min 1/4T_i where T_i = sum_{i <= j} ( s_ij^2 / B_j)
1344 : * where the B_i are the orthogonalized vectors squared norms.
1345 : * For p inert, S = Id, T_i = 1 / p^{2k} and Rmax = p^k / 2 */
1346 : static GEN
1347 19275 : max_radius(GEN PRK)
1348 : {
1349 19275 : GEN B, S, smax = gen_0;
1350 19275 : pari_sp av = avma;
1351 19275 : long i, j, d = lg(PRK)-1; /* number field degree > 1 */
1352 :
1353 19275 : S = RgM_inv_upper( get_R(PRK, &B) ); if (!S) pari_err_PREC("max_radius");
1354 92121 : for (i=1; i<=d; i++) /* S[i,i] = gen_1 */
1355 : {
1356 72846 : GEN s = ginv(gel(B,i)); /* t_INT, t_FRAC or t_REAL */
1357 274720 : for (j = i+1; j <= d; j++) /* s += a t_INT, t_FRAC or t_REAL */
1358 201874 : s = gadd(s, gdiv(mpsqr(gcoeff(S,i,j)), gel(B,j)));
1359 72846 : if (gcmp(s, smax) > 0) smax = s;
1360 : }
1361 19275 : return gc_upto(av, ginv(gmul2n(smax, 2)));
1362 : }
1363 :
1364 : static void
1365 19057 : bestlift_init(long a, GEN nf, GEN C, nflift_t *L)
1366 : {
1367 19057 : const double alpha = 0.99; /* LLL parameter */
1368 19057 : const long d = nf_get_degree(nf);
1369 19057 : pari_sp av = avma, av2;
1370 19057 : GEN prk, PRK, iPRK, GSmin, T = L->Tp, p = L->p;
1371 19057 : long f = degpol(T);
1372 : pari_timer ti;
1373 :
1374 19057 : if (f == d)
1375 : { /* inert p, much simpler */
1376 35 : long a0 = bestlift_bound(C, d, alpha, p, f);
1377 : GEN q;
1378 35 : if (a < a0) a = a0; /* guarantees GSmin >= C */
1379 35 : if (DEBUGLEVEL>2) err_printf("exponent %ld\n",a);
1380 35 : q = powiu(p,a);
1381 35 : PRK = prk = scalarmat_shallow(q, d);
1382 35 : GSmin = shiftr(itor(q, DEFAULTPREC), -1);
1383 35 : iPRK = matid(d); goto END;
1384 : }
1385 19022 : timer_start(&ti);
1386 19022 : if (!a) a = (long)bestlift_bound(C, d, alpha, p, f);
1387 253 : for (;; set_avma(av), a += (a==1)? 1: (a>>1)) /* roughly a *= 1.5 */
1388 253 : {
1389 19275 : GEN q = powiu(p,a), Tq = FpXQ_powu(T, a, FpX_red(nf_get_pol(nf), q), q);
1390 19275 : if (DEBUGLEVEL>2) err_printf("exponent %ld\n",a);
1391 19275 : prk = idealhnf_two(nf, mkvec2(q, Tq));
1392 19275 : av2 = avma;
1393 19275 : PRK = ZM_lll(prk, alpha, LLL_INPLACE);
1394 19275 : GSmin = max_radius(PRK);
1395 19275 : if (gcmp(GSmin, C) >= 0) break;
1396 : }
1397 19022 : (void)gc_all(av2, 2, &PRK, &GSmin);
1398 19022 : iPRK = ZM_inv(PRK, NULL);
1399 19022 : if (DEBUGLEVEL>2)
1400 0 : err_printf("for this exponent, GSmin = %Ps\nTime reduction: %ld\n",
1401 : GSmin, timer_delay(&ti));
1402 19022 : END:
1403 19057 : L->k = a;
1404 19057 : L->pk = gcoeff(prk,1,1);
1405 19057 : L->prk = PRK;
1406 19057 : L->iprk = iPRK;
1407 19057 : L->GSmin= GSmin;
1408 19057 : init_proj(L, prk, nf_get_pol(nf));
1409 19057 : }
1410 :
1411 : /* Let X = Tra * M_L, Y = bestlift(X) return V s.t Y = X - PRK V
1412 : * and set *eT2 = gexpo(Y) [cf nf_bestlift, but memory efficient] */
1413 : static GEN
1414 476 : get_V(GEN Tra, GEN M_L, GEN PRK, GEN PRKinv, GEN pk, long *eT2)
1415 : {
1416 476 : long i, e = 0, l = lg(M_L);
1417 476 : GEN V = cgetg(l, t_MAT);
1418 476 : *eT2 = 0;
1419 6013 : for (i = 1; i < l; i++)
1420 : { /* cf nf_bestlift(Tra * c) */
1421 5537 : pari_sp av = avma, av2;
1422 5537 : GEN v, T2 = ZM_ZC_mul(Tra, gel(M_L,i));
1423 :
1424 5537 : v = gdivround(ZM_ZC_mul(PRKinv, T2), pk); /* small */
1425 5537 : av2 = avma;
1426 5537 : T2 = ZC_sub(T2, ZM_ZC_mul(PRK, v));
1427 5537 : e = gexpo(T2); if (e > *eT2) *eT2 = e;
1428 5537 : set_avma(av2);
1429 5537 : gel(V,i) = gc_upto(av, v); /* small */
1430 : }
1431 476 : return V;
1432 : }
1433 :
1434 : static GEN
1435 91 : nf_LLL_cmbf(nfcmbf_t *T, long rec)
1436 : {
1437 91 : const double BitPerFactor = 0.4; /* nb bits / modular factor */
1438 91 : nflift_t *L = T->L;
1439 91 : GEN famod = T->fact, ZC = T->ZC, Br = T->Br, P = T->pol, dn = T->L->dn;
1440 91 : long dnf = nf_get_degree(T->nf), dP = degpol(P);
1441 : long i, C, tmax, n0;
1442 : GEN lP, Bnorm, Tra, T2, TT, CM_L, m, list, ZERO, Btra;
1443 : double Bhigh;
1444 : pari_sp av, av2;
1445 91 : long ti_LLL = 0, ti_CF = 0;
1446 : pari_timer ti2, TI;
1447 :
1448 91 : lP = absi_shallow(leading_coeff(P));
1449 91 : if (is_pm1(lP)) lP = NULL;
1450 :
1451 91 : n0 = lg(famod) - 1;
1452 : /* Lattice: (S PRK), small vector (vS vP). To find k bound for the image,
1453 : * write S = S1 q + S0, P = P1 q + P0
1454 : * |S1 vS + P1 vP|^2 <= Bhigh for all (vS,vP) assoc. to true factors */
1455 91 : Btra = mulrr(ZC, mulur(dP*dP, normTp(Br, 2, dnf)));
1456 91 : Bhigh = get_Bhigh(n0, dnf);
1457 91 : C = (long)ceil(sqrt(Bhigh/n0)) + 1; /* C^2 n0 ~ Bhigh */
1458 91 : Bnorm = dbltor( n0 * C * C + Bhigh );
1459 91 : ZERO = zeromat(n0, dnf);
1460 :
1461 91 : av = avma;
1462 91 : TT = const_vec(n0, NULL);
1463 91 : Tra = cgetg(n0+1, t_MAT);
1464 91 : CM_L = scalarmat_s(C, n0);
1465 : /* tmax = current number of traces used (and computed so far) */
1466 91 : for(tmax = 0;; tmax++)
1467 231 : {
1468 322 : long a, b, bmin, bgood, delta, tnew = tmax + 1, r = lg(CM_L)-1;
1469 : GEN M_L, q, CM_Lp, oldCM_L, S1, P1, VV;
1470 322 : int first = 1;
1471 :
1472 : /* bound for f . S_k(genuine factor) = ZC * bound for T_2(S_tnew) */
1473 322 : Btra = mulrr(ZC, mulur(dP*dP, normTp(Br, 2*tnew, dnf)));
1474 322 : bmin = logint(ceil_safe(sqrtr(Btra)), gen_2) + 1;
1475 322 : if (DEBUGLEVEL>2)
1476 0 : err_printf("\nLLL_cmbf: %ld potential factors (tmax = %ld, bmin = %ld)\n",
1477 : r, tmax, bmin);
1478 :
1479 : /* compute Newton sums (possibly relifting first) */
1480 322 : if (gcmp(L->GSmin, Btra) < 0)
1481 : {
1482 : GEN polred;
1483 :
1484 7 : bestlift_init((L->k)<<1, T->nf, Btra, L);
1485 7 : polred = ZqX_normalize(T->polbase, lP, L);
1486 7 : famod = ZqX_liftfact(polred, famod, L->Tpk, L->pk, L->p, L->k);
1487 133 : for (i=1; i<=n0; i++) gel(TT,i) = NULL;
1488 : }
1489 6776 : for (i=1; i<=n0; i++)
1490 : {
1491 6454 : GEN h, lPpow = lP? powiu(lP, tnew): NULL;
1492 6454 : GEN z = polsym_gen(gel(famod,i), gel(TT,i), tnew, L->Tpk, L->pk);
1493 6454 : gel(TT,i) = z;
1494 6454 : h = gel(z,tnew+1);
1495 : /* make Newton sums integral */
1496 6454 : lPpow = mul_content(lPpow, dn);
1497 6454 : if (lPpow)
1498 126 : h = (typ(h) == t_INT)? Fp_mul(h, lPpow, L->pk): FpX_Fp_mul(h, lPpow, L->pk);
1499 6454 : gel(Tra,i) = nf_bestlift(h, NULL, L); /* S_tnew(famod) */
1500 : }
1501 :
1502 : /* compute truncation parameter */
1503 322 : if (DEBUGLEVEL>2) { timer_start(&ti2); timer_start(&TI); }
1504 322 : oldCM_L = CM_L;
1505 322 : av2 = avma;
1506 322 : b = delta = 0; /* -Wall */
1507 476 : AGAIN:
1508 476 : M_L = Q_div_to_int(CM_L, utoipos(C));
1509 476 : VV = get_V(Tra, M_L, L->prk, L->iprk, L->pk, &a);
1510 476 : if (first)
1511 : { /* initialize lattice, using few p-adic digits for traces */
1512 322 : bgood = (long)(a - maxss(32, (long)(BitPerFactor * r)));
1513 322 : b = maxss(bmin, bgood);
1514 322 : delta = a - b;
1515 : }
1516 : else
1517 : { /* add more p-adic digits and continue reduction */
1518 154 : if (a < b) b = a;
1519 154 : b = maxss(b-delta, bmin);
1520 154 : if (b - delta/2 < bmin) b = bmin; /* near there. Go all the way */
1521 : }
1522 :
1523 : /* restart with truncated entries */
1524 476 : q = int2n(b);
1525 476 : P1 = gdivround(L->prk, q);
1526 476 : S1 = gdivround(Tra, q);
1527 476 : T2 = ZM_sub(ZM_mul(S1, M_L), ZM_mul(P1, VV));
1528 476 : m = vconcat( CM_L, T2 );
1529 476 : if (first)
1530 : {
1531 322 : first = 0;
1532 322 : m = shallowconcat( m, vconcat(ZERO, P1) );
1533 : /* [ C M_L 0 ]
1534 : * m = [ ] square matrix
1535 : * [ T2' PRK ] T2' = Tra * M_L truncated
1536 : */
1537 : }
1538 476 : CM_L = LLL_check_progress(Bnorm, n0, m, b == bmin, /*dbg:*/ &ti_LLL);
1539 476 : if (DEBUGLEVEL>2)
1540 0 : err_printf("LLL_cmbf: (a,b) =%4ld,%4ld; r =%3ld -->%3ld, time = %ld\n",
1541 0 : a,b, lg(m)-1, CM_L? lg(CM_L)-1: 1, timer_delay(&TI));
1542 476 : if (!CM_L) { list = mkcol(RgX_int_normalize(P)); break; }
1543 455 : if (b > bmin)
1544 : {
1545 154 : CM_L = gc_GEN(av2, CM_L);
1546 154 : goto AGAIN;
1547 : }
1548 301 : if (DEBUGLEVEL>2) timer_printf(&ti2, "for this trace");
1549 :
1550 301 : i = lg(CM_L) - 1;
1551 301 : if (i == r && ZM_equal(CM_L, oldCM_L))
1552 : {
1553 91 : CM_L = oldCM_L;
1554 91 : set_avma(av2); continue;
1555 : }
1556 :
1557 210 : CM_Lp = FpM_image(CM_L, utoipos(27449)); /* inexpensive test */
1558 210 : if (lg(CM_Lp) != lg(CM_L))
1559 : {
1560 0 : if (DEBUGLEVEL>2) err_printf("LLL_cmbf: rank decrease\n");
1561 0 : CM_L = ZM_hnf(CM_L);
1562 : }
1563 :
1564 210 : if (i <= r && i*rec < n0)
1565 : {
1566 : pari_timer ti;
1567 98 : if (DEBUGLEVEL>2) timer_start(&ti);
1568 98 : list = nf_chk_factors(T, P, Q_div_to_int(CM_L,utoipos(C)), famod, L->pk);
1569 98 : if (DEBUGLEVEL>2) ti_CF += timer_delay(&ti);
1570 98 : if (list) break;
1571 : }
1572 140 : if (gc_needed(av,1))
1573 : {
1574 0 : if(DEBUGMEM>1) pari_warn(warnmem,"nf_LLL_cmbf");
1575 0 : (void)gc_all(av, L->Tpk? 9: 8,
1576 : &CM_L,&TT,&Tra,&famod,&L->GSmin,&L->pk,&L->prk,&L->iprk,
1577 : &L->Tpk);
1578 : }
1579 140 : else CM_L = gc_GEN(av2, CM_L);
1580 : }
1581 91 : if (DEBUGLEVEL>2)
1582 0 : err_printf("* Time LLL: %ld\n* Time Check Factor: %ld\n",ti_LLL,ti_CF);
1583 91 : return list;
1584 : }
1585 :
1586 : static GEN
1587 12075 : nf_combine_factors(nfcmbf_t *T, GEN polred, long klim)
1588 : {
1589 12075 : nflift_t *L = T->L;
1590 : GEN res;
1591 : long maxK;
1592 : int done;
1593 : pari_timer ti;
1594 :
1595 12075 : if (DEBUGLEVEL>2) timer_start(&ti);
1596 12075 : T->fact = ZqX_liftfact(polred, T->fact, L->Tpk, L->pk, L->p, L->k);
1597 12075 : if (DEBUGLEVEL>2) timer_printf(&ti, "Hensel lift");
1598 12075 : res = nfcmbf(T, klim, &maxK, &done);
1599 12075 : if (DEBUGLEVEL>2) timer_printf(&ti, "Naive recombination");
1600 12075 : if (!done)
1601 : {
1602 91 : long l = lg(res)-1;
1603 : GEN v;
1604 91 : if (l > 1)
1605 : {
1606 : GEN den;
1607 49 : T->pol = gel(res,l);
1608 49 : T->polbase = Q_remove_denom(RgX_to_nfX(T->nf, T->pol), &den);
1609 49 : if (den) { T->Br = gmul(T->Br, den); T->pol = RgX_Rg_mul(T->pol, den); }
1610 : }
1611 91 : v = nf_LLL_cmbf(T, maxK);
1612 : /* remove last elt, possibly unfactored. Add all new ones. */
1613 91 : setlg(res, l); res = shallowconcat(res, v);
1614 : }
1615 12075 : return res;
1616 : }
1617 :
1618 : GEN
1619 13950 : nfroots_worker(GEN r, GEN C, GEN lt, GEN T, long v)
1620 : {
1621 13950 : pari_sp av = avma;
1622 13950 : retgc_const(av, ZXQX_dvd(C, deg1pol_shallow(lt, gneg(r), v), T) ? gen_1: gen_0);
1623 : }
1624 :
1625 : static GEN
1626 6940 : nf_DDF_roots(GEN pol, GEN polred, GEN nfpol, long fl, nflift_t *L)
1627 : {
1628 : GEN z, ltdn, dvd, worker;
1629 6940 : long i, m, lz, v = varn(pol);
1630 : div_data D;
1631 :
1632 6940 : init_div_data(&D, pol, L);
1633 6940 : ltdn = mul_content(D.lt, L->dn);
1634 6940 : z = ZqX_roots(polred, L->Tpk, L->p, L->k);
1635 6940 : lz = lg(z);
1636 6940 : if (DEBUGLEVEL > 3) err_printf("Checking %ld roots:",lz-1);
1637 20890 : for (i=1; i<lz; i++)
1638 : {
1639 13950 : GEN r = gel(z,i);
1640 13950 : if (DEBUGLEVEL > 3) err_printf(" %ld",i);
1641 : /* lt*dn*topowden * r = Clt * r */
1642 13950 : gel(z,i) = nf_bestlift_to_pol(ltdn? gmul(ltdn,r): r, NULL, L);
1643 : }
1644 6940 : worker = snm_closure(is_entry("_nfroots_worker"),mkvec4(D.C2ltpol,D.Clt,nfpol,stoi(v)));
1645 6940 : dvd = gen_parapply_percent(worker, z, DEBUGLEVEL > 3);
1646 20890 : for (m=1,i=1; i<lz; i++)
1647 : {
1648 13950 : GEN r = gel(z,i);
1649 : /* don't go on with q, usually much larger that C2ltpol */
1650 13950 : if (signe(gel(dvd,i))) {
1651 12841 : if (D.Clt) r = gdiv(r, D.Clt);
1652 12841 : gel(z,m++) = r;
1653 : }
1654 1109 : else if (fl == ROOTS_SPLIT) return cgetg(1, t_VEC);
1655 : }
1656 6940 : if (DEBUGLEVEL > 3) err_printf(" done\n");
1657 6940 : z[0] = evaltyp(t_VEC) | _evallg(m);
1658 6940 : return z;
1659 : }
1660 :
1661 : /* returns a few factors of T in Fp of degree <= maxf, NULL if none exist */
1662 : static GEN
1663 462892 : get_good_factor(GEN T, ulong p, long maxf)
1664 : {
1665 462892 : pari_sp av = avma;
1666 462892 : GEN r, R = gel(Flx_factor(ZX_to_Flx(T,p),p), 1);
1667 462892 : if (maxf == 1)
1668 : { /* degree 1 factors are best */
1669 452812 : r = gel(R,1);
1670 452812 : if (degpol(r) == 1) return mkvec(r);
1671 : }
1672 : else
1673 : { /* otherwise, pick factor of largish degree */
1674 10080 : long i, j, dr, d = 0, l = lg(R);
1675 : GEN v;
1676 10080 : if (l == 2) return mkvec(gel(R,1)); /* inert is fine */
1677 9198 : v = cgetg(l, t_VEC);
1678 77038 : for (i = j = 1; i < l; i++)
1679 : {
1680 69674 : r = gel(R,i); dr = degpol(r);
1681 69674 : if (dr > maxf) break;
1682 67840 : if (dr != d) { gel(v,j++) = r; d = dr; }
1683 : }
1684 9198 : setlg(v,j); if (j > 1) return v;
1685 : }
1686 300608 : return gc_NULL(av); /* failure */
1687 : }
1688 :
1689 : /* n = number of modular factors, f = residue degree; nold/fold current best
1690 : * return 1 if new values are "better" than old ones */
1691 : static int
1692 129420 : record(long nold, long n, long fold, long f)
1693 : {
1694 129420 : if (!nold) return 1; /* uninitialized */
1695 104668 : if (fold == f) return n < nold;
1696 : /* if f increases, allow increasing n a little */
1697 8042 : if (fold < f) return n <= 20 || n < 1.1*nold;
1698 : /* f decreases, only allow if decreasing n a lot */
1699 3935 : return n < 0.7*nold;
1700 : }
1701 : /* Optimization problem: factorization of polynomials over large Fq is slow,
1702 : * BUT bestlift correspondingly faster.
1703 : * Return maximal residue degree to be considered when picking a prime ideal */
1704 : static long
1705 31929 : get_maxf(long nfdeg)
1706 : {
1707 31929 : long maxf = 1;
1708 31929 : if (nfdeg >= 45) maxf =32;
1709 31922 : else if (nfdeg >= 30) maxf =16;
1710 31873 : else if (nfdeg >= 15) maxf = 8;
1711 31929 : return maxf;
1712 : }
1713 : /* number of maximal ideals to test before settling on best prime and number
1714 : * of factors; B = [K:Q]*deg(P) */
1715 : static long
1716 31929 : get_nbprimes(long B)
1717 : {
1718 31929 : if (B <= 128) return 5;
1719 1162 : if (B <= 1024) return 20;
1720 63 : if (B <= 2048) return 65;
1721 7 : return 100;
1722 : }
1723 : /* Select a prime ideal pr over which to factor pol.
1724 : * Return the number of factors (or roots, according to flag fl) mod pr.
1725 : * Set:
1726 : * lt: leading term of polbase (t_INT or NULL [ for 1 ])
1727 : * pr: a suitable maximal ideal
1728 : * Fa: factors found mod pr
1729 : * Tp: polynomial defining Fq/Fp */
1730 : static long
1731 31929 : nf_pick_prime(GEN nf, GEN pol, long fl, GEN *lt, GEN *Tp, ulong *pp)
1732 : {
1733 31929 : GEN nfpol = nf_get_pol(nf), bad = mulii(nf_get_disc(nf), nf_get_index(nf));
1734 31929 : long nfdeg = degpol(nfpol), dpol = degpol(pol), nold = 0, fold = 1;
1735 31929 : long maxf = get_maxf(nfdeg), ct = get_nbprimes(nfdeg * dpol);
1736 : ulong p;
1737 : forprime_t S;
1738 : pari_timer ti_pr;
1739 :
1740 31929 : if (DEBUGLEVEL>3) timer_start(&ti_pr);
1741 : /* if field degree is large, try hard to find an inert prime */
1742 31929 : if (nfdeg > 30) ct = maxss(ct, 2*nfdeg);
1743 31929 : *lt = leading_coeff(pol); /* t_INT */
1744 31929 : if (gequal1(*lt)) *lt = NULL;
1745 31929 : *pp = 0;
1746 31929 : *Tp = NULL;
1747 31929 : (void)u_forprime_init(&S, 2, ULONG_MAX);
1748 : /* select pr such that pol has the smallest number of factors, ct attempts */
1749 520056 : while ((p = u_forprime_next(&S)))
1750 : {
1751 : GEN vT;
1752 520056 : long n, i, l, ok = 0;
1753 520056 : ulong ltp = 0;
1754 :
1755 520056 : if (! umodiu(bad,p)) continue;
1756 467148 : if (*lt) { ltp = umodiu(*lt, p); if (!ltp) continue; }
1757 462892 : vT = get_good_factor(nfpol, p, maxf);
1758 462892 : if (!vT) continue;
1759 162284 : l = lg(vT);
1760 314544 : for (i = 1; i < l; i++)
1761 : {
1762 165188 : pari_sp av2 = avma;
1763 165188 : GEN T = gel(vT,i), red = RgX_to_FlxqX(pol, T, p);
1764 165188 : long f = degpol(T);
1765 165188 : if (f == 1)
1766 : { /* degree 1 */
1767 154340 : red = FlxX_to_Flx(red);
1768 154340 : if (ltp) red = Flx_normalize(red, p);
1769 154340 : if (!Flx_is_squarefree(red, p)) { set_avma(av2); continue; }
1770 131850 : ok = 1;
1771 131850 : n = (fl == FACTORS)? Flx_nbfact(red,p): Flx_nbroots(red,p);
1772 : }
1773 : else
1774 : {
1775 10848 : if (ltp) red = FlxqX_normalize(red, T, p);
1776 10848 : if (!FlxqX_is_squarefree(red, T, p)) { set_avma(av2); continue; }
1777 10498 : ok = 1;
1778 10498 : n = (fl == FACTORS)? FlxqX_nbfact(red,T,p): FlxqX_nbroots(red,T,p);
1779 : }
1780 142348 : if (fl == ROOTS_SPLIT && n < dpol) return n; /* not split */
1781 142327 : if (n <= 1)
1782 : {
1783 24233 : if (fl == FACTORS) return n; /* irreducible */
1784 23792 : if (!n) return 0; /* no root */
1785 : }
1786 129434 : if (DEBUGLEVEL>3)
1787 0 : err_printf("%3ld %s at prime (%ld,x^%ld+...)\n Time: %ld\n",
1788 : n, (fl == FACTORS)? "factors": "roots", p,f, timer_delay(&ti_pr));
1789 :
1790 129434 : if (fl == ROOTS && f==nfdeg) { *Tp = T; *pp = p; return n; }
1791 129420 : if (record(nold, n, fold, f)) { nold = n; fold = f; *Tp = T; *pp = p; }
1792 97122 : else set_avma(av2);
1793 : }
1794 149356 : if (ok && --ct <= 0) break;
1795 : }
1796 19001 : if (!nold) pari_err_OVERFLOW("nf_pick_prime [ran out of primes]");
1797 19001 : return nold;
1798 : }
1799 :
1800 : /* Assume lt(T) is a t_INT and T square free. Return t_VEC of irred. factors */
1801 : static GEN
1802 175 : nfsqff_trager(GEN u, GEN T, GEN dent, long fl)
1803 : {
1804 175 : long k = 0, i, lx;
1805 175 : GEN U, P, x0, mx0, fa, n = ZX_ZXY_rnfequation(T, u, &k);
1806 : int tmonic;
1807 175 : if (DEBUGLEVEL>4) err_printf("nfsqff_trager: choosing k = %ld\n",k);
1808 :
1809 : /* n guaranteed to be squarefree */
1810 175 : fa = ZX_DDF_max(Q_primpart(n),fl==ROOTS || fl==ROOTS_SPLIT ? degpol(T): 0); lx = lg(fa);
1811 175 : if (lx == 2) return mkvec(u);
1812 :
1813 133 : tmonic = is_pm1(leading_coeff(T));
1814 133 : P = cgetg(lx,t_VEC);
1815 133 : x0 = deg1pol_shallow(stoi(-k), gen_0, varn(T));
1816 133 : mx0 = deg1pol_shallow(stoi(k), gen_0, varn(T));
1817 133 : U = RgXQX_RgXQ_translate(u, mx0, T);
1818 133 : if (!tmonic) U = Q_primpart(U);
1819 511 : for (i=lx-1; i>0; i--)
1820 : {
1821 378 : GEN f = gel(fa,i), F = nfgcd(U, f, T, dent);
1822 378 : F = RgXQX_RgXQ_translate(F, x0, T);
1823 : /* F = gcd(f, u(t - x0)) [t + x0] = gcd(f(t + x0), u), more efficient */
1824 378 : if (typ(F) != t_POL || degpol(F) == 0)
1825 0 : pari_err_IRREDPOL("factornf [modulus]",T);
1826 378 : gel(P,i) = QXQX_normalize(F, T);
1827 : }
1828 133 : gen_sort_inplace(P, (void*)&cmp_RgX, &gen_cmp_RgX, NULL);
1829 133 : return P;
1830 : }
1831 :
1832 : /* Factor polynomial a on the number field defined by polynomial T, using
1833 : * Trager's trick */
1834 : GEN
1835 14 : polfnf(GEN a, GEN T)
1836 : {
1837 14 : GEN rep = cgetg(3, t_MAT), A, B, y, dent, bad;
1838 : long dA;
1839 : int tmonic;
1840 :
1841 14 : if (typ(a)!=t_POL) pari_err_TYPE("polfnf",a);
1842 14 : if (typ(T)!=t_POL) pari_err_TYPE("polfnf",T);
1843 14 : T = Q_primpart(T); tmonic = is_pm1(leading_coeff(T));
1844 14 : RgX_check_ZX(T,"polfnf");
1845 14 : A = Q_primpart( QXQX_normalize(RgX_nffix("polfnf",T,a,1), T) );
1846 14 : dA = degpol(A);
1847 14 : if (dA <= 0)
1848 : {
1849 0 : set_avma((pari_sp)(rep + 3));
1850 0 : return (dA == 0)? trivial_fact(): zerofact(varn(A));
1851 : }
1852 14 : bad = dent = absi_shallow(ZX_disc(T));
1853 14 : if (tmonic) dent = indexpartial(T, dent);
1854 14 : (void)nfgcd_all(A,RgX_deriv(A), T, dent, &B);
1855 14 : if (degpol(B) != dA) B = Q_primpart( QXQX_normalize(B, T) );
1856 14 : ensure_lt_INT(B);
1857 14 : y = nfsqff_trager(B, T, dent, FACTORS);
1858 14 : fact_from_sqff(rep, A, B, y, T, bad);
1859 14 : return sort_factor_pol(rep, cmp_RgX);
1860 : }
1861 :
1862 : static int
1863 68555 : nfsqff_use_Trager(long n, long dpol)
1864 : {
1865 68555 : return dpol*3<n;
1866 : }
1867 :
1868 : /* return the factorization of the square-free polynomial pol. Not memory-clean
1869 : The coeffs of pol are in Z_nf and its leading term is a rational integer.
1870 : deg(pol) > 0, deg(nfpol) > 1
1871 : fl is either FACTORS (return factors), or ROOTS / ROOTS_SPLIT (return roots):
1872 : - ROOTS, return only the roots of x in nf
1873 : - ROOTS_SPLIT, as ROOTS if pol splits, [] otherwise
1874 : den is usually 1, otherwise nf.zk is doubtful, and den bounds the
1875 : denominator of an arbitrary element of Z_nf on nf.zk */
1876 : static GEN
1877 42317 : nfsqff(GEN nf, GEN pol, long fl, GEN den)
1878 : {
1879 42317 : long n, nbf, dpol = degpol(pol);
1880 : GEN C0, polbase;
1881 42317 : GEN N2, res, polred, lt, nfpol = typ(nf)==t_POL?nf:nf_get_pol(nf);
1882 : ulong pp;
1883 : nfcmbf_t T;
1884 : nflift_t L;
1885 : pari_timer ti, ti_tot;
1886 :
1887 42317 : if (DEBUGLEVEL>2) { timer_start(&ti); timer_start(&ti_tot); }
1888 42317 : n = degpol(nfpol);
1889 : /* deg = 1 => irreducible */
1890 42317 : if (dpol == 1) {
1891 10227 : if (fl == FACTORS) return mkvec(QXQX_normalize(pol, nfpol));
1892 10171 : return mkvec(gneg(gdiv(gel(pol,2),gel(pol,3))));
1893 : }
1894 32090 : if (typ(nf)==t_POL || nfsqff_use_Trager(n,dpol))
1895 : {
1896 : GEN z;
1897 161 : if (DEBUGLEVEL>2) err_printf("Using Trager's method\n");
1898 161 : if (typ(nf) != t_POL) den = mulii(den, nf_get_index(nf));
1899 161 : z = nfsqff_trager(Q_primpart(pol), nfpol, den, fl);
1900 161 : if (fl != FACTORS) {
1901 119 : long i, l = lg(z);
1902 329 : for (i = 1; i < l; i++)
1903 : {
1904 252 : GEN LT, t = gel(z,i); if (degpol(t) > 1) break;
1905 210 : LT = gel(t,3);
1906 210 : if (typ(LT) == t_POL) LT = gel(LT,2); /* constant */
1907 210 : gel(z,i) = gdiv(gel(t,2), negi(LT));
1908 : }
1909 119 : setlg(z, i);
1910 119 : if (fl == ROOTS_SPLIT && i != l) return cgetg(1,t_VEC);
1911 : }
1912 161 : return z;
1913 : }
1914 :
1915 31929 : polbase = RgX_to_nfX(nf, pol);
1916 31929 : nbf = nf_pick_prime(nf, pol, fl, <, &L.Tp, &pp);
1917 31929 : if (L.Tp)
1918 : {
1919 24759 : L.Tp = Flx_to_ZX(L.Tp);
1920 24759 : L.p = utoi(pp);
1921 : }
1922 :
1923 31929 : if (fl == ROOTS_SPLIT && nbf < dpol) return cgetg(1,t_VEC);
1924 31908 : if (nbf <= 1)
1925 : {
1926 16029 : if (fl == FACTORS) return mkvec(QXQX_normalize(pol, nfpol)); /* irred. */
1927 15588 : if (!nbf) return cgetg(1,t_VEC); /* no root */
1928 : }
1929 :
1930 19015 : if (DEBUGLEVEL>2) {
1931 0 : timer_printf(&ti, "choice of a prime ideal");
1932 0 : err_printf("Prime ideal chosen: (%lu,x^%ld+...)\n", pp, degpol(L.Tp));
1933 : }
1934 19015 : L.tozk = nf_get_invzk(nf);
1935 19015 : L.topow= nf_get_zkprimpart(nf);
1936 19015 : L.topowden = nf_get_zkden(nf);
1937 19015 : if (is_pm1(den)) den = NULL;
1938 19015 : L.dn = den;
1939 19015 : T.ZC = nf_L2_bound(nf, den, NULL);
1940 19015 : T.Br = nf_root_bounds(nf, pol); if (lt) T.Br = gmul(T.Br, lt);
1941 :
1942 : /* C0 = bound for T_2(Q_i), Q | P */
1943 19015 : if (fl != FACTORS) C0 = normTp(T.Br, 2, n);
1944 12075 : else C0 = nf_factor_bound(nf, polbase);
1945 19015 : T.bound = mulrr(T.ZC, C0); /* bound for |Q_i|^2 in Z^n on chosen Z-basis */
1946 :
1947 19015 : N2 = mulur(dpol*dpol, normTp(T.Br, 4, n)); /* bound for T_2(lt * S_2) */
1948 19015 : T.BS_2 = mulrr(T.ZC, N2); /* bound for |S_2|^2 on chosen Z-basis */
1949 :
1950 19015 : if (DEBUGLEVEL>2) {
1951 0 : timer_printf(&ti, "bound computation");
1952 0 : err_printf(" 1) T_2 bound for %s: %Ps\n",
1953 : fl == FACTORS?"factor": "root", C0);
1954 0 : err_printf(" 2) Conversion from T_2 --> | |^2 bound : %Ps\n", T.ZC);
1955 0 : err_printf(" 3) Final bound: %Ps\n", T.bound);
1956 : }
1957 :
1958 19015 : bestlift_init(0, nf, T.bound, &L);
1959 19015 : if (DEBUGLEVEL>2) timer_start(&ti);
1960 19015 : polred = ZqX_normalize(polbase, lt, &L); /* monic */
1961 :
1962 19015 : if (fl != FACTORS) {
1963 6940 : GEN z = nf_DDF_roots(pol, polred, nfpol, fl, &L);
1964 6940 : if (lg(z) == 1) return cgetg(1, t_VEC);
1965 6550 : return z;
1966 : }
1967 :
1968 12075 : T.fact = gel(FqX_factor(polred, L.Tp, L.p), 1);
1969 12075 : if (DEBUGLEVEL>2)
1970 0 : timer_printf(&ti, "splitting mod %Ps^%ld", L.p, degpol(L.Tp));
1971 12075 : T.L = &L;
1972 12075 : T.polbase = polbase;
1973 12075 : T.pol = pol;
1974 12075 : T.nf = nf;
1975 12075 : res = nf_combine_factors(&T, polred, dpol-1);
1976 12075 : if (DEBUGLEVEL>2)
1977 0 : err_printf("Total Time: %ld\n===========\n", timer_delay(&ti_tot));
1978 12075 : return res;
1979 : }
1980 :
1981 : /* assume pol monic in nf.zk[X] */
1982 : GEN
1983 21 : nfroots_if_split(GEN *pnf, GEN pol)
1984 : {
1985 21 : GEN T = get_nfpol(*pnf,pnf), den = fix_nf(pnf, &T, &pol);
1986 21 : pari_sp av = avma;
1987 21 : GEN z = nfsqff(*pnf, pol, ROOTS_SPLIT, den);
1988 21 : if (lg(z) == 1) return gc_NULL(av);
1989 0 : return gc_GEN(av, z);
1990 : }
1991 :
1992 : /*******************************************************************/
1993 : /* */
1994 : /* Roots of unity in a number field */
1995 : /* (alternative to nfrootsof1 using factorization in K[X]) */
1996 : /* */
1997 : /*******************************************************************/
1998 : /* Code adapted from nffactor. Structure of the algorithm; only step 1 is
1999 : * specific to roots of unity.
2000 : *
2001 : * [Step 1]: guess roots via ramification. If trivial output this.
2002 : * [Step 2]: select prime [p] unramified and ideal [pr] above
2003 : * [Step 4]: factor the cyclotomic polynomial mod [pr],
2004 : */
2005 :
2006 : /* Guesses the number of roots of unity in number field [nf].
2007 : * Computes gcd of N(P)-1 for some primes. The value returned is a proven
2008 : * multiple of the correct value. */
2009 : static long
2010 25046 : guess_roots(GEN T, GEN D, GEN index)
2011 : {
2012 25046 : long c = 0, nfdegree = degpol(T), B = nfdegree + 20, l;
2013 25046 : ulong p = 2;
2014 25046 : GEN nbroots = NULL;
2015 : forprime_t S;
2016 : pari_sp av;
2017 :
2018 25046 : (void)u_forprime_init(&S, 3, ULONG_MAX);
2019 25046 : av = avma;
2020 : /* result must be stationary (counter c) for at least B loops */
2021 745395 : for (l=1; (p = u_forprime_next(&S)); l++)
2022 : {
2023 : GEN old, F, pf_1, Tp;
2024 745395 : ulong i, gcdf = 0;
2025 : long nb;
2026 :
2027 745395 : if (!umodiu(D,p) || !umodiu(index,p)) continue;
2028 708242 : Tp = ZX_to_Flx(T,p); /* squarefree */
2029 708242 : F = Flx_nbfact_by_degree(Tp, &nb, p);
2030 : /* the gcd of the p^f - 1 is p^(gcd of the f's) - 1 */
2031 2585623 : for (i = 1; i <= (ulong) nfdegree; i++)
2032 2145130 : if (F[i]) {
2033 709068 : gcdf = gcdf? ugcd(gcdf, i): i;
2034 709068 : if (gcdf == 1) break;
2035 : }
2036 708242 : pf_1 = subiu(powuu(p, gcdf), 1);
2037 708242 : old = nbroots;
2038 708242 : nbroots = nbroots? gcdii(pf_1, nbroots): pf_1;
2039 708242 : if (DEBUGLEVEL>5)
2040 0 : err_printf("p=%lu; gcf(f(P/p))=%ld; nbroots | %Ps",p, gcdf, nbroots);
2041 : /* if same result go on else reset the stop counter [c] */
2042 708242 : if (old && equalii(nbroots,old))
2043 653085 : { if (!is_bigint(nbroots) && ++c > B) break; }
2044 : else
2045 55157 : c = 0;
2046 : }
2047 25046 : if (!nbroots) pari_err_OVERFLOW("guess_roots [ran out of primes]");
2048 25046 : if (DEBUGLEVEL>5) err_printf("%ld loops\n",l);
2049 25046 : return gc_long(av, itos(nbroots));
2050 : }
2051 :
2052 : /* T(x) an irreducible ZX. Is it of the form Phi_n(c \pm x) ?
2053 : * Return NULL if not, and a root of 1 of maximal order in Z[x]/(T) otherwise
2054 : *
2055 : * N.B. Set n_squarefree = 1 if n is squarefree, and 0 otherwise.
2056 : * This last parameter is inconvenient, but it allows a cheap
2057 : * stringent test. (n guessed from guess_roots())*/
2058 : static GEN
2059 2989 : ZXirred_is_cyclo_translate(GEN T, long n_squarefree)
2060 : {
2061 2989 : long r, m, d = degpol(T);
2062 2989 : GEN T1, c = divis_rem(gel(T, d+1), d, &r); /* d-1 th coeff divisible by d ? */
2063 : /* The trace coefficient of polcyclo(n) is \pm1 if n is square free, and 0
2064 : * otherwise. */
2065 2989 : if (!n_squarefree)
2066 1162 : { if (r) return NULL; }
2067 : else
2068 : {
2069 1827 : if (r < -1)
2070 : {
2071 0 : r += d;
2072 0 : c = subiu(c, 1);
2073 : }
2074 1827 : else if (r == d-1)
2075 : {
2076 49 : r = -1;
2077 49 : c = addiu(c, 1);
2078 : }
2079 1827 : if (r != 1 && r != -1) return NULL;
2080 : }
2081 2723 : if (signe(c)) /* presumably Phi_guess(c \pm x) */
2082 49 : T = RgX_Rg_translate(T, negi(c));
2083 2723 : if (!n_squarefree) T = RgX_deflate_max(T, &m);
2084 : /* presumably Phi_core(guess)(\pm x), cyclotomic iff original T was */
2085 2723 : T1 = ZX_graeffe(T);
2086 2723 : if (ZX_equal(T, T1)) /* T = Phi_n, n odd */
2087 434 : return deg1pol_shallow(gen_m1, negi(c), varn(T));
2088 2289 : else if (ZX_equal(T1, ZX_z_unscale(T, -1))) /* T = Phi_{2n}, nodd */
2089 2135 : return deg1pol_shallow(gen_1, c, varn(T));
2090 154 : return NULL;
2091 : }
2092 :
2093 : static GEN
2094 55587 : trivroots(void) { return mkvec2(gen_2, gen_m1); }
2095 : /* Number of roots of unity in number field [nf]. */
2096 : GEN
2097 64771 : nfrootsof1(GEN nf)
2098 : {
2099 : GEN T, fa, LP, LE, z, disc, index;
2100 : pari_timer ti;
2101 : long i, l, nbguessed, nbroots, nfdegree;
2102 : pari_sp av;
2103 :
2104 64771 : T = get_nfpol(nf, &nf); nfdegree = degpol(T);
2105 64771 : RgX_check_ZX(T, "nfrootsof1");
2106 64771 : if (nfdegree <= 0) pari_err_CONSTPOL("nfrootsof1");
2107 64764 : if (nf)
2108 : {
2109 64750 : if (nf_get_r1(nf)) return trivroots();
2110 25046 : disc = nf_get_disc(nf);
2111 25046 : index = nf_get_index(nf);
2112 : }
2113 : else
2114 : {
2115 14 : if (odd(nfdegree) || signe(leading_coeff(T)) != signe(constant_coeff(T)))
2116 14 : return trivroots();
2117 0 : disc = ZX_disc(T);
2118 0 : index = gen_1;
2119 : }
2120 : /* Step 1 : guess number of roots and discard trivial case 2 */
2121 25046 : if (DEBUGLEVEL>2) timer_start(&ti);
2122 25046 : nbguessed = guess_roots(T, disc, index);
2123 25046 : if (DEBUGLEVEL>2)
2124 0 : timer_printf(&ti, "guessing roots of 1 [guess = %ld]", nbguessed);
2125 25046 : if (nbguessed == 2) return trivroots();
2126 :
2127 9184 : fa = factoru(nbguessed);
2128 9184 : LP = gel(fa,1); l = lg(LP);
2129 9184 : LE = gel(fa,2);
2130 25543 : for (i = 1; i < l; i++)
2131 : {
2132 16359 : long p = LP[i];
2133 : /* Degree and ramification test: find largest k such that Q(zeta_{p^k})
2134 : * may be a subfield of K. Q(zeta_p^k) has degree (p-1)p^(k-1)
2135 : * and v_p(discriminant) = ((p-1)k-1)p^(k-1); so we must have
2136 : * v_p(disc_K) >= ((p-1)k-1) * n / (p-1) = kn - q, where q = n/(p-1) */
2137 16359 : if (p == 2)
2138 : { /* the test simplifies a little in that case */
2139 : long v, vnf, k;
2140 9184 : if (LE[i] == 1) continue;
2141 2926 : vnf = vals(nfdegree);
2142 2926 : v = vali(disc);
2143 3038 : for (k = minss(LE[i], vnf+1); k >= 1; k--)
2144 3038 : if (v >= nfdegree*(k-1)) { nbguessed >>= LE[i]-k; LE[i] = k; break; }
2145 : /* N.B the test above always works for k = 1: LE[i] >= 1 */
2146 : }
2147 : else
2148 : {
2149 : long v, vnf, k;
2150 7175 : ulong r, q = udivuu_rem(nfdegree, p-1, &r);
2151 7175 : if (r) { nbguessed /= upowuu(p, LE[i]); LE[i] = 0; continue; }
2152 : /* q = n/(p-1) */
2153 7175 : vnf = u_lval(q, p);
2154 7175 : v = Z_lval(disc, p);
2155 7196 : for (k = minss(LE[i], vnf+1); k >= 0; k--)
2156 7196 : if (v >= nfdegree*k-(long)q)
2157 7175 : { nbguessed /= upowuu(p, LE[i]-k); LE[i] = k; break; }
2158 : /* N.B the test above always works for k = 0: LE[i] >= 0 */
2159 : }
2160 : }
2161 9184 : if (DEBUGLEVEL>2)
2162 0 : timer_printf(&ti, "after ramification conditions [guess = %ld]", nbguessed);
2163 9184 : if (nbguessed == 2) return trivroots();
2164 9177 : av = avma;
2165 :
2166 : /* Step 1.5 : test if nf.pol == subst(polcyclo(nbguessed), x, \pm x+c) */
2167 9177 : if (eulerphiu_fact(fa) == (ulong)nfdegree)
2168 : {
2169 2989 : z = ZXirred_is_cyclo_translate(T, uissquarefree_fact(fa));
2170 2989 : if (DEBUGLEVEL>2) timer_printf(&ti, "checking for cyclotomic polynomial");
2171 2989 : if (z)
2172 : {
2173 2569 : if (nf) z = nf_to_scalar_or_basis(nf,z);
2174 2569 : return gc_GEN(av, mkvec2(utoipos(nbguessed), z));
2175 : }
2176 420 : set_avma(av);
2177 : }
2178 :
2179 : /* Step 2 : actual computation of roots */
2180 6608 : nbroots = 2; z = scalarpol(gen_m1, varn(T));
2181 18375 : for (i = 1; i < l; i++)
2182 : { /* for all prime power factors of nbguessed, find a p^k-th root of unity */
2183 11767 : long k, p = LP[i];
2184 16555 : for (k = LE[i]; k > 0; k--)
2185 : { /* find p^k-th roots */
2186 11781 : pari_sp av = avma;
2187 11781 : long pk = upowuu(p,k);
2188 : GEN r;
2189 11781 : if (pk==2) continue; /* no need to test second roots ! */
2190 7070 : r = nfisincl0(polcyclo(pk, 0), T, 1);
2191 7070 : if (!isintzero(r))
2192 : {
2193 6993 : if (DEBUGLEVEL>2) err_printf(" %s root of unity found\n",uordinal(pk));
2194 6993 : if (p==2) { nbroots = pk; z = r; }
2195 : else
2196 : {
2197 5096 : nbroots *= pk;
2198 5096 : z = nf ? nfmul(nf, z, r): QXQ_mul(z, r, T);
2199 : }
2200 6993 : break;
2201 : }
2202 77 : set_avma(av);
2203 77 : if (DEBUGLEVEL) pari_warn(warner,"nfrootsof1: wrong guess");
2204 : }
2205 : }
2206 6608 : if (nf) z = nf_to_scalar_or_basis(nf,z);
2207 6608 : return gc_GEN(av, mkvec2(utoi(nbroots), z));
2208 : }
2209 :
2210 : /*******************************************************************/
2211 : /* */
2212 : /* rnfgaloisconj */
2213 : /* */
2214 : /*******************************************************************/
2215 :
2216 : static GEN
2217 35 : rnffrobeniuslift(GEN nf, GEN P, GEN Tp, GEN p, GEN bound, GEN den, GEN iden, GEN T)
2218 : {
2219 : GEN q, Tq, Pp, Pq, s, S;
2220 35 : long v = varn(P), e;
2221 : nflift_t L;
2222 35 : L.p = p; L.Tp = Tp;
2223 35 : L.tozk = nf_get_invzk(nf);
2224 35 : L.topow= nf_get_zkprimpart(nf);
2225 35 : L.topowden = nf_get_zkden(nf);
2226 35 : bestlift_init(0, nf, bound, &L);
2227 35 : q = L.pk; Tq = L.Tpk; e = L.k;
2228 35 : P = gmul(P,L.topowden);
2229 35 : Pq = ZqX(RgX_to_nfX(nf, P), q, Tq, L.ZqProj);
2230 35 : Pp = FqX_red(Pq, Tp, p);
2231 35 : s = FpXQXQ_pow(pol_x(v), powiu(p, degpol(Tp)), Pp, Tp, p);
2232 35 : S = ZqX_ZqXQ_liftroot(Pq, s, Pq, Tq, p, e);
2233 35 : den = nf_to_Zq(den, Tq, q, shifti(q,-1), L.ZqProj);
2234 35 : S = FqX_Fq_mul(S, den, Tq, q);
2235 35 : P = nf_pol_lift(S, bound, &L);
2236 35 : if (!P) return NULL;
2237 28 : return gdiv(QXQX_QXQ_mul(P, iden, T), L.topowden);
2238 : }
2239 :
2240 : static GEN
2241 63 : RgX_galconj_bound(GEN T, long prec)
2242 : {
2243 63 : GEN L = QX_complex_roots(T, prec);
2244 63 : GEN prep = vandermondeinverseinit(L);
2245 63 : GEN M = vandermondeinverse(L, RgX_gtofp(T, prec), gen_1, prep);
2246 63 : return gmul(matrixnorm(M, prec), gsupnorm(L, prec));
2247 : }
2248 :
2249 : static GEN
2250 42 : rnfgaloisconj_bound(GEN P, GEN z, GEN den, long prec)
2251 : {
2252 42 : pari_sp av = avma;
2253 42 : long i, l = lg(z);
2254 42 : GEN s = gen_0;
2255 105 : for (i = 1; i < l; i++)
2256 : {
2257 63 : GEN a = gabs(RgX_cxeval(den, gel(z,i), NULL),prec);
2258 63 : GEN b = RgX_galconj_bound(RgXY_cxevalx(P, gel(z,i), NULL), prec);
2259 63 : s = gc_upto(av, gadd(s, gsqr(gmul(a,b))));
2260 : }
2261 42 : return gc_upto(av, ceil_safe(s));
2262 : }
2263 :
2264 : static GEN
2265 0 : FlxqXQV_fixedalg(GEN aut, GEN S, GEN T, ulong p)
2266 : {
2267 0 : pari_sp av = avma;
2268 0 : long d = get_FlxqX_degree(S), sv = get_Flx_var(T), v = get_FlxqX_var(S);
2269 0 : long i, l = lg(aut);
2270 0 : GEN A = RgX_to_FlxqX(gel(aut,1), T, p);
2271 0 : GEN M = FlxXV_to_FlxM(FlxqXQ_powers(A, d-1, S, T, p), d, sv);
2272 0 : GEN M2 = FlxM_Flx_add_shallow(M, Fl_to_Flx(p-1, sv), p);
2273 0 : GEN C = FlxM_to_FlxXV(FlxqM_ker(M2, T, p), v);
2274 0 : for (i = 2; i < l; i++)
2275 : {
2276 0 : GEN A = RgX_to_FlxqX(gel(aut,i), T, p);
2277 0 : GEN M = FlxXC_sub(FlxqXC_FlxqXQ_eval(C, A, S, T, p), C, p);
2278 0 : GEN K = FlxqM_ker(FlxXV_to_FlxM(M, d, sv), T, p);
2279 0 : C = FlxM_to_FlxXV(FlxqM_mul(FlxXV_to_FlxM(C, d, sv), K, T, p), v);
2280 : }
2281 0 : return gc_GEN(av, C);
2282 : }
2283 :
2284 : static long
2285 0 : fixedalg_order(GEN aut, GEN S, GEN T, ulong p)
2286 : {
2287 0 : pari_sp av = avma;
2288 0 : long e, i, d = get_Flx_degree(T), l = lg(aut);
2289 0 : for (e = 1;; e++, set_avma(av))
2290 : {
2291 0 : for (i = 1; i < l; i++, set_avma(av))
2292 0 : if (!gequal(FlxqXQ_pow(gel(aut,i), powuu(p, d*e), S, T, p), gel(aut,i)))
2293 0 : break;
2294 0 : if (i == l) return e;
2295 : }
2296 : }
2297 :
2298 : static GEN
2299 35 : rnfgaloisanalysis(GEN nf, GEN P, GEN aut, long m, GEN d, long *pt_o)
2300 : {
2301 35 : GEN T = nf_get_pol(nf), R = NULL, den = nf_get_zkden(nf);
2302 35 : long n = degpol(P), omax = 0, ntry = 0;
2303 35 : GEN daut = lg(aut) > 1 ? Q_denom(aut) : NULL;
2304 : forprime_t S;
2305 35 : u_forprime_init(&S, 2, ULONG_MAX);
2306 182 : while (ntry <= n || omax < 2)
2307 : {
2308 182 : ulong p = u_forprime_next(&S);
2309 : GEN F, Tp;
2310 : long i, lF;
2311 182 : if ((d && dvdiu(d,p)) || dvdiu(den, p) || (daut && dvdiu(daut,p))) continue;
2312 161 : Tp = ZX_to_Flx(T, p);
2313 161 : if (!Flx_is_squarefree(Tp,p)) continue;
2314 140 : F = gel(Flx_factor(Tp, p), 1); lF = lg(F);
2315 294 : for (i = 1; i < lF; i++)
2316 : {
2317 189 : pari_sp av = avma;
2318 : long o, d;
2319 189 : GEN D, Fi = gel(F,i), Pp = RgX_to_FlxqX(P, Fi, p);
2320 189 : if (degpol(Pp) < n || !FlxqX_is_squarefree(Pp, Fi, p)) continue;
2321 147 : D = FlxqX_nbfact_by_degree(Pp, &d, Fi, p);
2322 147 : o = n / d; /* d factors, all should have degree o */
2323 147 : if (D[o] != d)
2324 : {
2325 0 : if(DEBUGLEVEL) err_printf("rnfisabelian: not Galois at %lu: %Ps\n",p,D);
2326 35 : return NULL;
2327 : }
2328 147 : if (lg(aut) > 1)
2329 : {
2330 0 : GEN U = FlxqXQV_fixedalg(aut, Pp, Fi, p);
2331 0 : o = fixedalg_order(U, Pp, Fi, p);
2332 0 : d = m / o;
2333 : }
2334 147 : if (d == 1) { *pt_o = o; return mkvec2(Flx_to_ZX(Fi), utoi(p)); }
2335 112 : if (o > omax) { omax = o; R = mkvec2(Flx_to_ZX(Fi), utoi(p)); }
2336 84 : else set_avma(av);
2337 : }
2338 105 : ntry++;
2339 : }
2340 0 : *pt_o = omax; return R;
2341 : }
2342 :
2343 : static GEN
2344 42 : rnfabelianconjgen_i(GEN nf, GEN P)
2345 : {
2346 : long prec, m, i, l;
2347 : GEN bnd, den, iden, T, Pr, gen, orders, d;
2348 42 : nf = checknf(nf);
2349 42 : T = nf_get_pol(nf);
2350 42 : P = RgX_nffix("rnfgaloisconj", T, P, 1);
2351 42 : P = Q_remove_denom(P, &d);
2352 42 : Pr = RgX_to_nfX(nf, P);
2353 42 : prec = DEFAULTPREC + nbits2extraprec(gexpo(Pr) * degpol(T));
2354 42 : nf = nfnewprec_shallow(nf, prec);
2355 42 : den = nfX_disc(nf, P);
2356 42 : bnd = rnfgaloisconj_bound(P, nf_get_roots(nf), den, prec);
2357 42 : iden = QXQ_inv(den, T);
2358 42 : den = nf_to_scalar_or_basis(nf, den);
2359 42 : m = degpol(P); l = expu(m) + 1;
2360 42 : gen = cgetg(l, t_VEC); orders = cgetg(l,t_VECSMALL);
2361 70 : for (i = 1; m > 1; i++)
2362 : {
2363 : long o;
2364 : GEN S, Tp;
2365 35 : setlg(gen, i);
2366 35 : Tp = rnfgaloisanalysis(nf, P, gen, m, d, &o);
2367 42 : if (!Tp) return NULL;
2368 35 : S = rnffrobeniuslift(nf, P, gel(Tp,1), gel(Tp,2), bnd, den, iden, T);
2369 35 : if (!S) return NULL;
2370 28 : gel(gen,i) = S; orders[i] = o; m /= o;
2371 : }
2372 35 : setlg(gen,i); setlg(orders,i); return mkvec2(gen, orders);
2373 : }
2374 : GEN
2375 0 : rnfabelianconjgen(GEN nf, GEN P)
2376 : {
2377 0 : pari_sp av = avma;
2378 0 : GEN G = rnfabelianconjgen_i(nf, P);
2379 0 : return G? gc_GEN(av, G): gc_const(av, gen_0);
2380 : }
2381 :
2382 : long
2383 42 : rnfisabelian(GEN nf, GEN pol)
2384 : {
2385 42 : pari_sp av = avma;
2386 42 : GEN G = rnfabelianconjgen_i(nf, pol);
2387 42 : return gc_long(av, G? 1: 0);
2388 : }
|
