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