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. 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 : /** ARITHMETIC OPERATIONS ON POLYNOMIALS **/
17 : /** (second part) **/
18 : /** **/
19 : /***********************************************************************/
20 : #include "pari.h"
21 : #include "paripriv.h"
22 :
23 : /* compute Newton sums S_1(P), ... , S_n(P). S_k(P) = sum a_j^k, a_j root of P
24 : * If N != NULL, assume p-adic roots and compute mod N [assume integer coeffs]
25 : * If T != NULL, compute mod (T,N) [assume integer coeffs if N != NULL]
26 : * If y0!= NULL, precomputed i-th powers, i=1..m, m = length(y0).
27 : * Not memory clean in the latter case */
28 : GEN
29 33292 : polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N)
30 : {
31 33292 : long dP=degpol(P), i, k, m;
32 : pari_sp av1, av2;
33 : GEN s,y,P_lead;
34 :
35 33292 : if (n<0) pari_err_IMPL("polsym of a negative n");
36 33292 : if (typ(P) != t_POL) pari_err_TYPE("polsym",P);
37 33292 : if (!signe(P)) pari_err_ROOTS0("polsym");
38 33292 : y = cgetg(n+2,t_COL);
39 33292 : if (y0)
40 : {
41 11375 : if (typ(y0) != t_COL) pari_err_TYPE("polsym_gen",y0);
42 11375 : m = lg(y0)-1;
43 11375 : for (i=1; i<=m; i++) gel(y,i) = gel(y0,i); /* not memory clean */
44 : }
45 : else
46 : {
47 21917 : m = 1;
48 21917 : gel(y,1) = stoi(dP);
49 : }
50 33292 : P += 2; /* strip codewords */
51 :
52 33292 : P_lead = gel(P,dP); if (gequal1(P_lead)) P_lead = NULL;
53 33292 : if (P_lead)
54 : {
55 7 : if (N) P_lead = Fq_inv(P_lead,T,N);
56 7 : else if (T) P_lead = QXQ_inv(P_lead,T);
57 : }
58 86317 : for (k=m; k<=n; k++)
59 : {
60 53025 : av1 = avma; s = (dP>=k)? gmulsg(k,gel(P,dP-k)): gen_0;
61 233772 : for (i=1; i<k && i<=dP; i++)
62 180747 : s = gadd(s, gmul(gel(y,k-i+1),gel(P,dP-i)));
63 53025 : if (N)
64 : {
65 15960 : s = Fq_red(s, T, N);
66 15960 : if (P_lead) s = Fq_mul(s, P_lead, T, N);
67 : }
68 37065 : else if (T)
69 : {
70 0 : s = grem(s, T);
71 0 : if (P_lead) s = grem(gmul(s, P_lead), T);
72 : }
73 : else
74 37065 : if (P_lead) s = gdiv(s, P_lead);
75 53025 : av2 = avma; gel(y,k+1) = gerepile(av1,av2, gneg(s));
76 : }
77 33292 : return y;
78 : }
79 :
80 : GEN
81 18676 : polsym(GEN x, long n)
82 : {
83 18676 : return polsym_gen(x, NULL, n, NULL,NULL);
84 : }
85 :
86 : /* centered residue x mod p. po2 = shifti(p, -1) or NULL (euclidean residue) */
87 : GEN
88 55653776 : centermodii(GEN x, GEN p, GEN po2)
89 : {
90 55653776 : GEN y = remii(x, p);
91 55668107 : switch(signe(y))
92 : {
93 2060103 : case 0: break;
94 33442583 : case 1: if (po2 && abscmpii(y,po2) > 0) y = subii(y, p);
95 33401680 : break;
96 20276234 : case -1: if (!po2 || abscmpii(y,po2) > 0) y = addii(y, p);
97 20255839 : break;
98 : }
99 55606809 : return y;
100 : }
101 :
102 : static long
103 0 : s_centermod(long x, ulong pp, ulong pps2)
104 : {
105 0 : long y = x % (long)pp;
106 0 : if (y < 0) y += pp;
107 0 : return Fl_center(y, pp,pps2);
108 : }
109 :
110 : /* for internal use */
111 : GEN
112 8375085 : centermod_i(GEN x, GEN p, GEN ps2)
113 : {
114 : long i, lx;
115 : pari_sp av;
116 : GEN y;
117 :
118 8375085 : if (!ps2) ps2 = shifti(p,-1);
119 8375752 : switch(typ(x))
120 : {
121 3672886 : case t_INT: return centermodii(x,p,ps2);
122 :
123 3904411 : case t_POL: lx = lg(x);
124 3904411 : y = cgetg(lx,t_POL); y[1] = x[1];
125 32099766 : for (i=2; i<lx; i++)
126 : {
127 28194744 : av = avma;
128 28194744 : gel(y,i) = gerepileuptoint(av, centermodii(gel(x,i),p,ps2));
129 : }
130 3905022 : return normalizepol_lg(y, lx);
131 :
132 797496 : case t_COL: lx = lg(x);
133 797496 : y = cgetg(lx,t_COL);
134 797497 : for (i=1; i<lx; i++) gel(y,i) = centermodii(gel(x,i),p,ps2);
135 797495 : return y;
136 :
137 959 : case t_MAT: lx = lg(x);
138 959 : y = cgetg(lx,t_MAT);
139 959 : for (i=1; i<lx; i++) gel(y,i) = centermod_i(gel(x,i),p,ps2);
140 959 : return y;
141 :
142 0 : case t_VECSMALL: lx = lg(x);
143 : {
144 0 : ulong pp = itou(p), pps2 = itou(ps2);
145 0 : y = cgetg(lx,t_VECSMALL);
146 0 : for (i=1; i<lx; i++) y[i] = s_centermod(x[i], pp, pps2);
147 0 : return y;
148 : }
149 : }
150 0 : return x;
151 : }
152 :
153 : GEN
154 4805745 : centermod(GEN x, GEN p) { return centermod_i(x,p,NULL); }
155 :
156 : static GEN
157 28 : RgX_Frobenius_deflate(GEN S, ulong p)
158 : {
159 28 : GEN F = RgX_deflate(S, p);
160 28 : long i, l = lg(F);
161 126 : for (i=2; i<l; i++)
162 : {
163 42 : GEN Fi = gel(F,i), R;
164 42 : if (typ(Fi)==t_POL)
165 : {
166 21 : if (signe(RgX_deriv(Fi))==0)
167 14 : gel(F,i) = RgX_Frobenius_deflate(gel(F, i), p);
168 14 : else return NULL;
169 : }
170 21 : else if (ispower(Fi, utoi(p), &R))
171 21 : gel(F,i) = R;
172 0 : else return NULL;
173 : }
174 21 : return F;
175 : }
176 :
177 :
178 : static GEN
179 84 : RgXY_squff(GEN f)
180 : {
181 84 : long i, q, n = degpol(f);
182 84 : ulong p = itos_or_0(characteristic(f));
183 84 : GEN u = const_vec(n+1, pol_1(varn(f)));
184 91 : for(q = 1;;q *= p)
185 : {
186 91 : GEN t, v, tv, r = RgX_gcd(f, RgX_deriv(f));
187 91 : if (degpol(r) == 0) { gel(u, q) = f; break; }
188 21 : t = RgX_div(f, r);
189 21 : if (degpol(t) > 0)
190 : {
191 : long j;
192 14 : for(j = 1;;j++)
193 : {
194 14 : v = RgX_gcd(r, t);
195 14 : tv = RgX_div(t, v);
196 14 : if (degpol(tv) > 0) gel(u, j*q) = tv;
197 14 : if (degpol(v) <= 0) break;
198 7 : r = RgX_div(r, v);
199 7 : t = v;
200 7 : }
201 7 : if (degpol(r) == 0) break;
202 : }
203 14 : if (!p) break;
204 14 : r = RgX_Frobenius_deflate(f, p);
205 14 : if (!r) { gel(u, q) = f; break; }
206 7 : f = r;
207 7 : }
208 413 : for (i = n; i; i--)
209 413 : if (degpol(gel(u,i))) break;
210 84 : setlg(u,i+1); return u;
211 : }
212 :
213 : /* Lmod contains modular factors of *F (NULL codes an empty slot: used factor)
214 : * Lfac accumulates irreducible factors as they are found.
215 : * p is a product of modular factors in Lmod[1..i-1] (NULL for p = 1), not
216 : * a rational factor of *F
217 : * Find an irreducible factor of *F divisible by p (by including
218 : * exhaustively further factors from Lmod[i..]); return 0 on failure, else 1.
219 : * Update Lmod, Lfac and *F */
220 : static int
221 4487 : RgX_cmbf(GEN p, long i, GEN BLOC, GEN Lmod, GEN Lfac, GEN *F)
222 : {
223 : GEN q;
224 4487 : if (i == lg(Lmod)) return 0;
225 2352 : if (RgX_cmbf(p, i+1, BLOC, Lmod, Lfac, F) && p) return 1;
226 2226 : if (!gel(Lmod,i)) return 0;
227 2219 : p = p? RgX_mul(p, gel(Lmod,i)): gel(Lmod,i);
228 2219 : q = RgV_to_RgX(RgX_digits(p, BLOC), varn(*F));
229 2219 : if (degpol(q))
230 : {
231 2023 : GEN R, Q = RgX_divrem(*F, q, &R);
232 2023 : if (signe(R)==0) { vectrunc_append(Lfac, q); *F = Q; return 1; }
233 : }
234 2044 : if (RgX_cmbf(p, i+1, BLOC, Lmod, Lfac, F)) { gel(Lmod,i) = NULL; return 1; }
235 1904 : return 0;
236 : }
237 :
238 : static GEN
239 91 : RgXY_factor_squarefree(GEN f)
240 : {
241 91 : pari_sp av = avma;
242 91 : GEN Lfac, Lmod, F = NULL, BLOC = NULL;
243 91 : long vy = gvar2(f), n = RgXY_degreex(f);
244 91 : ulong i, c = itou_or_0(residual_characteristic(f));
245 91 : long val = RgX_valrem(f, &f);
246 : for(;;)
247 : {
248 154 : for (i = 0; !c || i < c; i++)
249 : {
250 154 : BLOC = gpowgs(gaddgs(pol_x(vy), i), n+1);
251 154 : F = poleval(f, BLOC);
252 154 : if (issquarefree(c ? gmul(F,gmodulss(1,c)): F)) break;
253 : }
254 91 : if (!c || i < c) break;
255 0 : n++;
256 0 : }
257 91 : if (DEBUGLEVEL >= 2)
258 0 : err_printf("bifactor: bloc:(x+%ld)^%ld, deg f=%ld\n",i,n,RgXY_degreex(f));
259 91 : Lmod = gel(factor(F),1);
260 91 : if (DEBUGLEVEL >= 2)
261 0 : err_printf("bifactor: %ld local factors\n",lg(Lmod)-1);
262 91 : Lfac = vectrunc_init(lg(Lmod)+1);
263 91 : settyp(Lfac, t_COL);
264 91 : (void)RgX_cmbf(NULL, 1, BLOC, Lmod, Lfac, &f);
265 91 : if (degpol(f)) vectrunc_append(Lfac, f);
266 91 : if (val) vectrunc_append(Lfac, pol_x(varn(f)));
267 91 : return gerepilecopy(av, Lfac);
268 : }
269 :
270 : static GEN
271 84 : FE_matconcat(GEN F, GEN E, long l)
272 : {
273 84 : setlg(E,l); E = shallowconcat1(E);
274 84 : setlg(F,l); F = shallowconcat1(F); return mkmat2(F,E);
275 : }
276 :
277 : static GEN
278 84 : RgXY_factor(GEN f)
279 : {
280 84 : pari_sp av = avma;
281 : GEN F, E;
282 84 : GEN cnt = content(f);
283 84 : GEN V = RgXY_squff(gdiv(f, cnt));
284 84 : long i, j, l = lg(V);
285 84 : GEN C = factor(cnt);
286 84 : F = cgetg(l+1, t_VEC);
287 84 : E = cgetg(l+1, t_VEC);
288 84 : gel(F,1) = gel(C,1);
289 84 : gel(E,1) = gel(C,2);
290 217 : for (i=1, j=2; i < l; i++)
291 133 : if (degpol(gel(V,i)))
292 : {
293 91 : GEN Fj = RgXY_factor_squarefree(gel(V,i));
294 91 : gel(F, j) = Fj;
295 91 : gel(E, j) = const_vec(lg(Fj)-1, stoi(i));
296 91 : j++;
297 : }
298 84 : return gerepilecopy(av, sort_factor_pol(FE_matconcat(F,E,j), cmp_universal));
299 : }
300 :
301 : /***********************************************************************/
302 : /** **/
303 : /** FACTORIZATION **/
304 : /** **/
305 : /***********************************************************************/
306 : #define assign_or_fail(x,y) { GEN __x = x;\
307 : if (!*y) *y=__x; else if (!gequal(__x,*y)) return 0;\
308 : }
309 : #define update_prec(x,y) { long __x = x; if (__x < *y) *y=__x; }
310 :
311 : static const long tsh = 6;
312 : static long
313 222004 : code(long t1, long t2) { return (t1 << tsh) | t2; }
314 : void
315 401982 : RgX_type_decode(long x, long *t1, long *t2)
316 : {
317 401982 : *t1 = x >> tsh;
318 401982 : *t2 = (x & ((1L<<tsh)-1));
319 401982 : }
320 : int
321 2455071 : RgX_type_is_composite(long t) { return t >= tsh; }
322 :
323 : static int
324 544216103 : settype(GEN c, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
325 : {
326 : long j;
327 544216103 : switch(typ(c))
328 : {
329 : case t_INT:
330 456477584 : break;
331 : case t_FRAC:
332 13350597 : t[1]=1; break;
333 : break;
334 : case t_REAL:
335 26373772 : update_prec(precision(c), pa);
336 26373772 : t[2]=1; break;
337 : case t_INTMOD:
338 20197277 : assign_or_fail(gel(c,1),p);
339 20197277 : t[3]=1; break;
340 : case t_FFELT:
341 1358762 : if (!*ff) *ff=c; else if (!FF_samefield(c,*ff)) return 0;
342 1358762 : assign_or_fail(FF_p_i(c),p);
343 1358762 : t[5]=1; break;
344 : case t_COMPLEX:
345 23673673 : for (j=1; j<=2; j++)
346 : {
347 15782451 : GEN d = gel(c,j);
348 15782451 : switch(typ(d))
349 : {
350 : case t_INT: case t_FRAC:
351 644267 : t[4]=1; break;
352 : case t_REAL:
353 15138163 : update_prec(precision(d), pa);
354 15138163 : t[6]=1; break;
355 : case t_INTMOD:
356 14 : assign_or_fail(gel(d,1),p);
357 14 : if (!signe(*p) || mod4(*p) != 3) return 0;
358 7 : if (!*t2) *t2 = t_COMPLEX;
359 7 : t[3]=1; break;
360 : case t_PADIC:
361 7 : update_prec(precp(d)+valp(d), pa);
362 7 : assign_or_fail(gel(d,2),p);
363 7 : if (!*t2) *t2 = t_COMPLEX;
364 7 : t[7]=1; break;
365 0 : default: return 0;
366 : }
367 : }
368 7891222 : if (!t[6]) assign_or_fail(mkpoln(3, gen_1,gen_0,gen_1), pol); /*x^2+1*/
369 7891222 : break;
370 : case t_PADIC:
371 208411 : update_prec(precp(c)+valp(c), pa);
372 208411 : assign_or_fail(gel(c,2),p);
373 208411 : t[7]=1; break;
374 : case t_QUAD:
375 210 : assign_or_fail(gel(c,1),pol);
376 630 : for (j=2; j<=3; j++)
377 : {
378 420 : GEN d = gel(c,j);
379 420 : switch(typ(d))
380 : {
381 : case t_INT: case t_FRAC:
382 385 : t[8]=1; break;
383 : case t_INTMOD:
384 28 : assign_or_fail(gel(d,1),p);
385 28 : if (*t2 != t_POLMOD) *t2 = t_QUAD;
386 28 : t[3]=1; break;
387 : case t_PADIC:
388 7 : update_prec(precp(d)+valp(d), pa);
389 7 : assign_or_fail(gel(d,2),p);
390 7 : if (*t2 != t_POLMOD) *t2 = t_QUAD;
391 7 : t[7]=1; break;
392 0 : default: return 0;
393 : }
394 : }
395 210 : break;
396 : case t_POLMOD:
397 3132752 : assign_or_fail(gel(c,1),pol);
398 3132752 : if (typ(gel(c,2))==t_POL && varn(gel(c,2))!=varn(gel(c,1)))
399 182 : return 0;
400 18790646 : for (j=1; j<=2; j++)
401 : {
402 : GEN pbis, polbis;
403 : long pabis;
404 6264545 : *t2 = t_POLMOD;
405 6264545 : switch(Rg_type(gel(c,j),&pbis,&polbis,&pabis))
406 : {
407 3284912 : case t_INT: break;
408 486891 : case t_FRAC: t[1]=1; break;
409 2490943 : case t_INTMOD: t[3]=1; break;
410 7 : case t_PADIC: t[7]=1; update_prec(pabis,pa); break;
411 3584 : default: return 0;
412 : }
413 6262753 : if (pbis) assign_or_fail(pbis,p);
414 6262753 : if (polbis) assign_or_fail(polbis,pol);
415 : }
416 3130778 : break;
417 14007975 : case t_POL: t[10] = 1;
418 : {
419 : GEN pbis, polbis;
420 : long tbis, pabis;
421 14007975 : tbis = RgX_type(c,&pbis,&polbis,&pabis);
422 14007975 : if (tbis && tbis!=t_POL && !polbis)
423 : {
424 13527753 : if (*var == NO_VARIABLE) { *var = varn(c); break; }
425 8324622 : if (*var == varn(c)) break;
426 : }
427 670946 : *var = MAXVARN+1; break; /* ensure varn() == *var fails in choosetype */
428 : }
429 1217534 : default: return 0;
430 : }
431 542996588 : return 1;
432 : }
433 : /* t[0] unused. Other values, if set, indicate a coefficient of type
434 : * t[1] : t_FRAC
435 : * t[2] : t_REAL
436 : * t[3] : t_INTMOD
437 : * t[4] : t_COMPLEX of rationals (t_INT/t_FRAC)
438 : * t[5] : t_FFELT
439 : * t[6] : t_COMPLEX of t_REAL
440 : * t[7] : t_PADIC
441 : * t[8] : t_QUAD of rationals (t_INT/t_FRAC)
442 : * t[9]: Unused
443 : * t[10]: t_POL (recursive factorisation) */
444 : /* if t2 != 0: t_POLMOD/t_QUAD/t_COMPLEX of modular (t_INTMOD/t_PADIC,
445 : * given by t) */
446 : static long
447 75377088 : choosetype(long *t, long t2, GEN ff, GEN *pol, long var)
448 : {
449 75377088 : if (t[10] && (!*pol || var!=varn(*pol))) return t_POL;
450 70046761 : if (t[5]) /* ffelt */
451 : {
452 222832 : if (t2 ||t[2]||t[4]||t[6]||t[8]||t[9]) return 0;
453 222832 : *pol=ff; return t_FFELT;
454 : }
455 69823929 : if (t[6]) /* inexact, complex */
456 : {
457 1145525 : if (t2 ||t[3]||t[7]||t[9]) return 0;
458 1145525 : return t_COMPLEX;
459 : }
460 68678404 : if (t[2]) /* inexact, real */
461 : {
462 4040397 : if (t2 ||t[3]||t[7]||t[9]) return 0;
463 4040397 : return t[4]?t_COMPLEX:t_REAL;
464 : }
465 64638007 : if (t2) /* polmod/quad/complex of intmod/padic */
466 : {
467 216558 : if (t[3]) return code(t2,t_INTMOD);
468 189888 : if (t[7]) return code(t2,t_PADIC);
469 189867 : if (t[1]) return code(t2,t_FRAC);
470 123393 : return code(t2,t_INT);
471 : }
472 64421449 : if (t[10]) return t_POL;
473 64421449 : if (t[8]) return code(t_QUAD,t_INT);
474 64421337 : if (t[4]) return code(t_COMPLEX,t_INT);
475 64416003 : if (t[3]) return t_INTMOD;
476 60821979 : if (t[7]) return t_PADIC;
477 60803506 : if (t[1]) return t_FRAC;
478 57805988 : return t_INT;
479 : }
480 :
481 : static long
482 121412476 : RgX_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
483 : {
484 121412476 : long i, lx = lg(x);
485 568716321 : for (i=2; i<lx; i++)
486 448526035 : if (!settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
487 120190286 : return 1;
488 : }
489 :
490 : static long
491 15956883 : RgC_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
492 : {
493 15956883 : long i, l = lg(x);
494 111420426 : for (i = 1; i<l; i++)
495 95465216 : if (!settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
496 15955210 : return 1;
497 : }
498 :
499 : static long
500 5355522 : RgM_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
501 : {
502 5355522 : long i, l = lg(x);
503 21048127 : for (i = 1; i < l; i++)
504 15694243 : if (!RgC_settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
505 5353884 : return 1;
506 : }
507 :
508 : long
509 8719616 : Rg_type(GEN x, GEN *p, GEN *pol, long *pa)
510 : {
511 8719616 : long t[] = {0,0,0,0,0,0,0,0,0,0,0};
512 8719616 : long tx = typ(x), t2 = 0, var = NO_VARIABLE;
513 8719616 : GEN ff = NULL;
514 8719616 : *p = *pol = NULL; *pa = LONG_MAX;
515 8719616 : if (is_scalar_t(tx))
516 : {
517 226578 : if (tx == t_POLMOD) return 0;
518 226193 : if (!settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
519 : }
520 8493038 : else if (tx == t_MAT)
521 7 : { if(!RgM_settype(x, t, p, pol, pa, &ff, &t2, &var)) return 0; }
522 : else /* t_POL, t_SER */
523 8493031 : if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
524 8719231 : return choosetype(t,t2,ff,pol,var);
525 : }
526 :
527 : long
528 16667361 : RgX_type(GEN x, GEN *p, GEN *pol, long *pa)
529 : {
530 16667361 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
531 16667361 : long t2 = 0, var = NO_VARIABLE;
532 16667361 : GEN ff = NULL;
533 16667361 : *p = *pol = NULL; *pa = LONG_MAX;
534 16667361 : if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
535 16648535 : return choosetype(t,t2,ff,pol,var);
536 : }
537 :
538 : long
539 47640502 : RgX_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
540 : {
541 47640502 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
542 47640502 : long t2 = 0, var = NO_VARIABLE;
543 47640502 : GEN ff = NULL;
544 47640502 : *p = *pol = NULL; *pa = LONG_MAX;
545 94093533 : if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
546 47651509 : !RgX_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
547 46442703 : return choosetype(t,t2,ff,pol,var);
548 : }
549 :
550 : long
551 718390 : RgX_type3(GEN x, GEN y, GEN z, GEN *p, GEN *pol, long *pa)
552 : {
553 718390 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
554 718390 : long t2 = 0, var = NO_VARIABLE;
555 718390 : GEN ff = NULL;
556 718390 : *p = *pol = NULL; *pa = LONG_MAX;
557 1436784 : if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
558 1436787 : !RgX_settype(y,t,p,pol,pa,&ff,&t2,&var) ||
559 718394 : !RgX_settype(z,t,p,pol,pa,&ff,&t2,&var)) return 0;
560 718397 : return choosetype(t,t2,ff,pol,var);
561 : }
562 :
563 : long
564 83538 : RgM_type(GEN x, GEN *p, GEN *pol, long *pa)
565 : {
566 83538 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
567 83538 : long t2 = 0, var = NO_VARIABLE;
568 83538 : GEN ff = NULL;
569 83538 : *p = *pol = NULL; *pa = LONG_MAX;
570 83538 : if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
571 82838 : return choosetype(t,t2,ff,pol,var);
572 : }
573 :
574 : long
575 263032 : RgM_RgC_type(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
576 : {
577 263032 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
578 263032 : long t2 = 0, var = NO_VARIABLE;
579 263032 : GEN ff = NULL;
580 263032 : *p = *pol = NULL; *pa = LONG_MAX;
581 525672 : if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
582 263067 : !RgC_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
583 262605 : return choosetype(t,t2,ff,pol,var);
584 : }
585 :
586 : long
587 2504693 : RgM_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
588 : {
589 2504693 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
590 2504693 : long t2 = 0, var = NO_VARIABLE;
591 2504693 : GEN ff = NULL;
592 2504693 : *p = *pol = NULL; *pa = LONG_MAX;
593 5008945 : if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
594 2504798 : !RgM_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
595 2504147 : return choosetype(t,t2,ff,pol,var);
596 : }
597 :
598 : GEN
599 0 : factor0(GEN x,long flag) { return (flag<0)? factor(x): boundfact(x,flag); }
600 :
601 : /* only present for interface with GP */
602 : GEN
603 37164 : gp_factor0(GEN x, GEN flag)
604 : {
605 : ulong B;
606 37164 : if (!flag) return factor(x);
607 35 : if (typ(flag) != t_INT || signe(flag) < 0) pari_err_FLAG("factor");
608 35 : switch(lgefint(flag))
609 : {
610 7 : case 2: B = 0; break;
611 28 : case 3: B = flag[2]; break;
612 0 : default: pari_err_OVERFLOW("factor [large prime bound]");
613 : return NULL; /*LCOV_EXCL_LINE*/
614 : }
615 35 : return boundfact(x, B);
616 : }
617 :
618 : GEN
619 83403 : deg1_from_roots(GEN L, long v)
620 : {
621 83403 : long i, l = lg(L);
622 83403 : GEN z = cgetg(l,t_COL);
623 190457 : for (i=1; i<l; i++)
624 107054 : gel(z,i) = deg1pol_shallow(gen_1, gneg(gel(L,i)), v);
625 83403 : return z;
626 : }
627 : GEN
628 980 : roots_from_deg1(GEN x)
629 : {
630 980 : long i,l = lg(x);
631 980 : GEN r = cgetg(l,t_VEC);
632 980 : for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = gneg(gel(P,2)); }
633 980 : return r;
634 : }
635 :
636 : static GEN
637 35 : gauss_factor_p(GEN p)
638 : {
639 35 : GEN a, b; (void)cornacchia(gen_1, p, &a,&b);
640 35 : return mkcomplex(a, b);
641 : }
642 :
643 : static GEN
644 42 : gauss_primpart(GEN x, GEN *c)
645 : {
646 42 : GEN a = gel(x,1), b = gel(x,2), n = gcdii(a, b);
647 42 : *c = n; if (n == gen_1) return x;
648 42 : retmkcomplex(diviiexact(a,n), diviiexact(b,n));
649 : }
650 :
651 : static GEN
652 70 : gauss_primpart_try(GEN x, GEN c)
653 : {
654 : GEN r, y;
655 70 : if (typ(x) == t_INT)
656 : {
657 42 : y = dvmdii(x, c, &r); if (r != gen_0) return NULL;
658 : }
659 : else
660 : {
661 28 : GEN a = gel(x,1), b = gel(x,2); y = cgetg(3, t_COMPLEX);
662 28 : gel(y,1) = dvmdii(a, c, &r); if (r != gen_0) return NULL;
663 14 : gel(y,2) = dvmdii(b, c, &r); if (r != gen_0) return NULL;
664 : }
665 56 : return y;
666 : }
667 :
668 : static int
669 77 : gauss_cmp(GEN x, GEN y)
670 : {
671 : int v;
672 77 : if (typ(x) != t_COMPLEX)
673 0 : return (typ(y) == t_COMPLEX)? -1: gcmp(x, y);
674 77 : if (typ(y) != t_COMPLEX) return 1;
675 49 : v = cmpii(gel(x,2), gel(y,2));
676 49 : if (v) return v;
677 21 : return gcmp(gel(x,1), gel(y,1));
678 : }
679 :
680 : /* 0 or canonical representative in Z[i]^* / <i> (impose imag(x) >= 0) */
681 : static GEN
682 448 : gauss_normal(GEN x)
683 : {
684 448 : if (typ(x) != t_COMPLEX) return (signe(x) < 0)? absi(x): x;
685 378 : if (signe(gel(x,1)) < 0) x = gneg(x);
686 378 : if (signe(gel(x,2)) < 0) x = mulcxI(x);
687 378 : return x;
688 : }
689 :
690 : static GEN
691 42 : gauss_factor(GEN x)
692 : {
693 42 : pari_sp av = avma;
694 42 : GEN a = gel(x,1), b = gel(x,2), d = gen_1, n, y, fa, P, E, P2, E2;
695 42 : long t1 = typ(a);
696 42 : long t2 = typ(b), i, j, l, exp = 0;
697 42 : if (t1 == t_FRAC) d = gel(a,2);
698 42 : if (t2 == t_FRAC) d = lcmii(d, gel(b,2));
699 42 : if (d == gen_1) y = x;
700 : else
701 : {
702 14 : y = gmul(x, d);
703 14 : a = gel(y,1); t1 = typ(a);
704 14 : b = gel(y,2); t2 = typ(b);
705 : }
706 42 : if (t1 != t_INT || t2 != t_INT) return NULL;
707 42 : y = gauss_primpart(y, &n);
708 42 : fa = factor(cxnorm(y));
709 42 : P = gel(fa,1);
710 42 : E = gel(fa,2); l = lg(P);
711 42 : P2 = cgetg(l, t_COL);
712 42 : E2 = cgetg(l, t_COL);
713 98 : for (j = 1, i = l-1; i > 0; i--) /* remove largest factors first */
714 : { /* either p = 2 (ramified) or those factors split in Q(i) */
715 56 : GEN p = gel(P,i), w, w2, t, we, pe;
716 56 : long v, e = itos(gel(E,i));
717 56 : int is2 = absequaliu(p, 2);
718 56 : w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
719 56 : w2 = gauss_normal( gconj(w) );
720 : /* w * w2 * I^3 = p, w2 = gconj(w) * I */
721 56 : pe = powiu(p, e);
722 56 : we = gpowgs(w, e);
723 56 : t = gauss_primpart_try( gmul(y, gconj(we)), pe );
724 56 : if (t) y = t; /* y /= w^e */
725 : else {
726 : /* y /= conj(w)^e, should be y /= w2^e */
727 14 : y = gauss_primpart_try( gmul(y, we), pe );
728 14 : swap(w, w2); exp -= e; /* += 3*e mod 4 */
729 : }
730 56 : gel(P,i) = w;
731 56 : v = Z_pvalrem(n, p, &n);
732 56 : if (v) {
733 7 : exp -= v; /* += 3*v mod 4 */
734 7 : if (is2) v <<= 1; /* 2 = w^2 I^3 */
735 : else {
736 0 : gel(P2,j) = w2;
737 0 : gel(E2,j) = utoipos(v); j++;
738 : }
739 7 : gel(E,i) = stoi(e + v);
740 : }
741 56 : v = Z_pvalrem(d, p, &d);
742 56 : if (v) {
743 7 : exp += v; /* -= 3*v mod 4 */
744 7 : if (is2) v <<= 1; /* 2 is ramified */
745 : else {
746 7 : gel(P2,j) = w2;
747 7 : gel(E2,j) = utoineg(v); j++;
748 : }
749 7 : gel(E,i) = stoi(e - v);
750 : }
751 56 : exp &= 3;
752 : }
753 42 : if (j > 1) {
754 7 : long k = 1;
755 7 : GEN P1 = cgetg(l, t_COL);
756 7 : GEN E1 = cgetg(l, t_COL);
757 : /* remove factors with exponent 0 */
758 14 : for (i = 1; i < l; i++)
759 7 : if (signe(gel(E,i)))
760 : {
761 0 : gel(P1,k) = gel(P,i);
762 0 : gel(E1,k) = gel(E,i);
763 0 : k++;
764 : }
765 7 : setlg(P1, k); setlg(E1, k);
766 7 : setlg(P2, j); setlg(E2, j);
767 7 : fa = famat_mul_shallow(mkmat2(P1,E1), mkmat2(P2,E2));
768 : }
769 42 : if (!is_pm1(n) || !is_pm1(d))
770 : {
771 21 : GEN Fa = factor(gdiv(n, d));
772 21 : P = gel(Fa,1); l = lg(P);
773 21 : E = gel(Fa,2);
774 49 : for (i = 1; i < l; i++)
775 : {
776 28 : GEN w, p = gel(P,i);
777 : long e;
778 : int is2;
779 28 : switch(mod4(p))
780 : {
781 14 : case 3: continue;
782 7 : case 2: is2 = 1; break;
783 7 : default:is2 = 0; break;
784 : }
785 14 : e = itos(gel(E,i));
786 14 : w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
787 14 : gel(P,i) = w;
788 14 : if (is2)
789 7 : gel(E,i) = stoi(2*e);
790 : else
791 : {
792 7 : P = shallowconcat(P, gauss_normal( gconj(w) ));
793 7 : E = shallowconcat(E, gel(E,i));
794 : }
795 14 : exp -= e; /* += 3*e mod 4 */
796 14 : exp &= 3;
797 : }
798 21 : gel(Fa,1) = P;
799 21 : gel(Fa,2) = E;
800 21 : fa = famat_mul_shallow(fa, Fa);
801 : }
802 42 : fa = sort_factor(fa, (void*)&gauss_cmp, &cmp_nodata);
803 :
804 42 : y = gmul(y, powIs(exp));
805 42 : if (!gequal1(y)) {
806 35 : gel(fa,1) = shallowconcat(mkcol(y), gel(fa,1));
807 35 : gel(fa,2) = shallowconcat(gen_1, gel(fa,2));
808 : }
809 42 : return gerepilecopy(av, fa);
810 : }
811 :
812 : GEN
813 0 : Q_factor_limit(GEN x, ulong lim)
814 : {
815 0 : pari_sp av = avma;
816 : GEN a, b;
817 0 : if (typ(x) == t_INT) return Z_factor_limit(x, lim);
818 0 : a = Z_factor_limit(gel(x,1), lim);
819 0 : b = Z_factor_limit(gel(x,2), lim); gel(b,2) = ZC_neg(gel(b,2));
820 0 : return gerepilecopy(av, merge_factor(a,b,(void*)&cmpii,cmp_nodata));
821 : }
822 : GEN
823 2527 : Q_factor(GEN x)
824 : {
825 2527 : pari_sp av = avma;
826 : GEN a, b;
827 2527 : if (typ(x) == t_INT) return Z_factor(x);
828 182 : a = Z_factor(gel(x,1));
829 182 : b = Z_factor(gel(x,2)); gel(b,2) = ZC_neg(gel(b,2));
830 182 : return gerepilecopy(av, merge_factor(a,b,(void*)&cmpii,cmp_nodata));
831 : }
832 : /* replace t_QUAD/t_COMPLEX coeffs by t_POLMOD in T */
833 : static GEN
834 56 : RgX_fix_quad(GEN x, GEN T)
835 : {
836 56 : long i, l, v = varn(T);
837 56 : GEN y = cgetg_copy(x,&l);
838 294 : for (i = 2; i < l; i++)
839 : {
840 238 : GEN c = gel(x,i);
841 238 : switch(typ(c))
842 : {
843 28 : case t_QUAD: c++;/* fall through */
844 70 : case t_COMPLEX: c = deg1pol_shallow(gel(c,2),gel(c,1),v);
845 : }
846 238 : gel(y,i) = c;
847 : }
848 56 : y[1] = x[1]; return y;
849 : }
850 :
851 : GEN
852 40384 : factor(GEN x)
853 : {
854 40384 : long tx=typ(x), lx, i, pa, v, r1;
855 : pari_sp av, tetpil;
856 : GEN y, p, p1, p2, pol;
857 :
858 40384 : if (gequal0(x))
859 63 : switch(tx)
860 : {
861 : case t_INT:
862 : case t_COMPLEX:
863 : case t_POL:
864 63 : case t_RFRAC: return prime_fact(x);
865 0 : default: pari_err_TYPE("factor",x);
866 : }
867 40321 : av = avma;
868 40321 : switch(tx)
869 : {
870 38228 : case t_INT: return Z_factor(x);
871 168 : case t_FRAC: return Q_factor(x);
872 :
873 : case t_POL:
874 1876 : tx=RgX_type(x,&p,&pol,&pa);
875 1876 : switch(tx)
876 : {
877 7 : case 0: pari_err_IMPL("factor for general polynomials");
878 84 : case t_POL: return RgXY_factor(x);
879 1400 : case t_INT: return ZX_factor(x);
880 7 : case t_FRAC: return QX_factor(x);
881 70 : case t_INTMOD: return factmod(x,p);
882 28 : case t_PADIC: return factorpadic(x,p,pa);
883 98 : case t_FFELT: return FFX_factor(x,pol);
884 :
885 14 : case t_COMPLEX: y = cgetg(3,t_MAT);
886 14 : av = avma; p1 = deg1_from_roots(roots(x,pa), varn(x));
887 14 : gel(y,1) = p1 = gerepileupto(av, p1);
888 14 : gel(y,2) = const_col(lg(p1)-1, gen_1); return y;
889 :
890 21 : case t_REAL: y=cgetg(3,t_MAT); v=varn(x);
891 21 : av=avma; p1=cleanroots(x,pa); tetpil=avma;
892 21 : lx = lg(p1);
893 49 : for (r1 = 1; r1 < lx; r1++)
894 35 : if (typ(gel(p1,r1)) == t_COMPLEX) break;
895 21 : lx=(r1+lx)>>1; p2=cgetg(lx,t_COL);
896 49 : for (i = 1; i < r1; i++)
897 28 : gel(p2,i) = deg1pol_shallow(gen_1, negr(gel(p1,i)), v);
898 28 : for ( ; i < lx; i++)
899 : {
900 7 : GEN a = gel(p1,2*i-r1);
901 7 : p = cgetg(5, t_POL); gel(p2,i) = p;
902 7 : p[1] = x[1];
903 7 : gel(p,2) = gnorm(a);
904 7 : gel(p,3) = gmul2n(gel(a,1),1); togglesign(gel(p,3));
905 7 : gel(p,4) = gen_1;
906 : }
907 21 : gel(y,1) = gerepile(av,tetpil,p2);
908 21 : gel(y,2) = const_col(lx-1, gen_1); return y;
909 :
910 : default:
911 : {
912 147 : GEN w = NULL, T = pol;
913 : long t1, t2;
914 147 : RgX_type_decode(tx, &t1, &t2);
915 147 : if (t1 == t_COMPLEX) w = gen_I();
916 119 : else if (t1 == t_QUAD) w = mkquad(pol,gen_0,gen_1);
917 147 : if (w)
918 : { /* substitute I or w by t_POLMOD */
919 56 : T = leafcopy(pol); setvarn(T, fetch_var());
920 56 : x = RgX_fix_quad(x, T);
921 : }
922 147 : switch (t2)
923 : {
924 91 : case t_INT: case t_FRAC: p1 = nffactor(T,x); break;
925 : case t_INTMOD:
926 35 : T = RgX_to_FpX(T,p);
927 35 : if (FpX_is_irred(T,p)) { p1 = factmod(x,mkvec2(p,T)); break; }
928 : /*fall through*/
929 : default:
930 35 : if (w) (void)delete_var();
931 35 : pari_err_IMPL("factor for general polynomial");
932 : return NULL; /* LCOV_EXCL_LINE */
933 : }
934 112 : if (t1 == t_POLMOD) return gerepileupto(av, p1);
935 : /* substitute back I or w */
936 35 : gel(p1,1) = gsubst(liftpol_shallow(gel(p1,1)), varn(T), w);
937 35 : (void)delete_var(); return gerepilecopy(av, p1);
938 : }
939 : }
940 : case t_RFRAC: {
941 7 : GEN a = gel(x,1), b = gel(x,2);
942 7 : y = famat_inv_shallow(factor(b));
943 7 : if (typ(a)==t_POL && varn(a)==varn(b)) y = famat_mul_shallow(factor(a),y);
944 7 : return gerepilecopy(av, y);
945 : }
946 : case t_COMPLEX:
947 42 : y = gauss_factor(x); if (y) return y;
948 : /* fall through */
949 : }
950 0 : pari_err_TYPE("factor",x);
951 : return NULL; /* LCOV_EXCL_LINE */
952 : }
953 :
954 : /*******************************************************************/
955 : /* */
956 : /* ROOTS --> MONIC POLYNOMIAL */
957 : /* */
958 : /*******************************************************************/
959 : static GEN
960 634070 : normalized_mul(void *E, GEN x, GEN y)
961 : {
962 634070 : long a = gel(x,1)[1], b = gel(y,1)[1];
963 : (void) E;
964 1268140 : return mkvec2(mkvecsmall(a + b),
965 1268140 : RgX_mul_normalized(gel(x,2),a, gel(y,2),b));
966 : }
967 : /* L = [Vecsmall([a]), A], with a > 0, A an RgX, deg(A) < a; return X^a + A */
968 : static GEN
969 320499 : normalized_to_RgX(GEN L)
970 : {
971 320499 : long i, a = gel(L,1)[1];
972 320499 : GEN A = gel(L,2);
973 320499 : GEN z = cgetg(a + 3, t_POL);
974 320499 : z[1] = evalsigne(1) | evalvarn(varn(A));
975 320499 : for (i = 2; i < lg(A); i++) gel(z,i) = gcopy(gel(A,i));
976 320499 : for ( ; i < a+2; i++) gel(z,i) = gen_0;
977 320499 : gel(z,i) = gen_1; return z;
978 : }
979 :
980 : /* compute prod (x - a[i]) */
981 : GEN
982 276716 : roots_to_pol(GEN a, long v)
983 : {
984 276716 : pari_sp av = avma;
985 276716 : long i, k, lx = lg(a);
986 : GEN L;
987 276716 : if (lx == 1) return pol_1(v);
988 276681 : L = cgetg(lx, t_VEC);
989 589928 : for (k=1,i=1; i<lx-1; i+=2)
990 : {
991 313247 : GEN s = gel(a,i), t = gel(a,i+1);
992 313247 : GEN x0 = gmul(s,t);
993 313247 : GEN x1 = gneg(gadd(s,t));
994 313247 : gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
995 : }
996 536299 : if (i < lx) gel(L,k++) = mkvec2(mkvecsmall(1),
997 259618 : scalarpol_shallow(gneg(gel(a,i)), v));
998 276681 : setlg(L, k); L = gen_product(L, NULL, normalized_mul);
999 276681 : return gerepileupto(av, normalized_to_RgX(L));
1000 : }
1001 :
1002 : /* prod_{i=1..r1} (x - a[i]) prod_{i=1..r2} (x - a[i])(x - conj(a[i]))*/
1003 : GEN
1004 43818 : roots_to_pol_r1(GEN a, long v, long r1)
1005 : {
1006 43818 : pari_sp av = avma;
1007 43818 : long i, k, lx = lg(a);
1008 : GEN L;
1009 43818 : if (lx == 1) return pol_1(v);
1010 43818 : L = cgetg(lx, t_VEC);
1011 299873 : for (k=1,i=1; i<r1; i+=2)
1012 : {
1013 256055 : GEN s = gel(a,i), t = gel(a,i+1);
1014 256055 : GEN x0 = gmul(s,t);
1015 256055 : GEN x1 = gneg(gadd(s,t));
1016 256055 : gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
1017 : }
1018 64272 : if (i < r1+1) gel(L,k++) = mkvec2(mkvecsmall(1),
1019 20454 : scalarpol_shallow(gneg(gel(a,i)), v));
1020 149013 : for (i=r1+1; i<lx; i++)
1021 : {
1022 105195 : GEN s = gel(a,i);
1023 105195 : GEN x0 = gnorm(s);
1024 105195 : GEN x1 = gneg(gtrace(s));
1025 105195 : gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
1026 : }
1027 43818 : setlg(L, k); L = gen_product(L, NULL, normalized_mul);
1028 43818 : return gerepileupto(av, normalized_to_RgX(L));
1029 : }
1030 :
1031 : /*******************************************************************/
1032 : /* */
1033 : /* FACTORBACK */
1034 : /* */
1035 : /*******************************************************************/
1036 : static GEN
1037 7 : idmulred(void *nf, GEN x, GEN y) { return idealmulred((GEN) nf, x, y); }
1038 : static GEN
1039 609 : idpowred(void *nf, GEN x, GEN n) { return idealpowred((GEN) nf, x, n); }
1040 : static GEN
1041 2198 : idmul(void *nf, GEN x, GEN y) { return idealmul((GEN) nf, x, y); }
1042 : static GEN
1043 3423 : idpow(void *nf, GEN x, GEN n) { return idealpow((GEN) nf, x, n); }
1044 : static GEN
1045 4164 : eltmul(void *nf, GEN x, GEN y) { return nfmul((GEN) nf, x, y); }
1046 : static GEN
1047 16940 : eltpow(void *nf, GEN x, GEN n) { return nfpow((GEN) nf, x, n); }
1048 : static GEN
1049 54162803 : mul(void *a, GEN x, GEN y) { (void)a; return gmul(x,y);}
1050 : static GEN
1051 78776658 : powi(void *a, GEN x, GEN y) { (void)a; return powgi(x,y);}
1052 : static GEN
1053 16094057 : Fpmul(void *a, GEN x, GEN y) { return Fp_mul(x,y,(GEN)a); }
1054 : static GEN
1055 119504 : Fppow(void *a, GEN x, GEN n) { return Fp_pow(x,n,(GEN)a); }
1056 :
1057 : /* [L,e] = [fa, NULL] or [elts, NULL] or [elts, exponents] */
1058 : GEN
1059 28043739 : gen_factorback(GEN L, GEN e, GEN (*_mul)(void*,GEN,GEN),
1060 : GEN (*_pow)(void*,GEN,GEN), void *data)
1061 : {
1062 28043739 : pari_sp av = avma;
1063 : long k, l, lx;
1064 : GEN p,x;
1065 :
1066 28043739 : if (e) /* supplied vector of exponents */
1067 852894 : p = L;
1068 : else
1069 : {
1070 27190845 : switch(typ(L)) {
1071 : case t_VEC:
1072 : case t_COL: /* product of the L[i] */
1073 3274378 : return gerepileupto(av, gen_product(L, data, _mul));
1074 : case t_MAT: /* genuine factorization */
1075 23916467 : l = lg(L);
1076 23916467 : if (l == 1) return gen_1;
1077 23916460 : if (l == 3) break;
1078 : /*fall through*/
1079 : default:
1080 7 : pari_err_TYPE("factorback [not a factorization]", L);
1081 : }
1082 23916453 : p = gel(L,1);
1083 23916453 : e = gel(L,2);
1084 : }
1085 : /* p = elts, e = expo */
1086 24769347 : lx = lg(p);
1087 : /* check whether e is an integral vector of correct length */
1088 24769347 : switch(typ(e))
1089 : {
1090 : case t_VECSMALL:
1091 294 : if (lx != lg(e))
1092 0 : pari_err_TYPE("factorback [not an exponent vector]", e);
1093 294 : if (lx == 1) return gen_1;
1094 266 : x = cgetg(lx,t_VEC);
1095 882 : for (l=1,k=1; k<lx; k++)
1096 616 : if (e[k]) gel(x,l++) = _pow(data, gel(p,k), stoi(e[k]));
1097 266 : break;
1098 : case t_VEC: case t_COL:
1099 24769053 : if (lx != lg(e) || !RgV_is_ZV(e))
1100 7 : pari_err_TYPE("factorback [not an exponent vector]", e);
1101 24769046 : if (lx == 1) return gen_1;
1102 24749570 : x = cgetg(lx,t_VEC);
1103 103836160 : for (l=1,k=1; k<lx; k++)
1104 79086590 : if (signe(gel(e,k))) gel(x,l++) = _pow(data, gel(p,k), gel(e,k));
1105 24749570 : break;
1106 : default:
1107 0 : pari_err_TYPE("factorback [not an exponent vector]", e);
1108 0 : return NULL;
1109 : }
1110 24749836 : x[0] = evaltyp(t_VEC) | _evallg(l);
1111 24749836 : return gerepileupto(av, gen_product(x, data, _mul));
1112 : }
1113 :
1114 : GEN
1115 1911 : idealfactorback(GEN nf, GEN L, GEN e, int red)
1116 : {
1117 1911 : nf = checknf(nf);
1118 1911 : if (red) return gen_factorback(L, e, &idmulred, &idpowred, (void*)nf);
1119 1309 : else return gen_factorback(L, e, &idmul, &idpow, (void*)nf);
1120 : }
1121 :
1122 : GEN
1123 13286 : nffactorback(GEN nf, GEN L, GEN e)
1124 13286 : { return gen_factorback(L, e, &eltmul, &eltpow, (void*)checknf(nf)); }
1125 :
1126 : GEN
1127 3369663 : FpV_factorback(GEN L, GEN e, GEN p)
1128 3369663 : { return gen_factorback(L, e, &Fpmul, &Fppow, (void*)p); }
1129 :
1130 : GEN
1131 24643682 : factorback2(GEN L, GEN e) { return gen_factorback(L, e, &mul, &powi, NULL); }
1132 : GEN
1133 1395129 : factorback(GEN fa) { return factorback2(fa, NULL); }
1134 :
1135 : GEN
1136 14 : vecprod(GEN v)
1137 : {
1138 14 : pari_sp av = avma;
1139 14 : if (!is_vec_t(typ(v)))
1140 0 : pari_err_TYPE("vecprod", v);
1141 14 : if (lg(v) == 1) return gen_1;
1142 7 : return gerepilecopy(av, gen_product(v, NULL, mul));
1143 : }
1144 :
1145 : static int
1146 63 : RgX_is_irred_i(GEN x)
1147 : {
1148 : GEN y, p, pol;
1149 63 : long l = lg(x), pa;
1150 :
1151 63 : if (!signe(x) || l <= 3) return 0;
1152 63 : switch(RgX_type(x,&p,&pol,&pa))
1153 : {
1154 14 : case t_INTMOD: return FpX_is_irred(RgX_to_FpX(x,p), p);
1155 0 : case t_COMPLEX: return l == 4;
1156 : case t_REAL:
1157 0 : if (l == 4) return 1;
1158 0 : if (l > 5) return 0;
1159 0 : return gsigne(RgX_disc(x)) > 0;
1160 : }
1161 49 : y = factor(x);
1162 49 : return (lg(gcoeff(y,1,1))==l);
1163 : }
1164 : static int
1165 63 : RgX_is_irred(GEN x)
1166 : {
1167 63 : pari_sp av = avma;
1168 63 : int r = RgX_is_irred_i(x);
1169 63 : avma = av; return r;
1170 : }
1171 : long
1172 63 : isirreducible(GEN x)
1173 : {
1174 63 : switch(typ(x))
1175 : {
1176 0 : case t_INT: case t_REAL: case t_FRAC: return 0;
1177 63 : case t_POL: return RgX_is_irred(x);
1178 : }
1179 0 : pari_err_TYPE("isirreducible",x);
1180 0 : return 0;
1181 : }
1182 :
1183 : /*******************************************************************/
1184 : /* */
1185 : /* GENERIC GCD */
1186 : /* */
1187 : /*******************************************************************/
1188 : /* x is a COMPLEX or a QUAD */
1189 : static GEN
1190 2499 : triv_cont_gcd(GEN x, GEN y)
1191 : {
1192 2499 : pari_sp av = avma;
1193 : GEN c;
1194 2499 : if (typ(x)==t_COMPLEX)
1195 : {
1196 2177 : GEN a = gel(x,1), b = gel(x,2);
1197 2177 : if (typ(a) == t_REAL || typ(b) == t_REAL) return gen_1;
1198 21 : c = ggcd(a,b);
1199 : }
1200 : else
1201 322 : c = ggcd(gel(x,2),gel(x,3));
1202 343 : return gerepileupto(av, ggcd(c,y));
1203 : }
1204 :
1205 : /* y is a PADIC, x a rational number or an INTMOD */
1206 : static GEN
1207 2684 : padic_gcd(GEN x, GEN y)
1208 : {
1209 2684 : GEN p = gel(y,2);
1210 2684 : long v = gvaluation(x,p), w = valp(y);
1211 2684 : if (w < v) v = w;
1212 2684 : return powis(p, v);
1213 : }
1214 :
1215 : /* x,y in Z[i], at least one of which is t_COMPLEX */
1216 : static GEN
1217 385 : gauss_gcd(GEN x, GEN y)
1218 : {
1219 385 : pari_sp av = avma;
1220 : GEN dx, dy;
1221 385 : dx = denom(x); x = gmul(x, dx);
1222 385 : dy = denom(y); y = gmul(y, dy);
1223 1197 : while (!gequal0(y))
1224 : {
1225 427 : GEN z = gsub(x, gmul(ground(gdiv(x,y)), y));
1226 427 : x = y; y = z;
1227 : }
1228 385 : x = gauss_normal(x);
1229 385 : if (typ(x) == t_COMPLEX)
1230 : {
1231 196 : if (gequal0(gel(x,2))) x = gel(x,1);
1232 196 : else if (gequal0(gel(x,1))) x = gel(x,2);
1233 : }
1234 385 : return gerepileupto(av, gdiv(x, lcmii(dx, dy)));
1235 : }
1236 :
1237 : static int
1238 2520 : c_is_rational(GEN x)
1239 2520 : { return is_rational_t(typ(gel(x,1))) && is_rational_t(typ(gel(x,2))); }
1240 : static GEN
1241 1288 : c_zero_gcd(GEN c)
1242 : {
1243 1288 : GEN x = gel(c,1), y = gel(c,2);
1244 1288 : long tx = typ(x), ty = typ(y);
1245 1288 : if (tx == t_REAL || ty == t_REAL) return gen_1;
1246 42 : if (tx == t_PADIC || tx == t_INTMOD
1247 35 : || ty == t_PADIC || ty == t_INTMOD) return ggcd(x, y);
1248 35 : return gauss_gcd(c, gen_0);
1249 : }
1250 :
1251 : /* gcd(x, 0) */
1252 : static GEN
1253 8204365 : zero_gcd(GEN x)
1254 : {
1255 : pari_sp av;
1256 8204365 : switch(typ(x))
1257 : {
1258 24268 : case t_INT: return absi(x);
1259 1295 : case t_FRAC: return absfrac(x);
1260 1288 : case t_COMPLEX: return c_zero_gcd(x);
1261 7210 : case t_REAL: return gen_1;
1262 735 : case t_PADIC: return powis(gel(x,2), valp(x));
1263 210 : case t_SER: return pol_xn(valp(x), varn(x));
1264 : case t_POLMOD: {
1265 8338 : GEN d = gel(x,2);
1266 8338 : if (typ(d) == t_POL && varn(d) == varn(gel(x,1))) return content(d);
1267 28 : return isinexact(d)? zero_gcd(d): gcopy(d);
1268 : }
1269 : case t_POL:
1270 7921099 : if (!isinexact(x)) break;
1271 14 : av = avma;
1272 14 : return gerepileupto(av, monomialcopy(content(x), RgX_val(x), varn(x)));
1273 :
1274 : case t_RFRAC:
1275 217137 : if (!isinexact(x)) break;
1276 0 : av = avma;
1277 0 : return gerepileupto(av, gdiv(zero_gcd(gel(x,1)), gel(x,2)));
1278 : }
1279 8161007 : return gcopy(x);
1280 : }
1281 : /* z is an exact zero, t_INT, t_INTMOD or t_FFELT */
1282 : static GEN
1283 8749464 : zero_gcd2(GEN y, GEN z)
1284 : {
1285 : pari_sp av;
1286 8749464 : switch(typ(z))
1287 : {
1288 8181890 : case t_INT: return zero_gcd(y);
1289 : case t_INTMOD:
1290 566531 : av = avma;
1291 566531 : return gerepileupto(av, gmul(y, mkintmod(gen_1,gel(z,1))));
1292 : case t_FFELT:
1293 1043 : av = avma;
1294 1043 : return gerepileupto(av, gmul(y, FF_1(z)));
1295 : default:
1296 0 : pari_err_TYPE("zero_gcd", z);
1297 : }
1298 0 : return NULL;
1299 : }
1300 : /* tx = t_RFRAC, y considered as constant */
1301 : static GEN
1302 864876 : cont_gcd_rfrac(GEN x, GEN y)
1303 : {
1304 864876 : pari_sp av = avma;
1305 864876 : GEN cx; x = primitive_part(x, &cx);
1306 864876 : return gerepileupto(av, gred_rfrac_simple(ggcd(cx? cx: gen_1, y), gel(x,2)));
1307 : }
1308 : /* tx = t_POL, y considered as constant */
1309 : static GEN
1310 2256744 : cont_gcd_pol(GEN x, GEN y)
1311 : {
1312 2256744 : pari_sp av = avma;
1313 2256744 : return gerepileupto(av, scalarpol(ggcd(content(x),y), varn(x)));
1314 : }
1315 : /* !is_const_t(tx), tx != t_POL,t_RFRAC, y considered as constant */
1316 : static GEN
1317 6588 : cont_gcd_gen(GEN x, GEN y)
1318 : {
1319 6588 : pari_sp av = avma;
1320 6588 : return gerepileupto(av, ggcd(content(x),y));
1321 : }
1322 : /* !is_const(tx), y considered as constant */
1323 : static GEN
1324 3071427 : cont_gcd(GEN x, long tx, GEN y)
1325 : {
1326 3071427 : switch(tx)
1327 : {
1328 808109 : case t_RFRAC: return cont_gcd_rfrac(x,y);
1329 2256730 : case t_POL: return cont_gcd_pol(x,y);
1330 6588 : default: return cont_gcd_gen(x,y);
1331 : }
1332 : }
1333 : static GEN
1334 9780703 : gcdiq(GEN x, GEN y)
1335 : {
1336 : GEN z;
1337 9780703 : if (!signe(x)) return Q_abs(y);
1338 2458920 : z = cgetg(3,t_FRAC);
1339 2458920 : gel(z,1) = gcdii(x,gel(y,1));
1340 2458920 : gel(z,2) = icopy(gel(y,2));
1341 2458920 : return z;
1342 : }
1343 : static GEN
1344 10946197 : gcdqq(GEN x, GEN y)
1345 : {
1346 10946197 : GEN z = cgetg(3,t_FRAC);
1347 10946197 : gel(z,1) = gcdii(gel(x,1), gel(y,1));
1348 10946197 : gel(z,2) = lcmii(gel(x,2), gel(y,2));
1349 10946197 : return z;
1350 : }
1351 : /* assume x,y t_INT or t_FRAC */
1352 : GEN
1353 104851665 : Q_gcd(GEN x, GEN y)
1354 : {
1355 104851665 : long tx = typ(x), ty = typ(y);
1356 104851665 : if (tx == t_INT)
1357 84909721 : { return (ty == t_INT)? gcdii(x,y): gcdiq(x,y); }
1358 : else
1359 19941944 : { return (ty == t_INT)? gcdiq(y,x): gcdqq(x,y); }
1360 : }
1361 :
1362 : GEN
1363 24077207 : ggcd(GEN x, GEN y)
1364 : {
1365 24077207 : long vx, vy, tx = typ(x), ty = typ(y);
1366 : pari_sp av, tetpil;
1367 : GEN p1,z;
1368 :
1369 48154414 : if (is_noncalc_t(tx) || is_matvec_t(tx) ||
1370 48154414 : is_noncalc_t(ty) || is_matvec_t(ty)) pari_err_TYPE2("gcd",x,y);
1371 24077207 : if (tx>ty) { swap(x,y); lswap(tx,ty); }
1372 : /* tx <= ty */
1373 24077207 : z = gisexactzero(x); if (z) return zero_gcd2(y,z);
1374 21032365 : z = gisexactzero(y); if (z) return zero_gcd2(x,z);
1375 15327743 : if (is_const_t(tx))
1376 : {
1377 8391780 : if (ty == tx) switch(tx)
1378 : {
1379 : case t_INT:
1380 3210277 : return gcdii(x,y);
1381 :
1382 2074065 : case t_INTMOD: z=cgetg(3,t_INTMOD);
1383 2074065 : if (equalii(gel(x,1),gel(y,1)))
1384 2074065 : gel(z,1) = icopy(gel(x,1));
1385 : else
1386 0 : gel(z,1) = gcdii(gel(x,1),gel(y,1));
1387 2074065 : if (gequal1(gel(z,1))) gel(z,2) = gen_0;
1388 : else
1389 : {
1390 2074065 : av = avma; p1 = gcdii(gel(z,1),gel(x,2));
1391 2074065 : if (!is_pm1(p1))
1392 : {
1393 0 : p1 = gcdii(p1,gel(y,2));
1394 0 : if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
1395 0 : else p1 = gerepileuptoint(av, p1);
1396 : }
1397 2074065 : gel(z,2) = p1;
1398 : }
1399 2074065 : return z;
1400 :
1401 : case t_FRAC:
1402 32571 : return gcdqq(x,y);
1403 :
1404 : case t_FFELT:
1405 8841 : if (!FF_samefield(x,y)) pari_err_OP("gcd",x,y);
1406 8841 : return FF_equal0(x) && FF_equal0(y)? FF_zero(y): FF_1(y);
1407 :
1408 : case t_COMPLEX:
1409 21 : if (c_is_rational(x) && c_is_rational(y)) return gauss_gcd(x,y);
1410 7 : return triv_cont_gcd(y,x);
1411 :
1412 : case t_PADIC:
1413 14 : if (!equalii(gel(x,2),gel(y,2))) return gen_1;
1414 7 : return powis(gel(y,2), minss(valp(x), valp(y)));
1415 :
1416 : case t_QUAD:
1417 133 : av=avma; p1=gdiv(x,y);
1418 133 : if (gequal0(gel(p1,3)))
1419 : {
1420 14 : p1=gel(p1,2);
1421 14 : if (typ(p1)==t_INT) { avma=av; return gcopy(y); }
1422 7 : tetpil=avma; return gerepile(av,tetpil, gdiv(y,gel(p1,2)));
1423 : }
1424 119 : if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) {avma=av; return gcopy(y);}
1425 112 : p1 = ginv(p1); avma=av;
1426 112 : if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) return gcopy(x);
1427 105 : return triv_cont_gcd(y,x);
1428 :
1429 0 : default: return gen_1; /* t_REAL */
1430 : }
1431 3065858 : if (is_const_t(ty)) switch(tx)
1432 : {
1433 : case t_INT:
1434 21521 : switch(ty)
1435 : {
1436 1386 : case t_INTMOD: z = cgetg(3,t_INTMOD);
1437 1386 : gel(z,1) = icopy(gel(y,1)); av = avma;
1438 1386 : p1 = gcdii(gel(y,1),gel(y,2));
1439 1386 : if (!is_pm1(p1)) {
1440 14 : p1 = gcdii(x,p1);
1441 14 : if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
1442 : else
1443 14 : p1 = gerepileuptoint(av, p1);
1444 : }
1445 1386 : gel(z,2) = p1; return z;
1446 :
1447 6825 : case t_REAL: return gen_1;
1448 :
1449 : case t_FRAC:
1450 7889 : return gcdiq(x,y);
1451 :
1452 : case t_COMPLEX:
1453 2401 : if (c_is_rational(y)) return gauss_gcd(x,y);
1454 2072 : return triv_cont_gcd(y,x);
1455 :
1456 : case t_FFELT:
1457 224 : if (!FF_equal0(y)) return FF_1(y);
1458 0 : return dvdii(x, gel(y,4))? FF_zero(y): FF_1(y);
1459 :
1460 : case t_PADIC:
1461 2670 : return padic_gcd(x,y);
1462 :
1463 : case t_QUAD:
1464 126 : return triv_cont_gcd(y,x);
1465 : default:
1466 0 : pari_err_TYPE2("gcd",x,y);
1467 : }
1468 :
1469 : case t_REAL:
1470 14 : switch(ty)
1471 : {
1472 : case t_INTMOD:
1473 : case t_FFELT:
1474 14 : case t_PADIC: pari_err_TYPE2("gcd",x,y);
1475 0 : default: return gen_1;
1476 : }
1477 :
1478 : case t_INTMOD:
1479 49 : switch(ty)
1480 : {
1481 : case t_FRAC:
1482 14 : av = avma; p1=gcdii(gel(x,1),gel(y,2)); avma = av;
1483 14 : if (!is_pm1(p1)) pari_err_OP("gcd",x,y);
1484 7 : return ggcd(gel(y,1), x);
1485 :
1486 : case t_FFELT:
1487 : {
1488 14 : GEN p = gel(y,4);
1489 14 : if (!dvdii(gel(x,1), p)) pari_err_OP("gcd",x,y);
1490 7 : if (!FF_equal0(y)) return FF_1(y);
1491 0 : return dvdii(gel(x,2),p)? FF_zero(y): FF_1(y);
1492 : }
1493 :
1494 : case t_COMPLEX: case t_QUAD:
1495 14 : return triv_cont_gcd(y,x);
1496 :
1497 : case t_PADIC:
1498 7 : return padic_gcd(x,y);
1499 :
1500 0 : default: pari_err_TYPE2("gcd",x,y);
1501 : }
1502 :
1503 : case t_FRAC:
1504 210 : switch(ty)
1505 : {
1506 : case t_COMPLEX:
1507 84 : if (c_is_rational(y)) return gauss_gcd(x,y);
1508 : case t_QUAD:
1509 154 : return triv_cont_gcd(y,x);
1510 : case t_FFELT:
1511 : {
1512 42 : GEN p = gel(y,4);
1513 42 : if (dvdii(gel(x,2), p)) pari_err_OP("gcd",x,y);
1514 21 : if (!FF_equal0(y)) return FF_1(y);
1515 0 : return dvdii(gel(x,1),p)? FF_zero(y): FF_1(y);
1516 : }
1517 :
1518 : case t_PADIC:
1519 7 : return padic_gcd(x,y);
1520 :
1521 0 : default: pari_err_TYPE2("gcd",x,y);
1522 : }
1523 : case t_FFELT:
1524 70 : switch(ty)
1525 : {
1526 : case t_PADIC:
1527 : {
1528 42 : GEN p = gel(y,2);
1529 42 : long v = valp(y);
1530 42 : if (!equalii(p, gel(x,4)) || v < 0) pari_err_OP("gcd",x,y);
1531 14 : return (v && FF_equal0(x))? FF_zero(x): FF_1(x);
1532 : }
1533 28 : default: pari_err_TYPE2("gcd",x,y);
1534 : }
1535 :
1536 : case t_COMPLEX:
1537 14 : switch(ty)
1538 : {
1539 : case t_PADIC:
1540 14 : case t_QUAD: return triv_cont_gcd(x,y);
1541 0 : default: pari_err_TYPE2("gcd",x,y);
1542 : }
1543 :
1544 : case t_PADIC:
1545 7 : switch(ty)
1546 : {
1547 7 : case t_QUAD: return triv_cont_gcd(y,x);
1548 0 : default: pari_err_TYPE2("gcd",x,y);
1549 : }
1550 :
1551 0 : default: return gen_1; /* tx = t_REAL */
1552 : }
1553 3043973 : return cont_gcd(y,ty, x);
1554 : }
1555 :
1556 6935963 : if (tx == t_POLMOD)
1557 : {
1558 19358 : if (ty == t_POLMOD)
1559 : {
1560 19302 : GEN T = gel(x,1);
1561 19302 : z = cgetg(3,t_POLMOD);
1562 19302 : T = RgX_equal_var(T,gel(y,1))? RgX_copy(T): RgX_gcd(T, gel(y,1));
1563 19302 : gel(z,1) = T;
1564 19302 : if (degpol(T) <= 0) gel(z,2) = gen_0;
1565 : else
1566 : {
1567 : GEN X, Y, d;
1568 19302 : av = avma; X = gel(x,2); Y = gel(y,2);
1569 19302 : d = ggcd(content(X), content(Y));
1570 19302 : if (!gequal1(d)) { X = gdiv(X,d); Y = gdiv(Y,d); }
1571 19302 : p1 = ggcd(T, X);
1572 19302 : gel(z,2) = gerepileupto(av, gmul(d, ggcd(p1, Y)));
1573 : }
1574 19302 : return z;
1575 : }
1576 56 : vx = varn(gel(x,1));
1577 56 : switch(ty)
1578 : {
1579 : case t_POL:
1580 28 : vy = varn(y);
1581 28 : if (varncmp(vy,vx) < 0) return cont_gcd_pol(y, x);
1582 14 : z = cgetg(3,t_POLMOD);
1583 14 : gel(z,1) = RgX_copy(gel(x,1));
1584 14 : av = avma; p1 = ggcd(gel(x,1),gel(x,2));
1585 14 : gel(z,2) = gerepileupto(av, ggcd(p1,y));
1586 14 : return z;
1587 :
1588 : case t_RFRAC:
1589 28 : vy = varn(gel(y,2));
1590 28 : if (varncmp(vy,vx) < 0) return cont_gcd_rfrac(y, x);
1591 28 : av = avma;
1592 28 : p1 = ggcd(gel(x,1),gel(y,2));
1593 28 : if (degpol(p1)) pari_err_OP("gcd",x,y);
1594 21 : avma = av; return gdiv(ggcd(gel(y,1),x), content(gel(y,2)));
1595 : }
1596 : }
1597 :
1598 6916605 : vx = gvar(x);
1599 6916605 : vy = gvar(y);
1600 6916605 : if (varncmp(vy, vx) < 0) return cont_gcd(y,ty, x);
1601 6904600 : if (varncmp(vy, vx) > 0) return cont_gcd(x,tx, y);
1602 :
1603 : /* vx = vy: same main variable */
1604 6889151 : switch(tx)
1605 : {
1606 : case t_POL:
1607 6728245 : switch(ty)
1608 : {
1609 6671464 : case t_POL: return RgX_gcd(x,y);
1610 : case t_SER:
1611 14 : z = ggcd(content(x), content(y));
1612 14 : return monomialcopy(z, minss(valp(y),gval(x,vx)), vx);
1613 56767 : case t_RFRAC: return cont_gcd_rfrac(y, x);
1614 : }
1615 0 : break;
1616 :
1617 : case t_SER:
1618 14 : z = ggcd(content(x), content(y));
1619 14 : switch(ty)
1620 : {
1621 7 : case t_SER: return monomialcopy(z, minss(valp(x),valp(y)), vx);
1622 7 : case t_RFRAC: return monomialcopy(z, minss(valp(x),gval(y,vx)), vx);
1623 : }
1624 0 : break;
1625 :
1626 : case t_RFRAC:
1627 : {
1628 160892 : GEN xd = gel(x,2), yd = gel(y,2);
1629 160892 : if (ty != t_RFRAC) pari_err_TYPE2("gcd",x,y);
1630 160892 : z = cgetg(3,t_RFRAC); av = avma;
1631 160892 : gel(z,2) = gerepileupto(av, RgX_mul(xd, RgX_div(yd, RgX_gcd(xd, yd))));
1632 160892 : gel(z,1) = ggcd(gel(x,1), gel(y,1)); return z;
1633 : }
1634 : }
1635 0 : pari_err_TYPE2("gcd",x,y);
1636 : return NULL; /* LCOV_EXCL_LINE */
1637 : }
1638 : GEN
1639 3960 : ggcd0(GEN x, GEN y) { return y? ggcd(x,y): content(x); }
1640 :
1641 : static GEN
1642 133 : fix_lcm(GEN x)
1643 : {
1644 : GEN t;
1645 133 : switch(typ(x))
1646 : {
1647 35 : case t_INT: if (signe(x)<0) x = negi(x);
1648 35 : break;
1649 : case t_POL:
1650 98 : if (lg(x) <= 2) break;
1651 98 : t = leading_coeff(x);
1652 98 : if (typ(t) == t_INT && signe(t) < 0) x = gneg(x);
1653 : }
1654 133 : return x;
1655 : }
1656 : GEN
1657 2905 : glcm0(GEN x, GEN y)
1658 : {
1659 2905 : if (!y) return fix_lcm(gassoc_proto(glcm,x,y));
1660 2821 : return glcm(x,y);
1661 : }
1662 : GEN
1663 3199 : glcm(GEN x, GEN y)
1664 : {
1665 : pari_sp av;
1666 : GEN z;
1667 3199 : if (typ(x)==t_INT && typ(y)==t_INT) return lcmii(x,y);
1668 63 : av = avma; z = ggcd(x,y);
1669 63 : if (!gequal1(z))
1670 : {
1671 63 : if (gequal0(z)) { avma = av; return gmul(x,y); }
1672 49 : y = gdiv(y,z);
1673 : }
1674 49 : return gerepileupto(av, fix_lcm(gmul(x,y)));
1675 : }
1676 :
1677 : /* x + r ~ x ? Assume x,r are t_POL, deg(r) <= deg(x) */
1678 : static int
1679 302589 : pol_approx0(GEN r, GEN x, int exact)
1680 : {
1681 : long i, lx,lr;
1682 302589 : if (exact) return gequal0(r);
1683 0 : lx = lg(x);
1684 0 : lr = lg(r); if (lr < lx) lx = lr;
1685 0 : for (i=2; i<lx; i++)
1686 0 : if (!approx_0(gel(r,i), gel(x,i))) return 0;
1687 0 : return 1;
1688 : }
1689 :
1690 : GEN
1691 87164 : RgX_gcd_simple(GEN x, GEN y)
1692 : {
1693 87164 : pari_sp av1, av = avma;
1694 87164 : GEN r, yorig = y;
1695 87164 : int exact = !(isinexactreal(x) || isinexactreal(y));
1696 :
1697 : for(;;)
1698 : {
1699 302589 : av1 = avma; r = RgX_rem(x,y);
1700 302589 : if (pol_approx0(r, x, exact))
1701 : {
1702 87164 : avma = av1;
1703 87164 : if (y == yorig) return RgX_copy(y);
1704 60942 : y = normalizepol_approx(y, lg(y));
1705 60942 : if (lg(y) == 3) { avma = av; return pol_1(varn(x)); }
1706 4340 : return gerepileupto(av,y);
1707 : }
1708 215425 : x = y; y = r;
1709 215425 : if (gc_needed(av,1)) {
1710 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd_simple");
1711 0 : gerepileall(av,2, &x,&y);
1712 : }
1713 215425 : }
1714 : }
1715 : GEN
1716 0 : RgX_extgcd_simple(GEN a, GEN b, GEN *pu, GEN *pv)
1717 : {
1718 0 : pari_sp av = avma;
1719 : GEN q, r, d, d1, u, v, v1;
1720 0 : int exact = !(isinexactreal(a) || isinexactreal(b));
1721 :
1722 0 : d = a; d1 = b; v = gen_0; v1 = gen_1;
1723 : for(;;)
1724 : {
1725 0 : if (pol_approx0(d1, a, exact)) break;
1726 0 : q = poldivrem(d,d1, &r);
1727 0 : v = gsub(v, gmul(q,v1));
1728 0 : u=v; v=v1; v1=u;
1729 0 : u=r; d=d1; d1=u;
1730 0 : }
1731 0 : u = gsub(d, gmul(b,v));
1732 0 : u = RgX_div(u,a);
1733 :
1734 0 : gerepileall(av, 3, &u,&v,&d);
1735 0 : *pu = u;
1736 0 : *pv = v; return d;
1737 : }
1738 : /*******************************************************************/
1739 : /* */
1740 : /* CONTENT / PRIMITIVE PART */
1741 : /* */
1742 : /*******************************************************************/
1743 :
1744 : GEN
1745 77308784 : content(GEN x)
1746 : {
1747 77308784 : long lx, i, t, tx = typ(x);
1748 77308784 : pari_sp av = avma;
1749 : GEN c;
1750 :
1751 77308784 : if (is_scalar_t(tx)) return zero_gcd(x);
1752 77295703 : switch(tx)
1753 : {
1754 : case t_RFRAC:
1755 : {
1756 864890 : GEN n = gel(x,1), d = gel(x,2);
1757 : /* -- varncmp(vn, vd) < 0 can't happen
1758 : * -- if n is POLMOD, its main variable (in the sense of gvar2)
1759 : * has lower priority than denominator */
1760 864890 : if (typ(n) == t_POLMOD || varncmp(gvar(n), varn(d)) > 0)
1761 824566 : n = isinexact(n)? zero_gcd(n): gcopy(n);
1762 : else
1763 40324 : n = content(n);
1764 864890 : return gerepileupto(av, gdiv(n, content(d)));
1765 : }
1766 :
1767 : case t_VEC: case t_COL:
1768 7566135 : lx = lg(x); if (lx==1) return gen_0;
1769 7566128 : break;
1770 :
1771 : case t_MAT:
1772 : {
1773 : long hx, j;
1774 238 : lx = lg(x);
1775 238 : if (lx == 1) return gen_0;
1776 231 : hx = lgcols(x);
1777 231 : if (hx == 1) return gen_0;
1778 224 : if (lx == 2) { x = gel(x,1); lx = lg(x); break; }
1779 224 : if (hx == 2) { x = row_i(x, 1, 1, lx-1); break; }
1780 224 : c = content(gel(x,1));
1781 448 : for (j=2; j<lx; j++)
1782 224 : for (i=1; i<hx; i++) c = ggcd(c,gcoeff(x,i,j));
1783 224 : if (typ(c) == t_INTMOD || isinexact(c)) { avma=av; return gen_1; }
1784 224 : return gerepileupto(av,c);
1785 : }
1786 :
1787 : case t_POL: case t_SER:
1788 68864405 : lx = lg(x); if (lx == 2) return gen_0;
1789 68864391 : break;
1790 21 : case t_VECSMALL: return utoi(zv_content(x));
1791 : case t_QFR: case t_QFI:
1792 14 : lx = 4; break;
1793 :
1794 0 : default: pari_err_TYPE("content",x);
1795 : return NULL; /* LCOV_EXCL_LINE */
1796 : }
1797 224952391 : for (i=lontyp[tx]; i<lx; i++)
1798 161441355 : if (typ(gel(x,i)) != t_INT) break;
1799 76430533 : lx--; c = gel(x,lx);
1800 76430533 : t = typ(c); if (is_matvec_t(t)) c = content(c);
1801 76430533 : if (i > lx)
1802 : { /* integer coeffs */
1803 128282365 : while (lx-- > lontyp[tx])
1804 : {
1805 64013813 : c = gcdii(c, gel(x,lx));
1806 64013813 : if (is_pm1(c)) { avma=av; return gen_1; }
1807 : }
1808 : }
1809 : else
1810 : {
1811 12919497 : if (isinexact(c)) c = zero_gcd(c);
1812 46558949 : while (lx-- > lontyp[tx])
1813 : {
1814 20719955 : GEN d = gel(x,lx);
1815 20719955 : t = typ(d); if (is_matvec_t(t)) d = content(d);
1816 20719955 : c = ggcd(c, d);
1817 : }
1818 12919497 : if (isinexact(c)) { avma=av; return gen_1; }
1819 : }
1820 13677013 : switch(typ(c))
1821 : {
1822 : case t_INT:
1823 769374 : if (signe(c) < 0) c = negi(c);
1824 769374 : break;
1825 : case t_VEC: case t_COL: case t_MAT:
1826 0 : pari_err_TYPE("content",x);
1827 : }
1828 :
1829 13677013 : return av==avma? gcopy(c): gerepileupto(av,c);
1830 : }
1831 :
1832 : GEN
1833 1007796 : primitive_part(GEN x, GEN *ptc)
1834 : {
1835 1007796 : pari_sp av = avma;
1836 1007796 : GEN c = content(x);
1837 1007796 : if (gequal1(c)) { avma = av; c = NULL; }
1838 33135 : else if (!gequal0(c)) x = gdiv(x,c);
1839 1007796 : if (ptc) *ptc = c;
1840 1007796 : return x;
1841 : }
1842 : GEN
1843 2835 : primpart(GEN x) { return primitive_part(x, NULL); }
1844 :
1845 : static GEN
1846 24079382 : Q_content_v(GEN x, long i, long l)
1847 : {
1848 24079382 : pari_sp av = avma;
1849 24079382 : GEN d = Q_content_safe(gel(x,i));
1850 24079382 : if (!d) return NULL;
1851 125570630 : for (i++; i < l; i++)
1852 : {
1853 101491451 : GEN c = Q_content_safe(gel(x,i));
1854 101491451 : if (!c) return NULL;
1855 101491451 : d = Q_gcd(d, c);
1856 : }
1857 24079179 : return gerepileupto(av, d);
1858 : }
1859 : /* As content(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
1860 : * of Q(x2,...,xn)[x1] */
1861 : GEN
1862 140209274 : Q_content_safe(GEN x)
1863 : {
1864 : long l;
1865 140209274 : switch(typ(x))
1866 : {
1867 101892561 : case t_INT: return absi(x);
1868 13939592 : case t_FRAC: return absfrac(x);
1869 : case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
1870 12893790 : l = lg(x); return l==1? gen_1: Q_content_v(x, 1, l);
1871 : case t_POL:
1872 11461099 : l = lg(x); return l==2? gen_0: Q_content_v(x, 2, l);
1873 22092 : case t_POLMOD: return Q_content_safe(gel(x,2));
1874 : }
1875 140 : return NULL;
1876 : }
1877 : GEN
1878 99951 : Q_content(GEN x)
1879 : {
1880 99951 : GEN c = Q_content_safe(x);
1881 99951 : if (!c) pari_err_TYPE("Q_content",x);
1882 99951 : return c;
1883 : }
1884 :
1885 : GEN
1886 38678 : ZX_content(GEN x)
1887 : {
1888 38678 : long i, l = lg(x);
1889 : GEN d;
1890 : pari_sp av;
1891 :
1892 38678 : if (l == 2) return gen_0;
1893 38678 : d = gel(x,2);
1894 38678 : if (l == 3) return absi(d);
1895 32775 : av = avma;
1896 32775 : for (i=3; !is_pm1(d) && i<l; i++) d = gcdii(d, gel(x,i));
1897 32775 : if (signe(d) < 0) d = absi(d);
1898 32775 : return gerepileuptoint(av, d);
1899 : }
1900 :
1901 : static GEN
1902 82110 : Z_content_v(GEN x, long i, long l)
1903 : {
1904 82110 : pari_sp av = avma;
1905 82110 : GEN d = Z_content(gel(x,i));
1906 82110 : if (!d) return NULL;
1907 356867 : for (i++; i<l; i++)
1908 : {
1909 277984 : GEN c = Z_content(gel(x,i));
1910 277984 : if (!c) return NULL;
1911 277858 : d = gcdii(d, c); if (is_pm1(d)) return NULL;
1912 277816 : if ((i & 255) == 0) d = gerepileuptoint(av, d);
1913 : }
1914 78883 : return gerepileuptoint(av, d);
1915 : }
1916 : /* return NULL for 1 */
1917 : GEN
1918 361326 : Z_content(GEN x)
1919 : {
1920 : long l;
1921 361326 : switch(typ(x))
1922 : {
1923 : case t_INT:
1924 264110 : if (is_pm1(x)) return NULL;
1925 263340 : return absi(x);
1926 : case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
1927 7903 : l = lg(x); return l==1? NULL: Z_content_v(x, 1, l);
1928 : case t_POL:
1929 89313 : l = lg(x); return l==2? gen_0: Z_content_v(x, 2, l);
1930 0 : case t_POLMOD: return Z_content(gel(x,2));
1931 : }
1932 0 : pari_err_TYPE("Z_content", x);
1933 : return NULL; /* LCOV_EXCL_LINE */
1934 : }
1935 :
1936 : static GEN
1937 10761295 : Q_denom_v(GEN x, long i, long l)
1938 : {
1939 10761295 : pari_sp av = avma;
1940 10761295 : GEN d = Q_denom(gel(x,i));
1941 46591719 : for (i++; i<l; i++)
1942 : {
1943 35830424 : GEN D = Q_denom(gel(x,i));
1944 35830424 : if (D != gen_1) d = lcmii(d, D);
1945 35830424 : if ((i & 255) == 0) d = gerepileuptoint(av, d);
1946 : }
1947 10761295 : return gerepileuptoint(av, d);
1948 : }
1949 : /* NOT MEMORY CLEAN (because of t_FRAC).
1950 : * As denom(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
1951 : * of Q(x2,...,xn)[x1] */
1952 : GEN
1953 69170418 : Q_denom(GEN x)
1954 : {
1955 : long l;
1956 69170418 : switch(typ(x))
1957 : {
1958 45455912 : case t_INT: return gen_1;
1959 12806358 : case t_FRAC: return gel(x,2);
1960 : case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
1961 8191698 : l = lg(x); return l==1? gen_1: Q_denom_v(x, 1, l);
1962 : case t_POL:
1963 2640514 : l = lg(x); return l==2? gen_1: Q_denom_v(x, 2, l);
1964 75936 : case t_POLMOD: return Q_denom(gel(x,2));
1965 : }
1966 0 : pari_err_TYPE("Q_denom",x);
1967 : return NULL; /* LCOV_EXCL_LINE */
1968 : }
1969 :
1970 : GEN
1971 8184865 : Q_remove_denom(GEN x, GEN *ptd)
1972 : {
1973 8184865 : GEN d = Q_denom(x);
1974 8184865 : if (d == gen_1) d = NULL; else x = Q_muli_to_int(x,d);
1975 8184865 : if (ptd) *ptd = d;
1976 8184865 : return x;
1977 : }
1978 :
1979 : /* return y = x * d, assuming x rational, and d,y integral */
1980 : GEN
1981 44009800 : Q_muli_to_int(GEN x, GEN d)
1982 : {
1983 : long i, l;
1984 : GEN y, xn, xd;
1985 : pari_sp av;
1986 :
1987 44009800 : if (typ(d) != t_INT) pari_err_TYPE("Q_muli_to_int",d);
1988 44009800 : switch (typ(x))
1989 : {
1990 : case t_INT:
1991 18224557 : return mulii(x,d);
1992 :
1993 : case t_FRAC:
1994 18980597 : xn = gel(x,1);
1995 18980597 : xd = gel(x,2); av = avma;
1996 18980597 : y = mulii(xn, diviiexact(d, xd));
1997 18980597 : return gerepileuptoint(av, y);
1998 : case t_COMPLEX:
1999 0 : y = cgetg(3,t_COMPLEX);
2000 0 : gel(y,1) = Q_muli_to_int(gel(x,1),d);
2001 0 : gel(y,2) = Q_muli_to_int(gel(x,2),d);
2002 0 : return y;
2003 :
2004 : case t_VEC: case t_COL: case t_MAT:
2005 4918135 : y = cgetg_copy(x, &l);
2006 4918135 : for (i=1; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
2007 4918135 : return y;
2008 :
2009 : case t_POL:
2010 1844840 : y = cgetg_copy(x, &l); y[1] = x[1];
2011 1844840 : for (i=2; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
2012 1844840 : return y;
2013 :
2014 : case t_POLMOD:
2015 41671 : y = cgetg(3, t_POLMOD);
2016 41671 : gel(y,1) = RgX_copy(gel(x,1));
2017 41671 : gel(y,2) = Q_muli_to_int(gel(x,2), d);
2018 41671 : return y;
2019 : }
2020 0 : pari_err_TYPE("Q_muli_to_int",x);
2021 : return NULL; /* LCOV_EXCL_LINE */
2022 : }
2023 :
2024 : static void
2025 1840356 : rescale_init(GEN c, int *exact, long *emin, GEN *D)
2026 : {
2027 : long e;
2028 1840356 : switch(typ(c))
2029 : {
2030 : case t_REAL:
2031 1232321 : *exact = 0;
2032 1232321 : if (!signe(c)) return;
2033 1229942 : e = expo(c) - bit_prec(c);
2034 1229942 : break;
2035 : case t_INT:
2036 607937 : if (!signe(c)) return;
2037 248674 : e = expi(c) + 32;
2038 248674 : break;
2039 : case t_FRAC:
2040 70 : e = expi(gel(c,1)) - expi(gel(c,2)) + 32;
2041 70 : if (exact) *D = lcmii(*D, gel(c,2));
2042 70 : break;
2043 : default:
2044 28 : pari_err_TYPE("rescale_to_int",c);
2045 : return; /* LCOV_EXCL_LINE */
2046 : }
2047 1478686 : if (e < *emin) *emin = e;
2048 : }
2049 : GEN
2050 223019 : RgM_rescale_to_int(GEN x)
2051 : {
2052 223019 : long lx = lg(x), i,j, hx, emin;
2053 : GEN D;
2054 : int exact;
2055 :
2056 223019 : if (lx == 1) return cgetg(1,t_MAT);
2057 223019 : hx = lgcols(x);
2058 223019 : exact = 1;
2059 223019 : emin = HIGHEXPOBIT;
2060 223019 : D = gen_1;
2061 808750 : for (j = 1; j < lx; j++)
2062 585759 : for (i = 1; i < hx; i++) rescale_init(gcoeff(x,i,j), &exact, &emin, &D);
2063 222991 : if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
2064 222935 : return grndtoi(gmul2n(x, -emin), &i);
2065 : }
2066 : GEN
2067 238 : RgX_rescale_to_int(GEN x)
2068 : {
2069 238 : long lx = lg(x), i, emin;
2070 : GEN D;
2071 : int exact;
2072 238 : if (lx == 2) return gcopy(x); /* rare */
2073 238 : exact = 1;
2074 238 : emin = HIGHEXPOBIT;
2075 238 : D = gen_1;
2076 238 : for (i = 2; i < lx; i++) rescale_init(gel(x,i), &exact, &emin, &D);
2077 238 : if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
2078 238 : return grndtoi(gmul2n(x, -emin), &i);
2079 : }
2080 :
2081 : /* return x * n/d. x: rational; d,n,result: integral; d,n coprime */
2082 : static GEN
2083 8259014 : Q_divmuli_to_int(GEN x, GEN d, GEN n)
2084 : {
2085 : long i, l;
2086 : GEN y, xn, xd;
2087 : pari_sp av;
2088 :
2089 8259014 : switch(typ(x))
2090 : {
2091 : case t_INT:
2092 2493594 : av = avma; y = diviiexact(x,d);
2093 2493594 : return gerepileuptoint(av, mulii(y,n));
2094 :
2095 : case t_FRAC:
2096 3974819 : xn = gel(x,1);
2097 3974819 : xd = gel(x,2); av = avma;
2098 3974819 : y = mulii(diviiexact(xn, d), diviiexact(n, xd));
2099 3974819 : return gerepileuptoint(av, y);
2100 :
2101 : case t_VEC: case t_COL: case t_MAT:
2102 427444 : y = cgetg_copy(x, &l);
2103 427444 : for (i=1; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
2104 427444 : return y;
2105 :
2106 : case t_POL:
2107 1363157 : y = cgetg_copy(x, &l); y[1] = x[1];
2108 1363157 : for (i=2; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
2109 1363157 : return y;
2110 :
2111 : case t_POLMOD:
2112 0 : y = cgetg(3, t_POLMOD);
2113 0 : gel(y,1) = RgX_copy(gel(x,1));
2114 0 : gel(y,2) = Q_divmuli_to_int(gel(x,2), d,n);
2115 0 : return y;
2116 : }
2117 0 : pari_err_TYPE("Q_divmuli_to_int",x);
2118 : return NULL; /* LCOV_EXCL_LINE */
2119 : }
2120 :
2121 : /* return x / d. x: rational; d,result: integral. */
2122 : static GEN
2123 16975303 : Q_divi_to_int(GEN x, GEN d)
2124 : {
2125 : long i, l;
2126 : GEN y;
2127 :
2128 16975303 : switch(typ(x))
2129 : {
2130 : case t_INT:
2131 13835118 : return diviiexact(x,d);
2132 :
2133 : case t_VEC: case t_COL: case t_MAT:
2134 1308474 : y = cgetg_copy(x, &l);
2135 1308474 : for (i=1; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
2136 1308474 : return y;
2137 :
2138 : case t_POL:
2139 1828085 : y = cgetg_copy(x, &l); y[1] = x[1];
2140 1828085 : for (i=2; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
2141 1828085 : return y;
2142 :
2143 : case t_POLMOD:
2144 3626 : y = cgetg(3, t_POLMOD);
2145 3626 : gel(y,1) = RgX_copy(gel(x,1));
2146 3626 : gel(y,2) = Q_divi_to_int(gel(x,2), d);
2147 3626 : return y;
2148 : }
2149 0 : pari_err_TYPE("Q_divi_to_int",x);
2150 : return NULL; /* LCOV_EXCL_LINE */
2151 : }
2152 : /* c t_FRAC */
2153 : static GEN
2154 3025546 : Q_divq_to_int(GEN x, GEN c)
2155 : {
2156 3025546 : GEN n = gel(c,1), d = gel(c,2);
2157 3025546 : if (is_pm1(n)) {
2158 1608470 : GEN y = Q_muli_to_int(x,d);
2159 1608470 : if (signe(n) < 0) y = gneg(y);
2160 1608470 : return y;
2161 : }
2162 1417076 : return Q_divmuli_to_int(x, n,d);
2163 : }
2164 :
2165 : /* return y = x / c, assuming x,c rational, and y integral */
2166 : GEN
2167 7486 : Q_div_to_int(GEN x, GEN c)
2168 : {
2169 7486 : switch(typ(c))
2170 : {
2171 6275 : case t_INT: return Q_divi_to_int(x, c);
2172 1211 : case t_FRAC: return Q_divq_to_int(x, c);
2173 : }
2174 0 : pari_err_TYPE("Q_div_to_int",c);
2175 : return NULL; /* LCOV_EXCL_LINE */
2176 : }
2177 : /* return y = x * c, assuming x,c rational, and y integral */
2178 : GEN
2179 0 : Q_mul_to_int(GEN x, GEN c)
2180 : {
2181 : GEN d, n;
2182 0 : switch(typ(c))
2183 : {
2184 0 : case t_INT: return Q_muli_to_int(x, c);
2185 : case t_FRAC:
2186 0 : n = gel(c,1);
2187 0 : d = gel(c,2);
2188 0 : return Q_divmuli_to_int(x, d,n);
2189 : }
2190 0 : pari_err_TYPE("Q_mul_to_int",c);
2191 : return NULL; /* LCOV_EXCL_LINE */
2192 : }
2193 :
2194 : GEN
2195 14507483 : Q_primitive_part(GEN x, GEN *ptc)
2196 : {
2197 14507483 : pari_sp av = avma;
2198 14507483 : GEN c = Q_content_safe(x);
2199 14507483 : if (c)
2200 : {
2201 14507343 : if (typ(c) == t_INT)
2202 : {
2203 11483008 : if (is_pm1(c)) { avma = av; c = NULL; }
2204 2148312 : else if (signe(c)) x = Q_divi_to_int(x, c);
2205 : }
2206 3024335 : else x = Q_divq_to_int(x, c);
2207 : }
2208 14507483 : if (ptc) *ptc = c;
2209 14507483 : return x;
2210 : }
2211 : GEN
2212 1445084 : Q_primpart(GEN x) { return Q_primitive_part(x, NULL); }
2213 :
2214 : GEN
2215 68618 : vec_Q_primpart(GEN x)
2216 68618 : { pari_APPLY_same(Q_primpart(gel(x,i))) }
2217 :
2218 : /*******************************************************************/
2219 : /* */
2220 : /* SUBRESULTANT */
2221 : /* */
2222 : /*******************************************************************/
2223 : /* for internal use */
2224 : GEN
2225 1791241 : gdivexact(GEN x, GEN y)
2226 : {
2227 : long i,lx;
2228 : GEN z;
2229 1791241 : if (gequal1(y)) return x;
2230 1787384 : switch(typ(x))
2231 : {
2232 : case t_INT:
2233 1431206 : if (typ(y)==t_INT) return diviiexact(x,y);
2234 35 : if (!signe(x)) return gen_0;
2235 0 : break;
2236 : case t_INTMOD:
2237 : case t_FFELT:
2238 13696 : case t_POLMOD: return gmul(x,ginv(y));
2239 : case t_POL:
2240 338506 : switch(typ(y))
2241 : {
2242 : case t_INTMOD:
2243 : case t_FFELT:
2244 1344 : case t_POLMOD: return gmul(x,ginv(y));
2245 : case t_POL: { /* not stack-clean */
2246 : long v;
2247 20209 : if (varn(x)!=varn(y)) break;
2248 19243 : v = RgX_valrem(y,&y);
2249 19243 : if (v) x = RgX_shift_shallow(x,-v);
2250 19243 : if (!degpol(y)) { y = gel(y,2); break; }
2251 17164 : return RgX_div(x,y);
2252 : }
2253 : }
2254 319998 : return RgX_Rg_divexact(x, y);
2255 : case t_VEC: case t_COL: case t_MAT:
2256 3976 : lx = lg(x); z = new_chunk(lx);
2257 3976 : for (i=1; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
2258 3976 : z[0] = x[0]; return z;
2259 : }
2260 0 : if (DEBUGLEVEL) pari_warn(warner,"missing case in gdivexact");
2261 0 : return gdiv(x,y);
2262 : }
2263 :
2264 : static GEN
2265 59856 : init_resultant(GEN x, GEN y)
2266 : {
2267 59856 : long tx = typ(x), ty = typ(y), vx, vy;
2268 59856 : if (is_scalar_t(tx) || is_scalar_t(ty))
2269 : {
2270 210 : if (gequal0(x) || gequal0(y)) return gmul(x,y); /* keep type info */
2271 210 : if (tx==t_POL) return gpowgs(y, degpol(x));
2272 0 : if (ty==t_POL) return gpowgs(x, degpol(y));
2273 0 : return gen_1;
2274 : }
2275 59646 : if (tx!=t_POL) pari_err_TYPE("resultant_all",x);
2276 59646 : if (ty!=t_POL) pari_err_TYPE("resultant_all",y);
2277 59646 : if (!signe(x) || !signe(y)) return gmul(Rg_get_0(x),Rg_get_0(y)); /*type*/
2278 59569 : vx = varn(x);
2279 59569 : vy = varn(y); if (vx == vy) return NULL;
2280 0 : return (varncmp(vx,vy) < 0)? gpowgs(y,degpol(x)): gpowgs(x,degpol(y));
2281 : }
2282 :
2283 : static long
2284 96955 : RgX_simpletype(GEN x)
2285 : {
2286 96955 : long T = t_INT, i, lx = lg(x);
2287 574891 : for (i = 2; i < lx; i++)
2288 : {
2289 486609 : GEN c = gel(x,i);
2290 486609 : long tc = typ(c);
2291 486609 : switch(tc) {
2292 : case t_INT:
2293 415501 : break;
2294 : case t_FRAC:
2295 32693 : if (T == t_INT) T = t_FRAC;
2296 32693 : break;
2297 : default:
2298 38415 : if (isinexact(c)) return t_REAL;
2299 29742 : T = 0; break;
2300 : }
2301 : }
2302 88282 : return T;
2303 : }
2304 :
2305 : /* x an RgX, y a scalar */
2306 : static GEN
2307 0 : scalar_res(GEN x, GEN y, GEN *U, GEN *V)
2308 : {
2309 0 : *V = gpowgs(y,degpol(x)-1);
2310 0 : *U = gen_0; return gmul(y, *V);
2311 : }
2312 :
2313 : /* return 0 if the subresultant chain can be interrupted.
2314 : * Set u = NULL if the resultant is 0. */
2315 : static int
2316 5172 : subres_step(GEN *u, GEN *v, GEN *g, GEN *h, GEN *uze, GEN *um1, long *signh)
2317 : {
2318 5172 : GEN u0, c, r, q = RgX_pseudodivrem(*u,*v, &r);
2319 : long du, dv, dr, degq;
2320 :
2321 5172 : if (gequal0(leading_coeff(r))) r = RgX_renormalize(r);
2322 5172 : dr = lg(r); if (!signe(r)) { *u = NULL; return 0; }
2323 5046 : du = degpol(*u);
2324 5046 : dv = degpol(*v);
2325 5046 : degq = du - dv;
2326 5046 : if (*um1 == gen_1)
2327 1854 : u0 = gpowgs(gel(*v,dv+2),degq+1);
2328 3192 : else if (*um1 == gen_0)
2329 1596 : u0 = gen_0;
2330 : else /* except in those 2 cases, um1 is an RgX */
2331 1596 : u0 = RgX_Rg_mul(*um1, gpowgs(gel(*v,dv+2),degq+1));
2332 :
2333 5046 : if (*uze == gen_0) /* except in that case, uze is an RgX */
2334 1854 : u0 = scalarpol(u0, varn(*u)); /* now an RgX */
2335 : else
2336 3192 : u0 = gsub(u0, gmul(q,*uze));
2337 :
2338 5046 : *um1 = *uze;
2339 5046 : *uze = u0; /* uze <- lead(v)^(degq + 1) * um1 - q * uze */
2340 :
2341 5046 : *u = *v; c = *g; *g = leading_coeff(*u);
2342 5046 : switch(degq)
2343 : {
2344 21 : case 0: break;
2345 : case 1:
2346 3919 : c = gmul(*h,c); *h = *g; break;
2347 : default:
2348 1106 : c = gmul(gpowgs(*h,degq), c);
2349 1106 : *h = gdivexact(gpowgs(*g,degq), gpowgs(*h,degq-1));
2350 : }
2351 5046 : *v = RgX_Rg_divexact(r,c);
2352 5046 : *uze= RgX_Rg_divexact(*uze,c);
2353 5046 : if (both_odd(du, dv)) *signh = -*signh;
2354 5046 : return (dr > 3);
2355 : }
2356 :
2357 : /* compute U, V s.t Ux + Vy = resultant(x,y) */
2358 : static GEN
2359 1749 : subresext_i(GEN x, GEN y, GEN *U, GEN *V)
2360 : {
2361 : pari_sp av, av2;
2362 1749 : long dx, dy, du, signh, tx = typ(x), ty = typ(y);
2363 : GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze;
2364 :
2365 1749 : if (!is_extscalar_t(tx)) pari_err_TYPE("subresext",x);
2366 1749 : if (!is_extscalar_t(ty)) pari_err_TYPE("subresext",y);
2367 1749 : if (gequal0(x) || gequal0(y)) { *U = *V = gen_0; return gen_0; }
2368 1749 : if (tx != t_POL) {
2369 0 : if (ty != t_POL) { *U = ginv(x); *V = gen_0; return gen_1; }
2370 0 : return scalar_res(y,x,V,U);
2371 : }
2372 1749 : if (ty != t_POL) return scalar_res(x,y,U,V);
2373 1749 : if (varn(x) != varn(y))
2374 0 : return varncmp(varn(x), varn(y)) < 0? scalar_res(x,y,U,V)
2375 0 : : scalar_res(y,x,V,U);
2376 1749 : if (gequal0(leading_coeff(x))) x = RgX_renormalize(x);
2377 1749 : if (gequal0(leading_coeff(y))) y = RgX_renormalize(y);
2378 1749 : dx = degpol(x);
2379 1749 : dy = degpol(y);
2380 1749 : signh = 1;
2381 1749 : if (dx < dy)
2382 : {
2383 874 : pswap(U,V); lswap(dx,dy); swap(x,y);
2384 874 : if (both_odd(dx, dy)) signh = -signh;
2385 : }
2386 1749 : if (dy == 0)
2387 : {
2388 0 : *V = gpowgs(gel(y,2),dx-1);
2389 0 : *U = gen_0; return gmul(*V,gel(y,2));
2390 : }
2391 1749 : av = avma;
2392 1749 : u = x = primitive_part(x, &cu);
2393 1749 : v = y = primitive_part(y, &cv);
2394 1749 : g = h = gen_1; av2 = avma;
2395 1749 : um1 = gen_1; uze = gen_0;
2396 : for(;;)
2397 : {
2398 4703 : if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
2399 2954 : if (gc_needed(av2,1))
2400 : {
2401 0 : if(DEBUGMEM>1) pari_warn(warnmem,"subresext, dr = %ld", degpol(v));
2402 0 : gerepileall(av2,6, &u,&v,&g,&h,&uze,&um1);
2403 : }
2404 2954 : }
2405 : /* uze an RgX */
2406 1749 : if (!u) { *U = *V = gen_0; avma = av; return gen_0; }
2407 1742 : z = gel(v,2); du = degpol(u);
2408 1742 : if (du > 1)
2409 : { /* z = gdivexact(gpowgs(z,du), gpowgs(h,du-1)); */
2410 167 : p1 = gpowgs(gdiv(z,h),du-1);
2411 167 : z = gmul(z,p1);
2412 167 : uze = RgX_Rg_mul(uze, p1);
2413 : }
2414 1742 : if (signh < 0) { z = gneg_i(z); uze = RgX_neg(uze); }
2415 :
2416 1742 : vze = RgX_divrem(Rg_RgX_sub(z, RgX_mul(uze,x)), y, &r);
2417 1742 : if (signe(r)) pari_warn(warner,"inexact computation in subresext");
2418 : /* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */
2419 1742 : p1 = gen_1;
2420 1742 : if (cu) p1 = gmul(p1, gpowgs(cu,dy));
2421 1742 : if (cv) p1 = gmul(p1, gpowgs(cv,dx));
2422 1742 : cu = cu? gdiv(p1,cu): p1;
2423 1742 : cv = cv? gdiv(p1,cv): p1;
2424 1742 : z = gmul(z,p1);
2425 1742 : *U = RgX_Rg_mul(uze,cu);
2426 1742 : *V = RgX_Rg_mul(vze,cv);
2427 1742 : return z;
2428 : }
2429 : GEN
2430 0 : subresext(GEN x, GEN y, GEN *U, GEN *V)
2431 : {
2432 0 : pari_sp av = avma;
2433 0 : GEN z = subresext_i(x, y, U, V);
2434 0 : gerepileall(av, 3, &z, U, V);
2435 0 : return z;
2436 : }
2437 :
2438 : static GEN
2439 28 : zero_extgcd(GEN y, GEN *U, GEN *V, long vx)
2440 : {
2441 28 : GEN x=content(y);
2442 28 : *U=pol_0(vx); *V = scalarpol(ginv(x), vx); return gmul(y,*V);
2443 : }
2444 :
2445 : static int
2446 90188 : must_negate(GEN x)
2447 : {
2448 90188 : GEN t = leading_coeff(x);
2449 90188 : switch(typ(t))
2450 : {
2451 : case t_INT: case t_REAL:
2452 57274 : return (signe(t) < 0);
2453 : case t_FRAC:
2454 364 : return (signe(gel(t,1)) < 0);
2455 : }
2456 32550 : return 0;
2457 : }
2458 :
2459 : /* compute U, V s.t Ux + Vy = GCD(x,y) using subresultant */
2460 : GEN
2461 343 : RgX_extgcd(GEN x, GEN y, GEN *U, GEN *V)
2462 : {
2463 : pari_sp av, av2, tetpil;
2464 : long signh; /* junk */
2465 343 : long dx, dy, vx, tx = typ(x), ty = typ(y);
2466 : GEN z, g, h, p1, cu, cv, u, v, um1, uze, vze, *gptr[3];
2467 :
2468 343 : if (tx!=t_POL) pari_err_TYPE("RgX_extgcd",x);
2469 343 : if (ty!=t_POL) pari_err_TYPE("RgX_extgcd",y);
2470 343 : if ( varncmp(varn(x),varn(y))) pari_err_VAR("RgX_extgcd",x,y);
2471 343 : vx=varn(x);
2472 343 : if (!signe(x))
2473 : {
2474 14 : if (signe(y)) return zero_extgcd(y,U,V,vx);
2475 7 : *U = pol_0(vx); *V = pol_0(vx);
2476 7 : return pol_0(vx);
2477 : }
2478 329 : if (!signe(y)) return zero_extgcd(x,V,U,vx);
2479 308 : dx = degpol(x); dy = degpol(y);
2480 308 : if (dx < dy) { pswap(U,V); lswap(dx,dy); swap(x,y); }
2481 308 : if (dy==0) { *U=pol_0(vx); *V=ginv(y); return pol_1(vx); }
2482 :
2483 161 : av = avma;
2484 161 : u = x = primitive_part(x, &cu);
2485 161 : v = y = primitive_part(y, &cv);
2486 161 : g = h = gen_1; av2 = avma;
2487 161 : um1 = gen_1; uze = gen_0;
2488 : for(;;)
2489 : {
2490 182 : if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
2491 21 : if (gc_needed(av2,1))
2492 : {
2493 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_extgcd, dr = %ld",degpol(v));
2494 0 : gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
2495 : }
2496 21 : }
2497 161 : if (uze != gen_0) {
2498 : GEN r;
2499 42 : vze = RgX_divrem(RgX_sub(v, RgX_mul(uze,x)), y, &r);
2500 42 : if (signe(r)) pari_warn(warner,"inexact computation in RgX_extgcd");
2501 42 : if (cu) uze = RgX_Rg_div(uze,cu);
2502 42 : if (cv) vze = RgX_Rg_div(vze,cv);
2503 42 : p1 = ginv(content(v));
2504 : }
2505 : else /* y | x */
2506 : {
2507 119 : vze = cv ? RgX_Rg_div(pol_1(vx),cv): pol_1(vx);
2508 119 : uze = pol_0(vx);
2509 119 : p1 = gen_1;
2510 : }
2511 161 : if (must_negate(v)) p1 = gneg(p1);
2512 161 : tetpil = avma;
2513 161 : z = RgX_Rg_mul(v,p1);
2514 161 : *U = RgX_Rg_mul(uze,p1);
2515 161 : *V = RgX_Rg_mul(vze,p1);
2516 161 : gptr[0] = &z;
2517 161 : gptr[1] = U;
2518 161 : gptr[2] = V;
2519 161 : gerepilemanysp(av,tetpil,gptr,3); return z;
2520 : }
2521 :
2522 : int
2523 70 : RgXQ_ratlift(GEN x, GEN T, long amax, long bmax, GEN *P, GEN *Q)
2524 : {
2525 70 : pari_sp av = avma, av2, tetpil;
2526 : long signh; /* junk */
2527 : long vx;
2528 : GEN g, h, p1, cu, cv, u, v, um1, uze, *gptr[2];
2529 :
2530 70 : if (typ(x)!=t_POL) pari_err_TYPE("RgXQ_ratlift",x);
2531 70 : if (typ(T)!=t_POL) pari_err_TYPE("RgXQ_ratlift",T);
2532 70 : if ( varncmp(varn(x),varn(T)) ) pari_err_VAR("RgXQ_ratlift",x,T);
2533 70 : if (bmax < 0) pari_err_DOMAIN("ratlift", "bmax", "<", gen_0, stoi(bmax));
2534 70 : if (!signe(T)) {
2535 0 : if (degpol(x) <= amax) {
2536 0 : *P = RgX_copy(x);
2537 0 : *Q = pol_1(varn(x));
2538 0 : return 1;
2539 : }
2540 0 : return 0;
2541 : }
2542 70 : if (amax+bmax >= degpol(T))
2543 0 : pari_err_DOMAIN("ratlift", "amax+bmax", ">=", stoi(degpol(T)),
2544 : mkvec3(stoi(amax), stoi(bmax), T));
2545 70 : vx = varn(T);
2546 70 : u = x = primitive_part(x, &cu);
2547 70 : v = T = primitive_part(T, &cv);
2548 70 : g = h = gen_1; av2 = avma;
2549 70 : um1 = gen_1; uze = gen_0;
2550 : for(;;)
2551 : {
2552 287 : (void) subres_step(&u, &v, &g, &h, &uze, &um1, &signh);
2553 287 : if (!u || (typ(uze)==t_POL && degpol(uze)>bmax)) { avma=av; return 0; }
2554 287 : if (typ(v)!=t_POL || degpol(v)<=amax) break;
2555 217 : if (gc_needed(av2,1))
2556 : {
2557 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXQ_ratlift, dr = %ld", degpol(v));
2558 0 : gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
2559 : }
2560 217 : }
2561 70 : if (uze == gen_0)
2562 : {
2563 0 : avma = av; *P = pol_0(vx); *Q = pol_1(vx);
2564 0 : return 1;
2565 : }
2566 70 : if (cu) uze = RgX_Rg_div(uze,cu);
2567 70 : p1 = ginv(content(v));
2568 70 : if (must_negate(v)) p1 = gneg(p1);
2569 70 : tetpil = avma;
2570 70 : *P = RgX_Rg_mul(v,p1);
2571 70 : *Q = RgX_Rg_mul(uze,p1);
2572 70 : gptr[0] = P;
2573 70 : gptr[1] = Q;
2574 70 : gerepilemanysp(av,tetpil,gptr,2); return 1;
2575 : }
2576 :
2577 : /*******************************************************************/
2578 : /* */
2579 : /* RESULTANT USING DUCOS VARIANT */
2580 : /* */
2581 : /*******************************************************************/
2582 : /* x^n / y^(n-1), assume n > 0 */
2583 : static GEN
2584 19473 : Lazard(GEN x, GEN y, long n)
2585 : {
2586 : long a;
2587 : GEN c;
2588 :
2589 19473 : if (n == 1) return x;
2590 1397 : a = 1 << expu(n); /* a = 2^k <= n < 2^(k+1) */
2591 1397 : c=x; n-=a;
2592 4814 : while (a>1)
2593 : {
2594 2020 : a>>=1; c=gdivexact(gsqr(c),y);
2595 2020 : if (n>=a) { c=gdivexact(gmul(c,x),y); n -= a; }
2596 : }
2597 1397 : return c;
2598 : }
2599 :
2600 : /* F (x/y)^(n-1), assume n >= 1 */
2601 : static GEN
2602 19761 : Lazard2(GEN F, GEN x, GEN y, long n)
2603 : {
2604 19761 : if (n == 1) return F;
2605 735 : return RgX_Rg_divexact(RgX_Rg_mul(F, Lazard(x,y,n-1)), y);
2606 : }
2607 :
2608 : static GEN
2609 19761 : RgX_neg_i(GEN x, long lx)
2610 : {
2611 : long i;
2612 19761 : GEN y = cgetg(lx, t_POL); y[1] = x[1];
2613 19761 : for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
2614 19761 : return y;
2615 : }
2616 : static GEN
2617 59472 : RgX_Rg_mul_i(GEN y, GEN x, long ly)
2618 : {
2619 : long i;
2620 : GEN z;
2621 59472 : if (isrationalzero(x)) return pol_0(varn(y));
2622 59451 : z = cgetg(ly,t_POL); z[1] = y[1];
2623 59451 : for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
2624 59451 : return z;
2625 : }
2626 : static long
2627 56560 : reductum_lg(GEN x, long lx)
2628 : {
2629 56560 : long i = lx-2;
2630 56560 : while (i > 1 && gequal0(gel(x,i))) i--;
2631 56560 : return i+1;
2632 : }
2633 :
2634 : #define addshift(x,y) RgX_addmulXn_shallow((x),(y),1)
2635 : /* delta = deg(P) - deg(Q) > 0, deg(Q) > 0, P,Q,Z t_POL in the same variable,
2636 : * s "scalar". Return prem(P, -Q) / s^delta lc(P) */
2637 : static GEN
2638 19761 : nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
2639 : {
2640 19761 : GEN p0, q0, h0, TMP, H, A, z0 = leading_coeff(Z);
2641 : long p, q, j, lP, lQ;
2642 : pari_sp av;
2643 :
2644 19761 : p = degpol(P); p0 = gel(P,p+2); lP = reductum_lg(P,lg(P));
2645 19761 : q = degpol(Q); q0 = gel(Q,q+2); lQ = reductum_lg(Q,lg(Q));
2646 : /* p > q. Very often p - 1 = q */
2647 19761 : av = avma;
2648 : /* H = RgX_neg(reductum(Z)) optimized, using Q ~ Z */
2649 19761 : H = RgX_neg_i(Z, lQ); /* deg H < q */
2650 :
2651 19761 : A = (q+2 < lP)? RgX_Rg_mul_i(H, gel(P,q+2), lQ): NULL;
2652 23401 : for (j = q+1; j < p; j++)
2653 : {
2654 3640 : if (degpol(H) == q-1)
2655 : { /* h0 = coeff of degree q-1 = leading coeff */
2656 3206 : h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1);
2657 3206 : H = addshift(H, RgX_Rg_divexact(RgX_Rg_mul_i(Q, gneg(h0), lQ), q0));
2658 : }
2659 : else
2660 434 : H = RgX_shift_shallow(H, 1);
2661 3640 : if (j+2 < lP)
2662 : {
2663 3066 : TMP = RgX_Rg_mul(H, gel(P,j+2));
2664 3066 : A = A? RgX_add(A, TMP): TMP;
2665 : }
2666 3640 : if (gc_needed(av,1))
2667 : {
2668 147 : if(DEBUGMEM>1) pari_warn(warnmem,"nextSousResultant j = %ld/%ld",j,p);
2669 147 : gerepileall(av,A?2:1,&H,&A);
2670 : }
2671 : }
2672 19761 : if (q+2 < lP) lP = reductum_lg(P, q+3);
2673 19761 : TMP = RgX_Rg_mul_i(P, z0, lP);
2674 19761 : A = A? RgX_add(A, TMP): TMP;
2675 19761 : A = RgX_Rg_divexact(A, p0);
2676 19761 : if (degpol(H) == q-1)
2677 : {
2678 19467 : h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1); /* destroy old H */
2679 19467 : A = RgX_add(RgX_Rg_mul(addshift(H,A),q0), RgX_Rg_mul_i(Q, gneg(h0), lQ));
2680 : }
2681 : else
2682 294 : A = RgX_Rg_mul(addshift(H,A), q0);
2683 19761 : return RgX_Rg_divexact(A, s);
2684 : }
2685 : #undef addshift
2686 :
2687 : /* Ducos's subresultant */
2688 : GEN
2689 20432 : RgX_resultant_all(GEN P, GEN Q, GEN *sol)
2690 : {
2691 : pari_sp av, av2;
2692 20432 : long dP, dQ, delta, sig = 1;
2693 : GEN cP, cQ, Z, s;
2694 :
2695 20432 : dP = degpol(P);
2696 20432 : dQ = degpol(Q); delta = dP - dQ;
2697 20432 : if (delta < 0)
2698 : {
2699 1463 : if (both_odd(dP, dQ)) sig = -1;
2700 1463 : swap(P,Q); lswap(dP, dQ); delta = -delta;
2701 : }
2702 20432 : if (sol) *sol = gen_0;
2703 20432 : av = avma;
2704 20432 : if (dQ <= 0)
2705 : {
2706 1694 : if (dQ < 0) return Rg_get_0(P);
2707 1694 : s = gpowgs(gel(Q,2), dP);
2708 1694 : if (sig == -1) s = gerepileupto(av, gneg(s));
2709 1694 : return s;
2710 : }
2711 : /* primitive_part is also possible here, but possibly very costly,
2712 : * and hardly ever worth it */
2713 18738 : P = Q_primitive_part(P, &cP);
2714 18738 : Q = Q_primitive_part(Q, &cQ);
2715 18738 : av2 = avma;
2716 18738 : s = gpowgs(leading_coeff(Q),delta);
2717 18738 : if (both_odd(dP, dQ)) sig = -sig;
2718 18738 : Z = Q;
2719 18738 : Q = RgX_pseudorem(P, Q);
2720 18738 : P = Z;
2721 57237 : while(degpol(Q) > 0)
2722 : {
2723 19761 : delta = degpol(P) - degpol(Q); /* > 0 */
2724 19761 : Z = Lazard2(Q, leading_coeff(Q), s, delta);
2725 19761 : if (both_odd(degpol(P), degpol(Q))) sig = -sig;
2726 19761 : Q = nextSousResultant(P, Q, Z, s);
2727 19761 : P = Z;
2728 19761 : if (gc_needed(av,1))
2729 : {
2730 13 : if(DEBUGMEM>1) pari_warn(warnmem,"resultant_all, degpol Q = %ld",degpol(Q));
2731 13 : gerepileall(av2,2,&P,&Q);
2732 : }
2733 19761 : s = leading_coeff(P);
2734 : }
2735 18738 : if (!signe(Q)) { avma = av; return Rg_get_0(Q); }
2736 18738 : s = Lazard(leading_coeff(Q), s, degpol(P));
2737 18738 : if (sig == -1) s = gneg(s);
2738 18738 : if (cP) s = gmul(s, gpowgs(cP,dQ));
2739 18738 : if (cQ) s = gmul(s, gpowgs(cQ,dP));
2740 18738 : if (sol) { *sol = P; gerepileall(av, 2, &s, sol); return s; }
2741 16995 : return gerepilecopy(av, s);
2742 : }
2743 : /* Return resultant(P,Q).
2744 : * Uses Sylvester's matrix if P or Q inexact, a modular algorithm if they
2745 : * are in Q[X], and Ducos/Lazard optimization of the subresultant algorithm
2746 : * in the "generic" case. */
2747 : GEN
2748 51176 : resultant(GEN P, GEN Q)
2749 : {
2750 : long TP, TQ;
2751 51176 : GEN s, p = NULL;
2752 :
2753 51176 : if ((s = init_resultant(P,Q))) return s;
2754 50889 : if ((TP = RgX_simpletype(P)) == t_REAL || (TQ = RgX_simpletype(Q)) == t_REAL)
2755 8659 : return resultant2(P,Q); /* inexact */
2756 42230 : if (TP && TQ) /* rational */
2757 : {
2758 24682 : if (TP == t_INT && TQ == t_INT) return ZX_resultant(P,Q);
2759 10441 : return QX_resultant(P,Q);
2760 : }
2761 17548 : if (RgX_is_FpX(P, &p) && RgX_is_FpX(Q, &p) && p)
2762 : {
2763 14 : pari_sp av = avma;
2764 14 : GEN r = FpX_resultant(RgX_to_FpX(P, p), RgX_to_FpX(Q, p), p);
2765 14 : return gerepileupto(av, Fp_to_mod(r, p));
2766 : }
2767 17534 : return RgX_resultant_all(P, Q, NULL);
2768 : }
2769 :
2770 : /*******************************************************************/
2771 : /* */
2772 : /* RESULTANT USING SYLVESTER MATRIX */
2773 : /* */
2774 : /*******************************************************************/
2775 : static GEN
2776 72632 : syl_RgC(GEN x, long j, long d, long D, long cp)
2777 : {
2778 72632 : GEN c = cgetg(d+1,t_COL);
2779 : long i;
2780 72632 : for (i=1; i< j; i++) gel(c,i) = gen_0;
2781 72632 : for ( ; i<=D; i++) { GEN t = gel(x,D-i+2); gel(c,i) = cp? gcopy(t): t; }
2782 72632 : for ( ; i<=d; i++) gel(c,i) = gen_0;
2783 72632 : return c;
2784 : }
2785 : static GEN
2786 8687 : syl_RgM(GEN x, GEN y, long cp)
2787 : {
2788 8687 : long j, d, dx = degpol(x), dy = degpol(y);
2789 : GEN M;
2790 8687 : if (dx < 0) return dy < 0? cgetg(1,t_MAT): zeromat(dy,dy);
2791 8687 : if (dy < 0) return zeromat(dx,dx);
2792 8687 : d = dx+dy; M = cgetg(d+1,t_MAT);
2793 8687 : for (j=1; j<=dy; j++) gel(M,j) = syl_RgC(x,j,d,j+dx, cp);
2794 8687 : for (j=1; j<=dx; j++) gel(M,j+dy) = syl_RgC(y,j,d,j+dy, cp);
2795 8687 : return M;
2796 : }
2797 : GEN
2798 8680 : RgX_sylvestermatrix(GEN x, GEN y) { return syl_RgM(x,y,0); }
2799 : GEN
2800 7 : sylvestermatrix(GEN x, GEN y)
2801 : {
2802 7 : if (typ(x)!=t_POL) pari_err_TYPE("sylvestermatrix",x);
2803 7 : if (typ(y)!=t_POL) pari_err_TYPE("sylvestermatrix",y);
2804 7 : if (varn(x) != varn(y)) pari_err_VAR("sylvestermatrix",x,y);
2805 7 : return syl_RgM(x,y,1);
2806 : }
2807 :
2808 : GEN
2809 8680 : resultant2(GEN x, GEN y)
2810 : {
2811 8680 : pari_sp av = avma;
2812 8680 : GEN r = init_resultant(x,y);
2813 8680 : return r? r: gerepileupto(av, det(RgX_sylvestermatrix(x,y)));
2814 : }
2815 :
2816 : /* If x a t_POL, let vx = main variable of x; return a t_POL in variable v0:
2817 : * if vx <= v, return subst(x, v, pol_x(v0))
2818 : * if vx > v, return scalarpol(x, v0) */
2819 : static GEN
2820 3458 : fix_pol(GEN x, long v, long v0)
2821 : {
2822 : long vx;
2823 3458 : if (typ(x) != t_POL) return x;
2824 3458 : vx = varn(x);
2825 3458 : if (v == vx)
2826 : {
2827 3311 : if (v0 != v) { x = leafcopy(x); setvarn(x, v0); }
2828 3311 : return x;
2829 : }
2830 147 : if (varncmp(v, vx) > 0)
2831 : {
2832 140 : x = gsubst(x,v,pol_x(v0));
2833 140 : if (typ(x) == t_POL && varn(x) == v0) return x;
2834 : }
2835 35 : return scalarpol_shallow(x, v0);
2836 : }
2837 :
2838 : /* resultant of x and y with respect to variable v, or with respect to their
2839 : * main variable if v < 0. */
2840 : GEN
2841 1946 : polresultant0(GEN x, GEN y, long v, long flag)
2842 : {
2843 1946 : long v0 = 0;
2844 1946 : pari_sp av = avma;
2845 :
2846 1946 : if (v >= 0)
2847 : {
2848 1715 : v0 = fetch_var_higher();
2849 1715 : x = fix_pol(x,v, v0);
2850 1715 : y = fix_pol(y,v, v0);
2851 : }
2852 1946 : switch(flag)
2853 : {
2854 : case 2:
2855 1939 : case 0: x=resultant(x,y); break;
2856 7 : case 1: x=resultant2(x,y); break;
2857 0 : default: pari_err_FLAG("polresultant");
2858 : }
2859 1946 : if (v >= 0) (void)delete_var();
2860 1946 : return gerepileupto(av,x);
2861 : }
2862 : GEN
2863 875 : polresultantext0(GEN x, GEN y, long v)
2864 : {
2865 : GEN R, U, V;
2866 875 : long v0 = 0;
2867 875 : pari_sp av = avma;
2868 :
2869 875 : if (v >= 0)
2870 : {
2871 14 : v0 = fetch_var_higher();
2872 14 : x = fix_pol(x,v, v0);
2873 14 : y = fix_pol(y,v, v0);
2874 : }
2875 875 : R = subresext_i(x,y, &U,&V);
2876 875 : if (v >= 0)
2877 : {
2878 14 : (void)delete_var();
2879 14 : if (typ(U) == t_POL && varn(U) != v) U = poleval(U, pol_x(v));
2880 14 : if (typ(V) == t_POL && varn(V) != v) V = poleval(V, pol_x(v));
2881 : }
2882 875 : return gerepilecopy(av, mkvec3(U,V,R));
2883 : }
2884 : GEN
2885 847 : polresultantext(GEN x, GEN y) { return polresultantext0(x,y,-1); }
2886 :
2887 : /*******************************************************************/
2888 : /* */
2889 : /* CHARACTERISTIC POLYNOMIAL USING RESULTANT */
2890 : /* */
2891 : /*******************************************************************/
2892 :
2893 : /* (v - x)^d */
2894 : static GEN
2895 140 : caract_const(pari_sp av, GEN x, long v, long d)
2896 140 : { return gerepileupto(av, gpowgs(deg1pol_shallow(gen_1, gneg_i(x), v), d)); }
2897 :
2898 : /* return caract(Mod(x,T)) in variable v */
2899 : GEN
2900 14461 : RgXQ_charpoly(GEN x, GEN T, long v)
2901 : {
2902 14461 : pari_sp av = avma;
2903 14461 : long d = degpol(T), dx, vx, vp, v0;
2904 : GEN ch, L;
2905 :
2906 14461 : if (typ(x) != t_POL) return caract_const(av, x, v, d);
2907 14447 : vx = varn(x);
2908 14447 : vp = varn(T);
2909 14447 : if (varncmp(vx, vp) > 0) return caract_const(av, x, v, d);
2910 14447 : if (varncmp(vx, vp) < 0) pari_err_PRIORITY("RgXQ_charpoly", x, "<", vp);
2911 14447 : dx = degpol(x);
2912 14447 : if (dx >= degpol(T)) { x = RgX_rem(x, T); dx = degpol(x); }
2913 14447 : if (dx <= 0) return dx? pol_xn(d, v): caract_const(av, gel(x,2), v, d);
2914 :
2915 14391 : v0 = fetch_var_higher();
2916 14391 : x = RgX_neg(x);
2917 14391 : gel(x,2) = gadd(gel(x,2), pol_x(v));
2918 14391 : setvarn(x, v0);
2919 14391 : T = leafcopy(T); setvarn(T, v0);
2920 14391 : ch = resultant(T, x);
2921 14391 : (void)delete_var();
2922 : /* test for silly input: x mod (deg 0 polynomial) */
2923 14391 : if (typ(ch) != t_POL)
2924 7 : pari_err_PRIORITY("RgXQ_charpoly", pol_x(v), "<", gvar(ch));
2925 14384 : L = leading_coeff(ch);
2926 14384 : if (!gequal1(L)) ch = RgX_Rg_div(ch, L);
2927 14384 : return gerepileupto(av, ch);
2928 : }
2929 :
2930 : /* characteristic polynomial (in v) of x over nf, where x is an element of the
2931 : * algebra nf[t]/(Q(t)) */
2932 : GEN
2933 224 : rnfcharpoly(GEN nf, GEN Q, GEN x, long v)
2934 : {
2935 224 : const char *f = "rnfcharpoly";
2936 224 : long dQ = degpol(Q);
2937 224 : pari_sp av = avma;
2938 : GEN T;
2939 :
2940 224 : if (v < 0) v = 0;
2941 224 : nf = checknf(nf); T = nf_get_pol(nf);
2942 224 : Q = RgX_nffix(f, T,Q,0);
2943 224 : switch(typ(x))
2944 : {
2945 : case t_INT:
2946 28 : case t_FRAC: return caract_const(av, x, v, dQ);
2947 : case t_POLMOD:
2948 91 : x = polmod_nffix2(f,T,Q, x,0);
2949 49 : break;
2950 : case t_POL:
2951 56 : x = varn(x) == varn(T)? Rg_nffix(f,T,x,0): RgX_nffix(f, T,x,0);
2952 42 : break;
2953 49 : default: pari_err_TYPE(f,x);
2954 : }
2955 91 : if (typ(x) != t_POL) return caract_const(av, x, v, dQ);
2956 : /* x a t_POL in variable vQ */
2957 49 : if (degpol(x) >= dQ) x = RgX_rem(x, Q);
2958 49 : if (dQ <= 1) return caract_const(av, constant_coeff(x), v, 1);
2959 49 : return gerepilecopy(av, lift_if_rational( RgXQ_charpoly(x, Q, v) ));
2960 : }
2961 :
2962 : /*******************************************************************/
2963 : /* */
2964 : /* GCD USING SUBRESULTANT */
2965 : /* */
2966 : /*******************************************************************/
2967 : static int inexact(GEN x, int *simple, int *rational);
2968 : static int
2969 6944596 : isinexactall(GEN x, int *simple, int *rational)
2970 : {
2971 6944596 : long i, lx = lg(x);
2972 33458343 : for (i=2; i<lx; i++)
2973 26513761 : if (inexact(gel(x,i), simple, rational)) return 1;
2974 6944582 : return 0;
2975 : }
2976 : /* return 1 if coeff explosion is not possible */
2977 : static int
2978 26513971 : inexact(GEN x, int *simple, int *rational)
2979 : {
2980 26513971 : int junk = 0;
2981 26513971 : switch(typ(x))
2982 : {
2983 25088855 : case t_INT: case t_FRAC: return 0;
2984 :
2985 7 : case t_REAL: case t_PADIC: case t_SER: return 1;
2986 :
2987 : case t_INTMOD:
2988 : case t_FFELT:
2989 1383501 : *rational = 0;
2990 1383501 : if (!*simple) *simple = 1;
2991 1383501 : return 0;
2992 :
2993 : case t_COMPLEX:
2994 77 : *rational = 0;
2995 154 : return inexact(gel(x,1), simple, rational)
2996 77 : || inexact(gel(x,2), simple, rational);
2997 : case t_QUAD:
2998 0 : *rational = *simple = 0;
2999 0 : return inexact(gel(x,2), &junk, rational)
3000 0 : || inexact(gel(x,3), &junk, rational);
3001 :
3002 : case t_POLMOD:
3003 11718 : *rational = 0;
3004 11718 : return isinexactall(gel(x,1), simple, rational);
3005 : case t_POL:
3006 29785 : *rational = 0;
3007 29785 : *simple = -1;
3008 29785 : return isinexactall(x, &junk, rational);
3009 : case t_RFRAC:
3010 28 : *rational = 0;
3011 28 : *simple = -1;
3012 56 : return inexact(gel(x,1), &junk, rational)
3013 28 : || inexact(gel(x,2), &junk, rational);
3014 : }
3015 0 : *rational = 0;
3016 0 : *simple = -1; return 0;
3017 : }
3018 :
3019 : /* x monomial, y t_POL in the same variable */
3020 : static GEN
3021 7446453 : gcdmonome(GEN x, GEN y)
3022 : {
3023 7446453 : pari_sp av = avma;
3024 7446453 : long dx = degpol(x), e = RgX_valrem(y, &y);
3025 7446453 : long i, l = lg(y);
3026 7446453 : GEN t, v = cgetg(l, t_VEC);
3027 7446453 : gel(v,1) = gel(x,dx+2);
3028 7446453 : for (i = 2; i < l; i++) gel(v,i) = gel(y,i);
3029 7446453 : t = content(v); /* gcd(lc(x), cont(y)) */
3030 7446453 : t = simplify_shallow(t);
3031 7446453 : if (dx < e) e = dx;
3032 7446453 : return gerepileupto(av, monomialcopy(t, e, varn(x)));
3033 : }
3034 :
3035 : /* x, y are t_POL in the same variable */
3036 : GEN
3037 10898388 : RgX_gcd(GEN x, GEN y)
3038 : {
3039 : long dx, dy;
3040 : pari_sp av, av1;
3041 : GEN d, g, h, p1, p2, u, v;
3042 10898388 : int simple = 0, rational = 1;
3043 :
3044 10898388 : if (isexactzero(y)) return RgX_copy(x);
3045 10898010 : if (isexactzero(x)) return RgX_copy(y);
3046 10898003 : if (RgX_is_monomial(x)) return gcdmonome(x,y);
3047 4099077 : if (RgX_is_monomial(y)) return gcdmonome(y,x);
3048 3451550 : if (isinexactall(x,&simple,&rational) || isinexactall(y,&simple,&rational))
3049 : {
3050 7 : av = avma; u = ggcd(content(x), content(y));
3051 7 : return gerepileupto(av, scalarpol(u, varn(x)));
3052 : }
3053 3451543 : if (rational) return QX_gcd(x,y); /* Q[X] */
3054 :
3055 92099 : av = avma;
3056 92099 : if (simple > 0) x = RgX_gcd_simple(x,y);
3057 : else
3058 : {
3059 4935 : dx = lg(x); dy = lg(y);
3060 4935 : if (dx < dy) { swap(x,y); lswap(dx,dy); }
3061 4935 : if (dy==3)
3062 : {
3063 0 : d = ggcd(gel(y,2), content(x));
3064 0 : return gerepileupto(av, scalarpol(d, varn(x)));
3065 : }
3066 4935 : u = primitive_part(x, &p1); if (!p1) p1 = gen_1;
3067 4935 : v = primitive_part(y, &p2); if (!p2) p2 = gen_1;
3068 4935 : d = ggcd(p1,p2);
3069 4935 : av1 = avma;
3070 4935 : g = h = gen_1;
3071 : for(;;)
3072 : {
3073 6377 : GEN r = RgX_pseudorem(u,v);
3074 6377 : long degq, du, dv, dr = lg(r);
3075 :
3076 6377 : if (!signe(r)) break;
3077 3584 : if (dr <= 3)
3078 : {
3079 2142 : avma = av1; return gerepileupto(av, scalarpol(d, varn(x)));
3080 : }
3081 1442 : if (DEBUGLEVEL > 9) err_printf("RgX_gcd: dr = %ld\n", degpol(r));
3082 1442 : du = lg(u); dv = lg(v); degq = du-dv;
3083 1442 : u = v; p1 = g; g = leading_coeff(u);
3084 1442 : switch(degq)
3085 : {
3086 189 : case 0: break;
3087 : case 1:
3088 1029 : p1 = gmul(h,p1); h = g; break;
3089 : default:
3090 224 : p1 = gmul(gpowgs(h,degq), p1);
3091 224 : h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
3092 : }
3093 1442 : v = RgX_Rg_div(r,p1);
3094 1442 : if (gc_needed(av1,1))
3095 : {
3096 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd");
3097 0 : gerepileall(av1,4, &u,&v,&g,&h);
3098 : }
3099 1442 : }
3100 2793 : x = RgX_Rg_mul(primpart(v), d);
3101 : }
3102 89957 : if (must_negate(x)) x = RgX_neg(x);
3103 89957 : return gerepileupto(av,x);
3104 : }
3105 :
3106 : /* disc P = (-1)^(n(n-1)/2) lc(P)^(n - deg P' - 2) Res(P,P'), n = deg P */
3107 : static GEN
3108 5145 : RgX_disc_aux(GEN P)
3109 : {
3110 5145 : long n = degpol(P), TP, dd;
3111 : GEN D, L, y, p;
3112 5145 : if (!signe(P) || !n) return Rg_get_0(P);
3113 5145 : if (n == 1) return Rg_get_1(P);
3114 4970 : if (n == 2) {
3115 1134 : GEN a = gel(P,4), b = gel(P,3), c = gel(P,2);
3116 1134 : return gsub(gsqr(b), gmul2n(gmul(a,c),2));
3117 : }
3118 3836 : TP = RgX_simpletype(P);
3119 3836 : if (TP == t_INT) return ZX_disc(P);
3120 483 : if (TP == t_FRAC) return QX_disc(P);
3121 483 : p = NULL;
3122 483 : if (RgX_is_FpX(P, &p) && p)
3123 28 : return Fp_to_mod(FpX_disc(RgX_to_FpX(P,p), p), p);
3124 :
3125 455 : y = RgX_deriv(P);
3126 455 : if (!signe(y)) return Rg_get_0(y);
3127 455 : dd = degpol(P)-2 - degpol(y);
3128 455 : if (TP == t_REAL)
3129 14 : D = resultant2(P,y);
3130 : else
3131 : {
3132 441 : D = RgX_resultant_all(P, y, NULL);
3133 441 : if (D == gen_0) return Rg_get_0(y);
3134 : }
3135 455 : L = leading_coeff(P);
3136 455 : if (dd && !gequal1(L)) D = (dd == -1)? gdiv(D, L): gmul(D, gpowgs(L, dd));
3137 455 : if (n & 2) D = gneg(D);
3138 455 : return D;
3139 : }
3140 : GEN
3141 1953 : RgX_disc(GEN x) { pari_sp av = avma; return gerepileupto(av, RgX_disc_aux(x)); }
3142 :
3143 : GEN
3144 3206 : poldisc0(GEN x, long v)
3145 : {
3146 : pari_sp av;
3147 3206 : switch(typ(x))
3148 : {
3149 : case t_POL:
3150 : {
3151 : GEN D;
3152 3192 : long v0 = -1;
3153 3192 : av = avma;
3154 3192 : if (v >= 0 && v != varn(x))
3155 : {
3156 0 : v0 = fetch_var_higher();
3157 0 : x = fix_pol(x,v, v0);
3158 : }
3159 3192 : D = RgX_disc_aux(x);
3160 3192 : if (v0 >= 0) (void)delete_var();
3161 3192 : return gerepileupto(av, D);
3162 : }
3163 :
3164 : case t_COMPLEX:
3165 0 : return utoineg(4);
3166 :
3167 : case t_QUAD:
3168 0 : return quad_disc(x);
3169 : case t_POLMOD:
3170 0 : return poldisc0(gel(x,1), v);
3171 :
3172 : case t_QFR: case t_QFI:
3173 14 : av = avma; return gerepileuptoint(av, qfb_disc(x));
3174 :
3175 : case t_VEC: case t_COL: case t_MAT:
3176 : {
3177 : long i;
3178 0 : GEN z = cgetg_copy(x, &i);
3179 0 : for (i--; i; i--) gel(z,i) = poldisc0(gel(x,i), v);
3180 0 : return z;
3181 : }
3182 : }
3183 0 : pari_err_TYPE("poldisc",x);
3184 : return NULL; /* LCOV_EXCL_LINE */
3185 : }
3186 :
3187 : GEN
3188 7 : reduceddiscsmith(GEN x)
3189 : {
3190 7 : long j, n = degpol(x);
3191 7 : pari_sp av = avma;
3192 : GEN xp, M;
3193 :
3194 7 : if (typ(x) != t_POL) pari_err_TYPE("poldiscreduced",x);
3195 7 : if (n<=0) pari_err_CONSTPOL("poldiscreduced");
3196 7 : RgX_check_ZX(x,"poldiscreduced");
3197 7 : if (!gequal1(gel(x,n+2)))
3198 0 : pari_err_IMPL("non-monic polynomial in poldiscreduced");
3199 7 : M = cgetg(n+1,t_MAT);
3200 7 : xp = ZX_deriv(x);
3201 28 : for (j=1; j<=n; j++)
3202 : {
3203 21 : gel(M,j) = RgX_to_RgC(xp, n);
3204 21 : if (j<n) xp = RgX_rem(RgX_shift_shallow(xp, 1), x);
3205 : }
3206 7 : return gerepileupto(av, ZM_snf(M));
3207 : }
3208 :
3209 : /***********************************************************************/
3210 : /** **/
3211 : /** STURM ALGORITHM **/
3212 : /** (number of real roots of x in [a,b]) **/
3213 : /** **/
3214 : /***********************************************************************/
3215 : static GEN
3216 833 : R_to_Q_up(GEN x)
3217 : {
3218 : long e;
3219 833 : switch(typ(x))
3220 : {
3221 833 : case t_INT: case t_FRAC: case t_INFINITY: return x;
3222 : case t_REAL:
3223 0 : x = mantissa_real(x,&e);
3224 0 : return gmul2n(addiu(x,1), -e);
3225 0 : default: pari_err_TYPE("R_to_Q_up", x);
3226 : return NULL; /* LCOV_EXCL_LINE */
3227 : }
3228 : }
3229 : static GEN
3230 833 : R_to_Q_down(GEN x)
3231 : {
3232 : long e;
3233 833 : switch(typ(x))
3234 : {
3235 819 : case t_INT: case t_FRAC: case t_INFINITY: return x;
3236 : case t_REAL:
3237 14 : x = mantissa_real(x,&e);
3238 14 : return gmul2n(subiu(x,1), -e);
3239 0 : default: pari_err_TYPE("R_to_Q_down", x);
3240 : return NULL; /* LCOV_EXCL_LINE */
3241 : }
3242 : }
3243 :
3244 : static long
3245 847 : sturmpart_i(GEN x, GEN ab)
3246 : {
3247 847 : long tx = typ(x);
3248 847 : if (gequal0(x)) pari_err_ROOTS0("sturm");
3249 847 : if (tx != t_POL)
3250 : {
3251 0 : if (is_real_t(tx)) return 0;
3252 0 : pari_err_TYPE("sturm",x);
3253 : }
3254 847 : if (lg(x) == 3) return 0;
3255 847 : if (!RgX_is_ZX(x)) x = RgX_rescale_to_int(x);
3256 847 : if (!ZX_is_squarefree(x))
3257 14 : pari_err_DOMAIN("polsturm","issquarefree(pol)","=",gen_0,x);
3258 833 : if (ab)
3259 : {
3260 : GEN A, B;
3261 833 : if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("RgX_sturmpart", ab);
3262 833 : A = R_to_Q_down(gel(ab,1));
3263 833 : B = R_to_Q_up(gel(ab,2));
3264 833 : ab = mkvec2(A, B);
3265 : }
3266 833 : return ZX_sturmpart(x, ab);
3267 : }
3268 : /* Deprecated: RgX_sturmpart() should be preferred */
3269 : long
3270 847 : sturmpart(GEN x, GEN a, GEN b)
3271 : {
3272 847 : pari_sp av = avma;
3273 : long r;
3274 847 : if (!b && a && typ(a) == t_VEC) return RgX_sturmpart(x, a);
3275 714 : if (!a) a = mkmoo();
3276 714 : if (!b) b = mkoo();
3277 714 : r = sturmpart_i(x, mkvec2(a, b));
3278 700 : avma = av; return r;
3279 : }
3280 : long
3281 133 : RgX_sturmpart(GEN x, GEN ab)
3282 : {
3283 133 : pari_sp av = avma;
3284 133 : long r = sturmpart_i(x, ab);
3285 133 : avma = av; return r;
3286 : }
3287 :
3288 : /***********************************************************************/
3289 : /** **/
3290 : /** GENERIC EXTENDED GCD **/
3291 : /** **/
3292 : /***********************************************************************/
3293 : /* assume typ(x) = typ(y) = t_POL */
3294 : static GEN
3295 874 : RgXQ_inv_i(GEN x, GEN y)
3296 : {
3297 874 : long vx=varn(x), vy=varn(y);
3298 : pari_sp av;
3299 : GEN u, v, d;
3300 :
3301 1748 : while (vx != vy)
3302 : {
3303 0 : if (varncmp(vx,vy) > 0)
3304 : {
3305 0 : d = (vx == NO_VARIABLE)? ginv(x): gred_rfrac_simple(gen_1, x);
3306 0 : return scalarpol(d, vy);
3307 : }
3308 0 : if (lg(x)!=3) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
3309 0 : x = gel(x,2); vx = gvar(x);
3310 : }
3311 874 : av = avma; d = subresext_i(x,y,&u,&v/*junk*/);
3312 874 : if (gequal0(d)) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
3313 874 : d = gdiv(u,d);
3314 874 : if (typ(d) != t_POL || varn(d) != vy) d = scalarpol(d, vy);
3315 874 : return gerepileupto(av, d);
3316 : }
3317 :
3318 : /*Assume x is a polynomial and y is not */
3319 : static GEN
3320 112 : scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
3321 : {
3322 112 : long vx = varn(x);
3323 112 : int xis0 = signe(x)==0, yis0 = gequal0(y);
3324 112 : if (xis0 && yis0) { *U = *V = pol_0(vx); return pol_0(vx); }
3325 84 : if (yis0) { *U=pol_1(vx); *V = pol_0(vx); return RgX_copy(x);}
3326 56 : *U=pol_0(vx); *V= ginv(y); return pol_1(vx);
3327 : }
3328 : /* Assume x==0, y!=0 */
3329 : static GEN
3330 63 : zero_bezout(GEN y, GEN *U, GEN *V)
3331 : {
3332 63 : *U=gen_0; *V = ginv(y); return gen_1;
3333 : }
3334 :
3335 : GEN
3336 280 : gbezout(GEN x, GEN y, GEN *u, GEN *v)
3337 : {
3338 280 : long tx=typ(x), ty=typ(y), vx;
3339 280 : if (tx == t_INT && ty == t_INT) return bezout(x,y,u,v);
3340 245 : if (tx != t_POL)
3341 : {
3342 140 : if (ty == t_POL)
3343 56 : return scalar_bezout(y,x,v,u);
3344 : else
3345 : {
3346 84 : int xis0 = gequal0(x), yis0 = gequal0(y);
3347 84 : if (xis0 && yis0) { *u = *v = gen_0; return gen_0; }
3348 63 : if (xis0) return zero_bezout(y,u,v);
3349 42 : else return zero_bezout(x,v,u);
3350 : }
3351 : }
3352 105 : else if (ty != t_POL) return scalar_bezout(x,y,u,v);
3353 49 : vx = varn(x);
3354 49 : if (vx != varn(y))
3355 0 : return varncmp(vx, varn(y)) < 0? scalar_bezout(x,y,u,v)
3356 0 : : scalar_bezout(y,x,v,u);
3357 49 : return RgX_extgcd(x,y,u,v);
3358 : }
3359 :
3360 : GEN
3361 280 : gcdext0(GEN x, GEN y)
3362 : {
3363 280 : GEN z=cgetg(4,t_VEC);
3364 280 : gel(z,3) = gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
3365 280 : return z;
3366 : }
3367 :
3368 : /*******************************************************************/
3369 : /* */
3370 : /* GENERIC (modular) INVERSE */
3371 : /* */
3372 : /*******************************************************************/
3373 :
3374 : GEN
3375 14904 : ginvmod(GEN x, GEN y)
3376 : {
3377 14904 : long tx=typ(x);
3378 :
3379 14904 : switch(typ(y))
3380 : {
3381 : case t_POL:
3382 14904 : if (tx==t_POL) return RgXQ_inv(x,y);
3383 2023 : if (is_scalar_t(tx)) return ginv(x);
3384 0 : break;
3385 : case t_INT:
3386 0 : if (tx==t_INT) return Fp_inv(x,y);
3387 0 : if (tx==t_POL) return gen_0;
3388 : }
3389 0 : pari_err_TYPE2("ginvmod",x,y);
3390 : return NULL; /* LCOV_EXCL_LINE */
3391 : }
3392 :
3393 : /***********************************************************************/
3394 : /** **/
3395 : /** NEWTON POLYGON **/
3396 : /** **/
3397 : /***********************************************************************/
3398 :
3399 : /* assume leading coeff of x is non-zero */
3400 : GEN
3401 28 : newtonpoly(GEN x, GEN p)
3402 : {
3403 : GEN y;
3404 : long n,ind,a,b,c,u1,u2,r1,r2;
3405 28 : long *vval, num[] = {evaltyp(t_INT) | _evallg(3), 0, 0};
3406 :
3407 28 : if (typ(x)!=t_POL) pari_err_TYPE("newtonpoly",x);
3408 28 : n=degpol(x); if (n<=0) return cgetg(1,t_VEC);
3409 28 : y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
3410 28 : vval = (long *) pari_malloc(sizeof(long)*(n+1));
3411 28 : for (a=0; a<=n; a++) vval[a] = gvaluation(gel(x,a),p);
3412 42 : for (a=0, ind=1; a<n; a++)
3413 : {
3414 42 : if (vval[a] != LONG_MAX) break;
3415 14 : gel(y,ind++) = mkoo();
3416 : }
3417 84 : for (b=a+1; b<=n; a=b, b=a+1)
3418 : {
3419 56 : while (vval[b] == LONG_MAX) b++;
3420 56 : u1 = vval[a]-vval[b];
3421 56 : u2 = b-a;
3422 154 : for (c=b+1; c<=n; c++)
3423 : {
3424 98 : if (vval[c] == LONG_MAX) continue;
3425 70 : r1 = vval[a]-vval[c];
3426 70 : r2 = c-a;
3427 70 : if (u1*r2 <= u2*r1) { u1 = r1; u2 = r2; b = c; }
3428 : }
3429 56 : while (ind<=b) { affsi(u1,num); gel(y,ind++) = gdivgs(num,u2); }
3430 : }
3431 28 : pari_free(vval); return y;
3432 : }
3433 :
3434 : static GEN
3435 262591 : RgXQ_mul_FpXQ(GEN x, GEN y, GEN T, GEN p)
3436 : {
3437 262591 : pari_sp av = avma;
3438 : GEN r;
3439 262591 : if (lgefint(p) == 3)
3440 : {
3441 144256 : ulong pp = uel(p, 2);
3442 144256 : r = Flx_to_ZX_inplace(Flxq_mul(RgX_to_Flx(x, pp),
3443 : RgX_to_Flx(y, pp), RgX_to_Flx(T, pp), pp));
3444 : }
3445 : else
3446 118335 : r = FpXQ_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), RgX_to_FpX(T, p), p);
3447 262591 : return gerepileupto(av, FpX_to_mod(r, p));
3448 : }
3449 :
3450 : static GEN
3451 14 : RgXQ_sqr_FpXQ(GEN x, GEN y, GEN p)
3452 : {
3453 14 : pari_sp av = avma;
3454 : GEN r;
3455 14 : if (lgefint(p) == 3)
3456 : {
3457 7 : ulong pp = uel(p, 2);
3458 7 : r = Flx_to_ZX_inplace(Flxq_sqr(RgX_to_Flx(x, pp),
3459 : RgX_to_Flx(y, pp), pp));
3460 : }
3461 : else
3462 7 : r = FpXQ_sqr(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3463 14 : return gerepileupto(av, FpX_to_mod(r, p));
3464 : }
3465 :
3466 : static GEN
3467 11571 : RgXQ_inv_FpXQ(GEN x, GEN y, GEN p)
3468 : {
3469 11571 : pari_sp av = avma;
3470 : GEN r;
3471 11571 : if (lgefint(p) == 3)
3472 : {
3473 5789 : ulong pp = uel(p, 2);
3474 5789 : r = Flx_to_ZX_inplace(Flxq_inv(RgX_to_Flx(x, pp),
3475 : RgX_to_Flx(y, pp), pp));
3476 : }
3477 : else
3478 5782 : r = FpXQ_inv(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3479 11571 : return gerepileupto(av, FpX_to_mod(r, p));
3480 : }
3481 :
3482 : static GEN
3483 7 : RgXQ_mul_FpXQXQ(GEN x, GEN y, GEN S, GEN pol, GEN p)
3484 : {
3485 7 : pari_sp av = avma;
3486 : GEN r;
3487 7 : GEN T = RgX_to_FpX(pol, p);
3488 7 : if (lgefint(p) == 3)
3489 : {
3490 7 : ulong pp = uel(p, 2);
3491 7 : GEN Tp = ZX_to_Flx(T, pp);
3492 7 : r = FlxX_to_ZXX(FlxqXQ_mul(RgX_to_FlxqX(x, Tp, pp),
3493 : RgX_to_FlxqX(y, Tp, pp),
3494 : RgX_to_FlxqX(S, Tp, pp), Tp, pp));
3495 : }
3496 : else
3497 0 : r = FpXQXQ_mul(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p),
3498 : RgX_to_FpXQX(S, T, p), T, p);
3499 7 : return gerepileupto(av, FpXQX_to_mod(r, T, p));
3500 : }
3501 :
3502 : static GEN
3503 0 : RgXQ_sqr_FpXQXQ(GEN x, GEN y, GEN pol, GEN p)
3504 : {
3505 0 : pari_sp av = avma;
3506 : GEN r;
3507 0 : GEN T = RgX_to_FpX(pol, p);
3508 0 : if (lgefint(p) == 3)
3509 : {
3510 0 : ulong pp = uel(p, 2);
3511 0 : GEN Tp = ZX_to_Flx(T, pp);
3512 0 : r = FlxX_to_ZXX(FlxqXQ_sqr(RgX_to_FlxqX(x, Tp, pp),
3513 : RgX_to_FlxqX(y, Tp, pp), Tp, pp));
3514 : }
3515 : else
3516 0 : r = FpXQXQ_sqr(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
3517 0 : return gerepileupto(av, FpXQX_to_mod(r, T, p));
3518 : }
3519 :
3520 : static GEN
3521 7 : RgXQ_inv_FpXQXQ(GEN x, GEN y, GEN pol, GEN p)
3522 : {
3523 7 : pari_sp av = avma;
3524 : GEN r;
3525 7 : GEN T = RgX_to_FpX(pol, p);
3526 7 : if (lgefint(p) == 3)
3527 : {
3528 7 : ulong pp = uel(p, 2);
3529 7 : GEN Tp = ZX_to_Flx(T, pp);
3530 7 : r = FlxX_to_ZXX(FlxqXQ_inv(RgX_to_FlxqX(x, Tp, pp),
3531 : RgX_to_FlxqX(y, Tp, pp), Tp, pp));
3532 : }
3533 : else
3534 0 : r = FpXQXQ_inv(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
3535 7 : return gerepileupto(av, FpXQX_to_mod(r, T, p));
3536 : }
3537 :
3538 : #define code(t1,t2) ((t1 << 6) | t2)
3539 : static GEN
3540 718387 : RgXQ_mul_fast(GEN x, GEN y, GEN T)
3541 : {
3542 : GEN p, pol;
3543 : long pa;
3544 718387 : long t = RgX_type3(x,y,T, &p,&pol,&pa);
3545 718389 : switch(t)
3546 : {
3547 402378 : case t_INT: return ZX_is_monic(T) ? ZXQ_mul(x,y,T): NULL;
3548 42530 : case t_FRAC: return RgX_is_ZX(T) && ZX_is_monic(T) ? QXQ_mul(x,y,T): NULL;
3549 14 : case t_FFELT: return FFXQ_mul(x, y, T, pol);
3550 262591 : case t_INTMOD: return RgXQ_mul_FpXQ(x, y, T, p);
3551 : case code(t_POLMOD, t_INTMOD):
3552 7 : return RgXQ_mul_FpXQXQ(x, y, T, pol, p);
3553 10869 : default: return NULL;
3554 : }
3555 : }
3556 :
3557 : GEN
3558 718383 : RgXQ_mul(GEN x, GEN y, GEN T)
3559 : {
3560 718383 : GEN z = RgXQ_mul_fast(x, y, T);
3561 718389 : if (!z) z = RgX_rem(RgX_mul(x, y), T);
3562 718389 : return z;
3563 : }
3564 :
3565 : static GEN
3566 277409 : RgXQ_sqr_fast(GEN x, GEN T)
3567 : {
3568 : GEN p, pol;
3569 : long pa;
3570 277409 : long t = RgX_type2(x, T, &p,&pol,&pa);
3571 277412 : switch(t)
3572 : {
3573 265379 : case t_INT: return ZX_is_monic(T) ? ZXQ_sqr(x,T): NULL;
3574 10843 : case t_FRAC: return RgX_is_ZX(T) && ZX_is_monic(T) ? QXQ_sqr(x,T): NULL;
3575 0 : case t_FFELT: return FFXQ_sqr(x, T, pol);
3576 14 : case t_INTMOD: return RgXQ_sqr_FpXQ(x, T, p);
3577 : case code(t_POLMOD, t_INTMOD):
3578 0 : return RgXQ_sqr_FpXQXQ(x, T, pol, p);
3579 1176 : default: return NULL;
3580 : }
3581 : }
3582 :
3583 : GEN
3584 277411 : RgXQ_sqr(GEN x, GEN T)
3585 : {
3586 277411 : GEN z = RgXQ_sqr_fast(x, T);
3587 277400 : if (!z) z = RgX_rem(RgX_sqr(x), T);
3588 277400 : return z;
3589 : }
3590 :
3591 : static GEN
3592 29269 : RgXQ_inv_fast(GEN x, GEN y)
3593 : {
3594 : GEN p, pol;
3595 : long pa;
3596 29269 : long t = RgX_type2(x,y, &p,&pol,&pa);
3597 29269 : switch(t)
3598 : {
3599 16061 : case t_INT: return QXQ_inv(x,y);
3600 749 : case t_FRAC: return RgX_is_ZX(y)? QXQ_inv(x,y): NULL;
3601 14 : case t_FFELT: return FFXQ_inv(x, y, pol);
3602 11571 : case t_INTMOD: return RgXQ_inv_FpXQ(x, y, p);
3603 : case code(t_POLMOD, t_INTMOD):
3604 7 : return RgXQ_inv_FpXQXQ(x, y, pol, p);
3605 867 : default: return NULL;
3606 : }
3607 : }
3608 : #undef code
3609 :
3610 : GEN
3611 29269 : RgXQ_inv(GEN x, GEN y)
3612 : {
3613 29269 : GEN z = RgXQ_inv_fast(x, y);
3614 29255 : if (!z) z = RgXQ_inv_i(x, y);
3615 29255 : return z;
3616 : }
|
