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. It is distributed in the hope that it will be useful, but WITHOUT
8 : ANY WARRANTY WHATSOEVER.
9 :
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 : /* POLYNOMIAL FACTORIZATION IN A NUMBER FIELD */
17 : /* */
18 : /*******************************************************************/
19 : #include "pari.h"
20 : #include "paripriv.h"
21 :
22 : static GEN nfsqff(GEN nf,GEN pol,long fl,GEN den);
23 : static int nfsqff_use_Trager(long n, long dpol);
24 :
25 : enum { FACTORS = 0, ROOTS, ROOTS_SPLIT };
26 :
27 : /* for nf_bestlift: reconstruction of algebraic integers known mod P^k,
28 : * P maximal ideal above p */
29 : typedef struct {
30 : long k; /* input known mod P^k */
31 : GEN p, pk; /* p^k = denom(prk^-1) [ assume pr unramified ]*/
32 : GEN prk; /* |.|^2 LLL-reduced basis (b_i) of P^k (NOT T2-reduced) */
33 : GEN iprk; /* den * prk^-1 */
34 : GEN GSmin; /* min |b_i^*|^2 */
35 :
36 : GEN Tp; /* Tpk mod p */
37 : GEN Tpk;
38 : GEN ZqProj;/* projector to Zp / P^k = Z/p^k[X] / Tpk */
39 :
40 : GEN tozk;
41 : GEN topow;
42 : GEN topowden; /* topow x / topowden = basistoalg(x) */
43 : GEN dn; /* NULL (we trust nf.zk) or a t_INT > 1 (an alg. integer has
44 : denominator dividing dn, when expressed on nf.zk */
45 : } nflift_t;
46 :
47 : typedef struct
48 : {
49 : nflift_t *L;
50 : GEN nf;
51 : GEN pol, polbase; /* leading coeff is a t_INT */
52 : GEN fact;
53 : GEN Br, bound, ZC, BS_2;
54 : } nfcmbf_t;
55 :
56 : /*******************************************************************/
57 : /* RATIONAL RECONSTRUCTION (use ratlift) */
58 : /*******************************************************************/
59 : /* NOT stack clean. a, b stay on the stack */
60 : static GEN
61 12980640 : lift_to_frac(GEN t, GEN N, GEN amax, GEN bmax, GEN den, GEN tden)
62 : {
63 : GEN a, b;
64 12980640 : if (signe(t) < 0) t = addii(t, N); /* in case t is a centerlift */
65 12980637 : if (tden)
66 : {
67 3498453 : pari_sp av = avma;
68 3498453 : a = Fp_center_i(Fp_mul(t, tden, N), N, shifti(N,-1));
69 3497983 : if (abscmpii(a, amax) < 0) return gerepileupto(av, Qdivii(a, tden));
70 360932 : set_avma(av);
71 : }
72 9843116 : if (!Fp_ratlift(t, N, amax,bmax, &a,&b)
73 9771269 : || (den && !dvdii(den,b))
74 9847049 : || !is_pm1(gcdii(a,b))) return NULL;
75 9766724 : if (is_pm1(b)) { cgiv(b); return a; }
76 1000224 : return mkfrac(a, b);
77 : }
78 :
79 : /* Compute rational lifting for all the components of P modulo N. Assume
80 : * Fp_ratlift preconditions are met; we allow centerlifts. If one component
81 : * fails, return NULL. If den != NULL, check that the deninators divide den;
82 : * assume (N, den) = 1. */
83 : GEN
84 1565742 : FpC_ratlift(GEN P, GEN N, GEN amax, GEN bmax, GEN den)
85 : {
86 1565742 : pari_sp av = avma;
87 : long j, l;
88 1565742 : GEN tden = NULL, Q = cgetg_copy(P, &l);
89 1565743 : if (l==1) return Q;
90 1565743 : if (den && cmpii(bmax, den) > 0) bmax = den;
91 14143376 : for (j = 1; j < l; ++j)
92 : {
93 12629699 : GEN a = lift_to_frac(gel(P,j), N, amax, bmax, den, tden);
94 12630065 : if (!a) return gc_NULL(av);
95 12578117 : if (typ(a) == t_FRAC)
96 : {
97 2417090 : GEN d = gel(a,2);
98 2417090 : tden = tden? (cmpii(tden, d) < 0? d: tden): d;
99 : }
100 12577633 : gel(Q,j) = a;
101 : }
102 1513677 : return Q;
103 : }
104 : GEN
105 115018 : FpX_ratlift(GEN P, GEN N, GEN amax, GEN bmax, GEN den)
106 : {
107 115018 : pari_sp av = avma;
108 : long j, l;
109 115018 : GEN tden = NULL, Q = cgetg_copy(P, &l);
110 115018 : Q[1] = P[1];
111 115018 : if (den && cmpii(bmax, den) > 0) bmax = den;
112 438040 : for (j = 2; j < l; ++j)
113 : {
114 349899 : GEN a = lift_to_frac(gel(P,j), N, amax, bmax, den, tden);
115 349899 : if (!a) return gc_NULL(av);
116 323022 : if (typ(a) == t_FRAC)
117 : {
118 156189 : GEN d = gel(a,2);
119 156189 : tden = tden? (cmpii(tden, d) < 0? d: tden): d;
120 : }
121 323022 : gel(Q,j) = a;
122 : }
123 88141 : return Q;
124 : }
125 :
126 : GEN
127 666451 : FpM_ratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN den)
128 : {
129 666451 : pari_sp av = avma;
130 666451 : long j, l = lg(M);
131 666451 : GEN N = cgetg_copy(M, &l);
132 666451 : if (l == 1) return N;
133 2112380 : for (j = 1; j < l; ++j)
134 : {
135 1494060 : GEN a = FpC_ratlift(gel(M, j), mod, amax, bmax, den);
136 1494060 : if (!a) return gc_NULL(av);
137 1445929 : gel(N,j) = a;
138 : }
139 618320 : return N;
140 : }
141 :
142 : /*******************************************************************/
143 : /* GCD in K[X], K NUMBER FIELD */
144 : /*******************************************************************/
145 : /* P a non-zero ZXQX */
146 : static GEN
147 20155 : lead_simplify(GEN P)
148 : {
149 20155 : GEN x = gel(P, lg(P)-1); /* x a non-zero ZX or t_INT */
150 20155 : if (typ(x) == t_POL)
151 : {
152 2870 : if (degpol(x)) return x;
153 2590 : x = gel(x,2);
154 : }
155 19875 : return is_pm1(x)? NULL: x;
156 : }
157 : /* P,Q in Z[X,Y], T in Z[Y] irreducible. compute GCD in Q[Y]/(T)[X].
158 : *
159 : * M. Encarnacion "On a modular Algorithm for computing GCDs of polynomials
160 : * over number fields" (ISSAC'94).
161 : *
162 : * We procede as follows
163 : * 1:compute the gcd modulo primes discarding bad primes as they are detected.
164 : * 2:reconstruct the result via FpM_ratlift, stoping as soon as we get weird
165 : * denominators.
166 : * 3:if FpM_ratlift succeeds, try the full division.
167 : * Suppose accuracy is insufficient to get the result right: FpM_ratlift will
168 : * rarely succeed, and even if it does the polynomial we get has sensible
169 : * coefficients, so the full division will not be too costly.
170 : *
171 : * If not NULL, den must be a multiple of the denominator of the gcd,
172 : * for example the discriminant of T.
173 : *
174 : * NOTE: if T is not irreducible, nfgcd may loop forever, esp. if gcd | T */
175 : GEN
176 15604 : nfgcd_all(GEN P, GEN Q, GEN T, GEN den, GEN *Pnew)
177 : {
178 15604 : pari_sp btop, ltop = avma;
179 15604 : GEN lP, lQ, M, dsol, R, bo, sol, mod = NULL, lden = NULL;
180 15604 : long vP = varn(P), vT = varn(T), dT = degpol(T), dM = 0, dR;
181 : forprime_t S;
182 :
183 15604 : if (!signe(P)) { if (Pnew) *Pnew = pol_0(vT); return gcopy(Q); }
184 15604 : if (!signe(Q)) { if (Pnew) *Pnew = pol_1(vT); return gcopy(P); }
185 : /* Compute denominators */
186 15520 : if ((lP = lead_simplify(P)) && (lQ = lead_simplify(Q)))
187 : {
188 4103 : if (typ(lP) == t_INT && typ(lQ) == t_INT)
189 3963 : lden = powiu(gcdii(lP, lQ), dT);
190 140 : else if (typ(lP) == t_INT)
191 0 : lden = gcdii(powiu(lP, dT), ZX_resultant(lQ, T));
192 140 : else if (typ(lQ) == t_INT)
193 0 : lden = gcdii(powiu(lQ, dT), ZX_resultant(lP, T));
194 : else
195 140 : lden = gcdii(ZX_resultant(lP, T), ZX_resultant(lQ, T));
196 4103 : if (is_pm1(lden)) lden = NULL;
197 4103 : if (den && lden) den = mulii(den, lden);
198 : }
199 15520 : init_modular_small(&S);
200 15520 : btop = avma;
201 : for(;;)
202 4145 : {
203 19665 : ulong p = u_forprime_next(&S);
204 : GEN Tp;
205 19665 : if (!p) pari_err_OVERFLOW("nfgcd [ran out of primes]");
206 : /*Discard primes dividing disc(T) or lc(PQ) */
207 19665 : if (lden && !umodiu(lden, p)) continue;
208 19665 : Tp = ZX_to_Flx(T,p);
209 19665 : if (!Flx_is_squarefree(Tp, p)) continue;
210 : /*Discard primes when modular gcd does not exist*/
211 19665 : if ((R = FlxqX_safegcd(ZXX_to_FlxX(P,p,vT),
212 : ZXX_to_FlxX(Q,p,vT),
213 0 : Tp, p)) == NULL) continue;
214 19665 : dR = degpol(R);
215 19665 : if (dR == 0) { set_avma(ltop); if (Pnew) *Pnew = P; return pol_1(vP); }
216 5860 : if (mod && dR > dM) continue; /* p divides Res(P/gcd, Q/gcd). Discard. */
217 :
218 5860 : R = FlxX_to_Flm(R, dT);
219 : /* previous primes divided Res(P/gcd, Q/gcd)? Discard them. */
220 5860 : if (!mod || dR < dM) { M = ZM_init_CRT(R, p); mod = utoipos(p); dM = dR; continue; }
221 4145 : (void)ZM_incremental_CRT(&M,R, &mod,p);
222 4145 : if (gc_needed(btop, 1))
223 : {
224 0 : if (DEBUGMEM>1) pari_warn(warnmem,"nfgcd");
225 0 : gerepileall(btop, 2, &M, &mod);
226 : }
227 : /* I suspect it must be better to take amax > bmax*/
228 4145 : bo = sqrti(shifti(mod, -1));
229 4145 : if ((sol = FpM_ratlift(M, mod, bo, bo, den)) == NULL) continue;
230 1745 : sol = RgM_to_RgXX(sol,vP,vT);
231 1745 : dsol = Q_primpart(sol);
232 :
233 1745 : if (!ZXQX_dvd(Q, dsol, T)) continue;
234 1715 : if (Pnew)
235 : {
236 203 : *Pnew = RgXQX_pseudodivrem(P, dsol, T, &R);
237 203 : if (signe(R)) continue;
238 : }
239 : else
240 : {
241 1512 : if (!ZXQX_dvd(P, dsol, T)) continue;
242 : }
243 1715 : gerepileall(ltop, Pnew? 2: 1, &dsol, Pnew);
244 1715 : return dsol; /* both remainders are 0 */
245 : }
246 : }
247 : GEN
248 4278 : nfgcd(GEN P, GEN Q, GEN T, GEN den)
249 4278 : { return nfgcd_all(P, Q, T, den, NULL); }
250 :
251 : GEN
252 3556 : ZXQX_gcd(GEN P, GEN Q, GEN T)
253 3556 : { return nfgcd_all(P, Q, T, NULL, NULL); }
254 :
255 : GEN
256 1288 : QXQX_gcd(GEN P, GEN Q, GEN T)
257 : {
258 1288 : pari_sp av = avma;
259 1288 : GEN P1 = Q_remove_denom(P, NULL);
260 1288 : GEN Q1 = Q_remove_denom(Q, NULL);
261 1288 : return gerepileupto(av, ZXQX_gcd(P1, Q1, T));
262 : }
263 :
264 : int
265 4816 : nfissquarefree(GEN nf, GEN x)
266 : {
267 4816 : pari_sp av = avma;
268 4816 : GEN g, y = RgX_deriv(x);
269 4816 : if (RgX_is_rational(x)) g = QX_gcd(x, y);
270 : else
271 : {
272 3305 : GEN T = get_nfpol(nf,&nf);
273 3305 : x = Q_primpart( liftpol_shallow(x) );
274 3305 : y = Q_primpart( liftpol_shallow(y) );
275 3305 : g = nfgcd(x, y, T, nf? nf_get_index(nf): NULL);
276 : }
277 4816 : 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 gerepilecopy(av, rep);
306 : }
307 :
308 : /*******************************************************************/
309 : /* MAIN ROUTINES nfroots / nffactor */
310 : /*******************************************************************/
311 : static GEN
312 8631 : QXQX_normalize(GEN P, GEN T)
313 : {
314 8631 : GEN P0 = leading_coeff(P);
315 8631 : long t = typ(P0);
316 8631 : if (t == t_POL)
317 : {
318 1211 : if (degpol(P0)) return RgXQX_RgXQ_mul(P, QXQ_inv(P0,T), T);
319 455 : P0 = gel(P0,2); t = typ(P0);
320 : }
321 : /* t = t_INT/t_FRAC */
322 7875 : if (t == t_INT && is_pm1(P0) && signe(P0) > 0) return P; /* monic */
323 2681 : return RgX_Rg_div(P, P0);
324 : }
325 : /* assume leading term of P is an integer */
326 : static GEN
327 7623 : RgX_int_normalize(GEN P)
328 : {
329 7623 : GEN P0 = leading_coeff(P);
330 : /* cater for t_POL */
331 7623 : if (typ(P0) == t_POL)
332 : {
333 143 : P0 = gel(P0,2); /* non-0 constant */
334 143 : P = shallowcopy(P);
335 143 : gel(P,lg(P)-1) = P0; /* now leading term is a t_INT */
336 : }
337 7623 : if (typ(P0) != t_INT) pari_err_BUG("RgX_int_normalize");
338 7623 : if (is_pm1(P0)) return signe(P0) > 0? P: RgX_neg(P);
339 3878 : 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 non-monic polynomials */
344 : static GEN
345 2051 : proper_nf(GEN nf)
346 2051 : { 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 7567 : fix_nf(GEN *pnf, GEN *pT, GEN *pA)
353 : {
354 7567 : GEN nf, NF, fa, P, Q, q, D, T = *pT;
355 : nfmaxord_t S;
356 : long i, l;
357 :
358 7567 : if (*pnf) return gen_1;
359 2051 : nfmaxord(&S, T, nf_PARTIALFACT);
360 2051 : NF = nfinit_complete(&S, 0, DEFAULTPREC);
361 2051 : *pnf = nf = proper_nf(NF);
362 2051 : if (nf != NF) { /* t_POL defining base field changed (not monic) */
363 35 : GEN A = *pA, a = cgetg_copy(A, &l);
364 35 : GEN rev = gel(NF,2), pow, dpow;
365 :
366 35 : *pT = T = nf_get_pol(nf); /* need to update T */
367 35 : pow = QXQ_powers(lift_shallow(rev), degpol(T)-1, T);
368 35 : pow = Q_remove_denom(pow, &dpow);
369 35 : a[1] = A[1];
370 154 : for (i=2; i<l; i++) {
371 119 : GEN c = gel(A,i);
372 119 : if (typ(c) == t_POL) c = QX_ZXQV_eval(c, pow, dpow);
373 119 : gel(a,i) = c;
374 : }
375 35 : *pA = Q_primpart(a); /* need to update A */
376 : }
377 :
378 2051 : D = nf_get_disc(nf);
379 2051 : if (is_pm1(D)) return gen_1;
380 2044 : fa = absZ_factor_limit(D, 0);
381 2044 : P = gel(fa,1); q = gel(P, lg(P)-1);
382 2044 : if (BPSW_psp(q)) return gen_1;
383 : /* nf_get_disc(nf) may be incorrect */
384 14 : P = nf_get_ramified_primes(nf);
385 14 : l = lg(P);
386 14 : Q = q; q = gen_1;
387 70 : for (i = 1; i < l; i++)
388 : {
389 56 : GEN p = gel(P,i);
390 56 : if (Z_pvalrem(Q, p, &Q) && !BPSW_psp(p)) q = mulii(q, p);
391 : }
392 14 : return q;
393 : }
394 :
395 : /* lt(A) is an integer; ensure it is not a constant t_POL. In place */
396 : static void
397 7651 : ensure_lt_INT(GEN A)
398 : {
399 7651 : long n = lg(A)-1;
400 7651 : GEN lt = gel(A,n);
401 8484 : while (typ(lt) != t_INT) gel(A,n) = lt = gel(lt,2);
402 7651 : }
403 :
404 : /* set B = A/gcd(A,A'), squarefree */
405 : static GEN
406 7637 : get_nfsqff_data(GEN *pnf, GEN *pT, GEN *pA, GEN *pB, GEN *ptbad)
407 : {
408 7637 : GEN den, bad, D, B, A = *pA, T = *pT;
409 7637 : long n = degpol(T);
410 :
411 7637 : A = Q_primpart( QXQX_normalize(A, T) );
412 7637 : if (nfsqff_use_Trager(n, degpol(A)))
413 : {
414 154 : *pnf = T;
415 154 : bad = den = absi_shallow(ZX_disc(T));
416 154 : if (is_pm1(leading_coeff(T))) den = indexpartial(T, den);
417 : }
418 : else
419 : {
420 7483 : den = fix_nf(pnf, &T, &A);
421 7483 : bad = nf_get_index(*pnf);
422 7483 : if (den != gen_1) bad = mulii(bad, den);
423 : }
424 7637 : D = nfgcd_all(A, RgX_deriv(A), T, bad, &B);
425 7637 : if (degpol(D)) B = Q_primpart( QXQX_normalize(B, T) );
426 7637 : if (ptbad) *ptbad = bad;
427 7637 : *pA = A;
428 7637 : *pB = B; ensure_lt_INT(B);
429 7637 : *pT = T; return den;
430 : }
431 :
432 : /* return the roots of pol in nf */
433 : GEN
434 7756 : nfroots(GEN nf,GEN pol)
435 : {
436 7756 : pari_sp av = avma;
437 : GEN z, A, B, T, den;
438 : long d, dT;
439 :
440 7756 : if (!nf) return nfrootsQ(pol);
441 5334 : T = get_nfpol(nf, &nf);
442 5334 : RgX_check_ZX(T,"nfroots");
443 5334 : A = RgX_nffix("nfroots", T,pol,1);
444 5334 : d = degpol(A);
445 5334 : if (d < 0) pari_err_ROOTS0("nfroots");
446 5334 : if (d == 0) return cgetg(1,t_VEC);
447 5334 : if (d == 1)
448 : {
449 7 : A = QXQX_normalize(A,T);
450 7 : A = mkpolmod(gneg_i(gel(A,2)), T);
451 7 : return gerepilecopy(av, mkvec(A));
452 : }
453 5327 : dT = degpol(T);
454 5327 : if (dT == 1) return gerepileupto(av, nfrootsQ(simplify_shallow(A)));
455 :
456 5124 : den = get_nfsqff_data(&nf, &T, &A, &B, NULL);
457 5124 : if (RgX_is_ZX(B))
458 : {
459 1764 : GEN v = gel(ZX_factor(B), 1);
460 1764 : long i, l = lg(v), p = mael(factoru(dT),1,1); /* smallest prime divisor */
461 1764 : z = cgetg(1, t_VEC);
462 4690 : for (i = 1; i < l; i++)
463 : {
464 2926 : GEN b = gel(v,i); /* irreducible / Q */
465 2926 : long db = degpol(b);
466 2926 : if (db != 1 && degpol(b) < p) continue;
467 2926 : z = shallowconcat(z, nfsqff(nf, b, ROOTS, den));
468 : }
469 : }
470 : else
471 3360 : z = nfsqff(nf,B, ROOTS, den);
472 5124 : z = gerepileupto(av, QXQV_to_mod(z, T));
473 5124 : gen_sort_inplace(z, (void*)&cmp_RgX, &cmp_nodata, NULL);
474 5124 : return z;
475 : }
476 :
477 : static GEN
478 232165 : _norml2(GEN x) { return RgC_fpnorml2(x, DEFAULTPREC); }
479 :
480 : /* return a minimal lift of elt modulo id, as a ZC */
481 : static GEN
482 66002 : nf_bestlift(GEN elt, GEN bound, nflift_t *L)
483 : {
484 : GEN u;
485 66002 : long i,l = lg(L->prk), t = typ(elt);
486 66002 : if (t != t_INT)
487 : {
488 14736 : if (t == t_POL) elt = ZM_ZX_mul(L->tozk, elt);
489 14736 : u = ZM_ZC_mul(L->iprk,elt);
490 303841 : for (i=1; i<l; i++) gel(u,i) = diviiround(gel(u,i), L->pk);
491 : }
492 : else
493 : {
494 51266 : u = ZC_Z_mul(gel(L->iprk,1), elt);
495 556318 : for (i=1; i<l; i++) gel(u,i) = diviiround(gel(u,i), L->pk);
496 51267 : elt = scalarcol(elt, l-1);
497 : }
498 66002 : u = ZC_sub(elt, ZM_ZC_mul(L->prk, u));
499 66002 : if (bound && gcmp(_norml2(u), bound) > 0) u = NULL;
500 66002 : return u;
501 : }
502 :
503 : /* Warning: return L->topowden * (best lift). */
504 : static GEN
505 40313 : nf_bestlift_to_pol(GEN elt, GEN bound, nflift_t *L)
506 : {
507 40313 : pari_sp av = avma;
508 40313 : GEN u,v = nf_bestlift(elt,bound,L);
509 40313 : if (!v) return NULL;
510 34664 : if (ZV_isscalar(v))
511 : {
512 19194 : if (L->topowden)
513 19194 : u = mulii(L->topowden, gel(v,1));
514 : else
515 0 : u = icopy(gel(v,1));
516 19194 : u = gerepileuptoint(av, u);
517 : }
518 : else
519 : {
520 15470 : v = gclone(v); set_avma(av);
521 15470 : u = RgV_dotproduct(L->topow, v);
522 15470 : gunclone(v);
523 : }
524 34664 : return u;
525 : }
526 :
527 : /* return the T->powden * (lift of pol with coefficients of T2-norm <= C)
528 : * if it exists. */
529 : static GEN
530 11095 : nf_pol_lift(GEN pol, GEN bound, nflift_t *L)
531 : {
532 11095 : long i, l = lg(pol);
533 11095 : GEN x = cgetg(l,t_POL);
534 :
535 11095 : x[1] = pol[1];
536 11095 : gel(x,l-1) = mul_content(gel(pol,l-1), L->topowden);
537 40565 : for (i=l-2; i>1; i--)
538 : {
539 35119 : GEN t = nf_bestlift_to_pol(gel(pol,i), bound, L);
540 35119 : if (!t) return NULL;
541 29470 : gel(x,i) = t;
542 : }
543 5446 : return x;
544 : }
545 :
546 : static GEN
547 0 : zerofact(long v)
548 : {
549 0 : GEN z = cgetg(3, t_MAT);
550 0 : gel(z,1) = mkcol(pol_0(v));
551 0 : gel(z,2) = mkcol(gen_1); return z;
552 : }
553 :
554 : /* Return the factorization of A in Q[X]/(T) in rep [pre-allocated with
555 : * cgetg(3,t_MAT)], reclaiming all memory between avma and rep.
556 : * y is the vector of irreducible factors of B = Q_primpart( A/gcd(A,A') ).
557 : * Bad primes divide 'bad' */
558 : static void
559 2527 : fact_from_sqff(GEN rep, GEN A, GEN B, GEN y, GEN T, GEN bad)
560 : {
561 2527 : pari_sp av = (pari_sp)rep;
562 2527 : long n = lg(y)-1;
563 : GEN ex;
564 :
565 2527 : if (A != B)
566 : { /* not squarefree */
567 98 : if (n == 1)
568 : { /* perfect power, simple ! */
569 14 : long e = degpol(A) / degpol(gel(y,1));
570 14 : y = gerepileupto(av, QXQXV_to_mod(y, T));
571 14 : ex = mkcol(utoipos(e));
572 : }
573 : else
574 : { /* compute valuations mod a prime of degree 1 (avoid coeff explosion) */
575 84 : GEN quo, p, r, Bp, lb = leading_coeff(B), E = cgetalloc(t_VECSMALL,n+1);
576 84 : pari_sp av1 = avma;
577 : ulong pp;
578 : long j;
579 : forprime_t S;
580 84 : u_forprime_init(&S, degpol(T), ULONG_MAX);
581 224 : for (; ; set_avma(av1))
582 : {
583 308 : pp = u_forprime_next(&S);
584 308 : if (! umodiu(bad,pp) || !umodiu(lb, pp)) continue;
585 287 : p = utoipos(pp);
586 287 : r = FpX_oneroot(T, p);
587 287 : if (!r) continue;
588 147 : Bp = FpXY_evalx(B, r, p);
589 147 : if (FpX_is_squarefree(Bp, p)) break;
590 : }
591 :
592 84 : quo = FpXY_evalx(Q_primpart(A), r, p);
593 182 : for (j=n; j>=2; j--)
594 : {
595 98 : GEN junk, fact = Q_remove_denom(gel(y,j), &junk);
596 98 : long e = 0;
597 98 : fact = FpXY_evalx(fact, r, p);
598 217 : for(;; e++)
599 217 : {
600 315 : GEN q = FpX_divrem(quo,fact,p,ONLY_DIVIDES);
601 315 : if (!q) break;
602 217 : quo = q;
603 : }
604 98 : E[j] = e;
605 : }
606 84 : E[1] = degpol(quo) / degpol(gel(y,1));
607 84 : y = gerepileupto(av, QXQXV_to_mod(y, T));
608 84 : ex = zc_to_ZC(E); pari_free((void*)E);
609 : }
610 : }
611 : else
612 : {
613 2429 : y = gerepileupto(av, QXQXV_to_mod(y, T));
614 2429 : ex = const_col(n, gen_1);
615 : }
616 2527 : gel(rep,1) = y; settyp(y, t_COL);
617 2527 : gel(rep,2) = ex;
618 2527 : }
619 :
620 : /* return the factorization of polynomial pol in nf */
621 : static GEN
622 2709 : nffactor_i(GEN nf,GEN T,GEN pol)
623 : {
624 2709 : GEN bad, A, B, y, den, rep = cgetg(3, t_MAT);
625 2709 : pari_sp av = avma;
626 : long dA;
627 : pari_timer ti;
628 :
629 2709 : if (DEBUGLEVEL>2) { timer_start(&ti); err_printf("\nEntering nffactor:\n"); }
630 2709 : A = RgX_nffix("nffactor",T,pol,1);
631 2709 : dA = degpol(A);
632 2709 : if (dA <= 0) {
633 7 : set_avma((pari_sp)(rep + 3));
634 7 : return (dA == 0)? trivial_fact(): zerofact(varn(pol));
635 : }
636 2702 : if (dA == 1) {
637 : GEN c;
638 119 : A = Q_primpart( QXQX_normalize(A, T) );
639 119 : A = gerepilecopy(av, A); c = gel(A,2);
640 119 : if (typ(c) == t_POL && degpol(c) > 0) gel(A,2) = mkpolmod(c, ZX_copy(T));
641 119 : gel(rep,1) = mkcol(A);
642 119 : gel(rep,2) = mkcol(gen_1); return rep;
643 : }
644 2583 : if (degpol(T) == 1) return gerepileupto(av, QX_factor(simplify_shallow(A)));
645 :
646 2513 : den = get_nfsqff_data(&nf, &T, &A, &B, &bad);
647 2513 : if (DEBUGLEVEL>2) timer_printf(&ti, "squarefree test");
648 2513 : if (RgX_is_ZX(B))
649 : {
650 2107 : GEN v = gel(ZX_factor(B), 1);
651 2107 : long i, l = lg(v);
652 2107 : y = cgetg(1, t_VEC);
653 4263 : for (i = 1; i < l; i++)
654 : {
655 2156 : GEN b = gel(v,i); /* irreducible / Q */
656 2156 : y = shallowconcat(y, nfsqff(nf, b, 0, den));
657 : }
658 : }
659 : else
660 406 : y = nfsqff(nf,B, 0, den);
661 2513 : if (DEBUGLEVEL>3) err_printf("number of factor(s) found: %ld\n", lg(y)-1);
662 :
663 2513 : fact_from_sqff(rep, A, B, y, T, bad);
664 2513 : return rep;
665 : }
666 :
667 : /* return the factorization of P in nf */
668 : GEN
669 2702 : nffactor(GEN nf, GEN P)
670 : {
671 2702 : GEN y, T = get_nfpol(nf, &nf);
672 2702 : if (!nf) RgX_check_ZX(T,"nffactor");
673 2702 : if (typ(P) == t_RFRAC)
674 : {
675 14 : pari_sp av = avma;
676 14 : GEN a = gel(P, 1), b = gel(P, 2);
677 14 : y = famat_inv_shallow(nffactor_i(nf, T, b));
678 14 : if (typ(a) == t_POL && varn(a) == varn(b))
679 7 : y = famat_mul_shallow(nffactor_i(nf, T, a), y);
680 14 : y = gerepilecopy(av, y);
681 : }
682 : else
683 2688 : y = nffactor_i(nf, T, P);
684 2702 : return sort_factor_pol(y, cmp_RgX);
685 : }
686 :
687 : /* assume x scalar or t_COL, G t_MAT */
688 : static GEN
689 43384 : arch_for_T2(GEN G, GEN x)
690 : {
691 1806 : return (typ(x) == t_COL)? RgM_RgC_mul(G,x)
692 45190 : : RgC_Rg_mul(gel(G,1),x);
693 : }
694 :
695 : /* polbase a zkX with t_INT leading coeff; return a bound for T_2(P),
696 : * P | polbase in C[X]. NB: Mignotte bound: A | S ==>
697 : * |a_i| <= binom(d-1, i-1) || S ||_2 + binom(d-1, i) lc(S)
698 : *
699 : * Apply to sigma(S) for all embeddings sigma, then take the L_2 norm over
700 : * sigma, then take the sup over i */
701 : static GEN
702 2240 : nf_Mignotte_bound(GEN nf, GEN polbase)
703 2240 : { GEN lS = leading_coeff(polbase); /* t_INT */
704 : GEN p1, C, N2, binlS, bin;
705 2240 : long prec = nf_get_prec(nf), n = nf_get_degree(nf), r1 = nf_get_r1(nf);
706 2240 : long i, j, d = degpol(polbase);
707 :
708 2240 : binlS = bin = vecbinomial(d-1);
709 2240 : if (!isint1(lS)) binlS = ZC_Z_mul(bin,lS);
710 :
711 2240 : N2 = cgetg(n+1, t_VEC);
712 : for (;;)
713 0 : {
714 2240 : GEN G = nf_get_G(nf), matGS = cgetg(d+2, t_MAT);
715 :
716 27887 : for (j=0; j<=d; j++) gel(matGS,j+1) = arch_for_T2(G, gel(polbase,j+2));
717 2240 : matGS = shallowtrans(matGS);
718 4879 : for (j=1; j <= r1; j++) /* N2[j] = || sigma_j(S) ||_2 */
719 : {
720 2639 : GEN c = sqrtr( _norml2(gel(matGS,j)) );
721 2639 : gel(N2,j) = c; if (!signe(c)) goto PRECPB;
722 : }
723 6363 : for ( ; j <= n; j+=2)
724 : {
725 4123 : GEN q1 = _norml2(gel(matGS, j));
726 4123 : GEN q2 = _norml2(gel(matGS, j+1));
727 4123 : GEN c = sqrtr( gmul2n(addrr(q1, q2), -1) );
728 4123 : gel(N2,j) = gel(N2,j+1) = c; if (!signe(c)) goto PRECPB;
729 : }
730 2240 : break; /* done */
731 0 : PRECPB:
732 0 : prec = precdbl(prec);
733 0 : nf = nfnewprec_shallow(nf, prec);
734 0 : if (DEBUGLEVEL>1) pari_warn(warnprec, "nf_factor_bound", prec);
735 : }
736 :
737 : /* Take sup over 0 <= i <= d of
738 : * sum_j | binom(d-1, i-1) ||sigma_j(S)||_2 + binom(d-1,i) lc(S) |^2 */
739 :
740 : /* i = 0: n lc(S)^2 */
741 2240 : C = mului(n, sqri(lS));
742 : /* i = d: sum_sigma ||sigma(S)||_2^2 */
743 2240 : p1 = gnorml2(N2); if (gcmp(C, p1) < 0) C = p1;
744 23408 : for (i = 1; i < d; i++)
745 : {
746 21168 : GEN B = gel(bin,i), L = gel(binlS,i+1);
747 21168 : GEN s = sqrr(addri(mulir(B, gel(N2,1)), L)); /* j=1 */
748 194995 : for (j = 2; j <= n; j++) s = addrr(s, sqrr(addri(mulir(B, gel(N2,j)), L)));
749 21167 : if (mpcmp(C, s) < 0) C = s;
750 : }
751 2240 : return C;
752 : }
753 :
754 : /* return a bound for T_2(P), P | polbase
755 : * max |b_i|^2 <= 3^{3/2 + d} / (4 \pi d) [P]_2,
756 : * where [P]_2 is Bombieri's 2-norm
757 : * Sum over conjugates */
758 : static GEN
759 2240 : nf_Beauzamy_bound(GEN nf, GEN polbase)
760 : {
761 : GEN lt, C, s, POL, bin;
762 2240 : long d = degpol(polbase), n = nf_get_degree(nf), prec = nf_get_prec(nf);
763 2240 : bin = vecbinomial(d);
764 2240 : POL = polbase + 2;
765 : /* compute [POL]_2 */
766 : for (;;)
767 0 : {
768 2240 : GEN G = nf_get_G(nf);
769 : long i;
770 :
771 2240 : s = real_0(prec);
772 27884 : for (i=0; i<=d; i++)
773 : {
774 25644 : GEN c = gel(POL,i);
775 25644 : if (gequal0(c)) continue;
776 17736 : c = _norml2(arch_for_T2(G,c));
777 17738 : if (!signe(c)) goto PRECPB;
778 : /* s += T2(POL[i]) / binomial(d,i) */
779 17738 : s = addrr(s, divri(c, gel(bin,i+1)));
780 : }
781 2240 : break;
782 0 : PRECPB:
783 0 : prec = precdbl(prec);
784 0 : nf = nfnewprec_shallow(nf, prec);
785 0 : if (DEBUGLEVEL>1) pari_warn(warnprec, "nf_factor_bound", prec);
786 : }
787 2240 : lt = leading_coeff(polbase);
788 2240 : s = mulri(s, muliu(sqri(lt), n));
789 2240 : C = powruhalf(utor(3,DEFAULTPREC), 3 + 2*d); /* 3^{3/2 + d} */
790 2240 : return divrr(mulrr(C, s), mulur(d, mppi(DEFAULTPREC)));
791 : }
792 :
793 : static GEN
794 2240 : nf_factor_bound(GEN nf, GEN polbase)
795 : {
796 2240 : pari_sp av = avma;
797 2240 : GEN a = nf_Mignotte_bound(nf, polbase);
798 2240 : GEN b = nf_Beauzamy_bound(nf, polbase);
799 2240 : if (DEBUGLEVEL>2)
800 : {
801 0 : err_printf("Mignotte bound: %Ps\n",a);
802 0 : err_printf("Beauzamy bound: %Ps\n",b);
803 : }
804 2240 : return gerepileupto(av, gmin(a, b));
805 : }
806 :
807 : /* True nf; return Bs: if r a root of sigma_i(P), |r| < Bs[i] */
808 : static GEN
809 4669 : nf_root_bounds(GEN nf, GEN P)
810 : {
811 : long lR, i, j, l, prec, r1;
812 : GEN Ps, R, V;
813 :
814 4669 : if (RgX_is_rational(P)) return polrootsbound(P, NULL);
815 2289 : r1 = nf_get_r1(nf);
816 2289 : P = Q_primpart(P);
817 2289 : prec = ZXX_max_lg(P) + 1;
818 2289 : l = lg(P);
819 2289 : if (nf_get_prec(nf) >= prec)
820 1966 : R = nf_get_roots(nf);
821 : else
822 323 : R = QX_complex_roots(nf_get_pol(nf), prec);
823 2289 : lR = lg(R);
824 2289 : V = cgetg(lR, t_VEC);
825 2289 : Ps = cgetg(l, t_POL); /* sigma (P) */
826 2289 : Ps[1] = P[1];
827 6692 : for (j=1; j<lg(R); j++)
828 : {
829 4403 : GEN r = gel(R,j);
830 26985 : for (i=2; i<l; i++) gel(Ps,i) = poleval(gel(P,i), r);
831 4403 : gel(V,j) = polrootsbound(Ps, NULL);
832 : }
833 2289 : return mkvec2(vecslice(V,1,r1), vecslice(V,r1+1,lg(V)-1));
834 : }
835 :
836 : /* return B such that, if x = sum x_i K.zk[i] in O_K, then ||x||_2^2 <= B T_2(x)
837 : * den = multiplicative bound for denom(x) [usually NULL, for 1, but when we
838 : * use nf_PARTIALFACT K.zk may not generate O_K] */
839 : static GEN
840 4669 : L2_bound(GEN nf, GEN den)
841 : {
842 4669 : GEN M, L, prep, T = nf_get_pol(nf), tozk = nf_get_invzk(nf);
843 4669 : long prec = ZM_max_lg(tozk) + ZX_max_lg(T) + nbits2prec(degpol(T));
844 4669 : (void)initgaloisborne(nf, den? den: gen_1, prec, &L, &prep, NULL);
845 4669 : M = vandermondeinverse(L, RgX_gtofp(T,prec), den, prep);
846 4669 : return RgM_fpnorml2(RgM_mul(tozk,M), DEFAULTPREC);
847 : }
848 :
849 : /* sum_i L[i]^p */
850 : static GEN
851 8624 : normlp(GEN L, long p)
852 : {
853 8624 : long i, l = lg(L);
854 : GEN z;
855 8624 : if (l == 1) return gen_0;
856 4557 : z = gpowgs(gel(L,1), p);
857 8323 : for (i=2; i<l; i++) z = gadd(z, gpowgs(gel(L,i), p));
858 4557 : return z;
859 : }
860 : /* \sum_i deg(sigma_i) L[i]^p in dimension n (L may be a scalar
861 : * or [L1,L2], where Ld corresponds to the archimedean places of degree d) */
862 : static GEN
863 7413 : normTp(GEN L, long p, long n)
864 : {
865 7413 : if (typ(L) != t_VEC) return gmulsg(n, gpowgs(L, p));
866 4312 : return gadd(normlp(gel(L,1),p), gmul2n(normlp(gel(L,2),p), 1));
867 : }
868 :
869 : /* S = S0 + tS1, P = P0 + tP1 (Euclidean div. by t integer). For a true
870 : * factor (vS, vP), we have:
871 : * | S vS + P vP |^2 < Btra
872 : * This implies | S1 vS + P1 vP |^2 < Bhigh, assuming t > sqrt(Btra).
873 : * d = dimension of low part (= [nf:Q])
874 : * n0 = bound for |vS|^2
875 : * */
876 : static double
877 2310 : get_Bhigh(long n0, long d)
878 : {
879 2310 : double sqrtd = sqrt((double)d);
880 2310 : double z = n0*sqrtd + sqrtd/2 * (d * (n0+1));
881 2310 : z = 1. + 0.5 * z; return z * z;
882 : }
883 :
884 : typedef struct {
885 : GEN d;
886 : GEN dPinvS; /* d P^(-1) S [ integral ] */
887 : double **PinvSdbl; /* P^(-1) S as double */
888 : GEN S1, P1; /* S = S0 + S1 q, idem P */
889 : } trace_data;
890 :
891 : /* S1 * u - P1 * round(P^-1 S u). K non-zero coords in u given by ind */
892 : static GEN
893 168427 : get_trace(GEN ind, trace_data *T)
894 : {
895 168427 : long i, j, l, K = lg(ind)-1;
896 : GEN z, s, v;
897 :
898 168427 : s = gel(T->S1, ind[1]);
899 168427 : if (K == 1) return s;
900 :
901 : /* compute s = S1 u */
902 469217 : for (j=2; j<=K; j++) s = ZC_add(s, gel(T->S1, ind[j]));
903 :
904 : /* compute v := - round(P^1 S u) */
905 163485 : l = lg(s);
906 163485 : v = cgetg(l, t_VECSMALL);
907 2686404 : for (i=1; i<l; i++)
908 : {
909 2522919 : double r, t = 0.;
910 : /* quick approximate computation */
911 9741067 : for (j=1; j<=K; j++) t += T->PinvSdbl[ ind[j] ][i];
912 2522919 : r = floor(t + 0.5);
913 2522919 : if (fabs(t + 0.5 - r) < 0.0001)
914 : { /* dubious, compute exactly */
915 266 : z = gen_0;
916 1036 : for (j=1; j<=K; j++) z = addii(z, ((GEN**)T->dPinvS)[ ind[j] ][i]);
917 266 : v[i] = - itos( diviiround(z, T->d) );
918 : }
919 : else
920 2522653 : v[i] = - (long)r;
921 : }
922 163485 : return ZC_add(s, ZM_zc_mul(T->P1, v));
923 : }
924 :
925 : static trace_data *
926 4480 : init_trace(trace_data *T, GEN S, nflift_t *L, GEN q)
927 : {
928 4480 : long e = gexpo(S), i,j, l,h;
929 : GEN qgood, S1, invd;
930 :
931 4480 : if (e < 0) return NULL; /* S = 0 */
932 :
933 4109 : qgood = int2n(e - 32); /* single precision check */
934 4109 : if (cmpii(qgood, q) > 0) q = qgood;
935 :
936 4109 : S1 = gdivround(S, q);
937 4109 : if (gequal0(S1)) return NULL;
938 :
939 679 : invd = invr(itor(L->pk, DEFAULTPREC));
940 :
941 679 : T->dPinvS = ZM_mul(L->iprk, S);
942 679 : l = lg(S);
943 679 : h = lgcols(T->dPinvS);
944 679 : T->PinvSdbl = (double**)cgetg(l, t_MAT);
945 7714 : for (j = 1; j < l; j++)
946 : {
947 7035 : double *t = (double *) stack_malloc_align(h * sizeof(double), sizeof(double));
948 7035 : GEN c = gel(T->dPinvS,j);
949 7035 : pari_sp av = avma;
950 7035 : T->PinvSdbl[j] = t;
951 116886 : for (i=1; i < h; i++) t[i] = rtodbl(mulri(invd, gel(c,i)));
952 7035 : set_avma(av);
953 : }
954 :
955 679 : T->d = L->pk;
956 679 : T->P1 = gdivround(L->prk, q);
957 679 : T->S1 = S1; return T;
958 : }
959 :
960 : static void
961 55594 : update_trace(trace_data *T, long k, long i)
962 : {
963 55594 : if (!T) return;
964 25501 : gel(T->S1,k) = gel(T->S1,i);
965 25501 : gel(T->dPinvS,k) = gel(T->dPinvS,i);
966 25501 : T->PinvSdbl[k] = T->PinvSdbl[i];
967 : }
968 :
969 : /* reduce coeffs mod (T,pk), then center mod pk */
970 : static GEN
971 20909 : FqX_centermod(GEN z, GEN T, GEN pk, GEN pks2)
972 : {
973 : long i, l;
974 : GEN y;
975 20909 : if (!T) return centermod_i(z, pk, pks2);
976 15092 : y = FpXQX_red(z, T, pk); l = lg(y);
977 137305 : for (i = 2; i < l; i++)
978 : {
979 122213 : GEN c = gel(y,i);
980 122213 : if (typ(c) == t_INT)
981 81074 : c = Fp_center_i(c, pk, pks2);
982 : else
983 41139 : c = FpX_center_i(c, pk, pks2);
984 122213 : gel(y,i) = c;
985 : }
986 15092 : return y;
987 : }
988 :
989 : typedef struct {
990 : GEN lt, C, Clt, C2lt, C2ltpol;
991 : } div_data;
992 :
993 : static void
994 4732 : init_div_data(div_data *D, GEN pol, nflift_t *L)
995 : {
996 4732 : GEN C2lt, Clt, C = mul_content(L->topowden, L->dn);
997 4732 : GEN lc = leading_coeff(pol), lt = is_pm1(lc)? NULL: absi_shallow(lc);
998 4732 : if (C)
999 : {
1000 4732 : GEN C2 = sqri(C);
1001 4732 : if (lt) {
1002 1036 : C2lt = mulii(C2, lt);
1003 1036 : Clt = mulii(C,lt);
1004 : } else {
1005 3696 : C2lt = C2;
1006 3696 : Clt = C;
1007 : }
1008 : }
1009 : else
1010 0 : C2lt = Clt = lt;
1011 4732 : D->lt = lt;
1012 4732 : D->C = C;
1013 4732 : D->Clt = Clt;
1014 4732 : D->C2lt = C2lt;
1015 4732 : D->C2ltpol = C2lt? RgX_Rg_mul(pol, C2lt): pol;
1016 4732 : }
1017 : static void
1018 5264 : update_target(div_data *D, GEN pol)
1019 5264 : { D->C2ltpol = D->Clt? RgX_Rg_mul(pol, D->Clt): pol; }
1020 :
1021 : /* nb = number of modular factors; return a "good" K such that naive
1022 : * recombination of up to maxK modular factors is not too costly */
1023 : long
1024 20860 : cmbf_maxK(long nb)
1025 : {
1026 20860 : if (nb > 10) return 3;
1027 19607 : return nb-1;
1028 : }
1029 : /* Naive recombination of modular factors: combine up to maxK modular
1030 : * factors, degree <= klim
1031 : *
1032 : * target = polynomial we want to factor
1033 : * famod = array of modular factors. Product should be congruent to
1034 : * target/lc(target) modulo p^a
1035 : * For true factors: S1,S2 <= p^b, with b <= a and p^(b-a) < 2^31 */
1036 : /* set *done = 1 if factorisation is known to be complete */
1037 : static GEN
1038 2240 : nfcmbf(nfcmbf_t *T, long klim, long *pmaxK, int *done)
1039 : {
1040 2240 : GEN nf = T->nf, famod = T->fact, bound = T->bound;
1041 2240 : GEN ltdn, nfpol = nf_get_pol(nf);
1042 2240 : long K = 1, cnt = 1, i,j,k, curdeg, lfamod = lg(famod)-1, dnf = degpol(nfpol);
1043 2240 : pari_sp av0 = avma;
1044 2240 : GEN Tpk = T->L->Tpk, pk = T->L->pk, pks2 = shifti(pk,-1);
1045 2240 : GEN ind = cgetg(lfamod+1, t_VECSMALL);
1046 2240 : GEN deg = cgetg(lfamod+1, t_VECSMALL);
1047 2240 : GEN degsofar = cgetg(lfamod+1, t_VECSMALL);
1048 2240 : GEN fa = cgetg(lfamod+1, t_VEC);
1049 2240 : const double Bhigh = get_Bhigh(lfamod, dnf);
1050 : trace_data _T1, _T2, *T1, *T2;
1051 : div_data D;
1052 : pari_timer ti;
1053 :
1054 2240 : timer_start(&ti);
1055 :
1056 2240 : *pmaxK = cmbf_maxK(lfamod);
1057 2240 : init_div_data(&D, T->pol, T->L);
1058 2240 : ltdn = mul_content(D.lt, T->L->dn);
1059 : {
1060 2240 : GEN q = ceil_safe(sqrtr(T->BS_2));
1061 2240 : GEN t1,t2, lt2dn = mul_content(ltdn, D.lt);
1062 2240 : GEN trace1 = cgetg(lfamod+1, t_MAT);
1063 2240 : GEN trace2 = cgetg(lfamod+1, t_MAT);
1064 12362 : for (i=1; i <= lfamod; i++)
1065 : {
1066 10122 : pari_sp av = avma;
1067 10122 : GEN P = gel(famod,i);
1068 10122 : long d = degpol(P);
1069 :
1070 10122 : deg[i] = d; P += 2;
1071 10122 : t1 = gel(P,d-1);/* = - S_1 */
1072 10122 : t2 = Fq_sqr(t1, Tpk, pk);
1073 10121 : if (d > 1) t2 = Fq_sub(t2, gmul2n(gel(P,d-2), 1), Tpk, pk);
1074 : /* t2 = S_2 Newton sum */
1075 10121 : if (ltdn)
1076 : {
1077 294 : t1 = Fq_Fp_mul(t1, ltdn, Tpk, pk);
1078 294 : t2 = Fq_Fp_mul(t2, lt2dn, Tpk, pk);
1079 : }
1080 10121 : gel(trace1,i) = gclone( nf_bestlift(t1, NULL, T->L) );
1081 10122 : gel(trace2,i) = gclone( nf_bestlift(t2, NULL, T->L) ); set_avma(av);
1082 : }
1083 2240 : T1 = init_trace(&_T1, trace1, T->L, q);
1084 2240 : T2 = init_trace(&_T2, trace2, T->L, q);
1085 12362 : for (i=1; i <= lfamod; i++) {
1086 10122 : gunclone(gel(trace1,i));
1087 10122 : gunclone(gel(trace2,i));
1088 : }
1089 : }
1090 2240 : degsofar[0] = 0; /* sentinel */
1091 :
1092 : /* ind runs through strictly increasing sequences of length K,
1093 : * 1 <= ind[i] <= lfamod */
1094 3010 : nextK:
1095 3010 : if (K > *pmaxK || 2*K > lfamod) goto END;
1096 2667 : if (DEBUGLEVEL > 3)
1097 0 : err_printf("\n### K = %d, %Ps combinations\n", K,binomial(utoipos(lfamod), K));
1098 2667 : setlg(ind, K+1); ind[1] = 1;
1099 2667 : i = 1; curdeg = deg[ind[1]];
1100 : for(;;)
1101 : { /* try all combinations of K factors */
1102 197393 : for (j = i; j < K; j++)
1103 : {
1104 21539 : degsofar[j] = curdeg;
1105 21539 : ind[j+1] = ind[j]+1; curdeg += deg[ind[j+1]];
1106 : }
1107 175854 : if (curdeg <= klim) /* trial divide */
1108 : {
1109 : GEN t, y, q;
1110 : pari_sp av;
1111 :
1112 175854 : av = avma;
1113 175854 : if (T1)
1114 : { /* d-1 test */
1115 84308 : t = get_trace(ind, T1);
1116 84308 : if (rtodbl(_norml2(t)) > Bhigh)
1117 : {
1118 82341 : if (DEBUGLEVEL>6) err_printf(".");
1119 82341 : set_avma(av); goto NEXT;
1120 : }
1121 : }
1122 93513 : if (T2)
1123 : { /* d-2 test */
1124 84119 : t = get_trace(ind, T2);
1125 84119 : if (rtodbl(_norml2(t)) > Bhigh)
1126 : {
1127 82558 : if (DEBUGLEVEL>3) err_printf("|");
1128 82558 : set_avma(av); goto NEXT;
1129 : }
1130 : }
1131 10955 : set_avma(av);
1132 10955 : y = ltdn; /* full computation */
1133 31864 : for (i=1; i<=K; i++)
1134 : {
1135 20909 : GEN q = gel(famod, ind[i]);
1136 20909 : if (y) q = gmul(y, q);
1137 20909 : y = FqX_centermod(q, Tpk, pk, pks2);
1138 : }
1139 10955 : y = nf_pol_lift(y, bound, T->L);
1140 10955 : if (!y)
1141 : {
1142 5635 : if (DEBUGLEVEL>3) err_printf("@");
1143 5635 : set_avma(av); goto NEXT;
1144 : }
1145 : /* y = topowden*dn*lt*\prod_{i in ind} famod[i] is apparently in O_K[X],
1146 : * in fact in (Z[Y]/nf.pol)[X] due to multiplication by C = topowden*dn.
1147 : * Try out this candidate factor */
1148 5320 : q = RgXQX_divrem(D.C2ltpol, y, nfpol, ONLY_DIVIDES);
1149 5320 : if (!q)
1150 : {
1151 133 : if (DEBUGLEVEL>3) err_printf("*");
1152 133 : set_avma(av); goto NEXT;
1153 : }
1154 : /* Original T->pol in O_K[X] with leading coeff lt in Z,
1155 : * y = C*lt \prod famod[i] is in O_K[X] with leading coeff in Z
1156 : * q = C^2*lt*pol / y = C * (lt*pol) / (lt*\prod famod[i]) is a
1157 : * K-rational factor, in fact in Z[Y]/nf.pol)[X] as above, with
1158 : * leading term C*lt. */
1159 5187 : update_target(&D, q);
1160 5187 : gel(fa,cnt++) = D.C2lt? RgX_int_normalize(y): y; /* make monic */
1161 38843 : for (i=j=k=1; i <= lfamod; i++)
1162 : { /* remove used factors */
1163 33656 : if (j <= K && i == ind[j]) j++;
1164 : else
1165 : {
1166 27797 : gel(famod,k) = gel(famod,i);
1167 27797 : update_trace(T1, k, i);
1168 27797 : update_trace(T2, k, i);
1169 27797 : deg[k] = deg[i]; k++;
1170 : }
1171 : }
1172 5187 : lfamod -= K;
1173 5187 : *pmaxK = cmbf_maxK(lfamod);
1174 5187 : if (lfamod < 2*K) goto END;
1175 3290 : i = 1; curdeg = deg[ind[1]];
1176 3290 : if (DEBUGLEVEL > 2)
1177 : {
1178 0 : err_printf("\n"); timer_printf(&ti, "to find factor %Ps",gel(fa,cnt-1));
1179 0 : err_printf("remaining modular factor(s): %ld\n", lfamod);
1180 : }
1181 3290 : continue;
1182 : }
1183 :
1184 0 : NEXT:
1185 170667 : for (i = K+1;;)
1186 : {
1187 192304 : if (--i == 0) { K++; goto nextK; }
1188 191534 : if (++ind[i] <= lfamod - K + i)
1189 : {
1190 169897 : curdeg = degsofar[i-1] + deg[ind[i]];
1191 169897 : if (curdeg <= klim) break;
1192 : }
1193 : }
1194 : }
1195 2240 : END:
1196 2240 : *done = 1;
1197 2240 : if (degpol(D.C2ltpol) > 0)
1198 : { /* leftover factor */
1199 2240 : GEN q = D.C2ltpol;
1200 2240 : if (D.C2lt) q = RgX_int_normalize(q);
1201 2240 : if (lfamod >= 2*K)
1202 : { /* restore leading coefficient [#930] */
1203 70 : if (D.lt) q = RgX_Rg_mul(q, D.lt);
1204 70 : *done = 0; /* ... may still be reducible */
1205 : }
1206 2240 : setlg(famod, lfamod+1);
1207 2240 : gel(fa,cnt++) = q;
1208 : }
1209 2240 : if (DEBUGLEVEL>6) err_printf("\n");
1210 2240 : setlg(fa, cnt);
1211 2240 : return gerepilecopy(av0, fa);
1212 : }
1213 :
1214 : static GEN
1215 77 : nf_chk_factors(nfcmbf_t *T, GEN P, GEN M_L, GEN famod, GEN pk)
1216 : {
1217 77 : GEN nf = T->nf, bound = T->bound;
1218 77 : GEN nfT = nf_get_pol(nf);
1219 : long i, r;
1220 77 : GEN pol = P, list, piv, y;
1221 77 : GEN Tpk = T->L->Tpk;
1222 : div_data D;
1223 :
1224 77 : piv = ZM_hnf_knapsack(M_L);
1225 77 : if (!piv) return NULL;
1226 63 : if (DEBUGLEVEL>3) err_printf("ZM_hnf_knapsack output:\n%Ps\n",piv);
1227 :
1228 63 : r = lg(piv)-1;
1229 63 : list = cgetg(r+1, t_VEC);
1230 63 : init_div_data(&D, pol, T->L);
1231 63 : for (i = 1;;)
1232 77 : {
1233 140 : pari_sp av = avma;
1234 140 : if (DEBUGLEVEL) err_printf("nf_LLL_cmbf: checking factor %ld\n", i);
1235 140 : y = chk_factors_get(D.lt, famod, gel(piv,i), Tpk, pk);
1236 :
1237 140 : if (! (y = nf_pol_lift(y, bound, T->L)) ) return NULL;
1238 126 : y = gerepilecopy(av, y);
1239 : /* y is the candidate factor */
1240 126 : pol = RgXQX_divrem(D.C2ltpol, y, nfT, ONLY_DIVIDES);
1241 126 : if (!pol) return NULL;
1242 :
1243 126 : if (D.C2lt) y = RgX_int_normalize(y);
1244 126 : gel(list,i) = y;
1245 126 : if (++i >= r) break;
1246 :
1247 77 : update_target(&D, pol);
1248 : }
1249 49 : gel(list,i) = RgX_int_normalize(pol); return list;
1250 : }
1251 :
1252 : static GEN
1253 38129 : nf_to_Zq(GEN x, GEN T, GEN pk, GEN pks2, GEN proj)
1254 : {
1255 : GEN y;
1256 38129 : if (typ(x) != t_COL) return centermodii(x, pk, pks2);
1257 6258 : if (!T)
1258 : {
1259 6062 : y = ZV_dotproduct(proj, x);
1260 6062 : return centermodii(y, pk, pks2);
1261 : }
1262 196 : y = ZM_ZC_mul(proj, x);
1263 196 : y = RgV_to_RgX(y, varn(T));
1264 196 : return FpX_center_i(FpX_rem(y, T, pk), pk, pks2);
1265 : }
1266 :
1267 : /* Assume P in nfX form, lc(P) != 0 mod p. Reduce P to Zp[X]/(T) mod p^a */
1268 : static GEN
1269 4676 : ZqX(GEN P, GEN pk, GEN T, GEN proj)
1270 : {
1271 4676 : long i, l = lg(P);
1272 4676 : GEN z, pks2 = shifti(pk,-1);
1273 :
1274 4676 : z = cgetg(l,t_POL); z[1] = P[1];
1275 42805 : for (i=2; i<l; i++) gel(z,i) = nf_to_Zq(gel(P,i),T,pk,pks2,proj);
1276 4676 : return normalizepol_lg(z, l);
1277 : }
1278 :
1279 : static GEN
1280 4676 : ZqX_normalize(GEN P, GEN lt, nflift_t *L)
1281 : {
1282 4676 : GEN R = lt? RgX_Rg_mul(P, Fp_inv(lt, L->pk)): P;
1283 4676 : return ZqX(R, L->pk, L->Tpk, L->ZqProj);
1284 : }
1285 :
1286 : /* k allowing to reconstruct x, |x|^2 < C, from x mod pr^k */
1287 : /* return log [ 2sqrt(C/d) * ( (3/2)sqrt(gamma) )^(d-1) ] ^d / log N(pr)
1288 : * cf. Belabas relative van Hoeij algorithm, lemma 3.12 */
1289 : static double
1290 4676 : bestlift_bound(GEN C, long d, double alpha, GEN p, long f)
1291 : {
1292 4676 : const double g = 1 / (alpha - 0.25); /* = 2 if alpha = 3/4 */
1293 4676 : GEN C4 = shiftr(gtofp(C,DEFAULTPREC), 2);
1294 4676 : double t, logp = log(gtodouble(p));
1295 4676 : if (f == d)
1296 : { /* p inert, no LLL fudge factor: p^(2k) / 4 > C */
1297 28 : t = 0.5 * rtodbl(mplog(C4));
1298 28 : return ceil(t / logp);
1299 : }
1300 : /* (1/2)log (4C/d) + (d-1)(log 3/2 sqrt(gamma)) */
1301 4648 : t = 0.5 * rtodbl(mplog(divru(C4,d))) + (d-1) * log(1.5 * sqrt(g));
1302 4648 : return ceil((t * d) / (logp * f));
1303 : }
1304 :
1305 : static GEN
1306 4718 : get_R(GEN M)
1307 : {
1308 : GEN R;
1309 4718 : long i, l, prec = nbits2prec( gexpo(M) + 64 );
1310 :
1311 : for(;;)
1312 : {
1313 4718 : R = gaussred_from_QR(M, prec);
1314 4718 : if (R) break;
1315 0 : prec = precdbl(prec);
1316 : }
1317 4718 : l = lg(R);
1318 22918 : for (i=1; i<l; i++) gcoeff(R,i,i) = gen_1;
1319 4718 : return R;
1320 : }
1321 :
1322 : static void
1323 4676 : init_proj(nflift_t *L, GEN prkHNF, GEN nfT)
1324 : {
1325 4676 : if (degpol(L->Tp)>1)
1326 : {
1327 168 : GEN coTp = FpX_div(FpX_red(nfT, L->p), L->Tp, L->p); /* Tp's cofactor */
1328 : GEN z, proj;
1329 168 : z = ZpX_liftfact(nfT, mkvec2(L->Tp, coTp), L->pk, L->p, L->k);
1330 168 : L->Tpk = gel(z,1);
1331 168 : proj = QXQV_to_FpM(L->topow, L->Tpk, L->pk);
1332 168 : if (L->topowden)
1333 168 : proj = FpM_red(ZM_Z_mul(proj, Fp_inv(L->topowden, L->pk)), L->pk);
1334 168 : L->ZqProj = proj;
1335 : }
1336 : else
1337 : {
1338 4508 : L->Tpk = NULL;
1339 4508 : L->ZqProj = dim1proj(prkHNF);
1340 : }
1341 4676 : }
1342 :
1343 : /* Square of the radius of largest ball inscript in PRK's fundamental domain,
1344 : * whose orthogonalized vector's norms are the Bi
1345 : * Rmax ^2 = min 1/4T_i where T_i = sum_j ( s_ij^2 / B_j)
1346 : * For p inert, S = Id, T_i = 1 / p^{2k} and Rmax = p^k / 2 */
1347 : static GEN
1348 4718 : max_radius(GEN PRK, GEN B)
1349 : {
1350 4718 : GEN S, smax = gen_0;
1351 4718 : pari_sp av = avma;
1352 4718 : long i, j, d = lg(PRK)-1;
1353 :
1354 4718 : S = RgM_inv( get_R(PRK) ); if (!S) pari_err_PREC("max_radius");
1355 22918 : for (i=1; i<=d; i++)
1356 : {
1357 18200 : GEN s = gen_0;
1358 180362 : for (j=1; j<=d; j++)
1359 162162 : s = mpadd(s, mpdiv( mpsqr(gcoeff(S,i,j)), gel(B,j)));
1360 18200 : if (mpcmp(s, smax) > 0) smax = s;
1361 : }
1362 4718 : return gerepileupto(av, ginv(gmul2n(smax, 2)));
1363 : }
1364 :
1365 : static void
1366 4676 : bestlift_init(long a, GEN nf, GEN C, nflift_t *L)
1367 : {
1368 4676 : const double alpha = 0.99; /* LLL parameter */
1369 4676 : const long d = nf_get_degree(nf);
1370 4676 : pari_sp av = avma, av2;
1371 4676 : GEN prk, PRK, iPRK, GSmin, T = L->Tp, p = L->p;
1372 4676 : long f = degpol(T);
1373 : pari_timer ti;
1374 :
1375 4676 : if (f == d)
1376 : { /* inert p, much simpler */
1377 28 : long a0 = bestlift_bound(C, d, alpha, p, f);
1378 : GEN q;
1379 28 : if (a < a0) a = a0; /* guarantees GSmin >= C */
1380 28 : if (DEBUGLEVEL>2) err_printf("exponent %ld\n",a);
1381 28 : q = powiu(p,a);
1382 28 : PRK = prk = scalarmat_shallow(q, d);
1383 28 : GSmin = shiftr(itor(q, DEFAULTPREC), -1);
1384 28 : iPRK = matid(d); goto END;
1385 : }
1386 4648 : timer_start(&ti);
1387 4648 : if (!a) a = (long)bestlift_bound(C, d, alpha, p, f);
1388 70 : for (;; set_avma(av), a += (a==1)? 1: (a>>1)) /* roughly a *= 1.5 */
1389 70 : {
1390 4718 : GEN B, q = powiu(p,a), Tq = FpXQ_powu(T, a, FpX_red(nf_get_pol(nf), q), q);
1391 4718 : if (DEBUGLEVEL>2) err_printf("exponent %ld\n",a);
1392 4718 : prk = idealhnf_two(nf, mkvec2(q, Tq));
1393 4718 : av2 = avma;
1394 4718 : PRK = ZM_lll_norms(prk, alpha, LLL_INPLACE, &B);
1395 4718 : GSmin = max_radius(PRK, B);
1396 4718 : if (gcmp(GSmin, C) >= 0) break;
1397 : }
1398 4648 : gerepileall(av2, 2, &PRK, &GSmin);
1399 4648 : iPRK = ZM_inv(PRK, NULL);
1400 4648 : if (DEBUGLEVEL>2)
1401 0 : err_printf("for this exponent, GSmin = %Ps\nTime reduction: %ld\n",
1402 : GSmin, timer_delay(&ti));
1403 4648 : END:
1404 4676 : L->k = a;
1405 4676 : L->pk = gcoeff(prk,1,1);
1406 4676 : L->prk = PRK;
1407 4676 : L->iprk = iPRK;
1408 4676 : L->GSmin= GSmin;
1409 4676 : init_proj(L, prk, nf_get_pol(nf));
1410 4676 : }
1411 :
1412 : /* Let X = Tra * M_L, Y = bestlift(X) return V s.t Y = X - PRK V
1413 : * and set *eT2 = gexpo(Y) [cf nf_bestlift, but memory efficient] */
1414 : static GEN
1415 336 : get_V(GEN Tra, GEN M_L, GEN PRK, GEN PRKinv, GEN pk, long *eT2)
1416 : {
1417 336 : long i, e = 0, l = lg(M_L);
1418 336 : GEN V = cgetg(l, t_MAT);
1419 336 : *eT2 = 0;
1420 4529 : for (i = 1; i < l; i++)
1421 : { /* cf nf_bestlift(Tra * c) */
1422 4193 : pari_sp av = avma, av2;
1423 4193 : GEN v, T2 = ZM_ZC_mul(Tra, gel(M_L,i));
1424 :
1425 4193 : v = gdivround(ZM_ZC_mul(PRKinv, T2), pk); /* small */
1426 4193 : av2 = avma;
1427 4193 : T2 = ZC_sub(T2, ZM_ZC_mul(PRK, v));
1428 4193 : e = gexpo(T2); if (e > *eT2) *eT2 = e;
1429 4193 : set_avma(av2);
1430 4193 : gel(V,i) = gerepileupto(av, v); /* small */
1431 : }
1432 336 : return V;
1433 : }
1434 :
1435 : static GEN
1436 70 : nf_LLL_cmbf(nfcmbf_t *T, long rec)
1437 : {
1438 70 : const double BitPerFactor = 0.4; /* nb bits / modular factor */
1439 70 : nflift_t *L = T->L;
1440 70 : GEN famod = T->fact, ZC = T->ZC, Br = T->Br, P = T->pol, dn = T->L->dn;
1441 70 : long dnf = nf_get_degree(T->nf), dP = degpol(P);
1442 : long i, C, tmax, n0;
1443 : GEN lP, Bnorm, Tra, T2, TT, CM_L, m, list, ZERO, Btra;
1444 : double Bhigh;
1445 : pari_sp av, av2;
1446 70 : long ti_LLL = 0, ti_CF = 0;
1447 : pari_timer ti2, TI;
1448 :
1449 70 : lP = absi_shallow(leading_coeff(P));
1450 70 : if (is_pm1(lP)) lP = NULL;
1451 :
1452 70 : n0 = lg(famod) - 1;
1453 : /* Lattice: (S PRK), small vector (vS vP). To find k bound for the image,
1454 : * write S = S1 q + S0, P = P1 q + P0
1455 : * |S1 vS + P1 vP|^2 <= Bhigh for all (vS,vP) assoc. to true factors */
1456 70 : Btra = mulrr(ZC, mulur(dP*dP, normTp(Br, 2, dnf)));
1457 70 : Bhigh = get_Bhigh(n0, dnf);
1458 70 : C = (long)ceil(sqrt(Bhigh/n0)) + 1; /* C^2 n0 ~ Bhigh */
1459 70 : Bnorm = dbltor( n0 * C * C + Bhigh );
1460 70 : ZERO = zeromat(n0, dnf);
1461 :
1462 70 : av = avma;
1463 70 : TT = const_vec(n0, NULL);
1464 70 : Tra = cgetg(n0+1, t_MAT);
1465 70 : CM_L = scalarmat_s(C, n0);
1466 : /* tmax = current number of traces used (and computed so far) */
1467 70 : for(tmax = 0;; tmax++)
1468 175 : {
1469 245 : long a, b, bmin, bgood, delta, tnew = tmax + 1, r = lg(CM_L)-1;
1470 : GEN M_L, q, CM_Lp, oldCM_L, S1, P1, VV;
1471 245 : int first = 1;
1472 :
1473 : /* bound for f . S_k(genuine factor) = ZC * bound for T_2(S_tnew) */
1474 245 : Btra = mulrr(ZC, mulur(dP*dP, normTp(Br, 2*tnew, dnf)));
1475 245 : bmin = logint(ceil_safe(sqrtr(Btra)), gen_2) + 1;
1476 245 : if (DEBUGLEVEL>2)
1477 0 : err_printf("\nLLL_cmbf: %ld potential factors (tmax = %ld, bmin = %ld)\n",
1478 : r, tmax, bmin);
1479 :
1480 : /* compute Newton sums (possibly relifting first) */
1481 245 : if (gcmp(L->GSmin, Btra) < 0)
1482 : {
1483 : GEN polred;
1484 :
1485 7 : bestlift_init((L->k)<<1, T->nf, Btra, L);
1486 7 : polred = ZqX_normalize(T->polbase, lP, L);
1487 7 : famod = ZqX_liftfact(polred, famod, L->Tpk, L->pk, L->p, L->k);
1488 133 : for (i=1; i<=n0; i++) gel(TT,i) = NULL;
1489 : }
1490 5691 : for (i=1; i<=n0; i++)
1491 : {
1492 5446 : GEN h, lPpow = lP? powiu(lP, tnew): NULL;
1493 5446 : GEN z = polsym_gen(gel(famod,i), gel(TT,i), tnew, L->Tpk, L->pk);
1494 5446 : gel(TT,i) = z;
1495 5446 : h = gel(z,tnew+1);
1496 : /* make Newton sums integral */
1497 5446 : lPpow = mul_content(lPpow, dn);
1498 5446 : if (lPpow)
1499 126 : h = (typ(h) == t_INT)? Fp_mul(h, lPpow, L->pk): FpX_Fp_mul(h, lPpow, L->pk);
1500 5446 : gel(Tra,i) = nf_bestlift(h, NULL, L); /* S_tnew(famod) */
1501 : }
1502 :
1503 : /* compute truncation parameter */
1504 245 : if (DEBUGLEVEL>2) { timer_start(&ti2); timer_start(&TI); }
1505 245 : oldCM_L = CM_L;
1506 245 : av2 = avma;
1507 245 : b = delta = 0; /* -Wall */
1508 336 : AGAIN:
1509 336 : M_L = Q_div_to_int(CM_L, utoipos(C));
1510 336 : VV = get_V(Tra, M_L, L->prk, L->iprk, L->pk, &a);
1511 336 : if (first)
1512 : { /* initialize lattice, using few p-adic digits for traces */
1513 245 : bgood = (long)(a - maxss(32, (long)(BitPerFactor * r)));
1514 245 : b = maxss(bmin, bgood);
1515 245 : delta = a - b;
1516 : }
1517 : else
1518 : { /* add more p-adic digits and continue reduction */
1519 91 : if (a < b) b = a;
1520 91 : b = maxss(b-delta, bmin);
1521 91 : if (b - delta/2 < bmin) b = bmin; /* near there. Go all the way */
1522 : }
1523 :
1524 : /* restart with truncated entries */
1525 336 : q = int2n(b);
1526 336 : P1 = gdivround(L->prk, q);
1527 336 : S1 = gdivround(Tra, q);
1528 336 : T2 = ZM_sub(ZM_mul(S1, M_L), ZM_mul(P1, VV));
1529 336 : m = vconcat( CM_L, T2 );
1530 336 : if (first)
1531 : {
1532 245 : first = 0;
1533 245 : m = shallowconcat( m, vconcat(ZERO, P1) );
1534 : /* [ C M_L 0 ]
1535 : * m = [ ] square matrix
1536 : * [ T2' PRK ] T2' = Tra * M_L truncated
1537 : */
1538 : }
1539 336 : CM_L = LLL_check_progress(Bnorm, n0, m, b == bmin, /*dbg:*/ &ti_LLL);
1540 336 : if (DEBUGLEVEL>2)
1541 0 : err_printf("LLL_cmbf: (a,b) =%4ld,%4ld; r =%3ld -->%3ld, time = %ld\n",
1542 0 : a,b, lg(m)-1, CM_L? lg(CM_L)-1: 1, timer_delay(&TI));
1543 336 : if (!CM_L) { list = mkcol(RgX_int_normalize(P)); break; }
1544 315 : if (b > bmin)
1545 : {
1546 91 : CM_L = gerepilecopy(av2, CM_L);
1547 91 : goto AGAIN;
1548 : }
1549 224 : if (DEBUGLEVEL>2) timer_printf(&ti2, "for this trace");
1550 :
1551 224 : i = lg(CM_L) - 1;
1552 224 : if (i == r && ZM_equal(CM_L, oldCM_L))
1553 : {
1554 91 : CM_L = oldCM_L;
1555 91 : set_avma(av2); continue;
1556 : }
1557 :
1558 133 : CM_Lp = FpM_image(CM_L, utoipos(27449)); /* inexpensive test */
1559 133 : if (lg(CM_Lp) != lg(CM_L))
1560 : {
1561 0 : if (DEBUGLEVEL>2) err_printf("LLL_cmbf: rank decrease\n");
1562 0 : CM_L = ZM_hnf(CM_L);
1563 : }
1564 :
1565 133 : if (i <= r && i*rec < n0)
1566 : {
1567 : pari_timer ti;
1568 77 : if (DEBUGLEVEL>2) timer_start(&ti);
1569 77 : list = nf_chk_factors(T, P, Q_div_to_int(CM_L,utoipos(C)), famod, L->pk);
1570 77 : if (DEBUGLEVEL>2) ti_CF += timer_delay(&ti);
1571 77 : if (list) break;
1572 : }
1573 84 : if (gc_needed(av,1))
1574 : {
1575 6 : if(DEBUGMEM>1) pari_warn(warnmem,"nf_LLL_cmbf");
1576 6 : gerepileall(av, L->Tpk? 9: 8,
1577 : &CM_L,&TT,&Tra,&famod,&L->GSmin,&L->pk,&L->prk,&L->iprk,
1578 : &L->Tpk);
1579 : }
1580 78 : else CM_L = gerepilecopy(av2, CM_L);
1581 : }
1582 70 : if (DEBUGLEVEL>2)
1583 0 : err_printf("* Time LLL: %ld\n* Time Check Factor: %ld\n",ti_LLL,ti_CF);
1584 70 : return list;
1585 : }
1586 :
1587 : static GEN
1588 2240 : nf_combine_factors(nfcmbf_t *T, GEN polred, long klim)
1589 : {
1590 2240 : nflift_t *L = T->L;
1591 : GEN res;
1592 : long maxK;
1593 : int done;
1594 : pari_timer ti;
1595 :
1596 2240 : if (DEBUGLEVEL>2) timer_start(&ti);
1597 2240 : T->fact = ZqX_liftfact(polred, T->fact, L->Tpk, L->pk, L->p, L->k);
1598 2240 : if (DEBUGLEVEL>2) timer_printf(&ti, "Hensel lift");
1599 2240 : res = nfcmbf(T, klim, &maxK, &done);
1600 2240 : if (DEBUGLEVEL>2) timer_printf(&ti, "Naive recombination");
1601 2240 : if (!done)
1602 : {
1603 70 : long l = lg(res)-1;
1604 : GEN v;
1605 70 : if (l > 1)
1606 : {
1607 : GEN den;
1608 35 : T->pol = gel(res,l);
1609 35 : T->polbase = Q_remove_denom(RgX_to_nfX(T->nf, T->pol), &den);
1610 35 : if (den) { T->Br = gmul(T->Br, den); T->pol = RgX_Rg_mul(T->pol, den); }
1611 : }
1612 70 : v = nf_LLL_cmbf(T, maxK);
1613 : /* remove last elt, possibly unfactored. Add all new ones. */
1614 70 : setlg(res, l); res = shallowconcat(res, v);
1615 : }
1616 2240 : return res;
1617 : }
1618 :
1619 : static GEN
1620 2429 : nf_DDF_roots(GEN pol, GEN polred, GEN nfpol, long fl, nflift_t *L)
1621 : {
1622 : GEN z, Cltx_r, ltdn;
1623 : long i, m, lz;
1624 : div_data D;
1625 :
1626 2429 : init_div_data(&D, pol, L);
1627 2429 : ltdn = mul_content(D.lt, L->dn);
1628 2429 : z = ZqX_roots(polred, L->Tpk, L->p, L->k);
1629 2429 : Cltx_r = deg1pol_shallow(D.Clt? D.Clt: gen_1, NULL, varn(pol));
1630 2429 : lz = lg(z);
1631 2429 : if (DEBUGLEVEL > 3) err_printf("Checking %ld roots:",lz-1);
1632 7623 : for (m=1,i=1; i<lz; i++)
1633 : {
1634 5194 : GEN r = gel(z,i);
1635 : int dvd;
1636 : pari_sp av;
1637 5194 : if (DEBUGLEVEL > 3) err_printf(" %ld",i);
1638 : /* lt*dn*topowden * r = Clt * r */
1639 5194 : r = nf_bestlift_to_pol(ltdn? gmul(ltdn,r): r, NULL, L);
1640 5194 : av = avma;
1641 5194 : gel(Cltx_r,2) = gneg(r); /* check P(r) == 0 */
1642 5194 : dvd = ZXQX_dvd(D.C2ltpol, Cltx_r, nfpol); /* integral */
1643 5194 : set_avma(av);
1644 : /* don't go on with q, usually much larger that C2ltpol */
1645 5194 : if (dvd) {
1646 4767 : if (D.Clt) r = gdiv(r, D.Clt);
1647 4767 : gel(z,m++) = r;
1648 : }
1649 427 : else if (fl == ROOTS_SPLIT) return cgetg(1, t_VEC);
1650 : }
1651 2429 : if (DEBUGLEVEL > 3) err_printf(" done\n");
1652 2429 : z[0] = evaltyp(t_VEC) | evallg(m);
1653 2429 : return z;
1654 : }
1655 :
1656 : /* returns a few factors of T in Fp of degree <= maxf, NULL if none exist */
1657 : static GEN
1658 118769 : get_good_factor(GEN T, ulong p, long maxf)
1659 : {
1660 118769 : pari_sp av = avma;
1661 118769 : GEN r, R = gel(Flx_factor(ZX_to_Flx(T,p),p), 1);
1662 118768 : if (maxf == 1)
1663 : { /* degree 1 factors are best */
1664 113021 : r = gel(R,1);
1665 113021 : if (degpol(r) == 1) return mkvec(r);
1666 : }
1667 : else
1668 : { /* otherwise, pick factor of largish degree */
1669 5747 : long i, j, dr, d = 0, l = lg(R);
1670 : GEN v;
1671 5747 : if (l == 2) return mkvec(gel(R,1)); /* inert is fine */
1672 4949 : v = cgetg(l, t_VEC);
1673 34118 : for (i = j = 1; i < l; i++)
1674 : {
1675 30275 : r = gel(R,i); dr = degpol(r);
1676 30275 : if (dr > maxf) break;
1677 29169 : if (dr != d) { gel(v,j++) = r; d = dr; }
1678 : }
1679 4949 : setlg(v,j); if (j > 1) return v;
1680 : }
1681 77832 : return gc_NULL(av); /* failure */
1682 : }
1683 :
1684 : /* n = number of modular factors, f = residue degree; nold/fold current best
1685 : * return 1 if new values are "better" than old ones */
1686 : static int
1687 34111 : record(long nold, long n, long fold, long f)
1688 : {
1689 34111 : if (!nold) return 1; /* uninitialized */
1690 28189 : if (fold == f) return n < nold;
1691 : /* if f increases, allow increasing n a little */
1692 3500 : if (fold < f) return n <= 20 || n < 1.1*nold;
1693 : /* f decreases, only allow if decreasing n a lot */
1694 1512 : return n < 0.7*nold;
1695 : }
1696 : /* Optimization problem: factorization of polynomials over large Fq is slow,
1697 : * BUT bestlift correspondingly faster.
1698 : * Return maximal residue degree to be considered when picking a prime ideal */
1699 : static long
1700 6979 : get_maxf(long nfdeg)
1701 : {
1702 6979 : long maxf = 1;
1703 6979 : if (nfdeg >= 45) maxf =32;
1704 6979 : else if (nfdeg >= 30) maxf =16;
1705 6944 : else if (nfdeg >= 15) maxf = 8;
1706 6979 : return maxf;
1707 : }
1708 : /* number of maximal ideals to test before settling on best prime and number
1709 : * of factors; B = [K:Q]*deg(P) */
1710 : static long
1711 6979 : get_nbprimes(long B)
1712 : {
1713 6979 : if (B <= 128) return 5;
1714 399 : if (B <= 1024) return 20;
1715 42 : if (B <= 2048) return 65;
1716 0 : return 100;
1717 : }
1718 : /* Select a prime ideal pr over which to factor pol.
1719 : * Return the number of factors (or roots, according to flag fl) mod pr.
1720 : * Set:
1721 : * lt: leading term of polbase (t_INT or NULL [ for 1 ])
1722 : * pr: a suitable maximal ideal
1723 : * Fa: factors found mod pr
1724 : * Tp: polynomial defining Fq/Fp */
1725 : static long
1726 6979 : nf_pick_prime(GEN nf, GEN pol, long fl, GEN *lt, GEN *Tp, ulong *pp)
1727 : {
1728 6979 : GEN nfpol = nf_get_pol(nf), bad = mulii(nf_get_disc(nf), nf_get_index(nf));
1729 6979 : long nfdeg = degpol(nfpol), dpol = degpol(pol), nold = 0, fold = 1;
1730 6979 : long maxf = get_maxf(nfdeg), ct = get_nbprimes(nfdeg * dpol);
1731 : ulong p;
1732 : forprime_t S;
1733 : pari_timer ti_pr;
1734 :
1735 6979 : if (DEBUGLEVEL>3) timer_start(&ti_pr);
1736 6979 : *lt = leading_coeff(pol); /* t_INT */
1737 6979 : if (gequal1(*lt)) *lt = NULL;
1738 6979 : *pp = 0;
1739 6979 : *Tp = NULL;
1740 6979 : (void)u_forprime_init(&S, 2, ULONG_MAX);
1741 : /* select pr such that pol has the smallest number of factors, ct attempts */
1742 130220 : while ((p = u_forprime_next(&S)))
1743 : {
1744 : GEN vT;
1745 130221 : long n, i, l, ok = 0;
1746 130221 : ulong ltp = 0;
1747 :
1748 130221 : if (! umodiu(bad,p)) continue;
1749 119910 : if (*lt) { ltp = umodiu(*lt, p); if (!ltp) continue; }
1750 118769 : vT = get_good_factor(nfpol, p, maxf);
1751 118768 : if (!vT) continue;
1752 40935 : l = lg(vT);
1753 80751 : for (i = 1; i < l; i++)
1754 : {
1755 42133 : pari_sp av2 = avma;
1756 42133 : GEN T = gel(vT,i), red = RgX_to_FlxqX(pol, T, p);
1757 42133 : long f = degpol(T);
1758 42133 : if (f == 1)
1759 : { /* degree 1 */
1760 36652 : red = FlxX_to_Flx(red);
1761 36652 : if (ltp) red = Flx_normalize(red, p);
1762 36652 : if (!Flx_is_squarefree(red, p)) { set_avma(av2); continue; }
1763 31247 : ok = 1;
1764 31247 : n = (fl == FACTORS)? Flx_nbfact(red,p): Flx_nbroots(red,p);
1765 : }
1766 : else
1767 : {
1768 5481 : if (ltp) red = FlxqX_normalize(red, T, p);
1769 5481 : if (!FlxqX_is_squarefree(red, T, p)) { set_avma(av2); continue; }
1770 5180 : ok = 1;
1771 5180 : n = (fl == FACTORS)? FlxqX_nbfact(red,T,p): FlxqX_nbroots(red,T,p);
1772 : }
1773 36428 : if (fl == ROOTS_SPLIT && n < dpol) return n; /* not split */
1774 36386 : if (n <= 1)
1775 : {
1776 6391 : if (fl == FACTORS) return n; /* irreducible */
1777 6139 : if (!n) return 0; /* no root */
1778 : }
1779 34118 : if (DEBUGLEVEL>3)
1780 0 : err_printf("%3ld %s at prime (%ld,x^%ld+...)\n Time: %ld\n",
1781 : n, (fl == FACTORS)? "factors": "roots", p,f, timer_delay(&ti_pr));
1782 :
1783 34118 : if (fl == ROOTS && f==nfdeg) { *Tp = T; *pp = p; return n; }
1784 34111 : if (record(nold, n, fold, f)) { nold = n; fold = f; *Tp = T; *pp = p; }
1785 26208 : else set_avma(av2);
1786 : }
1787 38618 : if (ok && --ct <= 0) break;
1788 : }
1789 4662 : if (!nold) pari_err_OVERFLOW("nf_pick_prime [ran out of primes]");
1790 4662 : return nold;
1791 : }
1792 :
1793 : /* Assume lt(T) is a t_INT and T square free. Return t_VEC of irred. factors */
1794 : static GEN
1795 154 : nfsqff_trager(GEN u, GEN T, GEN dent)
1796 : {
1797 154 : long k = 0, i, lx;
1798 154 : GEN U, P, x0, mx0, fa, n = ZX_ZXY_rnfequation(T, u, &k);
1799 : int tmonic;
1800 154 : if (DEBUGLEVEL>4) err_printf("nfsqff_trager: choosing k = %ld\n",k);
1801 :
1802 : /* n guaranteed to be squarefree */
1803 154 : fa = ZX_DDF(Q_primpart(n)); lx = lg(fa);
1804 154 : if (lx == 2) return mkvec(u);
1805 :
1806 133 : tmonic = is_pm1(leading_coeff(T));
1807 133 : P = cgetg(lx,t_VEC);
1808 133 : x0 = deg1pol_shallow(stoi(-k), gen_0, varn(T));
1809 133 : mx0 = deg1pol_shallow(stoi(k), gen_0, varn(T));
1810 133 : U = RgXQX_translate(u, mx0, T);
1811 133 : if (!tmonic) U = Q_primpart(U);
1812 511 : for (i=lx-1; i>0; i--)
1813 : {
1814 378 : GEN f = gel(fa,i), F = nfgcd(U, f, T, dent);
1815 378 : F = RgXQX_translate(F, x0, T);
1816 : /* F = gcd(f, u(t - x0)) [t + x0] = gcd(f(t + x0), u), more efficient */
1817 378 : if (typ(F) != t_POL || degpol(F) == 0)
1818 0 : pari_err_IRREDPOL("factornf [modulus]",T);
1819 378 : gel(P,i) = QXQX_normalize(F, T);
1820 : }
1821 133 : gen_sort_inplace(P, (void*)&cmp_RgX, &gen_cmp_RgX, NULL);
1822 133 : return P;
1823 : }
1824 :
1825 : /* Factor polynomial a on the number field defined by polynomial T, using
1826 : * Trager's trick */
1827 : GEN
1828 14 : polfnf(GEN a, GEN T)
1829 : {
1830 14 : GEN rep = cgetg(3, t_MAT), A, B, y, dent, bad;
1831 : long dA;
1832 : int tmonic;
1833 :
1834 14 : if (typ(a)!=t_POL) pari_err_TYPE("polfnf",a);
1835 14 : if (typ(T)!=t_POL) pari_err_TYPE("polfnf",T);
1836 14 : T = Q_primpart(T); tmonic = is_pm1(leading_coeff(T));
1837 14 : RgX_check_ZX(T,"polfnf");
1838 14 : A = Q_primpart( QXQX_normalize(RgX_nffix("polfnf",T,a,1), T) );
1839 14 : dA = degpol(A);
1840 14 : if (dA <= 0)
1841 : {
1842 0 : set_avma((pari_sp)(rep + 3));
1843 0 : return (dA == 0)? trivial_fact(): zerofact(varn(A));
1844 : }
1845 14 : bad = dent = absi_shallow(ZX_disc(T));
1846 14 : if (tmonic) dent = indexpartial(T, dent);
1847 14 : (void)nfgcd_all(A,RgX_deriv(A), T, dent, &B);
1848 14 : if (degpol(B) != dA) B = Q_primpart( QXQX_normalize(B, T) );
1849 14 : ensure_lt_INT(B);
1850 14 : y = nfsqff_trager(B, T, dent);
1851 14 : fact_from_sqff(rep, A, B, y, T, bad);
1852 14 : return sort_factor_pol(rep, cmp_RgX);
1853 : }
1854 :
1855 : static int
1856 14637 : nfsqff_use_Trager(long n, long dpol)
1857 : {
1858 14637 : return dpol*3<n;
1859 : }
1860 :
1861 : /* return the factorization of the square-free polynomial pol. Not memory-clean
1862 : The coeffs of pol are in Z_nf and its leading term is a rational integer.
1863 : deg(pol) > 0, deg(nfpol) > 1
1864 : fl is either FACTORS (return factors), or ROOTS / ROOTS_SPLIT (return roots):
1865 : - ROOTS, return only the roots of x in nf
1866 : - ROOTS_SPLIT, as ROOTS if pol splits, [] otherwise
1867 : den is usually 1, otherwise nf.zk is doubtful, and den bounds the
1868 : denominator of an arbitrary element of Z_nf on nf.zk */
1869 : static GEN
1870 8932 : nfsqff(GEN nf, GEN pol, long fl, GEN den)
1871 : {
1872 8932 : long n, nbf, dpol = degpol(pol);
1873 : GEN C0, polbase;
1874 8932 : GEN N2, res, polred, lt, nfpol = typ(nf)==t_POL?nf:nf_get_pol(nf);
1875 : ulong pp;
1876 : nfcmbf_t T;
1877 : nflift_t L;
1878 : pari_timer ti, ti_tot;
1879 :
1880 8932 : if (DEBUGLEVEL>2) { timer_start(&ti); timer_start(&ti_tot); }
1881 8932 : n = degpol(nfpol);
1882 : /* deg = 1 => irreducible */
1883 8932 : if (dpol == 1) {
1884 1813 : if (fl == FACTORS) return mkvec(QXQX_normalize(pol, nfpol));
1885 1778 : return mkvec(gneg(gdiv(gel(pol,2),gel(pol,3))));
1886 : }
1887 7119 : if (typ(nf)==t_POL || nfsqff_use_Trager(n,dpol))
1888 : {
1889 : GEN z;
1890 140 : if (DEBUGLEVEL>2) err_printf("Using Trager's method\n");
1891 140 : if (typ(nf) != t_POL) den = mulii(den, nf_get_index(nf));
1892 140 : z = nfsqff_trager(Q_primpart(pol), nfpol, den);
1893 140 : if (fl != FACTORS) {
1894 105 : long i, l = lg(z);
1895 315 : for (i = 1; i < l; i++)
1896 : {
1897 238 : GEN LT, t = gel(z,i); if (degpol(t) > 1) break;
1898 210 : LT = gel(t,3);
1899 210 : if (typ(LT) == t_POL) LT = gel(LT,2); /* constant */
1900 210 : gel(z,i) = gdiv(gel(t,2), negi(LT));
1901 : }
1902 105 : setlg(z, i);
1903 105 : if (fl == ROOTS_SPLIT && i != l) return cgetg(1,t_VEC);
1904 : }
1905 140 : return z;
1906 : }
1907 :
1908 6979 : polbase = RgX_to_nfX(nf, pol);
1909 6979 : nbf = nf_pick_prime(nf, pol, fl, <, &L.Tp, &pp);
1910 6979 : if (L.Tp)
1911 : {
1912 5929 : L.Tp = Flx_to_ZX(L.Tp);
1913 5929 : L.p = utoi(pp);
1914 : }
1915 :
1916 6979 : if (fl == ROOTS_SPLIT && nbf < dpol) return cgetg(1,t_VEC);
1917 6937 : if (nbf <= 1)
1918 : {
1919 3073 : if (fl == FACTORS) return mkvec(QXQX_normalize(pol, nfpol)); /* irred. */
1920 2821 : if (!nbf) return cgetg(1,t_VEC); /* no root */
1921 : }
1922 :
1923 4669 : if (DEBUGLEVEL>2) {
1924 0 : timer_printf(&ti, "choice of a prime ideal");
1925 0 : err_printf("Prime ideal chosen: (%lu,x^%ld+...)\n", pp, degpol(L.Tp));
1926 : }
1927 4669 : L.tozk = nf_get_invzk(nf);
1928 4669 : L.topow= nf_get_zkprimpart(nf);
1929 4669 : L.topowden = nf_get_zkden(nf);
1930 4669 : if (is_pm1(den)) den = NULL;
1931 4669 : L.dn = den;
1932 4669 : T.ZC = L2_bound(nf, den);
1933 4669 : T.Br = nf_root_bounds(nf, pol); if (lt) T.Br = gmul(T.Br, lt);
1934 :
1935 : /* C0 = bound for T_2(Q_i), Q | P */
1936 4669 : if (fl != FACTORS) C0 = normTp(T.Br, 2, n);
1937 2240 : else C0 = nf_factor_bound(nf, polbase);
1938 4669 : T.bound = mulrr(T.ZC, C0); /* bound for |Q_i|^2 in Z^n on chosen Z-basis */
1939 :
1940 4669 : N2 = mulur(dpol*dpol, normTp(T.Br, 4, n)); /* bound for T_2(lt * S_2) */
1941 4669 : T.BS_2 = mulrr(T.ZC, N2); /* bound for |S_2|^2 on chosen Z-basis */
1942 :
1943 4669 : if (DEBUGLEVEL>2) {
1944 0 : timer_printf(&ti, "bound computation");
1945 0 : err_printf(" 1) T_2 bound for %s: %Ps\n",
1946 : fl == FACTORS?"factor": "root", C0);
1947 0 : err_printf(" 2) Conversion from T_2 --> | |^2 bound : %Ps\n", T.ZC);
1948 0 : err_printf(" 3) Final bound: %Ps\n", T.bound);
1949 : }
1950 :
1951 4669 : bestlift_init(0, nf, T.bound, &L);
1952 4669 : if (DEBUGLEVEL>2) timer_start(&ti);
1953 4669 : polred = ZqX_normalize(polbase, lt, &L); /* monic */
1954 :
1955 4669 : if (fl != FACTORS) {
1956 2429 : GEN z = nf_DDF_roots(pol, polred, nfpol, fl, &L);
1957 2429 : if (lg(z) == 1) return cgetg(1, t_VEC);
1958 2373 : return z;
1959 : }
1960 :
1961 2240 : T.fact = gel(FqX_factor(polred, L.Tp, L.p), 1);
1962 2240 : if (DEBUGLEVEL>2)
1963 0 : timer_printf(&ti, "splitting mod %Ps^%ld", L.p, degpol(L.Tp));
1964 2240 : T.L = &L;
1965 2240 : T.polbase = polbase;
1966 2240 : T.pol = pol;
1967 2240 : T.nf = nf;
1968 2240 : res = nf_combine_factors(&T, polred, dpol-1);
1969 2240 : if (DEBUGLEVEL>2)
1970 0 : err_printf("Total Time: %ld\n===========\n", timer_delay(&ti_tot));
1971 2240 : return res;
1972 : }
1973 :
1974 : /* assume pol monic in nf.zk[X] */
1975 : GEN
1976 84 : nfroots_if_split(GEN *pnf, GEN pol)
1977 : {
1978 84 : GEN T = get_nfpol(*pnf,pnf), den = fix_nf(pnf, &T, &pol);
1979 84 : pari_sp av = avma;
1980 84 : GEN z = nfsqff(*pnf, pol, ROOTS_SPLIT, den);
1981 84 : if (lg(z) == 1) return gc_NULL(av);
1982 42 : return gerepilecopy(av, z);
1983 : }
1984 :
1985 : /*******************************************************************/
1986 : /* */
1987 : /* Roots of unity in a number field */
1988 : /* (alternative to nfrootsof1 using factorization in K[X]) */
1989 : /* */
1990 : /*******************************************************************/
1991 : /* Code adapted from nffactor. Structure of the algorithm; only step 1 is
1992 : * specific to roots of unity.
1993 : *
1994 : * [Step 1]: guess roots via ramification. If trivial output this.
1995 : * [Step 2]: select prime [p] unramified and ideal [pr] above
1996 : * [Step 4]: factor the cyclotomic polynomial mod [pr],
1997 : */
1998 :
1999 : /* Guesses the number of roots of unity in number field [nf].
2000 : * Computes gcd of N(P)-1 for some primes. The value returned is a proven
2001 : * multiple of the correct value. */
2002 : static long
2003 8848 : guess_roots(GEN T, GEN D, GEN index)
2004 : {
2005 8848 : long c = 0, nfdegree = degpol(T), B = nfdegree + 20, l;
2006 8848 : ulong p = 2;
2007 8848 : GEN nbroots = NULL;
2008 : forprime_t S;
2009 : pari_sp av;
2010 :
2011 8848 : (void)u_forprime_init(&S, 3, ULONG_MAX);
2012 8848 : av = avma;
2013 : /* result must be stationary (counter c) for at least B loops */
2014 249760 : for (l=1; (p = u_forprime_next(&S)); l++)
2015 : {
2016 : GEN old, F, pf_1, Tp;
2017 249760 : ulong i, gcdf = 0;
2018 : long nb;
2019 :
2020 249760 : if (!umodiu(D,p) || !umodiu(index,p)) continue;
2021 238546 : Tp = ZX_to_Flx(T,p); /* squarefree */
2022 238546 : F = Flx_nbfact_by_degree(Tp, &nb, p);
2023 : /* the gcd of the p^f - 1 is p^(gcd of the f's) - 1 */
2024 766815 : for (i = 1; i <= (ulong) nfdegree; i++)
2025 633444 : if (F[i]) {
2026 239120 : gcdf = gcdf? ugcd(gcdf, i): i;
2027 239120 : if (gcdf == 1) break;
2028 : }
2029 238546 : pf_1 = subiu(powuu(p, gcdf), 1);
2030 238546 : old = nbroots;
2031 238546 : nbroots = nbroots? gcdii(pf_1, nbroots): pf_1;
2032 238546 : if (DEBUGLEVEL>5)
2033 0 : err_printf("p=%lu; gcf(f(P/p))=%ld; nbroots | %Ps",p, gcdf, nbroots);
2034 : /* if same result go on else reset the stop counter [c] */
2035 238546 : if (old && equalii(nbroots,old))
2036 221298 : { if (!is_bigint(nbroots) && ++c > B) break; }
2037 : else
2038 17248 : c = 0;
2039 : }
2040 8848 : if (!nbroots) pari_err_OVERFLOW("guess_roots [ran out of primes]");
2041 8848 : if (DEBUGLEVEL>5) err_printf("%ld loops\n",l);
2042 8848 : return gc_long(av, itos(nbroots));
2043 : }
2044 :
2045 : /* T(x) an irreducible ZX. Is it of the form Phi_n(c \pm x) ?
2046 : * Return NULL if not, and a root of 1 of maximal order in Z[x]/(T) otherwise
2047 : *
2048 : * N.B. Set n_squarefree = 1 if n is squarefree, and 0 otherwise.
2049 : * This last parameter is inconvenient, but it allows a cheap
2050 : * stringent test. (n guessed from guess_roots())*/
2051 : static GEN
2052 1267 : ZXirred_is_cyclo_translate(GEN T, long n_squarefree)
2053 : {
2054 1267 : long r, m, d = degpol(T);
2055 1267 : GEN T1, c = divis_rem(gel(T, d+1), d, &r); /* d-1 th coeff divisible by d ? */
2056 : /* The trace coefficient of polcyclo(n) is \pm1 if n is square free, and 0
2057 : * otherwise. */
2058 1267 : if (!n_squarefree)
2059 637 : { if (r) return NULL; }
2060 : else
2061 : {
2062 630 : if (r < -1)
2063 : {
2064 0 : r += d;
2065 0 : c = subiu(c, 1);
2066 : }
2067 630 : else if (r == d-1)
2068 : {
2069 35 : r = -1;
2070 35 : c = addiu(c, 1);
2071 : }
2072 630 : if (r != 1 && r != -1) return NULL;
2073 : }
2074 1232 : if (signe(c)) /* presumably Phi_guess(c \pm x) */
2075 35 : T = RgX_translate(T, negi(c));
2076 1232 : if (!n_squarefree) T = RgX_deflate_max(T, &m);
2077 : /* presumably Phi_core(guess)(\pm x), cyclotomic iff original T was */
2078 1232 : T1 = ZX_graeffe(T);
2079 1232 : if (ZX_equal(T, T1)) /* T = Phi_n, n odd */
2080 35 : return deg1pol_shallow(gen_m1, negi(c), varn(T));
2081 1197 : else if (ZX_equal(T1, ZX_z_unscale(T, -1))) /* T = Phi_{2n}, nodd */
2082 1176 : return deg1pol_shallow(gen_1, c, varn(T));
2083 21 : return NULL;
2084 : }
2085 :
2086 : static GEN
2087 10521 : trivroots(void) { return mkvec2(gen_2, gen_m1); }
2088 : /* Number of roots of unity in number field [nf]. */
2089 : GEN
2090 12411 : nfrootsof1(GEN nf)
2091 : {
2092 : GEN T, fa, LP, LE, z, disc, index;
2093 : pari_timer ti;
2094 : long i, l, nbguessed, nbroots, nfdegree;
2095 : pari_sp av;
2096 :
2097 12411 : T = get_nfpol(nf, &nf);
2098 12411 : RgX_check_ZX(T, "nfrootsof1");
2099 12411 : if (nf && nf_get_r1(nf)) return trivroots();
2100 8848 : disc = nf ? nf_get_disc(nf): ZX_disc(T);
2101 8848 : index = nf ? nf_get_index(nf): gen_1;
2102 : /* Step 1 : guess number of roots and discard trivial case 2 */
2103 8848 : if (DEBUGLEVEL>2) timer_start(&ti);
2104 8848 : nbguessed = guess_roots(T, disc, index);
2105 8848 : if (DEBUGLEVEL>2)
2106 0 : timer_printf(&ti, "guessing roots of 1 [guess = %ld]", nbguessed);
2107 8848 : if (nbguessed == 2) return trivroots();
2108 :
2109 1897 : nfdegree = degpol(T);
2110 1897 : fa = factoru(nbguessed);
2111 1897 : LP = gel(fa,1); l = lg(LP);
2112 1897 : LE = gel(fa,2);
2113 5019 : for (i = 1; i < l; i++)
2114 : {
2115 3122 : long p = LP[i];
2116 : /* Degree and ramification test: find largest k such that Q(zeta_{p^k})
2117 : * may be a subfield of K. Q(zeta_p^k) has degree (p-1)p^(k-1)
2118 : * and v_p(discriminant) = ((p-1)k-1)p^(k-1); so we must have
2119 : * v_p(disc_K) >= ((p-1)k-1) * n / (p-1) = kn - q, where q = n/(p-1) */
2120 3122 : if (p == 2)
2121 : { /* the test simplifies a little in that case */
2122 : long v, vnf, k;
2123 1897 : if (LE[i] == 1) continue;
2124 742 : vnf = vals(nfdegree);
2125 742 : v = vali(disc);
2126 784 : for (k = minss(LE[i], vnf+1); k >= 1; k--)
2127 784 : if (v >= nfdegree*(k-1)) { nbguessed >>= LE[i]-k; LE[i] = k; break; }
2128 : /* N.B the test above always works for k = 1: LE[i] >= 1 */
2129 : }
2130 : else
2131 : {
2132 : long v, vnf, k;
2133 1225 : ulong r, q = udivuu_rem(nfdegree, p-1, &r);
2134 1225 : if (r) { nbguessed /= upowuu(p, LE[i]); LE[i] = 0; continue; }
2135 : /* q = n/(p-1) */
2136 1225 : vnf = u_lval(q, p);
2137 1225 : v = Z_lval(disc, p);
2138 1232 : for (k = minss(LE[i], vnf+1); k >= 0; k--)
2139 1232 : if (v >= nfdegree*k-(long)q)
2140 1225 : { nbguessed /= upowuu(p, LE[i]-k); LE[i] = k; break; }
2141 : /* N.B the test above always works for k = 0: LE[i] >= 0 */
2142 : }
2143 : }
2144 1897 : if (DEBUGLEVEL>2)
2145 0 : timer_printf(&ti, "after ramification conditions [guess = %ld]", nbguessed);
2146 1897 : if (nbguessed == 2) return trivroots();
2147 1890 : av = avma;
2148 :
2149 : /* Step 1.5 : test if nf.pol == subst(polcyclo(nbguessed), x, \pm x+c) */
2150 1890 : if (eulerphiu_fact(fa) == (ulong)nfdegree)
2151 : {
2152 1267 : z = ZXirred_is_cyclo_translate(T, uissquarefree_fact(fa));
2153 1267 : if (DEBUGLEVEL>2) timer_printf(&ti, "checking for cyclotomic polynomial");
2154 1267 : if (z)
2155 : {
2156 1211 : if (nf) z = nf_to_scalar_or_basis(nf,z);
2157 1211 : return gerepilecopy(av, mkvec2(utoipos(nbguessed), z));
2158 : }
2159 56 : set_avma(av);
2160 : }
2161 :
2162 : /* Step 2 : actual computation of roots */
2163 679 : nbroots = 2; z = scalarpol(gen_m1, varn(T));
2164 1953 : for (i = 1; i < l; i++)
2165 : { /* for all prime power factors of nbguessed, find a p^k-th root of unity */
2166 1274 : long k, p = LP[i];
2167 1911 : for (k = LE[i]; k > 0; k--)
2168 : { /* find p^k-th roots */
2169 1302 : pari_sp av = avma;
2170 1302 : long pk = upowuu(p,k);
2171 : GEN r;
2172 1302 : if (pk==2) continue; /* no need to test second roots ! */
2173 714 : r = nfisincl0(polcyclo(pk, 0), T, 1);
2174 714 : if (!isintzero(r))
2175 : {
2176 665 : if (DEBUGLEVEL>2) err_printf(" %s root of unity found\n",uordinal(pk));
2177 665 : if (p==2) { nbroots = pk; z = r; }
2178 : else
2179 : {
2180 574 : nbroots *= pk;
2181 574 : z = nf ? nfmul(nf, z, r): QXQ_mul(z, r, T);
2182 : }
2183 665 : break;
2184 : }
2185 49 : set_avma(av);
2186 49 : if (DEBUGLEVEL) pari_warn(warner,"nfrootsof1: wrong guess");
2187 : }
2188 : }
2189 679 : if (nf) z = nf_to_scalar_or_basis(nf,z);
2190 679 : return gerepilecopy(av, mkvec2(utoi(nbroots), z));
2191 : }
|
