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