Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 : #include "pari.h"
15 : #include "paripriv.h"
16 :
17 : #define DEBUGLEVEL DEBUGLEVEL_factor
18 :
19 : /* x,y two ZX, y non constant. Return q = x/y if y divides x in Z[X] and NULL
20 : * otherwise. If not NULL, B is a t_INT upper bound for ||q||_oo. */
21 : static GEN
22 6588023 : ZX_divides_i(GEN x, GEN y, GEN B)
23 : {
24 : long dx, dy, dz, i, j;
25 : pari_sp av;
26 : GEN z,p1,y_lead;
27 :
28 6588023 : dy=degpol(y);
29 6588023 : dx=degpol(x);
30 6588025 : dz=dx-dy; if (dz<0) return NULL;
31 6586947 : z=cgetg(dz+3,t_POL); z[1] = x[1];
32 6586946 : x += 2; y += 2; z += 2;
33 6586946 : y_lead = gel(y,dy);
34 6586946 : if (equali1(y_lead)) y_lead = NULL;
35 :
36 6586946 : p1 = gel(x,dx);
37 6586946 : if (y_lead) {
38 : GEN r;
39 33929 : p1 = dvmdii(p1,y_lead, &r);
40 33929 : if (r != gen_0) return NULL;
41 : }
42 6553017 : else p1 = icopy(p1);
43 6582786 : gel(z,dz) = p1;
44 8284833 : for (i=dx-1; i>=dy; i--)
45 : {
46 1707322 : av = avma; p1 = gel(x,i);
47 5896468 : for (j=i-dy+1; j<=i && j<=dz; j++)
48 4189177 : p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
49 1707291 : if (y_lead) {
50 : GEN r;
51 61061 : p1 = dvmdii(p1,y_lead, &r);
52 61061 : if (r != gen_0) return NULL;
53 : }
54 1705792 : if (B && abscmpii(p1, B) > 0) return NULL;
55 1702019 : p1 = gc_INT(av, p1);
56 1702047 : gel(z,i-dy) = p1;
57 : }
58 6577511 : av = avma;
59 17151147 : for (; i >= 0; i--)
60 : {
61 10616800 : p1 = gel(x,i);
62 : /* we always enter this loop at least once */
63 23902082 : for (j=0; j<=i && j<=dz; j++)
64 13285306 : p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
65 10616776 : if (signe(p1)) return NULL;
66 10573634 : set_avma(av);
67 : }
68 6534347 : return z - 2;
69 : }
70 : static GEN
71 6558622 : ZX_divides(GEN x, GEN y) { return ZX_divides_i(x,y,NULL); }
72 :
73 : #if 0
74 : /* cf Beauzamy et al: upper bound for
75 : * lc(x) * [2^(5/8) / pi^(3/8)] e^(1/4n) 2^(n/2) sqrt([x]_2)/ n^(3/8)
76 : * where [x]_2 = sqrt(\sum_i=0^n x[i]^2 / binomial(n,i)). One factor has
77 : * all coeffs less than then bound */
78 : static GEN
79 : two_factor_bound(GEN x)
80 : {
81 : long i, j, n = lg(x) - 3;
82 : pari_sp av = avma;
83 : GEN *invbin, c, r = cgetr(LOWDEFAULTPREC), z;
84 :
85 : x += 2; invbin = (GEN*)new_chunk(n+1);
86 : z = real_1(LOWDEFAULTPREC); /* invbin[i] = 1 / binomial(n, i) */
87 : for (i=0,j=n; j >= i; i++,j--)
88 : {
89 : invbin[i] = invbin[j] = z;
90 : z = divru(mulru(z, i+1), n-i);
91 : }
92 : z = invbin[0]; /* = 1 */
93 : for (i=0; i<=n; i++)
94 : {
95 : c = gel(x,i); if (!signe(c)) continue;
96 : affir(c, r);
97 : z = addrr(z, mulrr(sqrr(r), invbin[i]));
98 : }
99 : z = shiftr(sqrtr(z), n);
100 : z = divrr(z, dbltor(pow((double)n, 0.75)));
101 : z = roundr_safe(sqrtr(z));
102 : z = mulii(z, absi_shallow(gel(x,n)));
103 : return gc_INT(av, shifti(z, 1));
104 : }
105 : #endif
106 :
107 : /* A | S ==> |a_i| <= binom(d-1, i-1) || S ||_2 + binom(d-1, i) lc(S) */
108 : static GEN
109 61278 : Mignotte_bound(GEN S)
110 : {
111 61278 : long i, d = degpol(S);
112 61278 : GEN C, N2, t, binlS, lS = leading_coeff(S), bin = vecbinomial(d-1);
113 :
114 61277 : N2 = sqrtr(RgX_fpnorml2(S,DEFAULTPREC));
115 61278 : binlS = is_pm1(lS)? bin: ZC_Z_mul(bin, lS);
116 :
117 : /* i = 0 */
118 61278 : C = gel(binlS,1);
119 : /* i = d */
120 61278 : t = N2; if (gcmp(C, t) < 0) C = t;
121 525408 : for (i = 1; i < d; i++)
122 : {
123 464130 : t = addri(mulir(gel(bin,i), N2), gel(binlS,i+1));
124 464126 : if (mpcmp(C, t) < 0) C = t;
125 : }
126 61278 : return C;
127 : }
128 : /* A | S ==> |a_i|^2 <= 3^{3/2 + d} / (4 \pi d) [P]_2^2,
129 : * where [P]_2 is Bombieri's 2-norm */
130 : static GEN
131 61278 : Beauzamy_bound(GEN S)
132 : {
133 61278 : const long prec = DEFAULTPREC;
134 61278 : long i, d = degpol(S);
135 : GEN bin, lS, s, C;
136 61278 : bin = vecbinomial(d);
137 :
138 61278 : s = real_0(prec);
139 647957 : for (i=0; i<=d; i++)
140 : {
141 586679 : GEN c = gel(S,i+2);
142 586679 : if (gequal0(c)) continue;
143 : /* s += P_i^2 / binomial(d,i) */
144 492494 : s = addrr(s, divri(itor(sqri(c), prec), gel(bin,i+1)));
145 : }
146 : /* s = [S]_2^2 */
147 61278 : C = powruhalf(utor(3,prec), 3 + 2*d); /* 3^{3/2 + d} */
148 61278 : C = divrr(mulrr(C, s), mulur(4*d, mppi(prec)));
149 61278 : lS = absi_shallow(leading_coeff(S));
150 61278 : return mulir(lS, sqrtr(C));
151 : }
152 :
153 : static GEN
154 61278 : factor_bound(GEN S)
155 : {
156 61278 : pari_sp av = avma;
157 61278 : GEN a = Mignotte_bound(S);
158 61278 : GEN b = Beauzamy_bound(S);
159 61277 : if (DEBUGLEVEL>2)
160 : {
161 0 : err_printf("Mignotte bound: %Ps\n",a);
162 0 : err_printf("Beauzamy bound: %Ps\n",b);
163 : }
164 61277 : return gc_upto(av, ceil_safe(gmin_shallow(a, b)));
165 : }
166 :
167 : /* Naive recombination of modular factors: combine up to maxK modular
168 : * factors, degree <= klim
169 : *
170 : * target = polynomial we want to factor
171 : * famod = array of modular factors. Product should be congruent to
172 : * target/lc(target) modulo p^a
173 : * For true factors: S1,S2 <= p^b, with b <= a and p^(b-a) < 2^31 */
174 : static GEN
175 50784 : cmbf(GEN pol, GEN famod, GEN bound, GEN p, long a, long b,
176 : long klim, long *pmaxK, int *done)
177 : {
178 50784 : long K = 1, cnt = 1, i,j,k, curdeg, lfamod = lg(famod)-1;
179 : ulong spa_b, spa_bs2, Sbound;
180 50784 : GEN lc, lcpol, pa = powiu(p,a), pas2 = shifti(pa,-1);
181 50784 : GEN trace1 = cgetg(lfamod+1, t_VECSMALL);
182 50784 : GEN trace2 = cgetg(lfamod+1, t_VECSMALL);
183 50784 : GEN ind = cgetg(lfamod+1, t_VECSMALL);
184 50784 : GEN deg = cgetg(lfamod+1, t_VECSMALL);
185 50784 : GEN degsofar = cgetg(lfamod+1, t_VECSMALL);
186 50784 : GEN listmod = cgetg(lfamod+1, t_VEC);
187 50784 : GEN fa = cgetg(lfamod+1, t_VEC);
188 :
189 50784 : *pmaxK = cmbf_maxK(lfamod);
190 50784 : lc = absi_shallow(leading_coeff(pol));
191 50784 : if (equali1(lc)) lc = NULL;
192 50784 : lcpol = lc? ZX_Z_mul(pol, lc): pol;
193 :
194 : {
195 50784 : GEN pa_b,pa_bs2,pb, lc2 = lc? sqri(lc): NULL;
196 :
197 50784 : pa_b = powiu(p, a-b); /* < 2^31 */
198 50784 : pa_bs2 = shifti(pa_b,-1);
199 50784 : pb= powiu(p, b);
200 179368 : for (i=1; i <= lfamod; i++)
201 : {
202 128583 : GEN T1,T2, P = gel(famod,i);
203 128583 : long d = degpol(P);
204 :
205 128583 : deg[i] = d; P += 2;
206 128583 : T1 = gel(P,d-1);/* = - S_1 */
207 128583 : T2 = sqri(T1);
208 128583 : if (d > 1) T2 = subii(T2, shifti(gel(P,d-2),1));
209 128583 : T2 = modii(T2, pa); /* = S_2 Newton sum */
210 128583 : if (lc)
211 : {
212 4193 : T1 = Fp_mul(lc, T1, pa);
213 4193 : T2 = Fp_mul(lc2,T2, pa);
214 : }
215 128583 : uel(trace1,i) = itou(diviiround(T1, pb));
216 128582 : uel(trace2,i) = itou(diviiround(T2, pb));
217 : }
218 50785 : spa_b = uel(pa_b,2); /* < 2^31 */
219 50785 : spa_bs2 = uel(pa_bs2,2); /* < 2^31 */
220 : }
221 50785 : degsofar[0] = 0; /* sentinel */
222 :
223 : /* ind runs through strictly increasing sequences of length K,
224 : * 1 <= ind[i] <= lfamod */
225 92295 : nextK:
226 92295 : if (K > *pmaxK || 2*K > lfamod) goto END;
227 56462 : if (DEBUGLEVEL > 3)
228 0 : err_printf("\n### K = %d, %Ps combinations\n", K,binomial(utoipos(lfamod), K));
229 56462 : setlg(ind, K+1); ind[1] = 1;
230 56462 : Sbound = (ulong) ((K+1)>>1);
231 56462 : i = 1; curdeg = deg[ind[1]];
232 : for(;;)
233 : { /* try all combinations of K factors */
234 629543 : for (j = i; j < K; j++)
235 : {
236 87710 : degsofar[j] = curdeg;
237 87710 : ind[j+1] = ind[j]+1; curdeg += deg[ind[j+1]];
238 : }
239 541833 : if (curdeg <= klim) /* trial divide */
240 10486 : {
241 : GEN y, q, list;
242 : pari_sp av;
243 : ulong t;
244 :
245 : /* d - 1 test */
246 1365131 : for (t=uel(trace1,ind[1]),i=2; i<=K; i++)
247 823298 : t = Fl_add(t, uel(trace1,ind[i]), spa_b);
248 541833 : if (t > spa_bs2) t = spa_b - t;
249 541833 : if (t > Sbound)
250 : {
251 431144 : if (DEBUGLEVEL>6) err_printf(".");
252 431144 : goto NEXT;
253 : }
254 : /* d - 2 test */
255 238824 : for (t=uel(trace2,ind[1]),i=2; i<=K; i++)
256 128135 : t = Fl_add(t, uel(trace2,ind[i]), spa_b);
257 110689 : if (t > spa_bs2) t = spa_b - t;
258 110689 : if (t > Sbound)
259 : {
260 58884 : if (DEBUGLEVEL>6) err_printf("|");
261 58884 : goto NEXT;
262 : }
263 :
264 51805 : av = avma;
265 : /* check trailing coeff */
266 51805 : y = lc;
267 152015 : for (i=1; i<=K; i++)
268 : {
269 100210 : GEN q = constant_coeff(gel(famod,ind[i]));
270 100210 : if (y) q = mulii(y, q);
271 100210 : y = centermodii(q, pa, pas2);
272 : }
273 51805 : if (!signe(y) || !dvdii(constant_coeff(lcpol), y))
274 : {
275 22566 : if (DEBUGLEVEL>3) err_printf("T");
276 22566 : set_avma(av); goto NEXT;
277 : }
278 29239 : y = lc; /* full computation */
279 65338 : for (i=1; i<=K; i++)
280 : {
281 36099 : GEN q = gel(famod,ind[i]);
282 36099 : if (y) q = gmul(y, q);
283 36099 : y = centermod_i(q, pa, pas2);
284 : }
285 :
286 : /* y is the candidate factor */
287 29239 : if (! (q = ZX_divides_i(lcpol,y,bound)) )
288 : {
289 3801 : if (DEBUGLEVEL>3) err_printf("*");
290 3801 : set_avma(av); goto NEXT;
291 : }
292 : /* found a factor */
293 25438 : list = cgetg(K+1, t_VEC);
294 25438 : gel(listmod,cnt) = list;
295 51590 : for (i=1; i<=K; i++) list[i] = famod[ind[i]];
296 :
297 25438 : y = Q_primpart(y);
298 25438 : gel(fa,cnt++) = y;
299 : /* fix up pol */
300 25438 : pol = q;
301 25438 : if (lc) pol = Q_div_to_int(pol, leading_coeff(y));
302 93877 : for (i=j=k=1; i <= lfamod; i++)
303 : { /* remove used factors */
304 68439 : if (j <= K && i == ind[j]) j++;
305 : else
306 : {
307 42287 : gel(famod,k) = gel(famod,i);
308 42287 : uel(trace1,k) = uel(trace1,i);
309 42287 : uel(trace2,k) = uel(trace2,i);
310 42287 : deg[k] = deg[i]; k++;
311 : }
312 : }
313 25438 : lfamod -= K;
314 25438 : *pmaxK = cmbf_maxK(lfamod);
315 25438 : if (lfamod < 2*K) goto END;
316 10486 : i = 1; curdeg = deg[ind[1]];
317 10486 : bound = factor_bound(pol);
318 10486 : if (lc) lc = absi_shallow(leading_coeff(pol));
319 10486 : lcpol = lc? ZX_Z_mul(pol, lc): pol;
320 10486 : if (DEBUGLEVEL>3)
321 0 : err_printf("\nfound factor %Ps\nremaining modular factor(s): %ld\n",
322 : y, lfamod);
323 10486 : continue;
324 : }
325 :
326 0 : NEXT:
327 516395 : for (i = K+1;;)
328 : {
329 644901 : if (--i == 0) { K++; goto nextK; }
330 603391 : if (++ind[i] <= lfamod - K + i)
331 : {
332 474885 : curdeg = degsofar[i-1] + deg[ind[i]];
333 474885 : if (curdeg <= klim) break;
334 : }
335 : }
336 : }
337 50785 : END:
338 50785 : *done = 1;
339 50785 : if (degpol(pol) > 0)
340 : { /* leftover factor */
341 50785 : if (signe(leading_coeff(pol)) < 0) pol = ZX_neg(pol);
342 50785 : if (lfamod >= 2*K) *done = 0;
343 :
344 50785 : setlg(famod, lfamod+1);
345 50785 : gel(listmod,cnt) = leafcopy(famod);
346 50785 : gel(fa,cnt++) = pol;
347 : }
348 50785 : if (DEBUGLEVEL>6) err_printf("\n");
349 50785 : setlg(listmod, cnt);
350 50785 : setlg(fa, cnt); return mkvec2(fa, listmod);
351 : }
352 :
353 : /* recombination of modular factors: van Hoeij's algorithm */
354 :
355 : /* Q in Z[X], return Q(2^n) */
356 : static GEN
357 184819 : shifteval(GEN Q, long n)
358 : {
359 184819 : pari_sp av = avma;
360 184819 : long i, l = lg(Q);
361 : GEN s;
362 :
363 184819 : if (!signe(Q)) return gen_0;
364 184819 : s = gel(Q,l-1);
365 976878 : for (i = l-2; i > 1; i--)
366 : {
367 792061 : s = addii(gel(Q,i), shifti(s, n));
368 792059 : if (gc_needed(av,1)) s = gc_INT(av, s);
369 : }
370 184817 : return s;
371 : }
372 :
373 : /* return integer y such that all |a| <= y if P(a) = 0 */
374 : static GEN
375 111029 : root_bound(GEN P0)
376 : {
377 111029 : GEN Q = leafcopy(P0), lP = absi_shallow(leading_coeff(Q)), x,y,z;
378 111029 : long k, d = degpol(Q);
379 :
380 : /* P0 = lP x^d + Q, deg Q < d */
381 111029 : Q = normalizepol_lg(Q, d+2);
382 696676 : for (k=lg(Q)-1; k>1; k--) gel(Q,k) = absi_shallow(gel(Q,k));
383 111029 : k = (long)(fujiwara_bound(P0));
384 185645 : for ( ; k >= 0; k--)
385 : {
386 184819 : pari_sp av = avma;
387 : /* y = 2^k; Q(y) >= lP y^d ? */
388 184819 : if (cmpii(shifteval(Q,k), shifti(lP, d*k)) >= 0) break;
389 74616 : set_avma(av);
390 : }
391 111029 : if (k < 0) k = 0;
392 111029 : y = int2n(k+1);
393 111029 : if (d > 2000) return y; /* likely to be expensive, don't bother */
394 111029 : x = int2n(k);
395 111029 : for(k=0; ; k++)
396 : {
397 600953 : z = shifti(addii(x,y), -1);
398 600953 : if (equalii(x,z) || k > 5) break;
399 489924 : if (cmpii(ZX_Z_eval(Q,z), mulii(lP, powiu(z, d))) < 0)
400 262588 : y = z;
401 : else
402 227336 : x = z;
403 : }
404 111029 : return y;
405 : }
406 :
407 : GEN
408 350 : chk_factors_get(GEN lt, GEN famod, GEN c, GEN T, GEN N)
409 : {
410 350 : long i = 1, j, l = lg(famod);
411 350 : GEN V = cgetg(l, t_VEC);
412 8414 : for (j = 1; j < l; j++)
413 8064 : if (signe(gel(c,j))) gel(V,i++) = gel(famod,j);
414 350 : if (lt && i > 1) gel(V,1) = RgX_Rg_mul(gel(V,1), lt);
415 350 : setlg(V, i);
416 350 : return T? FpXQXV_prod(V, T, N): FpXV_prod(V,N);
417 : }
418 :
419 : static GEN
420 140 : chk_factors(GEN P, GEN M_L, GEN bound, GEN famod, GEN pa)
421 : {
422 : long i, r;
423 140 : GEN pol = P, list, piv, y, ltpol, lt, paov2;
424 :
425 140 : piv = ZM_hnf_knapsack(M_L);
426 140 : if (!piv) return NULL;
427 70 : if (DEBUGLEVEL>7) err_printf("ZM_hnf_knapsack output:\n%Ps\n",piv);
428 :
429 70 : r = lg(piv)-1;
430 70 : list = cgetg(r+1, t_VEC);
431 70 : lt = absi_shallow(leading_coeff(pol));
432 70 : if (equali1(lt)) lt = NULL;
433 70 : ltpol = lt? ZX_Z_mul(pol, lt): pol;
434 70 : paov2 = shifti(pa,-1);
435 70 : for (i = 1;;)
436 : {
437 161 : if (DEBUGLEVEL) err_printf("LLL_cmbf: checking factor %ld\n",i);
438 161 : y = chk_factors_get(lt, famod, gel(piv,i), NULL, pa);
439 161 : y = FpX_center_i(y, pa, paov2);
440 161 : if (! (pol = ZX_divides_i(ltpol,y,bound)) ) return NULL;
441 133 : if (lt) y = Q_primpart(y);
442 133 : gel(list,i) = y;
443 133 : if (++i >= r) break;
444 :
445 91 : if (lt)
446 : {
447 35 : pol = ZX_Z_divexact(pol, leading_coeff(y));
448 35 : lt = absi_shallow(leading_coeff(pol));
449 35 : ltpol = ZX_Z_mul(pol, lt);
450 : }
451 : else
452 56 : ltpol = pol;
453 : }
454 42 : y = Q_primpart(pol);
455 42 : gel(list,i) = y; return list;
456 : }
457 :
458 : GEN
459 1547 : LLL_check_progress(GEN Bnorm, long n0, GEN m, int final, long *ti_LLL)
460 : {
461 : GEN norm, u;
462 : long i, R;
463 : pari_timer T;
464 :
465 1547 : if (DEBUGLEVEL>2) timer_start(&T);
466 1547 : u = ZM_lll_norms(m, final? 0.999: 0.75, LLL_INPLACE | LLL_NOFLATTER, &norm);
467 1547 : if (DEBUGLEVEL>2) *ti_LLL += timer_delay(&T);
468 10542 : for (R=lg(m)-1; R > 0; R--)
469 10542 : if (cmprr(gel(norm,R), Bnorm) < 0) break;
470 16156 : for (i=1; i<=R; i++) setlg(u[i], n0+1);
471 1547 : if (R <= 1)
472 : {
473 105 : if (!R) pari_err_BUG("LLL_cmbf [no factor]");
474 105 : return NULL; /* irreducible */
475 : }
476 1442 : setlg(u, R+1); return u;
477 : }
478 :
479 : static ulong
480 14 : next2pow(ulong a)
481 : {
482 14 : ulong b = 1;
483 112 : while (b < a) b <<= 1;
484 14 : return b;
485 : }
486 :
487 : /* Recombination phase of Berlekamp-Zassenhaus algorithm using a variant of
488 : * van Hoeij's knapsack
489 : *
490 : * P = squarefree in Z[X].
491 : * famod = array of (lifted) modular factors mod p^a
492 : * bound = Mignotte bound for the size of divisors of P (for the sup norm)
493 : * previously recombined all set of factors with less than rec elts */
494 : static GEN
495 126 : LLL_cmbf(GEN P, GEN famod, GEN p, GEN pa, GEN bound, long a, long rec)
496 : {
497 126 : const long N0 = 1; /* # of traces added at each step */
498 126 : double BitPerFactor = 0.4; /* nb bits in p^(a-b) / modular factor */
499 126 : long i,j,tmax,n0,C, dP = degpol(P);
500 126 : double logp = log((double)itos(p)), LOGp2 = M_LN2/logp;
501 126 : double b0 = log((double)dP*2) / logp, logBr;
502 : GEN lP, Br, Bnorm, Tra, T2, TT, CM_L, m, list, ZERO;
503 : pari_sp av, av2;
504 126 : long ti_LLL = 0, ti_CF = 0;
505 :
506 126 : lP = absi_shallow(leading_coeff(P));
507 126 : if (equali1(lP)) lP = NULL;
508 126 : Br = root_bound(P);
509 126 : if (lP) Br = mulii(lP, Br);
510 126 : logBr = gtodouble(glog(Br, DEFAULTPREC)) / logp;
511 :
512 126 : n0 = lg(famod) - 1;
513 126 : C = (long)ceil( sqrt(N0 * n0 / 4.) ); /* > 1 */
514 126 : Bnorm = dbltor(n0 * (C*C + N0*n0/4.) * 1.00001);
515 126 : ZERO = zeromat(n0, N0);
516 :
517 126 : av = avma;
518 126 : TT = cgetg(n0+1, t_VEC);
519 126 : Tra = cgetg(n0+1, t_MAT);
520 2401 : for (i=1; i<=n0; i++)
521 : {
522 2275 : TT[i] = 0;
523 2275 : gel(Tra,i) = cgetg(N0+1, t_COL);
524 : }
525 126 : CM_L = scalarmat_s(C, n0);
526 : /* tmax = current number of traces used (and computed so far) */
527 126 : for (tmax = 0;; tmax += N0)
528 469 : {
529 595 : long b, bmin, bgood, delta, tnew = tmax + N0, r = lg(CM_L)-1;
530 : GEN M_L, q, CM_Lp, oldCM_L;
531 595 : int first = 1;
532 : pari_timer ti2, TI;
533 :
534 595 : bmin = (long)ceil(b0 + tnew*logBr);
535 595 : if (DEBUGLEVEL>2)
536 0 : err_printf("\nLLL_cmbf: %ld potential factors (tmax = %ld, bmin = %ld)\n",
537 : r, tmax, bmin);
538 :
539 : /* compute Newton sums (possibly relifting first) */
540 595 : if (a <= bmin)
541 : {
542 14 : a = (long)ceil(bmin + 3*N0*logBr) + 1; /* enough for 3 more rounds */
543 14 : a = (long)next2pow((ulong)a);
544 :
545 14 : pa = powiu(p,a);
546 14 : famod = ZpX_liftfact(P, famod, pa, p, a);
547 266 : for (i=1; i<=n0; i++) TT[i] = 0;
548 : }
549 11557 : for (i=1; i<=n0; i++)
550 : {
551 10962 : GEN p1 = gel(Tra,i);
552 10962 : GEN p2 = polsym_gen(gel(famod,i), gel(TT,i), tnew, NULL, pa);
553 10962 : gel(TT,i) = p2;
554 10962 : p2 += 1+tmax; /* ignore traces number 0...tmax */
555 21924 : for (j=1; j<=N0; j++) gel(p1,j) = gel(p2,j);
556 10962 : if (lP)
557 : { /* make Newton sums integral */
558 1848 : GEN lPpow = powiu(lP, tmax);
559 3696 : for (j=1; j<=N0; j++)
560 : {
561 1848 : lPpow = mulii(lPpow,lP);
562 1848 : gel(p1,j) = mulii(gel(p1,j), lPpow);
563 : }
564 : }
565 : }
566 :
567 : /* compute truncation parameter */
568 595 : if (DEBUGLEVEL>2) { timer_start(&ti2); timer_start(&TI); }
569 595 : oldCM_L = CM_L;
570 595 : av2 = avma;
571 595 : delta = b = 0; /* -Wall */
572 1071 : AGAIN:
573 1071 : M_L = Q_div_to_int(CM_L, utoipos(C));
574 1071 : T2 = centermod( ZM_mul(Tra, M_L), pa );
575 1071 : if (first)
576 : { /* initialize lattice, using few p-adic digits for traces */
577 595 : double t = gexpo(T2) - maxdd(32.0, BitPerFactor*r);
578 595 : bgood = (long) (t * LOGp2);
579 595 : b = maxss(bmin, bgood);
580 595 : delta = a - b;
581 : }
582 : else
583 : { /* add more p-adic digits and continue reduction */
584 476 : long b0 = (long)(gexpo(T2) * LOGp2);
585 476 : if (b0 < b) b = b0;
586 476 : b = maxss(b-delta, bmin);
587 476 : if (b - delta/2 < bmin) b = bmin; /* near there. Go all the way */
588 : }
589 :
590 1071 : q = powiu(p, b);
591 1071 : m = vconcat( CM_L, gdivround(T2, q) );
592 1071 : if (first)
593 : {
594 595 : GEN P1 = scalarmat(powiu(p, a-b), N0);
595 595 : first = 0;
596 595 : m = shallowconcat( m, vconcat(ZERO, P1) );
597 : /* [ C M_L 0 ]
598 : * m = [ ] square matrix
599 : * [ T2' p^(a-b) I_N0 ] T2' = Tra * M_L truncated
600 : */
601 : }
602 :
603 1071 : CM_L = LLL_check_progress(Bnorm, n0, m, b == bmin, /*dbg:*/ &ti_LLL);
604 1071 : if (DEBUGLEVEL>2)
605 0 : err_printf("LLL_cmbf: (a,b) =%4ld,%4ld; r =%3ld -->%3ld, time = %ld\n",
606 0 : a,b, lg(m)-1, CM_L? lg(CM_L)-1: 1, timer_delay(&TI));
607 1071 : if (!CM_L) { list = mkvec(P); break; }
608 987 : if (b > bmin)
609 : {
610 476 : CM_L = gc_GEN(av2, CM_L);
611 476 : goto AGAIN;
612 : }
613 511 : if (DEBUGLEVEL>2) timer_printf(&ti2, "for this block of traces");
614 :
615 511 : i = lg(CM_L) - 1;
616 511 : if (i == r && ZM_equal(CM_L, oldCM_L))
617 : {
618 63 : CM_L = oldCM_L;
619 63 : set_avma(av2); continue;
620 : }
621 :
622 448 : CM_Lp = FpM_image(CM_L, utoipos(27449)); /* inexpensive test */
623 448 : if (lg(CM_Lp) != lg(CM_L))
624 : {
625 7 : if (DEBUGLEVEL>2) err_printf("LLL_cmbf: rank decrease\n");
626 7 : CM_L = ZM_hnf(CM_L);
627 : }
628 :
629 448 : if (i <= r && i*rec < n0)
630 : {
631 : pari_timer ti;
632 140 : if (DEBUGLEVEL>2) timer_start(&ti);
633 140 : list = chk_factors(P, Q_div_to_int(CM_L,utoipos(C)), bound, famod, pa);
634 140 : if (DEBUGLEVEL>2) ti_CF += timer_delay(&ti);
635 140 : if (list) break;
636 98 : if (DEBUGLEVEL>2) err_printf("LLL_cmbf: chk_factors failed");
637 : }
638 406 : CM_L = gc_GEN(av2, CM_L);
639 406 : if (gc_needed(av,1))
640 : {
641 0 : if(DEBUGMEM>1) pari_warn(warnmem,"LLL_cmbf");
642 0 : (void)gc_all(av, 5, &CM_L, &TT, &Tra, &famod, &pa);
643 : }
644 : }
645 126 : if (DEBUGLEVEL>2)
646 0 : err_printf("* Time LLL: %ld\n* Time Check Factor: %ld\n",ti_LLL,ti_CF);
647 126 : return list;
648 : }
649 :
650 : /* Find a,b minimal such that A < q^a, B < q^b, 1 << q^(a-b) < 2^31 */
651 : static int
652 50785 : cmbf_precs(GEN q, GEN A, GEN B, long *pta, long *ptb, GEN *qa, GEN *qb)
653 : {
654 50785 : long a,b,amin,d = (long)(31 * M_LN2/gtodouble(glog(q,DEFAULTPREC)) - 1e-5);
655 50785 : int fl = 0;
656 :
657 50785 : b = logintall(B, q, qb) + 1;
658 50785 : *qb = mulii(*qb, q);
659 50785 : amin = b + d;
660 50785 : if (gcmp(powiu(q, amin), A) <= 0)
661 : {
662 14260 : a = logintall(A, q, qa) + 1;
663 14260 : *qa = mulii(*qa, q);
664 14260 : b = a - d; *qb = powiu(q, b);
665 : }
666 : else
667 : { /* not enough room */
668 36525 : a = amin; *qa = powiu(q, a);
669 36525 : fl = 1;
670 : }
671 50785 : if (DEBUGLEVEL > 3) {
672 0 : err_printf("S_2 bound: %Ps^%ld\n", q,b);
673 0 : err_printf("coeff bound: %Ps^%ld\n", q,a);
674 : }
675 50785 : *pta = a;
676 50785 : *ptb = b; return fl;
677 : }
678 :
679 : /* use van Hoeij's knapsack algorithm */
680 : static GEN
681 50785 : combine_factors(GEN target, GEN famod, GEN p, long klim)
682 : {
683 : GEN la, B, A, res, L, pa, pb, listmod;
684 50785 : long a,b, l, maxK, n = degpol(target);
685 : int done;
686 : pari_timer T;
687 :
688 50785 : A = factor_bound(target);
689 :
690 50785 : la = absi_shallow(leading_coeff(target));
691 50785 : B = mului(n, sqri(mulii(la, root_bound(target)))); /* = bound for S_2 */
692 :
693 50785 : (void)cmbf_precs(p, A, B, &a, &b, &pa, &pb);
694 :
695 50785 : if (DEBUGLEVEL>2) timer_start(&T);
696 50785 : famod = ZpX_liftfact(target, famod, pa, p, a);
697 50784 : if (DEBUGLEVEL>2) timer_printf(&T, "Hensel lift (mod %Ps^%ld)", p,a);
698 50784 : L = cmbf(target, famod, A, p, a, b, klim, &maxK, &done);
699 50785 : if (DEBUGLEVEL>2) timer_printf(&T, "Naive recombination");
700 :
701 50785 : res = gel(L,1);
702 50785 : listmod = gel(L,2); l = lg(listmod)-1;
703 50785 : famod = gel(listmod,l);
704 50785 : if (maxK > 0 && lg(famod)-1 > 2*maxK)
705 : {
706 126 : if (l!=1) A = factor_bound(gel(res,l));
707 126 : if (DEBUGLEVEL > 4) err_printf("last factor still to be checked\n");
708 126 : L = LLL_cmbf(gel(res,l), famod, p, pa, A, a, maxK);
709 126 : if (DEBUGLEVEL>2) timer_printf(&T,"Knapsack");
710 : /* remove last elt, possibly unfactored. Add all new ones. */
711 126 : setlg(res, l); res = shallowconcat(res, L);
712 : }
713 50785 : return res;
714 : }
715 :
716 : /* Assume 'a' a squarefree ZX; return 0 if no root (fl=1) / irreducible (fl=0).
717 : * Otherwise return prime p such that a mod p has fewest roots / factors */
718 : static ulong
719 2083173 : pick_prime(GEN a, long fl, pari_timer *T)
720 : {
721 2083173 : pari_sp av = avma, av1;
722 2083173 : const long MAXNP = 7, da = degpol(a);
723 2083174 : long nmax = da+1, np;
724 2083174 : ulong chosenp = 0;
725 2083174 : GEN lead = gel(a,da+2);
726 : forprime_t S;
727 2083174 : if (equali1(lead)) lead = NULL;
728 2083174 : u_forprime_init(&S, 2, ULONG_MAX);
729 2083183 : av1 = avma;
730 8042252 : for (np = 0; np < MAXNP; set_avma(av1))
731 : {
732 7931336 : ulong p = u_forprime_next(&S);
733 : long nfacp;
734 : GEN z;
735 :
736 7931328 : if (!p) pari_err_OVERFLOW("DDF [out of small primes]");
737 7931328 : if (lead && !umodiu(lead,p)) continue;
738 7864876 : z = ZX_to_Flx(a, p);
739 7864870 : if (!Flx_is_squarefree(z, p)) continue;
740 :
741 4874429 : if (fl==1)
742 : {
743 4278446 : nfacp = Flx_nbroots(z, p);
744 4278447 : if (!nfacp) { chosenp = 0; break; } /* no root */
745 : }
746 595983 : else if(fl==0)
747 : {
748 595234 : nfacp = Flx_nbfact(z, p);
749 595238 : if (nfacp == 1) { chosenp = 0; break; } /* irreducible */
750 : } else
751 : {
752 749 : GEN f = gel(Flx_degfact(z, p),1);
753 749 : nfacp = lg(f)-1;
754 749 : if (f[1] > fl) { chosenp = 0; break; } /* no small factors */
755 : }
756 2902162 : if (DEBUGLEVEL>4)
757 0 : err_printf("...tried prime %3lu (%-3ld %s). Time = %ld\n",
758 : p, nfacp, fl==1? "roots": "factors", timer_delay(T));
759 2902171 : if (nfacp < nmax)
760 : {
761 1087461 : nmax = nfacp; chosenp = p;
762 1087461 : if (da > 100 && nmax < 5) break; /* large degree, few factors. Enough */
763 : }
764 2902164 : np++;
765 : }
766 2083189 : return gc_ulong(av, chosenp);
767 : }
768 :
769 : /* Assume A a squarefree ZX; return the vector of its rational roots */
770 : static GEN
771 1876570 : DDF_roots(GEN A)
772 : {
773 : GEN p, lc, lcpol, z, pe, pes2, bound;
774 : long i, m, e, lz;
775 : ulong pp;
776 : pari_sp av;
777 : pari_timer T;
778 :
779 1876570 : if (DEBUGLEVEL>2) timer_start(&T);
780 1876570 : pp = pick_prime(A, 1, &T);
781 1876571 : if (!pp) return cgetg(1,t_COL); /* no root */
782 60118 : p = utoipos(pp);
783 60118 : lc = leading_coeff(A);
784 60118 : if (is_pm1(lc))
785 52606 : { lc = NULL; lcpol = A; }
786 : else
787 7512 : { lc = absi_shallow(lc); lcpol = ZX_Z_mul(A, lc); }
788 60118 : bound = root_bound(A); if (lc) bound = mulii(lc, bound);
789 60118 : e = logintall(addiu(shifti(bound, 1), 1), p, &pe) + 1;
790 60118 : pe = mulii(pe, p);
791 60118 : pes2 = shifti(pe, -1);
792 60118 : if (DEBUGLEVEL>2) timer_printf(&T, "Root bound");
793 60118 : av = avma;
794 60118 : z = ZpX_roots(A, p, e); lz = lg(z);
795 60118 : z = deg1_from_roots(z, varn(A));
796 60118 : if (DEBUGLEVEL>2) timer_printf(&T, "Hensel lift (mod %lu^%ld)", pp,e);
797 130582 : for (m=1, i=1; i < lz; i++)
798 : {
799 70464 : GEN q, r, y = gel(z,i);
800 70464 : if (lc) y = ZX_Z_mul(y, lc);
801 70464 : y = centermod_i(y, pe, pes2);
802 70464 : if (! (q = ZX_divides(lcpol, y)) ) continue;
803 :
804 25974 : lcpol = q;
805 25974 : r = negi( constant_coeff(y) );
806 25974 : if (lc) {
807 8277 : r = gdiv(r,lc);
808 8277 : lcpol = Q_primpart(lcpol);
809 8277 : lc = absi_shallow( leading_coeff(lcpol) );
810 8277 : if (is_pm1(lc)) lc = NULL; else lcpol = ZX_Z_mul(lcpol, lc);
811 : }
812 25974 : gel(z,m++) = r;
813 25974 : if (gc_needed(av,2))
814 : {
815 0 : if (DEBUGMEM>1) pari_warn(warnmem,"DDF_roots, m = %ld", m);
816 0 : (void)gc_all(av, lc? 3:2, &z, &lcpol, &lc);
817 : }
818 : }
819 60118 : if (DEBUGLEVEL>2) timer_printf(&T, "Recombination");
820 60118 : setlg(z, m); return z;
821 : }
822 :
823 : /* Assume a squarefree ZX, deg(a) > 0, return rational factors.
824 : * In fact, a(0) != 0 but we don't use this
825 : * if dmax>0, Only look for factor of degree at most dmax */
826 : GEN
827 206603 : ZX_DDF_max(GEN a, long dmax)
828 : {
829 : GEN ap, prime, famod, z;
830 206603 : long ti = 0;
831 206603 : ulong p = 0;
832 206603 : pari_sp av = avma;
833 : pari_timer T, T2;
834 :
835 206603 : if (DEBUGLEVEL>2) { timer_start(&T); timer_start(&T2); }
836 206603 : p = pick_prime(a, dmax, &T2);
837 206605 : if (!p) return mkvec(a);
838 50785 : prime = utoipos(p);
839 50785 : ap = Flx_normalize(ZX_to_Flx(a, p), p);
840 50785 : famod = gel(Flx_factor(ap, p), 1);
841 50785 : if (DEBUGLEVEL>2)
842 : {
843 0 : if (DEBUGLEVEL>4) timer_printf(&T2, "splitting mod p = %lu", p);
844 0 : ti = timer_delay(&T);
845 0 : err_printf("Time setup: %ld\n", ti);
846 : }
847 50785 : z = combine_factors(a, FlxV_to_ZXV(famod), prime, degpol(a)-1);
848 50785 : if (DEBUGLEVEL>2)
849 0 : err_printf("Total Time: %ld\n===========\n", ti + timer_delay(&T));
850 50785 : return gc_GEN(av, z);
851 : }
852 :
853 : /* Distinct Degree Factorization (deflating first)
854 : * Assume x squarefree, degree(x) > 0, x(0) != 0 */
855 : GEN
856 153509 : ZX_DDF(GEN x)
857 : {
858 : GEN L;
859 : long m;
860 153509 : if (DEBUGLEVEL>2)
861 0 : err_printf("ZX_DDF: factoring pol of deg %ld, %ld bits\n",degpol(x),gexpo(x));
862 153509 : x = ZX_deflate_max(x, &m);
863 153507 : L = ZX_DDF_max(x,0);
864 153508 : if (m > 1)
865 : {
866 50735 : GEN e, v, fa = factoru(m);
867 : long i,j,k, l;
868 :
869 50735 : e = gel(fa,2); k = 0;
870 50735 : fa= gel(fa,1); l = lg(fa);
871 101792 : for (i=1; i<l; i++) k += e[i];
872 50735 : v = cgetg(k+1, t_VECSMALL); k = 1;
873 101794 : for (i=1; i<l; i++)
874 103617 : for (j=1; j<=e[i]; j++) v[k++] = fa[i];
875 103294 : for (k--; k; k--)
876 : {
877 52559 : GEN L2 = cgetg(1,t_VEC);
878 105481 : for (i=1; i < lg(L); i++)
879 52923 : L2 = shallowconcat(L2, ZX_DDF_max(RgX_inflate(gel(L,i), v[k]),0));
880 52558 : L = L2;
881 : }
882 : }
883 153508 : return L;
884 : }
885 :
886 : /* SquareFree Factorization in Z[X] (char 0 is enough, if ZX_gcd -> RgX_gcd)
887 : * f = prod Q[i]^E[i], E[1] < E[2] < ..., and Q[i] squarefree and coprime.
888 : * Return Q, set *pE = E. For efficiency, caller should have used ZX_valrem
889 : * so that f(0) != 0 */
890 : GEN
891 333726 : ZX_squff(GEN f, GEN *pE)
892 : {
893 : GEN T, V, P, E;
894 333726 : long i, k, n = 1 + degpol(f);
895 :
896 333726 : if (signe(leading_coeff(f)) < 0) f = ZX_neg(f);
897 333726 : E = cgetg(n, t_VECSMALL);
898 333726 : P = cgetg(n, t_COL);
899 333726 : T = ZX_gcd_all(f, ZX_deriv(f), &V);
900 333726 : for (k = i = 1;; k++)
901 3655 : {
902 337381 : GEN W = ZX_gcd_all(T,V, &T); /* V and W are squarefree */
903 337381 : long dW = degpol(W), dV = degpol(V);
904 : /* f = prod_i T_i^{e_i}
905 : * W = prod_{i: e_i > k} T_i,
906 : * V = prod_{i: e_i >= k} T_i,
907 : * T = prod_{i: e_i > k} T_i^{e_i - k} */
908 337381 : if (!dW)
909 : {
910 333726 : if (dV) { gel(P,i) = Q_primpart(V); E[i] = k; i++; }
911 333726 : break;
912 : }
913 3655 : if (dW == dV)
914 : {
915 : GEN U;
916 1659 : while ( (U = ZX_divides(T, V)) ) { k++; T = U; }
917 : }
918 : else
919 : {
920 2388 : gel(P,i) = Q_primpart(RgX_div(V,W));
921 2388 : E[i] = k; i++; V = W;
922 : }
923 : }
924 333726 : setlg(P,i);
925 333726 : setlg(E,i); *pE = E; return P;
926 : }
927 :
928 : static GEN
929 51555 : fact_from_DDF(GEN Q, GEN E, long n)
930 : {
931 51555 : GEN v,w, y = cgetg(3, t_MAT);
932 51555 : long i,j,k, l = lg(Q);
933 :
934 51555 : v = cgetg(n+1, t_COL); gel(y,1) = v;
935 51555 : w = cgetg(n+1, t_COL); gel(y,2) = w;
936 104909 : for (k = i = 1; i < l; i++)
937 : {
938 53354 : GEN L = gel(Q,i), e = utoipos(E[i]);
939 53354 : long J = lg(L);
940 131831 : for (j = 1; j < J; j++,k++)
941 : {
942 78477 : gel(v,k) = ZX_copy(gel(L,j));
943 78477 : gel(w,k) = e;
944 : }
945 : }
946 51555 : return y;
947 : }
948 :
949 : /* Factor T in Z[x] */
950 : static GEN
951 51562 : ZX_factor_i(GEN T)
952 : {
953 : GEN Q, E, y;
954 : long n, i, l, v;
955 :
956 51562 : if (!signe(T)) return prime_fact(T);
957 51555 : v = ZX_valrem(T, &T);
958 51555 : Q = ZX_squff(T, &E); l = lg(Q);
959 103355 : for (i = 1, n = 0; i < l; i++)
960 : {
961 51800 : gel(Q,i) = ZX_DDF(gel(Q,i));
962 51800 : n += lg(gel(Q,i)) - 1;
963 : }
964 51555 : if (v)
965 : {
966 1554 : Q = vec_append(Q, mkvec(pol_x(varn(T))));
967 1554 : E = vecsmall_append(E, v); n++;
968 : }
969 51555 : y = fact_from_DDF(Q, E, n);
970 51555 : return sort_factor_pol(y, cmpii);
971 : }
972 : GEN
973 50757 : ZX_factor(GEN x)
974 : {
975 50757 : pari_sp av = avma;
976 50757 : return gc_upto(av, ZX_factor_i(x));
977 : }
978 : GEN
979 805 : QX_factor(GEN x)
980 : {
981 805 : pari_sp av = avma;
982 805 : return gc_upto(av, ZX_factor_i(Q_primpart(x)));
983 : }
984 :
985 : long
986 102643 : ZX_is_irred(GEN x)
987 : {
988 102643 : pari_sp av = avma;
989 102643 : long l = lg(x);
990 : GEN y;
991 102643 : if (l <= 3) return 0; /* degree < 1 */
992 102643 : if (l == 4) return 1; /* degree 1 */
993 99056 : if (ZX_val(x)) return 0;
994 98832 : if (!ZX_is_squarefree(x)) return 0;
995 98678 : y = ZX_DDF(x); set_avma(av);
996 98677 : return (lg(y) == 2);
997 : }
998 :
999 : GEN
1000 1876569 : nfrootsQ(GEN x)
1001 : {
1002 1876569 : pari_sp av = avma;
1003 : GEN z;
1004 : long val;
1005 :
1006 1876569 : if (typ(x)!=t_POL) pari_err_TYPE("nfrootsQ",x);
1007 1876569 : if (!signe(x)) pari_err_ROOTS0("nfrootsQ");
1008 1876569 : x = Q_primpart(x);
1009 1876568 : RgX_check_ZX(x,"nfrootsQ");
1010 1876571 : val = ZX_valrem(x, &x);
1011 1876575 : z = DDF_roots( ZX_radical(x) );
1012 1876573 : if (val) z = vec_append(z, gen_0);
1013 1876573 : return gc_upto(av, sort(z));
1014 : }
1015 :
1016 : /************************************************************************
1017 : * GCD OVER Z[X] / Q[X] *
1018 : ************************************************************************/
1019 : int
1020 198734 : ZX_is_squarefree(GEN x)
1021 : {
1022 198734 : pari_sp av = avma;
1023 : GEN d;
1024 : long m;
1025 198734 : if (lg(x) == 2) return 0;
1026 198734 : m = ZX_deflate_order(x);
1027 198741 : if (m > 1)
1028 : {
1029 86194 : if (!signe(gel(x,2))) return 0;
1030 85956 : x = RgX_deflate(x, m);
1031 : }
1032 198506 : d = ZX_gcd(x,ZX_deriv(x));
1033 198510 : return gc_bool(av, lg(d) == 3);
1034 : }
1035 :
1036 : static int
1037 123595 : ZX_gcd_filter(GEN *pt_A, GEN *pt_P)
1038 : {
1039 123595 : GEN A = *pt_A, P = *pt_P;
1040 123595 : long i, j, l = lg(A), n = 1, d = degpol(gel(A,1));
1041 : GEN B, Q;
1042 253124 : for (i=2; i<l; i++)
1043 : {
1044 129529 : long di = degpol(gel(A,i));
1045 129529 : if (di==d) n++;
1046 36 : else if (d > di)
1047 36 : { n=1; d = di; }
1048 : }
1049 123595 : if (n == l-1)
1050 123559 : return 0;
1051 36 : B = cgetg(n+1, t_VEC);
1052 36 : Q = cgetg(n+1, typ(P));
1053 156 : for (i=1, j=1; i<l; i++)
1054 : {
1055 120 : if (degpol(gel(A,i))==d)
1056 : {
1057 84 : gel(B,j) = gel(A,i);
1058 84 : Q[j] = P[i];
1059 84 : j++;
1060 : }
1061 : }
1062 36 : *pt_A = B; *pt_P = Q; return 1;
1063 : }
1064 :
1065 : static GEN
1066 3370733 : ZX_gcd_Flx(GEN a, GEN b, ulong g, ulong p)
1067 : {
1068 3370733 : GEN H = Flx_gcd(a, b, p);
1069 3370733 : if (!g)
1070 3338339 : return Flx_normalize(H, p);
1071 : else
1072 : {
1073 32394 : ulong t = Fl_mul(g, Fl_inv(Flx_lead(H), p), p);
1074 32394 : return Flx_Fl_mul(H, t, p);
1075 : }
1076 : }
1077 :
1078 : static GEN
1079 3359087 : ZX_gcd_slice(GEN A, GEN B, GEN g, GEN P, GEN *mod)
1080 : {
1081 3359087 : pari_sp av = avma;
1082 3359087 : long i, n = lg(P)-1;
1083 : GEN H, T;
1084 3359087 : if (n == 1)
1085 : {
1086 3352933 : ulong p = uel(P,1), gp = g ? umodiu(g, p): 0;
1087 3352933 : GEN a = ZX_to_Flx(A, p), b = ZX_to_Flx(B, p);
1088 3352933 : GEN Hp = ZX_gcd_Flx(a, b, gp, p);
1089 3352933 : H = gc_upto(av, Flx_to_ZX(Hp));
1090 3352932 : *mod = utoi(p);
1091 3352932 : return H;
1092 : }
1093 6154 : T = ZV_producttree(P);
1094 6154 : A = ZX_nv_mod_tree(A, P, T);
1095 6154 : B = ZX_nv_mod_tree(B, P, T);
1096 6154 : g = g ? Z_ZV_mod_tree(g, P, T): NULL;
1097 6154 : H = cgetg(n+1, t_VEC);
1098 23954 : for(i=1; i <= n; i++)
1099 : {
1100 17800 : ulong p = P[i];
1101 17800 : GEN a = gel(A,i), b = gel(B,i);
1102 17800 : gel(H,i) = ZX_gcd_Flx(a, b, g? g[i]: 0, p);
1103 : }
1104 6154 : if (ZX_gcd_filter(&H, &P))
1105 12 : T = ZV_producttree(P);
1106 6154 : H = nxV_chinese_center_tree(H, P, T, ZV_chinesetree(P, T));
1107 6154 : *mod = gmael(T, lg(T)-1, 1); return gc_all(av, 2, &H, mod);
1108 : }
1109 :
1110 : GEN
1111 3359087 : ZX_gcd_worker(GEN P, GEN A, GEN B, GEN g)
1112 : {
1113 3359087 : GEN V = cgetg(3, t_VEC);
1114 3359087 : gel(V,1) = ZX_gcd_slice(A, B, equali1(g)? NULL: g , P, &gel(V,2));
1115 3359086 : return V;
1116 : }
1117 :
1118 : static GEN
1119 117441 : ZX_gcd_chinese(GEN A, GEN P, GEN *mod)
1120 : {
1121 117441 : ZX_gcd_filter(&A, &P);
1122 117441 : return nxV_chinese_center(A, P, mod);
1123 : }
1124 :
1125 : GEN
1126 14467765 : ZX_gcd_all(GEN A, GEN B, GEN *Anew)
1127 : {
1128 14467765 : pari_sp av = avma;
1129 14467765 : long k, valH, valA, valB, vA = varn(A), dA = degpol(A), dB = degpol(B);
1130 14467730 : GEN worker, c, cA, cB, g, Ag, Bg, H = NULL, mod = gen_1, R;
1131 : GEN Ap, Bp, Hp;
1132 : forprime_t S;
1133 : ulong pp;
1134 14467730 : if (dA < 0) { if (Anew) *Anew = pol_0(vA); return ZX_copy(B); }
1135 14467709 : if (dB < 0) { if (Anew) *Anew = pol_1(vA); return ZX_copy(A); }
1136 14466568 : A = Q_primitive_part(A, &cA);
1137 14466714 : B = Q_primitive_part(B, &cB);
1138 14466646 : valA = ZX_valrem(A, &A); dA -= valA;
1139 14466649 : valB = ZX_valrem(B, &B); dB -= valB;
1140 14466701 : valH = minss(valA, valB);
1141 14466722 : valA -= valH; /* valuation(Anew) */
1142 14466722 : c = (cA && cB)? gcdii(cA, cB): NULL; /* content(gcd) */
1143 14466731 : if (!dA || !dB)
1144 : {
1145 7338096 : if (Anew) *Anew = RgX_shift_shallow(A, valA);
1146 7338096 : return monomial(c? c: gen_1, valH, vA);
1147 : }
1148 7128635 : g = gcdii(leading_coeff(A), leading_coeff(B)); /* multiple of lead(gcd) */
1149 7128564 : if (is_pm1(g)) {
1150 6919221 : g = NULL;
1151 6919221 : Ag = A;
1152 6919221 : Bg = B;
1153 : } else {
1154 209338 : Ag = ZX_Z_mul(A,g);
1155 209338 : Bg = ZX_Z_mul(B,g);
1156 : }
1157 7128559 : init_modular_big(&S);
1158 : do {
1159 7128668 : pp = u_forprime_next(&S);
1160 7128665 : Ap = ZX_to_Flx(Ag, pp);
1161 7128691 : Bp = ZX_to_Flx(Bg, pp);
1162 7128680 : } while (degpol(Ap) != dA || degpol(Bp) != dB);
1163 7128665 : if (degpol(Flx_gcd(Ap, Bp, pp)) == 0)
1164 : {
1165 3887424 : if (Anew) *Anew = RgX_shift_shallow(A, valA);
1166 3887424 : return monomial(c? c: gen_1, valH, vA);
1167 : }
1168 3241204 : worker = snm_closure(is_entry("_ZX_gcd_worker"), mkvec3(A, B, g? g: gen_1));
1169 3241204 : av = avma;
1170 3356686 : for (k = 1; ;k *= 2)
1171 : {
1172 3356686 : gen_inccrt_i("ZX_gcd", worker, g, (k+1)>>1, 0, &S, &H, &mod, ZX_gcd_chinese, NULL);
1173 3356686 : (void)gc_all(av, 2, &H, &mod);
1174 3356686 : Hp = ZX_to_Flx(H, pp);
1175 3356686 : if (lgpol(Flx_rem(Ap, Hp, pp)) || lgpol(Flx_rem(Bp, Hp, pp))) continue;
1176 3245269 : if (!ZX_divides(Bg, H)) continue;
1177 3241231 : R = ZX_divides(Ag, H);
1178 3241231 : if (R) break;
1179 : }
1180 : /* lead(H) = g */
1181 3241201 : if (g) H = Q_primpart(H);
1182 3241201 : if (c) H = ZX_Z_mul(H,c);
1183 3241201 : if (DEBUGLEVEL>5) err_printf("done\n");
1184 3241201 : if (Anew) *Anew = RgX_shift_shallow(R, valA);
1185 3241201 : return valH? RgX_shift_shallow(H, valH): H;
1186 : }
1187 :
1188 : #if 0
1189 : /* ceil( || p ||_oo / lc(p) ) */
1190 : static GEN
1191 : maxnorm(GEN p)
1192 : {
1193 : long i, n = degpol(p), av = avma;
1194 : GEN x, m = gen_0;
1195 :
1196 : p += 2;
1197 : for (i=0; i<n; i++)
1198 : {
1199 : x = gel(p,i);
1200 : if (abscmpii(x,m) > 0) m = x;
1201 : }
1202 : m = divii(m, gel(p,n));
1203 : return gc_INT(av, addiu(absi_shallow(m),1));
1204 : }
1205 : #endif
1206 :
1207 : GEN
1208 10624527 : ZX_gcd(GEN A, GEN B)
1209 : {
1210 10624527 : pari_sp av = avma;
1211 10624527 : return gc_GEN(av, ZX_gcd_all(A,B,NULL));
1212 : }
1213 :
1214 : GEN
1215 3167884 : ZX_radical(GEN A) { GEN B; (void)ZX_gcd_all(A,ZX_deriv(A),&B); return B; }
1216 :
1217 : static GEN
1218 19183 : _gcd(GEN a, GEN b)
1219 : {
1220 19183 : if (!a) a = gen_1;
1221 19183 : if (!b) b = gen_1;
1222 19183 : return Q_gcd(a,b);
1223 : }
1224 : /* A0 and B0 in Q[X] */
1225 : GEN
1226 18994 : QX_gcd(GEN A0, GEN B0)
1227 : {
1228 : GEN a, b, D;
1229 18994 : pari_sp av = avma, av2;
1230 :
1231 18994 : D = ZX_gcd(Q_primitive_part(A0, &a), Q_primitive_part(B0, &b));
1232 18994 : av2 = avma; a = _gcd(a,b);
1233 18994 : if (isint1(a)) set_avma(av2); else D = ZX_Q_mul(D, a);
1234 18994 : return gc_upto(av, D);
1235 : }
1236 :
1237 : /***************************************************************************
1238 : *** ***
1239 : *** ZXk/QXk ***
1240 : *** ***
1241 : ***************************************************************************/
1242 :
1243 : /* ZXk/QXk: multivariate polynomials in Z[X_1,...,X_k] and Q[X_1,...,X_k] */
1244 :
1245 : INLINE GEN
1246 4961344 : ZXk_renormalize(GEN x, long lx) { return ZXX_renormalize(x,lx); }
1247 :
1248 : int
1249 0 : Rg_is_QXk(GEN z)
1250 : {
1251 0 : long i, t = typ(z), l = lg(z);
1252 0 : if (t==t_INT || t==t_FRAC) return 1;
1253 0 : if (t!=t_POL) return 0;
1254 0 : for (i = 2; i < l; i++)
1255 0 : if (!Rg_is_QXk(gel(z,i))) return 0;
1256 0 : return 1;
1257 : }
1258 :
1259 : static GEN
1260 3214540 : centeri2n(GEN z, long n, GEN N)
1261 : {
1262 3214540 : pari_sp av = avma;
1263 3214540 : z = remi2n(z, n);
1264 3214540 : if (expi(z)<n-1) return z;
1265 1253 : if (signe(z)<0) z = addii(z, N);
1266 763 : else z = subii(z, N);
1267 1253 : return gc_INT(av, z);
1268 : }
1269 :
1270 : static GEN
1271 4910915 : ZXk_center2n(GEN z, long n, GEN N)
1272 : {
1273 4910915 : if (typ(z) == t_INT)
1274 3214540 : return centeri2n(z, n, N);
1275 : else
1276 : {
1277 : long i,l;
1278 1696375 : GEN x = cgetg_copy(z, &l);
1279 1696375 : x[1] = z[1];
1280 3405042 : for (i = 2; i < l; i++)
1281 1708667 : gel(x,i) = ZXk_center2n(gel(z,i), n, N);
1282 1696375 : return ZXk_renormalize(x, l);
1283 : }
1284 : }
1285 :
1286 : static GEN ZXk_gcd_i(GEN A, GEN B);
1287 : static GEN
1288 3445402 : ZXk_content_shallow(GEN x)
1289 : {
1290 3445402 : long i, l = lg(x);
1291 : GEN c;
1292 3445402 : if (typ(x)==t_INT) return x;
1293 3445402 : if (!signe(x)) return gen_0;
1294 3445402 : c = gel(x, 2);
1295 3445402 : if (gequal1(c)) return gen_1;
1296 6098970 : for (i = 3; i < l; i++)
1297 : {
1298 3558036 : c = simplify_shallow(ZXk_gcd_i(c, gel(x,i)));
1299 3558036 : if (gequal1(c)) return gen_1;
1300 : }
1301 2540934 : return c;
1302 : }
1303 :
1304 : static GEN ZXk_divexact_s(GEN A, GEN B);
1305 :
1306 : static GEN
1307 2571335 : ZXkX_ZXk_divexact_s(GEN x, GEN B)
1308 7610944 : { pari_APPLY_ZX(ZXk_divexact_s(gel(x,i), B)); }
1309 :
1310 : static GEN
1311 2528103 : ZXkX_ZXk_divexact(GEN A, GEN B)
1312 : {
1313 2528103 : pari_sp av = avma;
1314 2528103 : return gc_upto(av, ZXkX_ZXk_divexact_s(A, simplify_shallow(B)));
1315 : }
1316 :
1317 : static GEN
1318 2529297 : ZXk_divexact_i(GEN x, GEN y)
1319 : {
1320 2529297 : long dx = degpol(x), dy = degpol(y), dz, i, j;
1321 2529297 : GEN z, y_lead = gel(y,dy+2);
1322 2529297 : if (dx < dy)
1323 0 : return gen_0;
1324 2529297 : dz = dx-dy;
1325 2529297 : z = cgetg(dz+3,t_POL); z[1] = x[1];
1326 2529297 : gel(z,dz+2) = ZXk_divexact_s(gel(x,dx+2), y_lead);
1327 2544858 : for (i=dx-1; i>=dy; i--)
1328 : {
1329 15561 : pari_sp btop = avma;
1330 15561 : GEN p1=gel(x,2+i);
1331 34034 : for (j=i-dy+1; j<=i && j<=dz; j++)
1332 18473 : p1 = gsub(p1, gmul(gel(z,2+j), gel(y,2+i-j)));
1333 15561 : gel(z,2+i-dy) = gc_upto(btop, ZXk_divexact_s(p1, y_lead));
1334 : }
1335 2529297 : return z;
1336 : }
1337 :
1338 : static GEN
1339 7584467 : ZXk_divexact_s(GEN A, GEN B)
1340 : {
1341 7584467 : if (!signe(A)) return gen_0;
1342 5137203 : if (typ(A)==t_INT && typ(B)==t_INT)
1343 2564674 : return diviiexact(A, B);
1344 2572529 : else if (typ(B)==t_INT || varn(A)!=varn(B))
1345 43232 : return ZXkX_ZXk_divexact_s(A, B);
1346 : else
1347 2529297 : return ZXk_divexact_i(A, B);
1348 : }
1349 :
1350 : GEN
1351 0 : ZXk_divexact(GEN A, GEN B)
1352 : {
1353 0 : pari_sp av = avma;
1354 0 : return gc_upto(av, ZXk_divexact_s(A, simplify_shallow(B)));
1355 : }
1356 :
1357 : static GEN ZXk_divides_s(GEN A, GEN B);
1358 :
1359 : static GEN
1360 3264969 : ZXkX_ZXk_divides_s(GEN x, GEN B)
1361 : {
1362 3264969 : pari_sp av = avma;
1363 : long i, l;
1364 3264969 : GEN y = cgetg_copy(x, &l); y[1] = x[1];
1365 3264969 : if (l == 2) return y;
1366 6916728 : for (i=2; i<l; i++)
1367 : {
1368 3651759 : GEN c = ZXk_divides_s(gel(x,i), B);
1369 3651759 : if (!c) return gc_NULL(av);
1370 3651759 : gel(y, i) = c;
1371 : }
1372 3264969 : return ZXk_renormalize(y, l);
1373 : }
1374 :
1375 : static GEN
1376 2789715 : ZXk_divides_i(GEN x, GEN y)
1377 : {
1378 2789715 : pari_sp av = avma, av2;
1379 2789715 : long dx = degpol(x), dy = degpol(y), dz, i, j, c;
1380 2789715 : GEN z, y_lead = gel(y,dy+2);
1381 2789715 : if (dx < dy)
1382 7 : return gen_0;
1383 2789708 : dz = dx-dy;
1384 2789708 : z = cgetg(dz+3,t_POL); z[1] = x[1];
1385 2789708 : gel(z,dz+2) = ZXk_divides(gel(x,dx+2), y_lead);
1386 2789708 : if (!gel(z,dz+2)) return gc_NULL(av);
1387 3042779 : for (i=dx-1; i>=dy; i--)
1388 : {
1389 253085 : pari_sp btop = avma;
1390 253085 : GEN p1 = gel(x,2+i), c;
1391 514178 : for (j=i-dy+1; j<=i && j<=dz; j++)
1392 261093 : p1 = gsub(p1, gmul(gel(z,2+j), gel(y,2+i-j)));
1393 253085 : c = ZXk_divides_s(p1, y_lead);
1394 253085 : if (!c) return gc_NULL(av);
1395 253078 : gel(z,2+i-dy) = gc_upto(btop, c);
1396 : }
1397 2789694 : av2 = avma;
1398 2789694 : c = gc_bool(av2, gequal(gmul(z,y),x));
1399 2789694 : return c ? z: gc_NULL(av);
1400 : }
1401 :
1402 : static GEN
1403 3481571 : dividesii(GEN A, GEN B)
1404 : {
1405 3481571 : GEN r, q = dvmdii(A, B, &r);
1406 3481571 : return signe(r) ? NULL: q;
1407 : }
1408 :
1409 : static GEN
1410 10123448 : ZXk_divides_s(GEN A, GEN B)
1411 : {
1412 10123448 : if (!signe(A)) return gen_0;
1413 9536255 : if (typ(B)==t_INT)
1414 3481571 : return typ(A)==t_INT ? dividesii(A, B)
1415 10227453 : : ZXkX_ZXk_divides_s(A, B);
1416 2790373 : else if (typ(A)==t_INT) return NULL;
1417 : else
1418 : {
1419 2790373 : long c = varncmp(varn(A),varn(B));
1420 2790373 : if (c < 0)
1421 658 : return ZXkX_ZXk_divides_s(A, B);
1422 2789715 : else if (c>0)
1423 0 : return NULL;
1424 : else
1425 2789715 : return ZXk_divides_i(A, B);
1426 : }
1427 : }
1428 :
1429 : GEN
1430 6218604 : ZXk_divides(GEN A, GEN B)
1431 : {
1432 6218604 : pari_sp av = avma;
1433 6218604 : GEN z = ZXk_divides_s(A, simplify_shallow(B));
1434 6218604 : return z ? gc_upto(av, z): z;
1435 : }
1436 :
1437 : static GEN
1438 1714455 : rec(GEN g, long e, GEN N, long v)
1439 : {
1440 1714455 : pari_sp av = avma;
1441 1714455 : long i, d = (gexpo(g)+2*e-1)/e;
1442 1714455 : GEN s = cgetg(d+3,t_POL);
1443 1714455 : s[1] = evalvarn(v);
1444 3202248 : for (i = 0; i <= d; i++)
1445 : {
1446 3202248 : GEN c = ZXk_center2n(g, e, N);
1447 3202248 : gel(s,i+2) = c;
1448 3202248 : g = gmul2n(gsub(g,c),-e);
1449 3202248 : if (!signe(g)) break;
1450 : }
1451 1714455 : s = RgX_renormalize_lg(s,i+3);
1452 1714455 : return gc_GEN(av, s);
1453 : }
1454 :
1455 : static GEN
1456 8689704 : ZXk_gcd_i(GEN A, GEN B)
1457 : {
1458 : pari_sp av;
1459 : long e, v, vc;
1460 : GEN c, cA, cB;
1461 8689704 : if (signe(A)==0) return gcopy(B);
1462 5162982 : if (signe(B)==0) return gcopy(A);
1463 5161141 : if (typ(A) == t_INT) return gcdii(A, typ(B)==t_INT ? B: Q_content(B));
1464 4292688 : if (typ(B) == t_INT) return gcdii(Q_content(A), B);
1465 3439928 : v = varn(A); vc = varncmp(v, varn(B));
1466 3439928 : if (vc < 0) return ZXk_gcd_i(ZXk_content_shallow(A), B);
1467 3428392 : if (vc > 0) return ZXk_gcd_i(A, ZXk_content_shallow(B));
1468 3423408 : if (RgX_is_ZX(A) && RgX_is_ZX(B)) return ZX_gcd(A,B);
1469 1714441 : cA = ZXk_content_shallow(A); if (!gequal1(cA)) A = ZXkX_ZXk_divexact(A, cA);
1470 1714441 : cB = ZXk_content_shallow(B); if (!gequal1(cB)) B = ZXkX_ZXk_divexact(B, cB);
1471 1714441 : c = ZXk_gcd_i(cA, cB); av = avma;
1472 1714441 : e = maxss(3, minss(gexpo(A), gexpo(B)) + 2);
1473 14 : for ( ; ; e++, set_avma(av))
1474 14 : {
1475 1714455 : GEN N = int2n(e), G = ZXk_gcd_i(poleval(A,N), poleval(B,N));
1476 1714455 : GEN g = Q_primpart(rec(G, e, N, v));
1477 1714455 : if (ZXk_divides(A,g) && ZXk_divides(B,g))
1478 1714441 : return gmul(c,g);
1479 : }
1480 : }
1481 : GEN
1482 1686063 : ZXk_gcd(GEN A, GEN B)
1483 1686063 : { pari_sp av = avma; return gc_upto(av, ZXk_gcd_i(A, B)); }
1484 :
1485 : GEN
1486 189 : QXk_gcd(GEN A, GEN B)
1487 : {
1488 : GEN a, b, D;
1489 189 : pari_sp av = avma, av2;
1490 189 : D = ZXk_gcd_i(Q_primitive_part(A, &a), Q_primitive_part(B, &b));
1491 189 : av2 = avma; a = _gcd(a,b);
1492 189 : if (isint1(a)) set_avma(av2); else D = gmul(D, a);
1493 189 : return gc_upto(av, D);
1494 : }
1495 :
1496 : /*****************************************************************************
1497 : * Variants of the Bradford-Davenport algorithm: look for cyclotomic *
1498 : * factors, and decide whether a ZX is cyclotomic or a product of cyclotomic *
1499 : *****************************************************************************/
1500 : /* f of degree 1, return a cyclotomic factor (Phi_1 or Phi_2) or NULL */
1501 : static GEN
1502 0 : BD_deg1(GEN f)
1503 : {
1504 0 : GEN a = gel(f,3), b = gel(f,2); /* f = ax + b */
1505 0 : if (!absequalii(a,b)) return NULL;
1506 0 : return polcyclo((signe(a) == signe(b))? 2: 1, varn(f));
1507 : }
1508 :
1509 : /* f a squarefree ZX; not divisible by any Phi_n, n even */
1510 : static GEN
1511 406 : BD_odd(GEN f)
1512 : {
1513 406 : while(degpol(f) > 1)
1514 : {
1515 406 : GEN f1 = ZX_graeffe(f); /* contain all cyclotomic divisors of f */
1516 406 : if (ZX_equal(f1, f)) return f; /* product of cyclotomics */
1517 0 : f = ZX_gcd(f, f1);
1518 : }
1519 0 : if (degpol(f) == 1) return BD_deg1(f);
1520 0 : return NULL; /* no cyclotomic divisor */
1521 : }
1522 :
1523 : static GEN
1524 2310 : myconcat(GEN v, GEN x)
1525 : {
1526 2310 : if (typ(x) != t_VEC) x = mkvec(x);
1527 2310 : if (!v) return x;
1528 1470 : return shallowconcat(v, x);
1529 : }
1530 :
1531 : /* Bradford-Davenport algorithm.
1532 : * f a squarefree ZX of degree > 0, return NULL or a vector of coprime
1533 : * cyclotomic factors of f [ possibly reducible ] */
1534 : static GEN
1535 2359 : BD(GEN f)
1536 : {
1537 2359 : GEN G = NULL, Gs = NULL, Gp = NULL, Gi = NULL;
1538 : GEN fs2, fp, f2, f1, fe, fo, fe1, fo1;
1539 2359 : RgX_even_odd(f, &fe, &fo);
1540 2359 : fe1 = ZX_eval1(fe);
1541 2359 : fo1 = ZX_eval1(fo);
1542 2359 : if (absequalii(fe1, fo1)) /* f(1) = 0 or f(-1) = 0 */
1543 : {
1544 1519 : long i, v = varn(f);
1545 1519 : if (!signe(fe1))
1546 371 : G = mkvec2(polcyclo(1, v), polcyclo(2, v)); /* both 0 */
1547 1148 : else if (signe(fe1) == signe(fo1))
1548 693 : G = mkvec(polcyclo(2, v)); /*f(-1) = 0*/
1549 : else
1550 455 : G = mkvec(polcyclo(1, v)); /*f(1) = 0*/
1551 3409 : for (i = lg(G)-1; i; i--) f = RgX_div(f, gel(G,i));
1552 : }
1553 : /* f no longer divisible by Phi_1 or Phi_2 */
1554 2359 : if (degpol(f) <= 1) return G;
1555 2058 : f1 = ZX_graeffe(f); /* has at most square factors */
1556 2058 : if (ZX_equal(f1, f)) return myconcat(G,f); /* f = product of Phi_n, n odd */
1557 :
1558 1183 : fs2 = ZX_gcd_all(f1, ZX_deriv(f1), &f2); /* fs2 squarefree */
1559 1183 : if (degpol(fs2))
1560 : { /* fs contains all Phi_n | f, 4 | n; and only those */
1561 : /* In that case, Graeffe(Phi_n) = Phi_{n/2}^2, and Phi_n = Phi_{n/2}(x^2) */
1562 1029 : GEN fs = RgX_inflate(fs2, 2);
1563 1029 : (void)ZX_gcd_all(f, fs, &f); /* remove those Phi_n | f, 4 | n */
1564 1029 : Gs = BD(fs2);
1565 1029 : if (Gs)
1566 : {
1567 : long i;
1568 2555 : for (i = lg(Gs)-1; i; i--) gel(Gs,i) = RgX_inflate(gel(Gs,i), 2);
1569 : /* prod Gs[i] is the product of all Phi_n | f, 4 | n */
1570 1029 : G = myconcat(G, Gs);
1571 : }
1572 : /* f2 = f1 / fs2 */
1573 1029 : f1 = RgX_div(f2, fs2); /* f1 / fs2^2 */
1574 : }
1575 1183 : fp = ZX_gcd(f, f1); /* contains all Phi_n | f, n > 1 odd; and only those */
1576 1183 : if (degpol(fp))
1577 : {
1578 196 : Gp = BD_odd(fp);
1579 : /* Gp is the product of all Phi_n | f, n odd */
1580 196 : if (Gp) G = myconcat(G, Gp);
1581 196 : f = RgX_div(f, fp);
1582 : }
1583 1183 : if (degpol(f))
1584 : { /* contains all Phi_n originally dividing f, n = 2 mod 4, n > 2;
1585 : * and only those
1586 : * In that case, Graeffe(Phi_n) = Phi_{n/2}, and Phi_n = Phi_{n/2}(-x) */
1587 210 : Gi = BD_odd(ZX_z_unscale(f, -1));
1588 210 : if (Gi)
1589 : { /* N.B. Phi_2 does not divide f */
1590 210 : Gi = ZX_z_unscale(Gi, -1);
1591 : /* Gi is the product of all Phi_n | f, n = 2 mod 4 */
1592 210 : G = myconcat(G, Gi);
1593 : }
1594 : }
1595 1183 : return G;
1596 : }
1597 :
1598 : /* Let f be a nonzero QX, return the (squarefree) product of cyclotomic
1599 : * divisors of f */
1600 : GEN
1601 315 : polcyclofactors(GEN f)
1602 : {
1603 315 : pari_sp av = avma;
1604 315 : if (typ(f) != t_POL || !signe(f)) pari_err_TYPE("polcyclofactors",f);
1605 315 : (void)RgX_valrem(f, &f);
1606 315 : f = Q_primpart(f);
1607 315 : RgX_check_ZX(f,"polcyclofactors");
1608 315 : if (degpol(f))
1609 : {
1610 315 : f = BD(ZX_radical(f));
1611 315 : if (f) return gc_GEN(av, f);
1612 : }
1613 0 : retgc_const(av, cgetg(1, t_VEC));
1614 : }
1615 :
1616 : /* list of all squarefree odd x such that phi(x) = n, P^-(x) > m. Unsorted */
1617 : static GEN
1618 19622 : invphi(ulong n, ulong m)
1619 : {
1620 : GEN C, D;
1621 : long l, i;
1622 19622 : if (n == 1) return mkvecsmall(1);
1623 14331 : D = divisorsu(n); l = lg(D);
1624 14331 : C = cgetg(1, t_VECSMALL);
1625 39791 : for (i = 2; i < l; i++) /* skip 1 */
1626 : {
1627 25460 : ulong d = D[i], p;
1628 25460 : if (d < m) continue;
1629 20290 : p = d + 1; if (!uisprime(p)) continue;
1630 10544 : C = vecsmall_concat(C, zv_z_mul(invphi(D[l-i], p), p));
1631 : }
1632 14331 : return C;
1633 : }
1634 :
1635 : long
1636 99203 : poliscyclo(GEN f)
1637 : {
1638 99203 : const ulong p = 2147483647; /* prime */
1639 : pari_sp av;
1640 : long i, n, e, l;
1641 : ulong f3, fm3;
1642 : GEN D, fp, _3;
1643 99203 : if (typ(f) != t_POL) pari_err_TYPE("poliscyclo", f);
1644 99196 : n = degpol(f);
1645 99196 : if (n <= 0 || !RgX_is_ZX(f)) return 0;
1646 99188 : if (!equali1(gel(f,n+2)) || !is_pm1(gel(f,2))) return 0;
1647 9183 : if (n == 1) return signe(gel(f,2)) > 0? 2: 1;
1648 9078 : av = avma;
1649 9078 : f = ZX_deflate_max(f, &e); if (e != 1) n = degpol(f);
1650 9078 : D = invphi(n, 1); /* squareefree odd d s.t. phi(d) = n */
1651 9078 : l = lg(D); _3 = gmodulss(3, p);
1652 9078 : fp = ZX_to_Flx(f, p);
1653 9078 : f3 = Flx_eval(fp, 3, p);
1654 9078 : fm3 = Flx_eval(fp, p-3, p);
1655 : /* f(x^e) is cyclotomic (= Phi_{de}) iff f = Phi_d, where all prime dividing
1656 : * e also divide d. */
1657 11926 : for (i = 1; i < l; i++)
1658 : {
1659 5158 : long d = D[i]; /* squarefree odd */
1660 5158 : if (odd(e))
1661 : {
1662 4081 : if (e == 1 || u_ppo(e, d) == 1)
1663 : { /* early abort: check whether f(3) = Phi_d(3) or Phi_2d(3) = Phi_d(-3)
1664 : * mod p before checking in Z. N.B. phi(d) and value at 3 mod p
1665 : * determine Phi_d for all d <= 10^7 */
1666 3840 : ulong F3 = Rg_to_Fl(polcyclo_eval(d, _3), p);
1667 3840 : if (F3 == f3 && ZX_equal(f, polcyclo(d, varn(f))))
1668 980 : return gc_long(av, d * e);
1669 2860 : if (F3 == fm3 && ZX_equal(f, polcyclo(2*d, varn(f))))
1670 749 : return gc_long(av, 2* d * e);
1671 : }
1672 : }
1673 : else
1674 : {
1675 1077 : if (u_ppo(e, 2*d) == 1)
1676 : { /* early abort: check whether f(3) = Phi_2d(3) mod p */
1677 1070 : ulong F3 = Rg_to_Fl(polcyclo_eval(2*d, _3), p);
1678 1070 : if (F3 == f3 && ZX_equal(f, polcyclo(2*d, varn(f))))
1679 581 : return gc_long(av, 2* d * e);
1680 : }
1681 : }
1682 : }
1683 6768 : return gc_long(av, 0);
1684 : }
1685 :
1686 : long
1687 1029 : poliscycloprod(GEN f)
1688 : {
1689 1029 : pari_sp av = avma;
1690 1029 : long i, d = degpol(f);
1691 1029 : if (typ(f) != t_POL) pari_err_TYPE("poliscycloprod",f);
1692 1029 : if (!RgX_is_ZX(f)) return 0;
1693 1029 : if (!ZX_is_monic(f) || !is_pm1(constant_coeff(f))) return 0;
1694 1029 : if (d < 2) return (d == 1);
1695 1022 : if ( degpol(ZX_gcd_all(f, ZX_deriv(f), &f)) )
1696 : {
1697 14 : d = degpol(f);
1698 14 : if (d == 1) return 1;
1699 : }
1700 1015 : f = BD(f); if (!f) return 0;
1701 3619 : for (i = lg(f)-1; i; i--) d -= degpol(gel(f,i));
1702 1015 : return gc_long(av, d == 0);
1703 : }
|
