Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /***********************************************************************/
16 : /** **/
17 : /** ARITHMETIC OPERATIONS ON POLYNOMIALS **/
18 : /** (second part) **/
19 : /** **/
20 : /***********************************************************************/
21 : #include "pari.h"
22 : #include "paripriv.h"
23 :
24 : #define DEBUGLEVEL DEBUGLEVEL_pol
25 :
26 : /* compute Newton sums S_1(P), ... , S_n(P). S_k(P) = sum a_j^k, a_j root of P
27 : * If N != NULL, assume p-adic roots and compute mod N [assume integer coeffs]
28 : * If T != NULL, compute mod (T,N) [assume integer coeffs if N != NULL]
29 : * If y0!= NULL, precomputed i-th powers, i=1..m, m = length(y0).
30 : * Not memory clean in the latter case */
31 : GEN
32 111943 : polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N)
33 : {
34 111943 : long dP = degpol(P), i, k, m;
35 : GEN y, P_lead;
36 :
37 111943 : if (n<0) pari_err_IMPL("polsym of a negative n");
38 111943 : if (typ(P) != t_POL) pari_err_TYPE("polsym",P);
39 111943 : if (!signe(P)) pari_err_ROOTS0("polsym");
40 111943 : y = cgetg(n+2,t_COL);
41 111943 : if (y0)
42 : {
43 11205 : if (typ(y0) != t_COL) pari_err_TYPE("polsym_gen",y0);
44 11205 : m = lg(y0)-1;
45 54289 : for (i=1; i<=m; i++) gel(y,i) = gel(y0,i); /* not memory clean */
46 : }
47 : else
48 : {
49 100738 : m = 1;
50 100738 : gel(y,1) = stoi(dP);
51 : }
52 111943 : P += 2; /* strip codewords */
53 :
54 111943 : P_lead = gel(P,dP); if (gequal1(P_lead)) P_lead = NULL;
55 111943 : if (P_lead)
56 : {
57 5 : if (N) P_lead = Fq_inv(P_lead,T,N);
58 5 : else if (T) P_lead = QXQ_inv(P_lead,T);
59 : }
60 332203 : for (k=m; k<=n; k++)
61 : {
62 220260 : pari_sp av1 = avma;
63 220260 : GEN s = (dP>=k)? gmulsg(k,gel(P,dP-k)): gen_0;
64 635726 : for (i=1; i<k && i<=dP; i++)
65 415466 : s = gadd(s, gmul(gel(y,k-i+1),gel(P,dP-i)));
66 220260 : if (N)
67 : {
68 15866 : s = Fq_red(s, T, N);
69 15866 : if (P_lead) s = Fq_mul(s, P_lead, T, N);
70 : }
71 204394 : else if (T)
72 : {
73 0 : s = grem(s, T);
74 0 : if (P_lead) s = grem(gmul(s, P_lead), T);
75 : }
76 : else
77 204394 : if (P_lead) s = gdiv(s, P_lead);
78 220260 : gel(y,k+1) = gc_upto(av1, gneg(s));
79 : }
80 111943 : return y;
81 : }
82 :
83 : GEN
84 97229 : polsym(GEN x, long n)
85 : {
86 97229 : return polsym_gen(x, NULL, n, NULL,NULL);
87 : }
88 :
89 : /* centered residue x mod p. po2 = shifti(p, -1) or NULL (euclidean residue) */
90 : GEN
91 59520710 : centermodii(GEN x, GEN p, GEN po2)
92 : {
93 59520710 : GEN y = remii(x, p);
94 59520710 : switch(signe(y))
95 : {
96 9622449 : case 0: break;
97 37929700 : case 1: if (po2 && abscmpii(y,po2) > 0) y = subii(y, p);
98 37929700 : break;
99 11968561 : case -1: if (!po2 || abscmpii(y,po2) > 0) y = addii(y, p);
100 11968561 : break;
101 : }
102 59520710 : return y;
103 : }
104 :
105 : static long
106 0 : s_centermod(long x, ulong pp, ulong pps2)
107 : {
108 0 : long y = x % (long)pp;
109 0 : if (y < 0) y += pp;
110 0 : return Fl_center(y, pp,pps2);
111 : }
112 :
113 : /* for internal use */
114 : GEN
115 11120833 : centermod_i(GEN x, GEN p, GEN ps2)
116 : {
117 : long i, lx;
118 : pari_sp av;
119 : GEN y;
120 :
121 11120833 : if (!ps2) ps2 = shifti(p,-1);
122 11120833 : switch(typ(x))
123 : {
124 1212470 : case t_INT: return centermodii(x,p,ps2);
125 :
126 3776009 : case t_POL: lx = lg(x);
127 3776009 : y = cgetg(lx,t_POL); y[1] = x[1];
128 21666289 : for (i=2; i<lx; i++)
129 : {
130 17890280 : av = avma;
131 17890280 : gel(y,i) = gc_INT(av, centermodii(gel(x,i),p,ps2));
132 : }
133 3776009 : return normalizepol_lg(y, lx);
134 :
135 25921100 : case t_COL: pari_APPLY_same(centermodii(gel(x,i),p,ps2));
136 11522 : case t_MAT: pari_APPLY_same(centermod_i(gel(x,i),p,ps2));
137 :
138 0 : case t_VECSMALL: lx = lg(x);
139 : {
140 0 : ulong pp = itou(p), pps2 = itou(ps2);
141 0 : pari_APPLY_long(s_centermod(x[i], pp, pps2));
142 : }
143 : }
144 0 : return x;
145 : }
146 :
147 : GEN
148 9523710 : centermod(GEN x, GEN p) { return centermod_i(x,p,NULL); }
149 :
150 : static GEN
151 282 : RgX_Frobenius_deflate(GEN S, ulong p)
152 : {
153 282 : if (degpol(S)%p)
154 0 : return NULL;
155 : else
156 : {
157 282 : GEN F = RgX_deflate(S, p);
158 282 : long i, l = lg(F);
159 870 : for (i=2; i<l; i++)
160 : {
161 606 : GEN Fi = gel(F,i), R;
162 606 : if (typ(Fi)==t_POL)
163 : {
164 210 : if (signe(RgX_deriv(Fi))==0)
165 192 : gel(F,i) = RgX_Frobenius_deflate(gel(F, i), p);
166 18 : else return NULL;
167 : }
168 396 : else if (ispower(Fi, utoi(p), &R))
169 396 : gel(F,i) = R;
170 0 : else return NULL;
171 : }
172 264 : return F;
173 : }
174 : }
175 :
176 : static GEN
177 228 : RgXY_squff(GEN f)
178 : {
179 228 : long i, q, n = degpol(f);
180 228 : ulong p = itos_or_0(characteristic(f));
181 228 : GEN u = const_vec(n+1, pol_1(varn(f)));
182 228 : for(q = 1;;q *= p)
183 72 : {
184 300 : GEN t, v, tv, r = RgX_gcd(f, RgX_deriv(f));
185 300 : if (degpol(r) == 0) { gel(u, q) = f; break; }
186 108 : t = RgX_div(f, r);
187 108 : if (degpol(t) > 0)
188 : {
189 : long j;
190 24 : for(j = 1;;j++)
191 : {
192 120 : v = RgX_gcd(r, t);
193 120 : tv = RgX_div(t, v);
194 120 : if (degpol(tv) > 0) gel(u, j*q) = tv;
195 120 : if (degpol(v) <= 0) break;
196 96 : r = RgX_div(r, v);
197 96 : t = v;
198 : }
199 24 : if (degpol(r) == 0) break;
200 : }
201 90 : if (!p) break;
202 90 : f = RgX_Frobenius_deflate(r, p);
203 90 : if (!f) { gel(u, q) = r; break; }
204 : }
205 822 : for (i = n; i; i--)
206 810 : if (degpol(gel(u,i))) break;
207 228 : setlg(u,i+1); return u;
208 : }
209 :
210 : /* Lmod contains modular factors of *F (NULL codes an empty slot: used factor)
211 : * Lfac accumulates irreducible factors as they are found.
212 : * p is a product of modular factors in Lmod[1..i-1] (NULL for p = 1), not
213 : * a rational factor of *F
214 : * Find an irreducible factor of *F divisible by p (by including
215 : * exhaustively further factors from Lmod[i..]); return 0 on failure, else 1.
216 : * Update Lmod, Lfac and *F */
217 : static int
218 5694 : RgX_cmbf(GEN p, long i, GEN BLOC, GEN Lmod, GEN Lfac, GEN *F)
219 : {
220 : pari_sp av;
221 : GEN q;
222 5694 : if (i == lg(Lmod)) return 0;
223 2982 : if (RgX_cmbf(p, i+1, BLOC, Lmod, Lfac, F) && p) return 1;
224 2814 : if (!gel(Lmod,i)) return 0;
225 2772 : p = p? RgX_mul(p, gel(Lmod,i)): gel(Lmod,i);
226 2772 : av = avma;
227 2772 : q = RgV_to_RgX(RgX_digits(p, BLOC), varn(*F));
228 2772 : if (degpol(q))
229 : {
230 : GEN R, Q;
231 2454 : q = simplify_shallow(q);
232 2454 : Q = RgX_divrem(*F, q, &R);
233 2454 : if (signe(R)==0) { vectrunc_append(Lfac, q); *F = Q; return 1; }
234 : }
235 2484 : set_avma(av);
236 2484 : if (RgX_cmbf(p, i+1, BLOC, Lmod, Lfac, F)) { gel(Lmod,i) = NULL; return 1; }
237 2250 : return 0;
238 : }
239 :
240 : static GEN factor_domain(GEN x, GEN flag);
241 :
242 : static GEN
243 390 : ok_bloc(GEN f, GEN BLOC, ulong c)
244 : {
245 390 : GEN F = poleval(f, BLOC);
246 390 : return issquarefree(c ? gmul(F,mkintmodu(1,c)): F)? F: NULL;
247 : }
248 : static GEN
249 102 : random_FpX_monic(long n, long v, GEN p)
250 : {
251 102 : long i, d = n + 2;
252 102 : GEN y = cgetg(d + 1, t_POL); y[1] = evalsigne(1) | evalvarn(v);
253 336 : for (i = 2; i < d; i++) gel(y,i) = randomi(p);
254 102 : gel(y,i) = gen_1; return y;
255 : }
256 : static GEN
257 240 : RgXY_factor_squarefree(GEN f, GEN dom)
258 : {
259 240 : pari_sp av = avma;
260 240 : ulong i, c = itou_or_0(residual_characteristic(f));
261 240 : long vy = gvar2(f), val = RgX_valrem(f, &f), n = RgXY_degreex(f);
262 240 : GEN y, Lmod, F = NULL, BLOC = NULL, Lfac = coltrunc_init(degpol(f)+2);
263 240 : GEN gc = c? utoipos(c): NULL;
264 240 : if (val)
265 : {
266 30 : GEN x = pol_x(varn(f));
267 30 : if (dom)
268 : {
269 12 : GEN one = Rg_get_1(dom);
270 12 : if (typ(one) != t_INT) x = RgX_Rg_mul(x, one);
271 : }
272 30 : vectrunc_append(Lfac, x); if (!degpol(f)) return Lfac;
273 : }
274 228 : y = pol_x(vy);
275 : for(;;)
276 : {
277 288 : for (i = 0; !c || i < c; i++)
278 : {
279 288 : BLOC = gpowgs(gaddgs(y, i), n+1);
280 288 : if ((F = ok_bloc(f, BLOC, c))) break;
281 132 : if (c)
282 : {
283 102 : BLOC = random_FpX_monic(n, vy, gc);
284 102 : if ((F = ok_bloc(f, BLOC, c))) break;
285 : }
286 : }
287 228 : if (!c || i < c) break;
288 0 : n++;
289 : }
290 228 : if (DEBUGLEVEL >= 2)
291 0 : err_printf("bifactor: bloc:(x+%ld)^%ld, deg f=%ld\n",i,n,RgXY_degreex(f));
292 228 : Lmod = gel(factor_domain(F,dom),1);
293 228 : if (DEBUGLEVEL >= 2)
294 0 : err_printf("bifactor: %ld local factors\n",lg(Lmod)-1);
295 228 : (void)RgX_cmbf(NULL, 1, BLOC, Lmod, Lfac, &f);
296 228 : if (degpol(f)) vectrunc_append(Lfac, f);
297 228 : return gc_GEN(av, Lfac);
298 : }
299 :
300 : static GEN
301 228 : FE_matconcat(GEN F, GEN E, long l)
302 : {
303 228 : setlg(E,l); E = shallowconcat1(E);
304 228 : setlg(F,l); F = shallowconcat1(F); return mkmat2(F,E);
305 : }
306 :
307 : static int
308 342 : gen_cmp_RgXY(void *data, GEN x, GEN y)
309 : {
310 342 : long vx = varn(x), vy = varn(y);
311 342 : return (vx == vy)? gen_cmp_RgX(data, x, y): -varncmp(vx, vy);
312 : }
313 : static GEN
314 228 : RgXY_factor(GEN f, GEN dom)
315 : {
316 228 : pari_sp av = avma;
317 : GEN C, F, E, cf, V;
318 : long i, j, l;
319 228 : if (dom) { GEN c = Rg_get_1(dom); if (typ(c) != t_INT) f = RgX_Rg_mul(f,c); }
320 228 : cf = content(f);
321 228 : V = RgXY_squff(gdiv(f, simplify_shallow(cf))); l = lg(V);
322 228 : C = factor_domain(cf, dom);
323 228 : F = cgetg(l+1, t_VEC); gel(F,1) = gel(C,1);
324 228 : E = cgetg(l+1, t_VEC); gel(E,1) = gel(C,2);
325 672 : for (i=1, j=2; i < l; i++)
326 : {
327 444 : GEN v = gel(V,i);
328 444 : if (degpol(v))
329 : {
330 240 : gel(F,j) = v = RgXY_factor_squarefree(v, dom);
331 240 : gel(E,j) = const_col(lg(v)-1, utoipos(i));
332 240 : j++;
333 : }
334 : }
335 228 : f = FE_matconcat(F,E,j);
336 228 : (void)sort_factor(f,(void*)cmp_universal, &gen_cmp_RgXY);
337 228 : return gc_GEN(av, f);
338 : }
339 :
340 : /***********************************************************************/
341 : /** **/
342 : /** FACTORIZATION **/
343 : /** **/
344 : /***********************************************************************/
345 : static long RgX_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var);
346 : #define assign_or_fail(x,y) { GEN __x = x;\
347 : if (!*y) *y=__x; else if (!gequal(__x,*y)) return 0;\
348 : }
349 : #define update_prec(x,y) { long __x = x; if (__x < *y) *y=__x; }
350 :
351 : static const long RgX_type_shift = 6;
352 : void
353 9573697 : RgX_type_decode(long x, long *t1, long *t2)
354 : {
355 9573697 : *t1 = x >> RgX_type_shift;
356 9573697 : *t2 = (x & ((1L<<RgX_type_shift)-1));
357 9573697 : }
358 : int
359 138012340 : RgX_type_is_composite(long t) { return t >= RgX_type_shift; }
360 :
361 : static int
362 3219517138 : settype(GEN c, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
363 : {
364 : long j;
365 3219517138 : switch(typ(c))
366 : {
367 2420228217 : case t_INT:
368 2420228217 : break;
369 28526207 : case t_FRAC:
370 28526207 : t[1]=1; break;
371 : break;
372 259144439 : case t_REAL:
373 259144439 : update_prec(precision(c), pa);
374 259144439 : t[2]=1; break;
375 25616763 : case t_INTMOD:
376 25616763 : assign_or_fail(gel(c,1),p);
377 25616763 : t[3]=1; break;
378 1574626 : case t_FFELT:
379 1574626 : if (!*ff) *ff=c; else if (!FF_samefield(c,*ff)) return 0;
380 1574626 : assign_or_fail(FF_p_i(c),p);
381 1574626 : t[5]=1; break;
382 278066242 : case t_COMPLEX:
383 834198714 : for (j=1; j<=2; j++)
384 : {
385 556132478 : GEN d = gel(c,j);
386 556132478 : switch(typ(d))
387 : {
388 2710348 : case t_INT: case t_FRAC:
389 2710348 : if (!*t2) *t2 = t_COMPLEX;
390 2710348 : t[1]=1; break;
391 553422102 : case t_REAL:
392 553422102 : update_prec(precision(d), pa);
393 553422102 : if (!*t2) *t2 = t_COMPLEX;
394 553422102 : t[2]=1; break;
395 12 : case t_INTMOD:
396 12 : assign_or_fail(gel(d,1),p);
397 12 : if (!signe(*p) || mod4(*p) != 3) return 0;
398 6 : if (!*t2) *t2 = t_COMPLEX;
399 6 : t[3]=1; break;
400 16 : case t_PADIC:
401 16 : update_prec(precp(d)+valp(d), pa);
402 16 : assign_or_fail(padic_p(d), p);
403 16 : if (!*t2) *t2 = t_COMPLEX;
404 16 : t[7]=1; break;
405 0 : default: return 0;
406 : }
407 : }
408 278066236 : if (!t[2]) assign_or_fail(mkpoln(3, gen_1,gen_0,gen_1), pol); /*x^2+1*/
409 278066230 : break;
410 2001647 : case t_PADIC:
411 2001647 : update_prec(precp(c)+valp(c), pa);
412 2001647 : assign_or_fail(padic_p(c), p);
413 2001647 : t[7]=1; break;
414 1722 : case t_QUAD:
415 1722 : assign_or_fail(gel(c,1),pol);
416 5166 : for (j=2; j<=3; j++)
417 : {
418 3444 : GEN d = gel(c,j);
419 3444 : switch(typ(d))
420 : {
421 3414 : case t_INT: case t_FRAC:
422 3414 : t[8]=1; break;
423 24 : case t_INTMOD:
424 24 : assign_or_fail(gel(d,1),p);
425 24 : if (*t2 != t_POLMOD) *t2 = t_QUAD;
426 24 : t[3]=1; break;
427 6 : case t_PADIC:
428 6 : update_prec(precp(d)+valp(d), pa);
429 6 : assign_or_fail(padic_p(d), p);
430 6 : if (*t2 != t_POLMOD) *t2 = t_QUAD;
431 6 : t[7]=1; break;
432 0 : default: return 0;
433 : }
434 : }
435 1722 : break;
436 3443821 : case t_POLMOD:
437 3443821 : assign_or_fail(gel(c,1),pol);
438 3443647 : if (typ(gel(c,2))==t_POL && varn(gel(c,2))!=varn(gel(c,1))) return 0;
439 10325279 : for (j=1; j<=2; j++)
440 : {
441 : GEN pbis, polbis;
442 : long pabis;
443 6884998 : *t2 = t_POLMOD;
444 6884998 : switch(Rg_type(gel(c,j),&pbis,&polbis,&pabis))
445 : {
446 3853618 : case t_INT: break;
447 819709 : case t_FRAC: t[1]=1; break;
448 2208304 : case t_INTMOD: t[3]=1; break;
449 6 : case t_PADIC: t[7]=1; update_prec(pabis,pa); break;
450 3361 : default: return 0;
451 : }
452 6881637 : if (pbis) assign_or_fail(pbis,p);
453 6881637 : if (polbis) assign_or_fail(polbis,pol);
454 : }
455 3440281 : break;
456 5806132 : case t_RFRAC: t[11] = 1;
457 5806132 : if (!settype(gel(c,1),t,p,pol,pa,ff,t2,var)) return 0;
458 5806132 : c = gel(c,2); /* fall through */
459 200911780 : case t_POL: t[10] = 1;
460 200911780 : if (!RgX_settype(c,t,p,pol,pa,ff,t2,var)) return 0;
461 200911775 : if (*var == NO_VARIABLE) { *var = varn(c); break; }
462 : /* if more than one free var, ensure varn() == *var fails. FIXME: should
463 : * keep the list of all variables, later t_POLMOD may cancel them */
464 168835053 : if (*var != varn(c)) *var = MAXVARN+1;
465 168835053 : break;
466 1674 : default: return 0;
467 : }
468 3219511907 : return 1;
469 : }
470 : /* t[0] unused. Other values, if set, indicate a coefficient of type
471 : * t[1] : t_FRAC
472 : * t[2] : t_REAL
473 : * t[3] : t_INTMOD
474 : * t[4] : Unused
475 : * t[5] : t_FFELT
476 : * t[6] : Unused
477 : * t[7] : t_PADIC
478 : * t[8] : t_QUAD of rationals (t_INT/t_FRAC)
479 : * t[9]: Unused
480 : * t[10]: t_POL (multivariate polynomials)
481 : * t[11]: t_RFRAC (recursive factorisation) */
482 : /* if t2 != 0: t_POLMOD/t_QUAD/t_COMPLEX of modular (t_INTMOD/t_PADIC,
483 : * given by t) */
484 :
485 : static long
486 32061719 : choosesubtype(long *t, long t2)
487 : {
488 32061719 : if (t2 || t[11]) return 0;
489 29961555 : if (t[2] && (t[3]||t[7]||t[5])) return 0;
490 29961555 : if (t[8]) return t_QUAD;
491 29961531 : if (t[7]) return t_PADIC;
492 29961519 : if (t[5]) return t_FFELT;
493 29959190 : if (t[3]) return t_INTMOD;
494 29956651 : if (t[2]) return t_REAL;
495 29956603 : if (t[1]) return t_FRAC;
496 29730627 : return t_INT;
497 : }
498 :
499 : static long
500 312072317 : choosetype(long *t, long t2, GEN ff, GEN *pol, long var)
501 : {
502 312072317 : if (t[10] && (!*pol || var!=varn(*pol)))
503 : {
504 32061719 : long ts = choosesubtype(t, t2);
505 32061719 : if (ts==t_FFELT) *pol=ff;
506 32061719 : return ts == 0 ? t_POL: RgX_type_code(t_POL,ts);
507 : }
508 280010598 : if (t2) /* polmod/quad/complex of intmod/padic */
509 : {
510 19924128 : if (t[2] && (t[3]||t[7])) return 0;
511 19924128 : if (t[3]) return RgX_type_code(t2,t_INTMOD);
512 19898732 : if (t[7]) return RgX_type_code(t2,t_PADIC);
513 19898692 : if (t[2]) return t_COMPLEX;
514 572156 : if (t[1]) return RgX_type_code(t2,t_FRAC);
515 202741 : return RgX_type_code(t2,t_INT);
516 : }
517 260086470 : if (t[5]) /* ffelt */
518 : {
519 192542 : if (t[2]||t[8]||t[9]) return 0;
520 192542 : *pol=ff; return t_FFELT;
521 : }
522 259893928 : if (t[2]) /* inexact, real */
523 : {
524 40839353 : if (t[3]||t[7]||t[9]) return 0;
525 40839353 : return t_REAL;
526 : }
527 219054575 : if (t[10])
528 : {
529 0 : long ts = choosesubtype(t, t2);
530 0 : if (ts==t_FFELT) *pol=ff;
531 0 : return ts == 0 ? t_POL: RgX_type_code(t_POL,ts);
532 : }
533 219054575 : if (t[8]) return RgX_type_code(t_QUAD,t_INT);
534 219053849 : if (t[3]) return t_INTMOD;
535 215527051 : if (t[7]) return t_PADIC;
536 215203074 : if (t[1]) return t_FRAC;
537 208533541 : return t_INT;
538 : }
539 :
540 : static long
541 555072842 : RgX_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
542 : {
543 555072842 : long i, lx = lg(x);
544 1940482370 : for (i=2; i<lx; i++)
545 1385412318 : if (!settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
546 555070052 : return 1;
547 : }
548 :
549 : static long
550 231414096 : RgC_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
551 : {
552 231414096 : long i, l = lg(x);
553 2006992007 : for (i = 1; i<l; i++)
554 1775580346 : if (!settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
555 231411661 : return 1;
556 : }
557 :
558 : static long
559 41754647 : RgM_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
560 : {
561 41754647 : long i, l = lg(x);
562 244249207 : for (i = 1; i < l; i++)
563 202496311 : if (!RgC_settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
564 41752896 : return 1;
565 : }
566 :
567 : long
568 144897355 : Rg_type(GEN x, GEN *p, GEN *pol, long *pa)
569 : {
570 144897355 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
571 144897355 : long t2 = 0, var = NO_VARIABLE;
572 144897355 : GEN ff = NULL;
573 144897355 : *p = *pol = NULL; *pa = LONG_MAX;
574 144897355 : switch(typ(x))
575 : {
576 47539178 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_FFELT:
577 : case t_COMPLEX: case t_PADIC: case t_QUAD:
578 47539178 : if (!settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
579 47539178 : break;
580 96896595 : case t_POL: case t_SER:
581 96896595 : if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
582 96896389 : break;
583 17 : case t_VEC: case t_COL:
584 17 : if(!RgC_settype(x, t, p, pol, pa, &ff, &t2, &var)) return 0;
585 17 : break;
586 90 : case t_MAT:
587 90 : if(!RgM_settype(x, t, p, pol, pa, &ff, &t2, &var)) return 0;
588 90 : break;
589 461475 : default: return 0;
590 : }
591 144435674 : return choosetype(t,t2,ff,pol,var);
592 : }
593 :
594 : long
595 2114518 : RgX_type(GEN x, GEN *p, GEN *pol, long *pa)
596 : {
597 2114518 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
598 2114518 : long t2 = 0, var = NO_VARIABLE;
599 2114518 : GEN ff = NULL;
600 2114518 : *p = *pol = NULL; *pa = LONG_MAX;
601 2114518 : if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
602 2114475 : return choosetype(t,t2,ff,pol,var);
603 : }
604 :
605 : long
606 5179164 : RgX_Rg_type(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
607 : {
608 5179164 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
609 5179164 : long t2 = 0, var = NO_VARIABLE;
610 5179164 : GEN ff = NULL;
611 5179164 : *p = *pol = NULL; *pa = LONG_MAX;
612 5179164 : if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
613 5179164 : if (!settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
614 5179158 : return choosetype(t,t2,ff,pol,var);
615 : }
616 :
617 : long
618 123038234 : RgX_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
619 : {
620 123038234 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
621 123038234 : long t2 = 0, var = NO_VARIABLE;
622 123038234 : GEN ff = NULL;
623 123038234 : *p = *pol = NULL; *pa = LONG_MAX;
624 246074034 : if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
625 123038336 : !RgX_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
626 123035698 : return choosetype(t,t2,ff,pol,var);
627 : }
628 :
629 : long
630 1298917 : RgX_type3(GEN x, GEN y, GEN z, GEN *p, GEN *pol, long *pa)
631 : {
632 1298917 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
633 1298917 : long t2 = 0, var = NO_VARIABLE;
634 1298917 : GEN ff = NULL;
635 1298917 : *p = *pol = NULL; *pa = LONG_MAX;
636 2597834 : if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
637 2597834 : !RgX_settype(y,t,p,pol,pa,&ff,&t2,&var) ||
638 1298917 : !RgX_settype(z,t,p,pol,pa,&ff,&t2,&var)) return 0;
639 1298917 : return choosetype(t,t2,ff,pol,var);
640 : }
641 :
642 : long
643 683256 : RgM_type(GEN x, GEN *p, GEN *pol, long *pa)
644 : {
645 683256 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
646 683256 : long t2 = 0, var = NO_VARIABLE;
647 683256 : GEN ff = NULL;
648 683256 : *p = *pol = NULL; *pa = LONG_MAX;
649 683256 : if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
650 682369 : return choosetype(t,t2,ff,pol,var);
651 : }
652 :
653 : long
654 665305 : RgV_type(GEN x, GEN *p, GEN *pol, long *pa)
655 : {
656 665305 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
657 665305 : long t2 = 0, var = NO_VARIABLE;
658 665305 : GEN ff = NULL;
659 665305 : *p = *pol = NULL; *pa = LONG_MAX;
660 665305 : if (!RgC_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
661 665305 : return choosetype(t,t2,ff,pol,var);
662 : }
663 :
664 : long
665 174 : RgV_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
666 : {
667 174 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
668 174 : long t2 = 0, var = NO_VARIABLE;
669 174 : GEN ff = NULL;
670 174 : *p = *pol = NULL; *pa = LONG_MAX;
671 348 : if (!RgC_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
672 174 : !RgC_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
673 174 : return choosetype(t,t2,ff,pol,var);
674 : }
675 :
676 : long
677 28252469 : RgM_RgC_type(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
678 : {
679 28252469 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
680 28252469 : long t2 = 0, var = NO_VARIABLE;
681 28252469 : GEN ff = NULL;
682 28252469 : *p = *pol = NULL; *pa = LONG_MAX;
683 56504584 : if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
684 28253153 : !RgC_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
685 28251431 : return choosetype(t,t2,ff,pol,var);
686 : }
687 :
688 : long
689 6409626 : RgM_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
690 : {
691 6409626 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0,0};
692 6409626 : long t2 = 0, var = NO_VARIABLE;
693 6409626 : GEN ff = NULL;
694 6409626 : *p = *pol = NULL; *pa = LONG_MAX;
695 12818832 : if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
696 6409716 : !RgM_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
697 6409116 : return choosetype(t,t2,ff,pol,var);
698 : }
699 :
700 : GEN
701 50760 : factor0(GEN x, GEN flag)
702 : {
703 : ulong B;
704 50760 : long tx = typ(x);
705 50760 : if (!flag) return factor(x);
706 226 : if ((tx != t_INT && tx!=t_FRAC) || typ(flag) != t_INT)
707 156 : return factor_domain(x, flag);
708 70 : if (signe(flag) < 0) pari_err_FLAG("factor");
709 70 : switch(lgefint(flag))
710 : {
711 11 : case 2: B = 0; break;
712 59 : case 3: B = flag[2]; break;
713 0 : default: pari_err_OVERFLOW("factor [large prime bound]");
714 : return NULL; /*LCOV_EXCL_LINE*/
715 : }
716 70 : return boundfact(x, B);
717 : }
718 :
719 : GEN
720 122340 : deg1_from_roots(GEN L, long v)
721 : {
722 122340 : long i, l = lg(L);
723 122340 : GEN z = cgetg(l,t_COL);
724 371891 : for (i=1; i<l; i++)
725 249551 : gel(z,i) = deg1pol_shallow(gen_1, gneg(gel(L,i)), v);
726 122340 : return z;
727 : }
728 : GEN
729 54548 : roots_from_deg1(GEN x)
730 : {
731 54548 : long i,l = lg(x);
732 54548 : GEN r = cgetg(l,t_VEC);
733 333528 : for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = gneg(gel(P,2)); }
734 54548 : return r;
735 : }
736 :
737 : static GEN
738 36 : Qi_factor_p(GEN p)
739 : {
740 36 : GEN a, b; (void)cornacchia(gen_1, p, &a,&b);
741 36 : return mkcomplex(a, b);
742 : }
743 :
744 : static GEN
745 42 : Qi_primpart(GEN x, GEN *c)
746 : {
747 42 : GEN a = real_i(x), b = imag_i(x), n = gcdii(a, b);
748 42 : *c = n; if (n == gen_1) return x;
749 42 : retmkcomplex(diviiexact(a,n), diviiexact(b,n));
750 : }
751 :
752 : static GEN
753 60 : Qi_primpart_try(GEN x, GEN c)
754 : {
755 : GEN r, y;
756 60 : if (typ(x) == t_INT)
757 : {
758 36 : y = dvmdii(x, c, &r); if (r != gen_0) return NULL;
759 : }
760 : else
761 : {
762 24 : GEN a = gel(x,1), b = gel(x,2); y = cgetg(3, t_COMPLEX);
763 24 : gel(y,1) = dvmdii(a, c, &r); if (r != gen_0) return NULL;
764 12 : gel(y,2) = dvmdii(b, c, &r); if (r != gen_0) return NULL;
765 : }
766 48 : return y;
767 : }
768 :
769 : static int
770 78 : Qi_cmp(GEN x, GEN y)
771 : {
772 : int v;
773 78 : if (typ(x) != t_COMPLEX)
774 0 : return (typ(y) == t_COMPLEX)? -1: gcmp(x, y);
775 78 : if (typ(y) != t_COMPLEX) return 1;
776 54 : v = cmpii(gel(x,2), gel(y,2));
777 54 : if (v) return v;
778 24 : return gcmp(gel(x,1), gel(y,1));
779 : }
780 :
781 : /* 0 or canonical representative in Z[i]^* / <i> (impose imag(x) >= 0) */
782 : static GEN
783 402 : Qi_normal(GEN x)
784 : {
785 402 : if (typ(x) != t_COMPLEX) return absi_shallow(x);
786 402 : if (signe(gel(x,1)) < 0) x = gneg(x);
787 402 : if (signe(gel(x,2)) < 0) x = mulcxI(x);
788 402 : return x;
789 : }
790 :
791 : static GEN
792 42 : Qi_factor(GEN x)
793 : {
794 42 : pari_sp av = avma;
795 42 : GEN a = real_i(x), b = imag_i(x), d = gen_1, n, y, fa, P, E, P2, E2;
796 42 : long t1 = typ(a);
797 42 : long t2 = typ(b), i, j, l, exp = 0;
798 42 : if (t1 == t_FRAC) d = gel(a,2);
799 42 : if (t2 == t_FRAC) d = lcmii(d, gel(b,2));
800 42 : if (d == gen_1) y = x;
801 : else
802 : {
803 18 : y = gmul(x, d);
804 18 : a = real_i(y); t1 = typ(a);
805 18 : b = imag_i(y); t2 = typ(b);
806 : }
807 42 : if (t1 != t_INT || t2 != t_INT) return NULL;
808 42 : y = Qi_primpart(y, &n);
809 42 : fa = factor(cxnorm(y));
810 42 : P = gel(fa,1);
811 42 : E = gel(fa,2); l = lg(P);
812 42 : P2 = cgetg(l, t_COL);
813 42 : E2 = cgetg(l, t_COL);
814 90 : for (j = 1, i = l-1; i > 0; i--) /* remove largest factors first */
815 : { /* either p = 2 (ramified) or those factors split in Q(i) */
816 48 : GEN p = gel(P,i), w, w2, t, we, pe;
817 48 : long v, e = itos(gel(E,i));
818 48 : int is2 = absequaliu(p, 2);
819 48 : w = is2? mkcomplex(gen_1,gen_1): Qi_factor_p(p);
820 48 : w2 = Qi_normal( conj_i(w) );
821 : /* w * w2 * I^3 = p, w2 = conj(w) * I */
822 48 : pe = powiu(p, e);
823 48 : we = gpowgs(w, e);
824 48 : t = Qi_primpart_try( gmul(y, conj_i(we)), pe );
825 48 : if (t) y = t; /* y /= w^e */
826 : else {
827 : /* y /= conj(w)^e, should be y /= w2^e */
828 12 : y = Qi_primpart_try( gmul(y, we), pe );
829 12 : swap(w, w2); exp -= e; /* += 3*e mod 4 */
830 : }
831 48 : gel(P,i) = w;
832 48 : v = Z_pvalrem(n, p, &n);
833 48 : if (v) {
834 6 : exp -= v; /* += 3*v mod 4 */
835 6 : if (is2) v <<= 1; /* 2 = w^2 I^3 */
836 : else {
837 0 : gel(P2,j) = w2;
838 0 : gel(E2,j) = utoipos(v); j++;
839 : }
840 6 : gel(E,i) = stoi(e + v);
841 : }
842 48 : v = Z_pvalrem(d, p, &d);
843 48 : if (v) {
844 6 : exp += v; /* -= 3*v mod 4 */
845 6 : if (is2) v <<= 1; /* 2 is ramified */
846 : else {
847 6 : gel(P2,j) = w2;
848 6 : gel(E2,j) = utoineg(v); j++;
849 : }
850 6 : gel(E,i) = stoi(e - v);
851 : }
852 48 : exp &= 3;
853 : }
854 42 : if (j > 1) {
855 6 : long k = 1;
856 6 : GEN P1 = cgetg(l, t_COL);
857 6 : GEN E1 = cgetg(l, t_COL);
858 : /* remove factors with exponent 0 */
859 12 : for (i = 1; i < l; i++)
860 6 : if (signe(gel(E,i)))
861 : {
862 0 : gel(P1,k) = gel(P,i);
863 0 : gel(E1,k) = gel(E,i);
864 0 : k++;
865 : }
866 6 : setlg(P1, k); setlg(E1, k);
867 6 : setlg(P2, j); setlg(E2, j);
868 6 : fa = famat_mul_shallow(mkmat2(P1,E1), mkmat2(P2,E2));
869 : }
870 42 : if (!equali1(n) || !equali1(d))
871 : {
872 24 : GEN Fa = factor(Qdivii(n, d));
873 24 : P = gel(Fa,1); l = lg(P);
874 24 : E = gel(Fa,2);
875 60 : for (i = 1; i < l; i++)
876 : {
877 36 : GEN w, p = gel(P,i);
878 : long e;
879 : int is2;
880 36 : switch(mod4(p))
881 : {
882 12 : case 3: continue;
883 12 : case 2: is2 = 1; break;
884 12 : default:is2 = 0; break;
885 : }
886 24 : e = itos(gel(E,i));
887 24 : w = is2? mkcomplex(gen_1,gen_1): Qi_factor_p(p);
888 24 : gel(P,i) = w;
889 24 : if (is2)
890 12 : gel(E,i) = stoi(2*e);
891 : else
892 : {
893 12 : P = vec_append(P, Qi_normal( conj_i(w) ));
894 12 : E = vec_append(E, gel(E,i));
895 : }
896 24 : exp -= e; /* += 3*e mod 4 */
897 24 : exp &= 3;
898 : }
899 24 : gel(Fa,1) = P;
900 24 : gel(Fa,2) = E;
901 24 : fa = famat_mul_shallow(fa, Fa);
902 : }
903 42 : fa = sort_factor(fa, (void*)&Qi_cmp, &cmp_nodata);
904 :
905 42 : y = gmul(y, powIs(exp));
906 42 : if (!gequal1(y)) {
907 30 : gel(fa,1) = vec_prepend(gel(fa,1), y);
908 30 : gel(fa,2) = vec_prepend(gel(fa,2), gen_1);
909 : }
910 42 : return gc_GEN(av, fa);
911 : }
912 :
913 : GEN
914 11908 : Q_factor_limit(GEN x, ulong lim)
915 : {
916 11908 : pari_sp av = avma;
917 : GEN a, b;
918 11908 : if (typ(x) == t_INT) return Z_factor_limit(x, lim);
919 4291 : a = Z_factor_limit(gel(x,1), lim);
920 4291 : b = Z_factor_limit(gel(x,2), lim); gel(b,2) = ZC_neg(gel(b,2));
921 4291 : return gc_GEN(av, ZM_merge_factor(a,b));
922 : }
923 : GEN
924 18413 : Q_factor(GEN x)
925 : {
926 18413 : pari_sp av = avma;
927 : GEN a, b;
928 18413 : if (typ(x) == t_INT) return Z_factor(x);
929 30 : a = Z_factor(gel(x,1));
930 30 : b = Z_factor(gel(x,2)); gel(b,2) = ZC_neg(gel(b,2));
931 30 : return gc_GEN(av, ZM_merge_factor(a,b));
932 : }
933 :
934 : /* replace quadratic number over Fp or Q by t_POL in v */
935 : static GEN
936 360 : quadratic_to_RgX(GEN z, long v)
937 : {
938 : GEN a, b;
939 360 : switch(typ(z))
940 : {
941 294 : case t_INT: case t_FRAC: case t_INTMOD: return z;
942 30 : case t_COMPLEX: a = gel(z,2); b = gel(z,1); break;
943 36 : case t_QUAD: a = gel(z,3); b = gel(z,2); break;
944 0 : default: pari_err_IMPL("factor for general polynomials"); /* paranoia */
945 : return NULL; /* LCOV_EXCL_LINE */
946 : }
947 66 : return deg1pol_shallow(a, b, v);
948 : }
949 : /* replace t_QUAD/t_COMPLEX [of rationals] coeffs by t_POL in v */
950 : static GEN
951 84 : RgX_fix_quadratic(GEN x, long v)
952 444 : { pari_APPLY_pol_normalized(quadratic_to_RgX(gel(x,i), v)); }
953 : static GEN
954 216 : RgXQ_factor_i(GEN x, GEN T, GEN p, long t1, long t2, long *pv)
955 : {
956 216 : *pv = -1;
957 216 : if (t2 == t_PADIC) return NULL;
958 186 : if (t2 == t_INTMOD)
959 : {
960 48 : T = RgX_to_FpX(T,p);
961 48 : if (!FpX_is_irred(T,p)) return NULL;
962 : }
963 168 : if (t1 != t_POLMOD)
964 : { /* replace w in x by t_POL */
965 84 : if (t2 != t_INTMOD) T = leafcopy(T);
966 84 : *pv = fetch_var(); setvarn(T, *pv);
967 84 : x = RgX_fix_quadratic(x, *pv);
968 : }
969 168 : if (t2 == t_INTMOD) return factmod(x, mkvec2(p,T));
970 138 : return nffactor(T, x);
971 : }
972 : static GEN
973 216 : RgXQ_factor(GEN x, GEN T, GEN p, long tx)
974 : {
975 216 : pari_sp av = avma;
976 : long t1, t2, v;
977 : GEN w, y;
978 216 : RgX_type_decode(tx, &t1, &t2);
979 216 : y = RgXQ_factor_i(x, T, p, t1, t2, &v);
980 216 : if (!y) pari_err_IMPL("factor for general polynomials");
981 168 : if (v < 0) return gc_upto(av, y);
982 : /* substitute back w */
983 84 : w = (t1 == t_COMPLEX)? gen_I(): mkquad(T,gen_0,gen_1);
984 84 : gel(y,1) = gsubst(liftpol_shallow(gel(y,1)), v, w);
985 84 : (void)delete_var(); return gc_GEN(av, y);
986 : }
987 :
988 : static GEN
989 24 : RX_factor(GEN x, long prec)
990 : {
991 24 : GEN y = cgetg(3,t_MAT), R, P;
992 24 : pari_sp av = avma;
993 24 : long v = varn(x), i, l, r1;
994 :
995 24 : R = cleanroots(x, prec); l = lg(R);
996 60 : for (r1 = 1; r1 < l; r1++)
997 42 : if (typ(gel(R,r1)) == t_COMPLEX) break;
998 24 : l = (r1+l)>>1; P = cgetg(l,t_COL);
999 60 : for (i = 1; i < r1; i++)
1000 36 : gel(P,i) = deg1pol_shallow(gen_1, negr(gel(R,i)), v);
1001 30 : for ( ; i < l; i++)
1002 : {
1003 6 : GEN a = gel(R,2*i-r1), t;
1004 6 : t = gmul2n(gel(a,1), 1); togglesign(t);
1005 6 : gel(P,i) = deg2pol_shallow(gen_1, t, gnorm(a), v);
1006 : }
1007 24 : gel(y,1) = gc_upto(av, P);
1008 24 : gel(y,2) = const_col(l-1, gen_1); return y;
1009 : }
1010 : static GEN
1011 18 : CX_factor(GEN x, long prec)
1012 : {
1013 18 : GEN y = cgetg(3,t_MAT), R;
1014 18 : pari_sp av = avma;
1015 18 : long v = varn(x);
1016 :
1017 18 : R = roots(x, prec);
1018 18 : gel(y,1) = gc_upto(av, deg1_from_roots(R, v));
1019 18 : gel(y,2) = const_col(degpol(x), gen_1); return y;
1020 : }
1021 :
1022 : static GEN
1023 11869 : RgX_factor(GEN x, GEN dom)
1024 : {
1025 : GEN p, T;
1026 11869 : long pa, tx = dom ? RgX_Rg_type(x,dom,&p,&T,&pa): RgX_type(x,&p,&T,&pa);
1027 11869 : if (tx>>RgX_type_shift==t_POL) tx = t_POL;
1028 11869 : switch(tx)
1029 : {
1030 6 : case 0: pari_err_IMPL("factor for general polynomials");
1031 228 : case t_POL: return RgXY_factor(x, dom);
1032 10958 : case t_INT: return ZX_factor(x);
1033 6 : case t_FRAC: return QX_factor(x);
1034 294 : case t_INTMOD: return factmod(x, p);
1035 35 : case t_PADIC: return factorpadic(x, p, pa);
1036 84 : case t_FFELT: return FFX_factor(x, T);
1037 18 : case t_COMPLEX: return CX_factor(x, pa);
1038 24 : case t_REAL: return RX_factor(x, pa);
1039 : }
1040 216 : return RgXQ_factor(x, T, p, tx);
1041 : }
1042 :
1043 : static GEN
1044 53708 : factor_domain(GEN x, GEN dom)
1045 : {
1046 53708 : long tx = typ(x), tdom = dom ? typ(dom): 0;
1047 : pari_sp av;
1048 :
1049 53708 : if (gequal0(x))
1050 45 : switch(tx)
1051 : {
1052 45 : case t_INT:
1053 : case t_COMPLEX:
1054 : case t_POL:
1055 45 : case t_RFRAC: return prime_fact(x);
1056 0 : default: pari_err_TYPE("factor",x);
1057 : }
1058 53663 : av = avma;
1059 53663 : switch(tx)
1060 : {
1061 2270 : case t_POL: return RgX_factor(x, dom);
1062 30 : case t_RFRAC: {
1063 30 : GEN a = gel(x,1), b = gel(x,2);
1064 30 : GEN y = famat_inv_shallow(RgX_factor(b, dom));
1065 30 : if (typ(a)==t_POL) y = famat_mul_shallow(RgX_factor(a, dom), y);
1066 30 : return gc_GEN(av, sort_factor_pol(y, cmp_universal));
1067 : }
1068 51303 : case t_INT: if (tdom==0 || tdom==t_INT) return Z_factor(x);
1069 24 : case t_FRAC: if (tdom==0 || tdom==t_INT) return Q_factor(x);
1070 : case t_COMPLEX: /* fall through */
1071 42 : if (tdom==0 || tdom==t_COMPLEX)
1072 42 : { GEN y = Qi_factor(x); if (y) return y; }
1073 : /* fall through */
1074 : }
1075 0 : pari_err_TYPE("factor",x);
1076 : return NULL; /* LCOV_EXCL_LINE */
1077 : }
1078 :
1079 : GEN
1080 53096 : factor(GEN x) { return factor_domain(x, NULL); }
1081 :
1082 : /*******************************************************************/
1083 : /* */
1084 : /* ROOTS --> MONIC POLYNOMIAL */
1085 : /* */
1086 : /*******************************************************************/
1087 : static GEN
1088 1431899 : normalized_mul(void *E, GEN x, GEN y)
1089 : {
1090 1431899 : long a = gel(x,1)[1], b = gel(y,1)[1];
1091 : (void) E;
1092 1431899 : return mkvec2(mkvecsmall(a + b),
1093 1431899 : RgX_mul_normalized(gel(x,2),a, gel(y,2),b));
1094 : }
1095 : /* L = [Vecsmall([a]), A], with a > 0, A an RgX, deg(A) < a; return X^a + A */
1096 : static GEN
1097 885942 : normalized_to_RgX(GEN L)
1098 : {
1099 885942 : long i, a = gel(L,1)[1];
1100 885942 : GEN A = gel(L,2);
1101 885942 : GEN z = cgetg(a + 3, t_POL);
1102 885942 : z[1] = evalsigne(1) | evalvarn(varn(A));
1103 4863856 : for (i = 2; i < lg(A); i++) gel(z,i) = gcopy(gel(A,i));
1104 890360 : for ( ; i < a+2; i++) gel(z,i) = gen_0;
1105 885942 : gel(z,i) = gen_1; return z;
1106 : }
1107 :
1108 : static GEN
1109 12 : roots_to_pol_FpV(GEN x, long v, GEN p)
1110 : {
1111 12 : pari_sp av = avma;
1112 : GEN r;
1113 12 : if (lgefint(p) == 3)
1114 : {
1115 12 : ulong pp = uel(p, 2);
1116 12 : r = Flx_to_ZX_inplace(Flv_roots_to_pol(RgV_to_Flv(x, pp), pp, v<<VARNSHIFT));
1117 : }
1118 : else
1119 0 : r = FpV_roots_to_pol(RgV_to_FpV(x, p), p, v);
1120 12 : return gc_upto(av, FpX_to_mod(r, p));
1121 : }
1122 :
1123 : static GEN
1124 6 : roots_to_pol_FqV(GEN x, long v, GEN pol, GEN p)
1125 : {
1126 6 : pari_sp av = avma;
1127 6 : GEN r, T = RgX_to_FpX(pol, p);
1128 6 : if (signe(T)==0) pari_err_OP("/", x, pol);
1129 6 : r = FqV_roots_to_pol(RgC_to_FqC(x, T, p), T, p, v);
1130 6 : return gc_upto(av, FpXQX_to_mod(r, T, p));
1131 : }
1132 :
1133 : static GEN
1134 665233 : roots_to_pol_fast(GEN x, long v)
1135 : {
1136 : GEN p, pol;
1137 : long pa;
1138 665233 : long t = RgV_type(x, &p,&pol,&pa);
1139 665233 : switch(t)
1140 : {
1141 12 : case t_INTMOD: return roots_to_pol_FpV(x, v, p);
1142 12 : case t_FFELT: return FFV_roots_to_pol(x, pol, v);
1143 6 : case RgX_type_code(t_POLMOD, t_INTMOD):
1144 6 : return roots_to_pol_FqV(x, v, pol, p);
1145 665203 : default: return NULL;
1146 : }
1147 : }
1148 :
1149 : /* compute prod (x - a[i]) */
1150 : GEN
1151 665296 : roots_to_pol(GEN a, long v)
1152 : {
1153 665296 : pari_sp av = avma;
1154 665296 : long i, k, lx = lg(a);
1155 : GEN L;
1156 665296 : if (lx == 1) return pol_1(v);
1157 665233 : L = roots_to_pol_fast(a, v);
1158 665233 : if (L) return L;
1159 665203 : L = cgetg(lx, t_VEC);
1160 1426077 : for (k=1,i=1; i<lx-1; i+=2)
1161 : {
1162 760874 : GEN s = gel(a,i), t = gel(a,i+1);
1163 760874 : GEN x0 = gmul(s,t);
1164 760874 : GEN x1 = gneg(gadd(s,t));
1165 760874 : gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
1166 : }
1167 1259697 : if (i < lx) gel(L,k++) = mkvec2(mkvecsmall(1),
1168 594494 : scalarpol_shallow(gneg(gel(a,i)), v));
1169 665203 : setlg(L, k); L = gen_product(L, NULL, normalized_mul);
1170 665203 : return gc_upto(av, normalized_to_RgX(L));
1171 : }
1172 :
1173 : /* prod_{i=1..r1} (x - a[i]) prod_{i=1..r2} (x - a[i])(x - conj(a[i]))*/
1174 : GEN
1175 220739 : roots_to_pol_r1(GEN a, long v, long r1)
1176 : {
1177 220739 : pari_sp av = avma;
1178 220739 : long i, k, lx = lg(a);
1179 : GEN L;
1180 220739 : if (lx == 1) return pol_1(v);
1181 220739 : L = cgetg(lx, t_VEC);
1182 577145 : for (k=1,i=1; i<r1; i+=2)
1183 : {
1184 356406 : GEN s = gel(a,i), t = gel(a,i+1);
1185 356406 : GEN x0 = gmul(s,t);
1186 356406 : GEN x1 = gneg(gadd(s,t));
1187 356406 : gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
1188 : }
1189 279595 : if (i < r1+1) gel(L,k++) = mkvec2(mkvecsmall(1),
1190 58856 : scalarpol_shallow(gneg(gel(a,i)), v));
1191 767950 : for (i=r1+1; i<lx; i++)
1192 : {
1193 547211 : GEN s = gel(a,i);
1194 547211 : GEN x0 = gnorm(s);
1195 547211 : GEN x1 = gneg(gtrace(s));
1196 547211 : gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
1197 : }
1198 220739 : setlg(L, k); L = gen_product(L, NULL, normalized_mul);
1199 220739 : return gc_upto(av, normalized_to_RgX(L));
1200 : }
1201 :
1202 : GEN
1203 48 : polfromroots(GEN a, long v)
1204 : {
1205 48 : if (!is_vec_t(typ(a)))
1206 0 : pari_err_TYPE("polfromroots",a);
1207 48 : if (v < 0) v = 0;
1208 48 : if (varncmp(gvar(a), v) <= 0) pari_err_PRIORITY("polfromroots",a,"<=",v);
1209 42 : return roots_to_pol(a, v);
1210 : }
1211 :
1212 : /*******************************************************************/
1213 : /* */
1214 : /* FACTORBACK */
1215 : /* */
1216 : /*******************************************************************/
1217 : static GEN
1218 47558551 : mul(void *a, GEN x, GEN y) { (void)a; return gmul(x,y);}
1219 : static GEN
1220 68848897 : powi(void *a, GEN x, GEN y) { (void)a; return powgi(x,y);}
1221 : static GEN
1222 24465988 : Fpmul(void *a, GEN x, GEN y) { return Fp_mul(x,y,(GEN)a); }
1223 : static GEN
1224 209647 : Fppow(void *a, GEN x, GEN n) { return Fp_pow(x,n,(GEN)a); }
1225 :
1226 : /* [L,e] = [fa, NULL] or [elts, NULL] or [elts, exponents] */
1227 : GEN
1228 29217186 : gen_factorback(GEN L, GEN e, void *data, GEN (*_mul)(void*,GEN,GEN),
1229 : GEN (*_pow)(void*,GEN,GEN), GEN (*_one)(void*))
1230 : {
1231 29217186 : pari_sp av = avma;
1232 : long k, l, lx;
1233 : GEN p,x;
1234 :
1235 29217186 : if (e) /* supplied vector of exponents */
1236 1606313 : p = L;
1237 : else
1238 : {
1239 27610873 : switch(typ(L)) {
1240 7140971 : case t_VEC:
1241 : case t_COL: /* product of the L[i] */
1242 7140971 : if (lg(L)==1) return _one? _one(data): gen_1;
1243 7075403 : return gc_upto(av, gen_product(L, data, _mul));
1244 20469902 : case t_MAT: /* genuine factorization */
1245 20469902 : l = lg(L);
1246 20469902 : if (l == 3) break;
1247 : /*fall through*/
1248 : default:
1249 6 : pari_err_TYPE("factorback [not a factorization]", L);
1250 : }
1251 20469896 : p = gel(L,1);
1252 20469896 : e = gel(L,2);
1253 : }
1254 22076209 : if (!is_vec_t(typ(p))) pari_err_TYPE("factorback [not a vector]", p);
1255 : /* p = elts, e = expo */
1256 22076197 : lx = lg(p);
1257 : /* check whether e is an integral vector of correct length */
1258 22076197 : switch(typ(e))
1259 : {
1260 105073 : case t_VECSMALL:
1261 105073 : if (lx != lg(e))
1262 0 : pari_err_TYPE("factorback [not an exponent vector]", e);
1263 105073 : if (lx == 1) return _one? _one(data): gen_1;
1264 104886 : x = cgetg(lx,t_VEC);
1265 1066365 : for (l=1,k=1; k<lx; k++)
1266 961479 : if (e[k]) gel(x,l++) = _pow(data, gel(p,k), stoi(e[k]));
1267 104886 : break;
1268 21971118 : case t_VEC: case t_COL:
1269 21971118 : if (lx != lg(e) || !RgV_is_ZV(e))
1270 12 : pari_err_TYPE("factorback [not an exponent vector]", e);
1271 21971106 : if (lx == 1) return _one? _one(data): gen_1;
1272 21874824 : x = cgetg(lx,t_VEC);
1273 91545700 : for (l=1,k=1; k<lx; k++)
1274 69670876 : if (signe(gel(e,k))) gel(x,l++) = _pow(data, gel(p,k), gel(e,k));
1275 21874824 : break;
1276 6 : default:
1277 6 : pari_err_TYPE("factorback [not an exponent vector]", e);
1278 : return NULL;/*LCOV_EXCL_LINE*/
1279 : }
1280 21979710 : if (l==1) return gc_upto(av, _one? _one(data): gen_1);
1281 21924744 : x[0] = evaltyp(t_VEC) | _evallg(l);
1282 21924744 : return gc_upto(av, gen_product(x, data, _mul));
1283 : }
1284 :
1285 : GEN
1286 7233747 : FpV_factorback(GEN L, GEN e, GEN p)
1287 7233747 : { return gen_factorback(L, e, (void*)p, &Fpmul, &Fppow, NULL); }
1288 :
1289 : ulong
1290 80348 : Flv_factorback(GEN L, GEN e, ulong p)
1291 : {
1292 80348 : long i, l = lg(e);
1293 80348 : ulong r = 1UL, ri = 1UL;
1294 387786 : for (i = 1; i < l; i++)
1295 : {
1296 307438 : long c = e[i];
1297 307438 : if (!c) continue;
1298 131432 : if (c < 0)
1299 0 : ri = Fl_mul(ri, Fl_powu(L[i],-c,p), p);
1300 : else
1301 131432 : r = Fl_mul(r, Fl_powu(L[i],c,p), p);
1302 : }
1303 80348 : if (ri != 1UL) r = Fl_div(r, ri, p);
1304 80348 : return r;
1305 : }
1306 : GEN
1307 2022 : FlxqV_factorback(GEN L, GEN e, GEN Tp, ulong p)
1308 : {
1309 2022 : pari_sp av = avma;
1310 2022 : GEN Hi = NULL, H = NULL;
1311 2022 : long i, l = lg(L), v = get_Flx_var(Tp);
1312 140120 : for (i = 1; i < l; i++)
1313 : {
1314 138098 : GEN x, ei = gel(e,i);
1315 138098 : long s = signe(ei);
1316 138098 : if (!s) continue;
1317 131465 : x = Flxq_pow(gel(L,i), s > 0? ei: negi(ei), Tp, p);
1318 131465 : if (s > 0)
1319 66201 : H = H? Flxq_mul(H, x, Tp, p): x;
1320 : else
1321 65264 : Hi = Hi? Flxq_mul(Hi, x, Tp, p): x;
1322 : }
1323 2022 : if (!Hi)
1324 : {
1325 0 : if (!H) { set_avma(av); return pol1_Flx(v); }
1326 0 : return gc_leaf(av, H);
1327 : }
1328 2022 : Hi = Flxq_inv(Hi, Tp, p);
1329 2022 : return gc_leaf(av, H? Flxq_mul(H,Hi,Tp,p): Hi);
1330 : }
1331 : GEN
1332 12 : FqV_factorback(GEN L, GEN e, GEN Tp, GEN p)
1333 : {
1334 12 : pari_sp av = avma;
1335 12 : GEN Hi = NULL, H = NULL;
1336 12 : long i, l = lg(L), small = typ(e) == t_VECSMALL;
1337 1332 : for (i = 1; i < l; i++)
1338 : {
1339 : GEN x;
1340 : long s;
1341 1320 : if (small)
1342 : {
1343 0 : s = e[i]; if (!s) continue;
1344 0 : x = Fq_powu(gel(L,i), labs(s), Tp, p);
1345 : }
1346 : else
1347 : {
1348 1320 : GEN ei = gel(e,i);
1349 1320 : s = signe(ei); if (!s) continue;
1350 1320 : x = Fq_pow(gel(L,i), s > 0? ei: negi(ei), Tp, p);
1351 : }
1352 1320 : if (s > 0)
1353 702 : H = H? Fq_mul(H, x, Tp, p): x;
1354 : else
1355 618 : Hi = Hi? Fq_mul(Hi, x, Tp, p): x;
1356 : }
1357 12 : if (Hi)
1358 : {
1359 6 : Hi = Fq_inv(Hi, Tp, p);
1360 6 : H = H? Fq_mul(H,Hi,Tp,p): Hi;
1361 : }
1362 6 : else if (!H) return gc_const(av, gen_1);
1363 12 : return gc_upto(av, H);
1364 : }
1365 :
1366 : GEN
1367 21686812 : factorback2(GEN L, GEN e) { return gen_factorback(L, e, NULL, &mul, &powi, NULL); }
1368 : GEN
1369 1126366 : factorback(GEN fa) { return factorback2(fa, NULL); }
1370 :
1371 : GEN
1372 8574 : vecprod(GEN v)
1373 : {
1374 8574 : pari_sp av = avma;
1375 8574 : if (!is_vec_t(typ(v)))
1376 0 : pari_err_TYPE("vecprod", v);
1377 8574 : if (lg(v) == 1) return gen_1;
1378 7908 : return gc_GEN(av, gen_product(v, NULL, mul));
1379 : }
1380 :
1381 : static int
1382 9569 : RgX_is_irred_i(GEN x)
1383 : {
1384 : GEN y, p, pol;
1385 9569 : long l = lg(x), pa;
1386 :
1387 9569 : if (!signe(x) || l <= 3) return 0;
1388 9569 : switch(RgX_type(x,&p,&pol,&pa))
1389 : {
1390 18 : case t_INTMOD: return FpX_is_irred(RgX_to_FpX(x,p), p);
1391 0 : case t_COMPLEX: return l == 4;
1392 0 : case t_REAL:
1393 0 : if (l == 4) return 1;
1394 0 : if (l > 5) return 0;
1395 0 : return gsigne(RgX_disc(x)) > 0;
1396 : }
1397 9551 : y = RgX_factor(x, NULL);
1398 9551 : return (lg(gcoeff(y,1,1))==l);
1399 : }
1400 : static int
1401 9569 : RgX_is_irred(GEN x)
1402 9569 : { pari_sp av = avma; return gc_bool(av, RgX_is_irred_i(x)); }
1403 : long
1404 9569 : polisirreducible(GEN x)
1405 : {
1406 9569 : long tx = typ(x);
1407 9569 : if (tx == t_POL) return RgX_is_irred(x);
1408 0 : if (!is_scalar_t(tx)) pari_err_TYPE("polisirreducible",x);
1409 0 : return 0;
1410 : }
1411 :
1412 : /*******************************************************************/
1413 : /* */
1414 : /* GENERIC GCD */
1415 : /* */
1416 : /*******************************************************************/
1417 : static GEN
1418 5365 : gcd3(GEN x, GEN y, GEN z) { return ggcd(ggcd(x, y), z); }
1419 :
1420 : /* x is a COMPLEX or a QUAD */
1421 : static GEN
1422 2808 : triv_cont_gcd(GEN x, GEN y)
1423 : {
1424 2808 : pari_sp av = avma;
1425 : GEN a, b;
1426 2808 : if (typ(x)==t_COMPLEX)
1427 : {
1428 2496 : a = gel(x,1); b = gel(x,2);
1429 2496 : if (typ(a) == t_REAL || typ(b) == t_REAL) return gen_1;
1430 : }
1431 : else
1432 : {
1433 312 : a = gel(x,2); b = gel(x,3);
1434 : }
1435 330 : return gc_upto(av, gcd3(a, b, y));
1436 : }
1437 :
1438 : /* y is a PADIC, x a rational number or an INTMOD */
1439 : static GEN
1440 2320 : padic_gcd(GEN x, GEN y)
1441 : {
1442 2320 : GEN p = padic_p(y);
1443 2320 : long v = gvaluation(x,p), w = valp(y);
1444 2320 : if (w < v) v = w;
1445 2320 : return powis(p, v);
1446 : }
1447 :
1448 : static void
1449 732 : Zi_mul3(GEN xr, GEN xi, GEN yr, GEN yi, GEN *zr, GEN *zi)
1450 : {
1451 732 : GEN p3 = addii(xr,xi);
1452 732 : GEN p4 = addii(yr,yi);
1453 732 : GEN p1 = mulii(xr,yr);
1454 732 : GEN p2 = mulii(xi,yi);
1455 732 : p3 = mulii(p3,p4);
1456 732 : p4 = addii(p2,p1);
1457 732 : *zr = subii(p1,p2); *zi = subii(p3,p4);
1458 732 : }
1459 :
1460 : static GEN
1461 366 : Zi_rem(GEN x, GEN y)
1462 : {
1463 366 : GEN xr = real_i(x), xi = imag_i(x);
1464 366 : GEN yr = real_i(y), yi = imag_i(y);
1465 366 : GEN n = addii(sqri(yr), sqri(yi));
1466 : GEN ur, ui, zr, zi;
1467 366 : Zi_mul3(xr, xi, yr, negi(yi), &ur, &ui);
1468 366 : Zi_mul3(yr, yi, diviiround(ur, n), diviiround(ui, n), &zr, &zi);
1469 366 : return mkcomplex(subii(xr,zr), subii(xi,zi));
1470 : }
1471 :
1472 : static GEN
1473 342 : Qi_gcd(GEN x, GEN y)
1474 : {
1475 342 : pari_sp av = avma, btop;
1476 : GEN dx, dy;
1477 342 : x = Q_remove_denom(x, &dx);
1478 342 : y = Q_remove_denom(y, &dy);
1479 342 : btop = avma;
1480 708 : while (!gequal0(y))
1481 : {
1482 366 : GEN z = Zi_rem(x,y);
1483 366 : x = y; y = z;
1484 366 : if (gc_needed(btop,1)) {
1485 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Qi_gcd");
1486 0 : (void)gc_all(btop,2, &x,&y);
1487 : }
1488 : }
1489 342 : x = Qi_normal(x);
1490 342 : if (typ(x) == t_COMPLEX)
1491 : {
1492 240 : if (gequal0(gel(x,2))) x = gel(x,1);
1493 180 : else if (gequal0(gel(x,1))) x = gel(x,2);
1494 : }
1495 342 : if (!dx && !dy) return gc_GEN(av, x);
1496 30 : return gc_upto(av, gdiv(x, dx? (dy? lcmii(dx, dy): dx): dy));
1497 : }
1498 :
1499 : static int
1500 2790 : c_is_rational(GEN x)
1501 2790 : { return is_rational_t(typ(gel(x,1))) && is_rational_t(typ(gel(x,2))); }
1502 : static GEN
1503 1200 : c_zero_gcd(GEN c)
1504 : {
1505 1200 : GEN x = gel(c,1), y = gel(c,2);
1506 1200 : long tx = typ(x), ty = typ(y);
1507 1200 : if (tx == t_REAL || ty == t_REAL) return gen_1;
1508 48 : if (tx == t_PADIC || tx == t_INTMOD
1509 48 : || ty == t_PADIC || ty == t_INTMOD) return ggcd(x, y);
1510 42 : return Qi_gcd(c, gen_0);
1511 : }
1512 :
1513 : /* gcd(x, 0) */
1514 : static GEN
1515 6771561 : zero_gcd(GEN x)
1516 : {
1517 : pari_sp av;
1518 6771561 : switch(typ(x))
1519 : {
1520 37868 : case t_INT: return absi(x);
1521 39718 : case t_FRAC: return absfrac(x);
1522 1200 : case t_COMPLEX: return c_zero_gcd(x);
1523 587 : case t_REAL: return gen_1;
1524 642 : case t_PADIC: return powis(padic_p(x), valp(x));
1525 184 : case t_SER: return pol_xnall(valser(x), varn(x));
1526 2807 : case t_POLMOD: {
1527 2807 : GEN d = gel(x,2);
1528 2807 : if (typ(d) == t_POL && varn(d) == varn(gel(x,1))) return content(d);
1529 366 : return isinexact(d)? zero_gcd(d): gcopy(d);
1530 : }
1531 6488136 : case t_POL:
1532 6488136 : if (!isinexact(x)) break;
1533 0 : av = avma;
1534 0 : return gc_upto(av, monomialcopy(content(x), RgX_val(x), varn(x)));
1535 :
1536 180470 : case t_RFRAC:
1537 180470 : if (!isinexact(x)) break;
1538 0 : av = avma;
1539 0 : return gc_upto(av, gdiv(zero_gcd(gel(x,1)), gel(x,2)));
1540 : }
1541 6688555 : return gcopy(x);
1542 : }
1543 : /* z is an exact zero, t_INT, t_INTMOD or t_FFELT */
1544 : static GEN
1545 7294611 : zero_gcd2(GEN y, GEN z)
1546 : {
1547 : pari_sp av;
1548 7294611 : switch(typ(z))
1549 : {
1550 6754946 : case t_INT: return zero_gcd(y);
1551 536884 : case t_INTMOD:
1552 536884 : av = avma;
1553 536884 : return gc_upto(av, gmul(y, mkintmod(gen_1,gel(z,1))));
1554 2781 : case t_FFELT:
1555 2781 : av = avma;
1556 2781 : return gc_upto(av, gmul(y, FF_1(z)));
1557 0 : default:
1558 0 : pari_err_TYPE("zero_gcd", z);
1559 : return NULL;/*LCOV_EXCL_LINE*/
1560 : }
1561 : }
1562 : static GEN
1563 1539917 : cont_gcd_pol_i(GEN x, GEN y) { return scalarpol(ggcd(content(x),y), varn(x));}
1564 : /* tx = t_POL, y considered as constant */
1565 : static GEN
1566 1539917 : cont_gcd_pol(GEN x, GEN y)
1567 1539917 : { pari_sp av = avma; return gc_upto(av, cont_gcd_pol_i(x,y)); }
1568 : /* tx = t_RFRAC, y considered as constant */
1569 : static GEN
1570 7707 : cont_gcd_rfrac(GEN x, GEN y)
1571 : {
1572 7707 : pari_sp av = avma;
1573 7707 : GEN cx; x = primitive_part(x, &cx);
1574 : /* e.g. Mod(1,2) / (2*y+1) => primitive_part = Mod(1,2)*y^0 */
1575 7707 : if (typ(x) != t_RFRAC) x = cont_gcd_pol_i(x, y);
1576 7707 : else x = gred_rfrac_simple(ggcd(cx? cx: gen_1, y), gel(x,2));
1577 7707 : return gc_upto(av, x);
1578 : }
1579 : /* !is_const_t(tx), tx != t_POL,t_RFRAC, y considered as constant */
1580 : static GEN
1581 2228 : cont_gcd_gen(GEN x, GEN y)
1582 : {
1583 2228 : pari_sp av = avma;
1584 2228 : return gc_upto(av, ggcd(content(x),y));
1585 : }
1586 : /* !is_const(tx), y considered as constant */
1587 : static GEN
1588 1549834 : cont_gcd(GEN x, long tx, GEN y)
1589 : {
1590 1549834 : switch(tx)
1591 : {
1592 7707 : case t_RFRAC: return cont_gcd_rfrac(x,y);
1593 1539899 : case t_POL: return cont_gcd_pol(x,y);
1594 2228 : default: return cont_gcd_gen(x,y);
1595 : }
1596 : }
1597 : static GEN
1598 10734445 : gcdiq(GEN x, GEN y)
1599 : {
1600 : GEN z;
1601 10734445 : if (!signe(x)) return Q_abs(y);
1602 3992028 : z = cgetg(3,t_FRAC);
1603 3992028 : gel(z,1) = gcdii(x,gel(y,1));
1604 3992028 : gel(z,2) = icopy(gel(y,2));
1605 3992028 : return z;
1606 : }
1607 : static GEN
1608 23586736 : gcdqq(GEN x, GEN y)
1609 : {
1610 23586736 : GEN z = cgetg(3,t_FRAC);
1611 23586736 : gel(z,1) = gcdii(gel(x,1), gel(y,1));
1612 23586736 : gel(z,2) = lcmii(gel(x,2), gel(y,2));
1613 23586736 : return z;
1614 : }
1615 : /* assume x,y t_INT or t_FRAC */
1616 : GEN
1617 846153333 : Q_gcd(GEN x, GEN y)
1618 : {
1619 846153333 : long tx = typ(x), ty = typ(y);
1620 846153333 : if (tx == t_INT)
1621 814458818 : { return (ty == t_INT)? gcdii(x,y): gcdiq(x,y); }
1622 : else
1623 31694515 : { return (ty == t_INT)? gcdiq(y,x): gcdqq(x,y); }
1624 : }
1625 :
1626 : /* t_QUADs */
1627 : static GEN
1628 132 : qgcd(GEN x, GEN y)
1629 : {
1630 132 : pari_sp av = avma;
1631 132 : GEN q = gdiv(x,y), u, v;
1632 : /* e.g. x = y with t_PADIC components */
1633 132 : if (typ(q) != t_QUAD) { set_avma(av); return triv_cont_gcd(x,y); }
1634 120 : u = gel(q,2); v = gel(q,3);
1635 120 : if (gequal0(v))
1636 : {
1637 18 : if (typ(u)==t_INT) { set_avma(av); return gcopy(y); }
1638 12 : if (typ(u)==t_FRAC) return gc_upto(av, gdiv(y, gel(u,2)));
1639 6 : set_avma(av); return triv_cont_gcd(x,y);
1640 : }
1641 102 : if (typ(u)==t_INT && typ(v)==t_INT) { set_avma(av); return gcopy(y); }
1642 96 : q = ginv(q); u = gel(q,2); v = gel(q,3); set_avma(av);
1643 96 : if (typ(u)==t_INT && typ(v)==t_INT) return gcopy(x);
1644 90 : return triv_cont_gcd(y,x);
1645 : }
1646 :
1647 : GEN
1648 21741349 : ggcd(GEN x, GEN y)
1649 : {
1650 21741349 : long vx, vy, tx = typ(x), ty = typ(y);
1651 : pari_sp av;
1652 : GEN p1,z;
1653 :
1654 43482698 : if (is_noncalc_t(tx) || is_matvec_t(tx) ||
1655 43482698 : is_noncalc_t(ty) || is_matvec_t(ty)) pari_err_TYPE2("gcd",x,y);
1656 21741349 : if (tx>ty) { swap(x,y); lswap(tx,ty); }
1657 : /* tx <= ty */
1658 21741349 : z = gisexactzero(x); if (z) return zero_gcd2(y,z);
1659 19273274 : z = gisexactzero(y); if (z) return zero_gcd2(x,z);
1660 14446738 : if (is_const_t(tx))
1661 : {
1662 9066076 : if (ty == tx) switch(tx)
1663 : {
1664 5397368 : case t_INT:
1665 5397368 : return gcdii(x,y);
1666 :
1667 1839546 : case t_INTMOD: z=cgetg(3,t_INTMOD);
1668 1839546 : if (equalii(gel(x,1),gel(y,1)))
1669 1839541 : gel(z,1) = icopy(gel(x,1));
1670 : else
1671 5 : gel(z,1) = gcdii(gel(x,1),gel(y,1));
1672 1839546 : if (gequal1(gel(z,1))) gel(z,2) = gen_0;
1673 : else
1674 : {
1675 1839546 : av = avma; p1 = gcdii(gel(z,1),gel(x,2));
1676 1839546 : if (!equali1(p1))
1677 : {
1678 5 : p1 = gcdii(p1,gel(y,2));
1679 5 : if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
1680 5 : else p1 = gc_INT(av, p1);
1681 : }
1682 1839546 : gel(z,2) = p1;
1683 : }
1684 1839546 : return z;
1685 :
1686 221798 : case t_FRAC:
1687 221798 : return gcdqq(x,y);
1688 :
1689 4995 : case t_FFELT:
1690 4995 : if (!FF_samefield(x,y)) pari_err_OP("gcd",x,y);
1691 4995 : return FF_equal0(x) && FF_equal0(y)? FF_zero(y): FF_1(y);
1692 :
1693 18 : case t_COMPLEX:
1694 18 : if (c_is_rational(x) && c_is_rational(y)) return Qi_gcd(x,y);
1695 6 : return triv_cont_gcd(y,x);
1696 :
1697 12 : case t_PADIC:
1698 12 : if (!equalii(padic_p(x), padic_p(y))) return gen_1;
1699 6 : return powis(padic_p(x), minss(valp(x), valp(y)));
1700 :
1701 132 : case t_QUAD: return qgcd(x, y);
1702 :
1703 0 : default: return gen_1; /* t_REAL */
1704 : }
1705 1602207 : if (is_const_t(ty)) switch(tx)
1706 : {
1707 64525 : case t_INT:
1708 64525 : switch(ty)
1709 : {
1710 66 : case t_INTMOD: z = cgetg(3,t_INTMOD);
1711 66 : gel(z,1) = icopy(gel(y,1)); av = avma;
1712 66 : p1 = gcdii(gel(y,1),gel(y,2));
1713 66 : if (!equali1(p1)) {
1714 12 : p1 = gcdii(x,p1);
1715 12 : if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
1716 : else
1717 12 : p1 = gc_INT(av, p1);
1718 : }
1719 66 : gel(z,2) = p1; return z;
1720 :
1721 6413 : case t_REAL: return gen_1;
1722 :
1723 52882 : case t_FRAC:
1724 52882 : return gcdiq(x,y);
1725 :
1726 2682 : case t_COMPLEX:
1727 2682 : if (c_is_rational(y)) return Qi_gcd(x,y);
1728 2400 : return triv_cont_gcd(y,x);
1729 :
1730 60 : case t_FFELT:
1731 60 : if (!FF_equal0(y)) return FF_1(y);
1732 0 : return dvdii(x, gel(y,4))? FF_zero(y): FF_1(y);
1733 :
1734 2302 : case t_PADIC:
1735 2302 : return padic_gcd(x,y);
1736 :
1737 120 : case t_QUAD:
1738 120 : return triv_cont_gcd(y,x);
1739 0 : default:
1740 0 : pari_err_TYPE2("gcd",x,y);
1741 : }
1742 :
1743 12 : case t_REAL:
1744 12 : switch(ty)
1745 : {
1746 12 : case t_INTMOD:
1747 : case t_FFELT:
1748 12 : case t_PADIC: pari_err_TYPE2("gcd",x,y);
1749 0 : default: return gen_1;
1750 : }
1751 :
1752 48 : case t_INTMOD:
1753 48 : switch(ty)
1754 : {
1755 12 : case t_FRAC:
1756 12 : av = avma;
1757 12 : if (!equali1(gcdii(gel(x,1),gel(y,2)))) pari_err_OP("gcd",x,y);
1758 6 : set_avma(av); return ggcd(gel(y,1), x);
1759 :
1760 12 : case t_FFELT:
1761 : {
1762 12 : GEN p = gel(y,4);
1763 12 : if (!dvdii(gel(x,1), p)) pari_err_OP("gcd",x,y);
1764 6 : if (!FF_equal0(y)) return FF_1(y);
1765 0 : return dvdii(gel(x,2),p)? FF_zero(y): FF_1(y);
1766 : }
1767 :
1768 18 : case t_COMPLEX: case t_QUAD:
1769 18 : return triv_cont_gcd(y,x);
1770 :
1771 6 : case t_PADIC:
1772 6 : return padic_gcd(x,y);
1773 :
1774 0 : default: pari_err_TYPE2("gcd",x,y);
1775 : }
1776 :
1777 192 : case t_FRAC:
1778 192 : switch(ty)
1779 : {
1780 78 : case t_COMPLEX:
1781 78 : if (c_is_rational(y)) return Qi_gcd(x,y);
1782 : case t_QUAD:
1783 138 : return triv_cont_gcd(y,x);
1784 36 : case t_FFELT:
1785 : {
1786 36 : GEN p = gel(y,4);
1787 36 : if (dvdii(gel(x,2), p)) pari_err_OP("gcd",x,y);
1788 18 : if (!FF_equal0(y)) return FF_1(y);
1789 0 : return dvdii(gel(x,1),p)? FF_zero(y): FF_1(y);
1790 : }
1791 :
1792 12 : case t_PADIC:
1793 12 : return padic_gcd(x,y);
1794 :
1795 0 : default: pari_err_TYPE2("gcd",x,y);
1796 : }
1797 60 : case t_FFELT:
1798 60 : switch(ty)
1799 : {
1800 36 : case t_PADIC:
1801 : {
1802 36 : GEN p = padic_p(y);
1803 36 : long v = valp(y);
1804 36 : if (!equalii(p, gel(x,4)) || v < 0) pari_err_OP("gcd",x,y);
1805 12 : return (v && FF_equal0(x))? FF_zero(x): FF_1(x);
1806 : }
1807 24 : default: pari_err_TYPE2("gcd",x,y);
1808 : }
1809 :
1810 12 : case t_COMPLEX:
1811 12 : switch(ty)
1812 : {
1813 12 : case t_PADIC:
1814 12 : case t_QUAD: return triv_cont_gcd(x,y);
1815 0 : default: pari_err_TYPE2("gcd",x,y);
1816 : }
1817 :
1818 6 : case t_PADIC:
1819 6 : switch(ty)
1820 : {
1821 6 : case t_QUAD: return triv_cont_gcd(y,x);
1822 0 : default: pari_err_TYPE2("gcd",x,y);
1823 : }
1824 :
1825 0 : default: return gen_1; /* tx = t_REAL */
1826 : }
1827 1537352 : return cont_gcd(y,ty, x);
1828 : }
1829 :
1830 5380662 : if (tx == t_POLMOD)
1831 : {
1832 5077 : if (ty == t_POLMOD)
1833 : {
1834 5011 : GEN T = gel(x,1), Ty = gel(y,1);
1835 5011 : vx = varn(T); vy = varn(Ty);
1836 5011 : z = cgetg(3,t_POLMOD);
1837 5011 : if (vx == vy)
1838 4999 : T = RgX_equal(T,Ty)? RgX_copy(T): RgX_gcd(T, Ty);
1839 : else
1840 12 : T = RgX_copy(varncmp(vx,vy) < 0? T: Ty);
1841 5011 : gel(z,1) = T;
1842 5011 : if (degpol(T) <= 0) gel(z,2) = gen_0;
1843 : else
1844 : {
1845 5011 : GEN X = gel(x,2), Y = gel(y,2), d;
1846 5011 : av = avma; d = ggcd(content(X), content(Y));
1847 5011 : if (!gequal1(d)) { X = gdiv(X,d); Y = gdiv(Y,d); }
1848 5011 : gel(z,2) = gc_upto(av, gmul(d, gcd3(T, X, Y)));
1849 : }
1850 5011 : return z;
1851 : }
1852 66 : vx = varn(gel(x,1));
1853 66 : switch(ty)
1854 : {
1855 42 : case t_POL:
1856 42 : vy = varn(y);
1857 42 : if (varncmp(vy,vx) < 0) return cont_gcd_pol(y, x);
1858 24 : z = cgetg(3,t_POLMOD);
1859 24 : gel(z,1) = RgX_copy(gel(x,1)); av = avma;
1860 24 : gel(z,2) = gc_upto(av, gcd3(gel(x,1), gel(x,2), y));
1861 24 : return z;
1862 :
1863 24 : case t_RFRAC:
1864 24 : vy = varn(gel(y,2));
1865 24 : if (varncmp(vy,vx) < 0) return cont_gcd_rfrac(y, x);
1866 24 : av = avma;
1867 24 : if (degpol(ggcd(gel(x,1),gel(y,2)))) pari_err_OP("gcd",x,y);
1868 18 : set_avma(av); return gdiv(ggcd(gel(y,1),x), content(gel(y,2)));
1869 : }
1870 : }
1871 :
1872 5375585 : vx = gvar(x);
1873 5375585 : vy = gvar(y);
1874 5375585 : if (varncmp(vy, vx) < 0) return cont_gcd(y,ty, x);
1875 5370650 : if (varncmp(vy, vx) > 0) return cont_gcd(x,tx, y);
1876 :
1877 : /* vx = vy: same main variable */
1878 5363103 : switch(tx)
1879 : {
1880 5221241 : case t_POL:
1881 : switch(ty)
1882 : {
1883 : GEN cz, cx, cy;
1884 4486109 : case t_POL: return RgX_gcd(x,y);
1885 22 : case t_SER:
1886 22 : z = ggcd(content(x), content(y));
1887 22 : return monomialcopy(z, minss(valser(y),gval(x,vx)), vx);
1888 735110 : case t_RFRAC:
1889 735110 : av = avma;
1890 735110 : x = primitive_part(x, &cx);
1891 735110 : y = primitive_part(y, &cy);
1892 735110 : z = gred_rfrac_simple(ggcd(gel(y,1), x), gel(y,2));
1893 735110 : if (cx) cz = cy? ggcd(cx, cy): cx; else cz = cy? cy: NULL;
1894 735110 : if (cz) z = gmul(z, cz);
1895 735110 : return gc_upto(av, z);
1896 : }
1897 0 : break;
1898 :
1899 12 : case t_SER:
1900 12 : z = ggcd(content(x), content(y));
1901 : switch(ty)
1902 : {
1903 6 : case t_SER: return monomialcopy(z, minss(valser(x),valser(y)), vx);
1904 6 : case t_RFRAC: return monomialcopy(z, minss(valser(x),gval(y,vx)), vx);
1905 : }
1906 0 : break;
1907 :
1908 141850 : case t_RFRAC:
1909 : {
1910 141850 : GEN xd = gel(x,2), yd = gel(y,2);
1911 141850 : if (ty != t_RFRAC) pari_err_TYPE2("gcd",x,y);
1912 141850 : z = cgetg(3,t_RFRAC); av = avma;
1913 141850 : gel(z,2) = gc_upto(av, RgX_mul(xd, RgX_div(yd, RgX_gcd(xd, yd))));
1914 141850 : gel(z,1) = ggcd(gel(x,1), gel(y,1)); return z;
1915 : }
1916 : }
1917 0 : pari_err_TYPE2("gcd",x,y);
1918 : return NULL; /* LCOV_EXCL_LINE */
1919 : }
1920 : GEN
1921 3717122 : ggcd0(GEN x, GEN y) { return y? ggcd(x,y): content(x); }
1922 :
1923 : static GEN
1924 90 : fix_lcm(GEN x)
1925 : {
1926 : GEN t;
1927 90 : switch(typ(x))
1928 : {
1929 0 : case t_INT:
1930 0 : x = absi_shallow(x); break;
1931 84 : case t_POL:
1932 84 : if (lg(x) <= 2) break;
1933 84 : t = leading_coeff(x);
1934 84 : if (typ(t) == t_INT && signe(t) < 0) x = gneg(x);
1935 : }
1936 90 : return x;
1937 : }
1938 : GEN
1939 2081 : glcm0(GEN x, GEN y)
1940 : {
1941 2081 : if (!y) return fix_lcm(gassoc_proto(glcm,x,y));
1942 2039 : return glcm(x,y);
1943 : }
1944 :
1945 : static GEN
1946 6985 : _lcmii(void*E, GEN x, GEN y)
1947 6985 : { (void) E; return lcmii(x,y); }
1948 :
1949 : GEN
1950 3067 : ZV_lcm(GEN x) { return gen_product(x, NULL, _lcmii); }
1951 :
1952 : GEN
1953 2423 : glcm(GEN x, GEN y)
1954 : {
1955 : pari_sp av;
1956 : GEN z;
1957 2423 : if (typ(x)==t_INT && typ(y)==t_INT) return lcmii(x,y);
1958 60 : av = avma; z = ggcd(x,y);
1959 60 : if (!gequal1(z))
1960 : {
1961 54 : if (gequal0(z)) { set_avma(av); return gmul(x,y); }
1962 42 : y = gdiv(y,z);
1963 : }
1964 48 : return gc_upto(av, fix_lcm(gmul(x,y)));
1965 : }
1966 :
1967 : /* x + r ~ x ? Assume x,r are t_POL, deg(r) <= deg(x) */
1968 : static int
1969 0 : pol_approx0(GEN r, GEN x, int exact)
1970 : {
1971 : long i, l;
1972 0 : if (exact) return !signe(r);
1973 0 : l = minss(lg(x), lg(r));
1974 0 : for (i = 2; i < l; i++)
1975 0 : if (!cx_approx0(gel(r,i), gel(x,i))) return 0;
1976 0 : return 1;
1977 : }
1978 :
1979 : GEN
1980 0 : RgX_gcd_simple(GEN x, GEN y)
1981 : {
1982 0 : pari_sp av1, av = avma;
1983 0 : GEN r, yorig = y;
1984 0 : int exact = !(isinexactreal(x) || isinexactreal(y));
1985 :
1986 : for(;;)
1987 : {
1988 0 : av1 = avma; r = RgX_rem(x,y);
1989 0 : if (pol_approx0(r, x, exact))
1990 : {
1991 0 : set_avma(av1);
1992 0 : if (y == yorig) return RgX_copy(y);
1993 0 : y = normalizepol_approx(y, lg(y));
1994 0 : if (lg(y) == 3) { set_avma(av); return pol_1(varn(x)); }
1995 0 : return gc_upto(av,y);
1996 : }
1997 0 : x = y; y = r;
1998 0 : if (gc_needed(av,1)) {
1999 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd_simple");
2000 0 : (void)gc_all(av,2, &x,&y);
2001 : }
2002 : }
2003 : }
2004 : GEN
2005 0 : RgX_extgcd_simple(GEN a, GEN b, GEN *pu, GEN *pv)
2006 : {
2007 0 : pari_sp av = avma;
2008 : GEN q, r, d, d1, u, v, v1;
2009 0 : int exact = !(isinexactreal(a) || isinexactreal(b));
2010 :
2011 0 : d = a; d1 = b; v = gen_0; v1 = gen_1;
2012 : for(;;)
2013 : {
2014 0 : if (pol_approx0(d1, a, exact)) break;
2015 0 : q = poldivrem(d,d1, &r);
2016 0 : v = gsub(v, gmul(q,v1));
2017 0 : u=v; v=v1; v1=u;
2018 0 : u=r; d=d1; d1=u;
2019 : }
2020 0 : u = gsub(d, gmul(b,v));
2021 0 : u = RgX_div(u,a);
2022 :
2023 0 : (void)gc_all(av, 3, &u,&v,&d);
2024 0 : *pu = u;
2025 0 : *pv = v; return d;
2026 : }
2027 :
2028 : GEN
2029 78 : ghalfgcd(GEN x, GEN y)
2030 : {
2031 78 : long tx = typ(x), ty = typ(y);
2032 78 : if (tx==t_INT && ty==t_INT) return halfgcdii(x, y);
2033 54 : if (tx==t_POL && ty==t_POL && varn(x)==varn(y))
2034 : {
2035 54 : pari_sp av = avma;
2036 54 : GEN a, b, M = RgX_halfgcd_all(x, y, &a, &b);
2037 54 : return gc_GEN(av, mkvec2(M, mkcol2(a,b)));
2038 : }
2039 0 : pari_err_OP("halfgcd", x, y);
2040 : return NULL; /* LCOV_EXCL_LINE */
2041 : }
2042 :
2043 : /*******************************************************************/
2044 : /* */
2045 : /* CONTENT / PRIMITIVE PART */
2046 : /* */
2047 : /*******************************************************************/
2048 :
2049 : GEN
2050 61135641 : content(GEN x)
2051 : {
2052 61135641 : long lx, i, t, tx = typ(x);
2053 61135641 : pari_sp av = avma;
2054 : GEN c;
2055 :
2056 61135641 : if (is_scalar_t(tx)) return zero_gcd(x);
2057 61121567 : switch(tx)
2058 : {
2059 742847 : case t_RFRAC:
2060 : {
2061 742847 : GEN n = gel(x,1), d = gel(x,2);
2062 : /* -- varncmp(vn, vd) < 0 can't happen
2063 : * -- if n is POLMOD, its main variable (in the sense of gvar2)
2064 : * has lower priority than denominator */
2065 742847 : if (typ(n) == t_POLMOD || varncmp(gvar(n), varn(d)) > 0)
2066 707341 : n = isinexact(n)? zero_gcd(n): gcopy(n);
2067 : else
2068 35506 : n = content(n);
2069 742847 : return gc_upto(av, gdiv(n, content(d)));
2070 : }
2071 :
2072 1030101 : case t_VEC: case t_COL:
2073 1030101 : lx = lg(x); if (lx==1) return gen_0;
2074 1030095 : break;
2075 :
2076 18 : case t_MAT:
2077 : {
2078 : long hx, j;
2079 18 : lx = lg(x);
2080 18 : if (lx == 1) return gen_0;
2081 12 : hx = lgcols(x);
2082 12 : if (hx == 1) return gen_0;
2083 6 : if (lx == 2) { x = gel(x,1); lx = lg(x); break; }
2084 6 : if (hx == 2) { x = row_i(x, 1, 1, lx-1); break; }
2085 6 : c = content(gel(x,1));
2086 12 : for (j=2; j<lx; j++)
2087 18 : for (i=1; i<hx; i++) c = ggcd(c,gcoeff(x,i,j));
2088 6 : if (typ(c) == t_INTMOD || isinexact(c)) return gc_const(av, gen_1);
2089 6 : return gc_upto(av,c);
2090 : }
2091 :
2092 59348443 : case t_POL: case t_SER:
2093 59348443 : lx = lg(x); if (lx == 2) return gen_0;
2094 59328565 : break;
2095 18 : case t_VECSMALL: return utoi(zv_content(x));
2096 140 : case t_QFB:
2097 140 : lx = 4; break;
2098 :
2099 0 : default: pari_err_TYPE("content",x);
2100 : return NULL; /* LCOV_EXCL_LINE */
2101 : }
2102 180275705 : for (i=lontyp[tx]; i<lx; i++)
2103 129209375 : if (typ(gel(x,i)) != t_INT) break;
2104 60358800 : lx--; c = gel(x,lx);
2105 60358800 : t = typ(c); if (is_matvec_t(t)) c = content(c);
2106 60358800 : if (i > lx)
2107 : { /* integer coeffs */
2108 53474738 : while (lx-- > lontyp[tx])
2109 : {
2110 51527828 : c = gcdii(c, gel(x,lx));
2111 51527828 : if (equali1(c)) return gc_const(av, gen_1);
2112 : }
2113 : }
2114 : else
2115 : {
2116 9292470 : if (isinexact(c)) c = zero_gcd(c);
2117 24870911 : while (lx-- > lontyp[tx])
2118 : {
2119 15578441 : GEN d = gel(x,lx);
2120 15578441 : t = typ(d); if (is_matvec_t(t)) d = content(d);
2121 15578441 : c = ggcd(c, d);
2122 : }
2123 9292470 : if (isinexact(c)) return gc_const(av, gen_1);
2124 : }
2125 11239380 : switch(typ(c))
2126 : {
2127 1951117 : case t_INT:
2128 1951117 : c = absi_shallow(c); break;
2129 0 : case t_VEC: case t_COL: case t_MAT:
2130 0 : pari_err_TYPE("content",x);
2131 : }
2132 :
2133 11239380 : return av==avma? gcopy(c): gc_upto(av,c);
2134 : }
2135 :
2136 : GEN
2137 2478314 : primitive_part(GEN x, GEN *ptc)
2138 : {
2139 2478314 : pari_sp av = avma;
2140 2478314 : GEN c = content(x);
2141 2478314 : if (gequal1(c)) { set_avma(av); c = NULL; }
2142 173153 : else if (!gequal0(c)) x = gdiv(x,c);
2143 2478314 : if (ptc) *ptc = c;
2144 2478314 : return x;
2145 : }
2146 : GEN
2147 136 : primpart(GEN x) { return primitive_part(x, NULL); }
2148 :
2149 : static GEN
2150 151379547 : Q_content_v(GEN x, long imin, long l)
2151 : {
2152 151379547 : pari_sp av = avma;
2153 151379547 : long i = l-1;
2154 151379547 : GEN d = Q_content_safe(gel(x,i));
2155 151379547 : if (!d) return NULL;
2156 997515863 : for (i--; i >= imin; i--)
2157 : {
2158 846136732 : GEN c = Q_content_safe(gel(x,i));
2159 846136732 : if (!c) return NULL;
2160 846136702 : d = Q_gcd(d, c);
2161 846136702 : if (gc_needed(av,1)) d = gc_upto(av, d);
2162 : }
2163 151379131 : return gc_upto(av, d);
2164 : }
2165 : /* As content(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
2166 : * of Q(x2,...,xn)[x1] */
2167 : GEN
2168 1073322017 : Q_content_safe(GEN x)
2169 : {
2170 : long l;
2171 1073322017 : switch(typ(x))
2172 : {
2173 888951165 : case t_INT: return absi(x);
2174 32025306 : case t_FRAC: return absfrac(x);
2175 108589720 : case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
2176 108589720 : l = lg(x); return l==1? gen_1: Q_content_v(x, 1, l);
2177 43728287 : case t_POL:
2178 43728287 : l = lg(x); return l==2? gen_0: Q_content_v(x, 2, l);
2179 27215 : case t_POLMOD: return Q_content_safe(gel(x,2));
2180 18 : case t_RFRAC:
2181 : {
2182 : GEN a, b;
2183 18 : a = Q_content_safe(gel(x,1)); if (!a) return NULL;
2184 18 : b = Q_content_safe(gel(x,2)); if (!b) return NULL;
2185 18 : return gdiv(a, b);
2186 : }
2187 : }
2188 306 : return NULL;
2189 : }
2190 : GEN
2191 1584217 : Q_content(GEN x)
2192 : {
2193 1584217 : GEN c = Q_content_safe(x);
2194 1584217 : if (!c) pari_err_TYPE("Q_content",x);
2195 1584217 : return c;
2196 : }
2197 :
2198 : GEN
2199 11268 : ZX_content(GEN x)
2200 : {
2201 11268 : long i, l = lg(x);
2202 : GEN d;
2203 : pari_sp av;
2204 :
2205 11268 : if (l == 2) return gen_0;
2206 11268 : d = gel(x,2);
2207 11268 : if (l == 3) return absi(d);
2208 7860 : av = avma;
2209 16254 : for (i=3; !is_pm1(d) && i<l; i++) d = gcdii(d, gel(x,i));
2210 7860 : if (signe(d) < 0) d = negi(d);
2211 7860 : return gc_INT(av, d);
2212 : }
2213 :
2214 : static GEN
2215 2003818 : Z_content_v(GEN x, long i, long l)
2216 : {
2217 2003818 : pari_sp av = avma;
2218 2003818 : GEN d = Z_content(gel(x,i));
2219 2003818 : if (!d) return NULL;
2220 5133148 : for (i++; i<l; i++)
2221 : {
2222 4661803 : GEN c = Z_content(gel(x,i));
2223 4661803 : if (!c) return NULL;
2224 4139087 : d = gcdii(d, c); if (equali1(d)) return NULL;
2225 3393822 : if ((i & 255) == 0) d = gc_INT(av, d);
2226 : }
2227 471345 : return gc_INT(av, d);
2228 : }
2229 : /* return NULL for 1 */
2230 : GEN
2231 8624763 : Z_content(GEN x)
2232 : {
2233 : long l;
2234 8624763 : switch(typ(x))
2235 : {
2236 6603449 : case t_INT:
2237 6603449 : if (is_pm1(x)) return NULL;
2238 5819625 : return absi(x);
2239 1965856 : case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
2240 1965856 : l = lg(x); return l==1? NULL: Z_content_v(x, 1, l);
2241 55458 : case t_POL:
2242 55458 : l = lg(x); return l==2? gen_0: Z_content_v(x, 2, l);
2243 0 : case t_POLMOD: return Z_content(gel(x,2));
2244 : }
2245 0 : pari_err_TYPE("Z_content", x);
2246 : return NULL; /* LCOV_EXCL_LINE */
2247 : }
2248 :
2249 : static GEN
2250 47035842 : Q_denom_v(GEN x, long i, long l)
2251 : {
2252 47035842 : pari_sp av = avma;
2253 47035842 : GEN d = Q_denom_safe(gel(x,i));
2254 47035842 : if (!d) return NULL;
2255 160387572 : for (i++; i<l; i++)
2256 : {
2257 113351735 : GEN D = Q_denom_safe(gel(x,i));
2258 113351735 : if (!D) return NULL;
2259 113351735 : if (D != gen_1) d = lcmii(d, D);
2260 113351735 : if ((i & 255) == 0) d = gc_INT(av, d);
2261 : }
2262 47035837 : return gc_INT(av, d);
2263 : }
2264 : /* NOT MEMORY CLEAN (because of t_FRAC).
2265 : * As denom(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
2266 : * of Q(x2,...,xn)[x1] */
2267 : GEN
2268 212437602 : Q_denom_safe(GEN x)
2269 : {
2270 : long l;
2271 212437602 : switch(typ(x))
2272 : {
2273 135567415 : case t_INT: return gen_1;
2274 24 : case t_PADIC: l = valp(x); return l < 0? powiu(padic_p(x), -l): gen_1;
2275 29596011 : case t_FRAC: return gel(x,2);
2276 410 : case t_QUAD: return Q_denom_v(x, 2, 4);
2277 36296002 : case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
2278 36296002 : l = lg(x); return l==1? gen_1: Q_denom_v(x, 1, l);
2279 10885625 : case t_POL: case t_SER:
2280 10885625 : l = lg(x); return l==2? gen_1: Q_denom_v(x, 2, l);
2281 85078 : case t_POLMOD: return Q_denom(gel(x,2));
2282 6972 : case t_RFRAC:
2283 : {
2284 : GEN a, b;
2285 6972 : a = Q_content(gel(x,1)); if (!a) return NULL;
2286 6972 : b = Q_content(gel(x,2)); if (!b) return NULL;
2287 6972 : return Q_denom(gdiv(a, b));
2288 : }
2289 : }
2290 65 : return NULL;
2291 : }
2292 : GEN
2293 3485848 : Q_denom(GEN x)
2294 : {
2295 3485848 : GEN d = Q_denom_safe(x);
2296 3485848 : if (!d) pari_err_TYPE("Q_denom",x);
2297 3485848 : return d;
2298 : }
2299 :
2300 : GEN
2301 48564057 : Q_remove_denom(GEN x, GEN *ptd)
2302 : {
2303 48564057 : GEN d = Q_denom_safe(x);
2304 48564057 : if (d) { if (d == gen_1) d = NULL; else x = Q_muli_to_int(x,d); }
2305 48564057 : if (ptd) *ptd = d;
2306 48564057 : return x;
2307 : }
2308 :
2309 : /* return y = x * d, assuming x rational, and d,y integral */
2310 : GEN
2311 120240579 : Q_muli_to_int(GEN x, GEN d)
2312 : {
2313 : GEN y, xn, xd;
2314 : pari_sp av;
2315 :
2316 120240579 : if (typ(d) != t_INT) pari_err_TYPE("Q_muli_to_int",d);
2317 120240579 : switch (typ(x))
2318 : {
2319 38262356 : case t_INT:
2320 38262356 : return mulii(x,d);
2321 :
2322 55745519 : case t_FRAC:
2323 55745519 : xn = gel(x,1);
2324 55745519 : xd = gel(x,2); av = avma;
2325 55745519 : y = mulii(xn, diviiexact(d, xd));
2326 55745519 : return gc_INT(av, y);
2327 36 : case t_COMPLEX:
2328 36 : y = cgetg(3,t_COMPLEX);
2329 36 : gel(y,1) = Q_muli_to_int(gel(x,1),d);
2330 36 : gel(y,2) = Q_muli_to_int(gel(x,2),d);
2331 36 : return y;
2332 12 : case t_PADIC:
2333 12 : y = gcopy(x); if (!isint1(d)) setvalp(y, 0);
2334 12 : return y;
2335 142 : case t_QUAD:
2336 142 : y = cgetg(4,t_QUAD);
2337 142 : gel(y,1) = ZX_copy(gel(x,1));
2338 142 : gel(y,2) = Q_muli_to_int(gel(x,2),d);
2339 142 : gel(y,3) = Q_muli_to_int(gel(x,3),d); return y;
2340 :
2341 17186543 : case t_VEC:
2342 : case t_COL:
2343 86532585 : case t_MAT: pari_APPLY_same(Q_muli_to_int(gel(x,i), d));
2344 40886188 : case t_POL: pari_APPLY_pol_normalized(Q_muli_to_int(gel(x,i), d));
2345 18 : case t_SER: pari_APPLY_ser_normalized(Q_muli_to_int(gel(x,i), d));
2346 :
2347 43837 : case t_POLMOD:
2348 43837 : retmkpolmod(Q_muli_to_int(gel(x,2), d), RgX_copy(gel(x,1)));
2349 18 : case t_RFRAC:
2350 18 : return gmul(x, d);
2351 : }
2352 0 : pari_err_TYPE("Q_muli_to_int",x);
2353 : return NULL; /* LCOV_EXCL_LINE */
2354 : }
2355 :
2356 : static void
2357 24691705 : rescale_init(GEN c, int *exact, long *emin, GEN *D)
2358 : {
2359 : long e, i;
2360 24691705 : switch(typ(c))
2361 : {
2362 16725814 : case t_REAL:
2363 16725814 : *exact = 0;
2364 16725814 : if (!signe(c)) return;
2365 16316748 : e = expo(c) + 1 - bit_prec(c);
2366 18452462 : for (i = lg(c)-1; i > 2; i--, e += BITS_IN_LONG)
2367 14092961 : if (c[i]) break;
2368 16316748 : e += vals(c[i]); break; /* e[2] != 0 */
2369 7961956 : case t_INT:
2370 7961956 : if (!signe(c)) return;
2371 1184646 : e = expi(c);
2372 1184646 : break;
2373 3887 : case t_FRAC:
2374 3887 : e = expi(gel(c,1)) - expi(gel(c,2));
2375 3887 : if (*exact) *D = lcmii(*D, gel(c,2));
2376 3887 : break;
2377 48 : default:
2378 48 : pari_err_TYPE("rescale_to_int",c);
2379 : return; /* LCOV_EXCL_LINE */
2380 : }
2381 17505281 : if (e < *emin) *emin = e;
2382 : }
2383 : GEN
2384 3897876 : RgM_rescale_to_int(GEN x)
2385 : {
2386 3897876 : long lx = lg(x), i,j, hx, emin;
2387 : GEN D;
2388 : int exact;
2389 :
2390 3897876 : if (lx == 1) return cgetg(1,t_MAT);
2391 3897876 : hx = lgcols(x);
2392 3897876 : exact = 1;
2393 3897876 : emin = HIGHEXPOBIT;
2394 3897876 : D = gen_1;
2395 12915892 : for (j = 1; j < lx; j++)
2396 33545808 : for (i = 1; i < hx; i++) rescale_init(gcoeff(x,i,j), &exact, &emin, &D);
2397 3897828 : if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
2398 3897747 : return grndtoi(gmul2n(x, -emin), NULL);
2399 : }
2400 : GEN
2401 32087 : RgX_rescale_to_int(GEN x)
2402 : {
2403 32087 : long lx = lg(x), i, emin;
2404 : GEN D;
2405 : int exact;
2406 32087 : if (lx == 2) return gcopy(x); /* rare */
2407 32087 : exact = 1;
2408 32087 : emin = HIGHEXPOBIT;
2409 32087 : D = gen_1;
2410 196000 : for (i = 2; i < lx; i++) rescale_init(gel(x,i), &exact, &emin, &D);
2411 32087 : if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
2412 31123 : return grndtoi(gmul2n(x, -emin), NULL);
2413 : }
2414 :
2415 : /* return x * n/d. x: rational; d,n,result: integral; d,n coprime */
2416 : static GEN
2417 9872580 : Q_divmuli_to_int(GEN x, GEN d, GEN n)
2418 : {
2419 : GEN y, xn, xd;
2420 : pari_sp av;
2421 :
2422 9872580 : switch(typ(x))
2423 : {
2424 2484046 : case t_INT:
2425 2484046 : av = avma; y = diviiexact(x,d);
2426 2484046 : return gc_INT(av, mulii(y,n));
2427 :
2428 4955391 : case t_FRAC:
2429 4955391 : xn = gel(x,1);
2430 4955391 : xd = gel(x,2); av = avma;
2431 4955391 : y = mulii(diviiexact(xn, d), diviiexact(n, xd));
2432 4955391 : return gc_INT(av, y);
2433 :
2434 386522 : case t_VEC:
2435 : case t_COL:
2436 3465983 : case t_MAT: pari_APPLY_same(Q_divmuli_to_int(gel(x,i), d,n));
2437 6721412 : case t_POL: pari_APPLY_pol_normalized(Q_divmuli_to_int(gel(x,i), d,n));
2438 :
2439 0 : case t_RFRAC:
2440 0 : av = avma;
2441 0 : return gc_upto(av, gmul(x,mkfrac(n,d)));
2442 :
2443 0 : case t_POLMOD:
2444 0 : retmkpolmod(Q_divmuli_to_int(gel(x,2), d,n), RgX_copy(gel(x,1)));
2445 : }
2446 0 : pari_err_TYPE("Q_divmuli_to_int",x);
2447 : return NULL; /* LCOV_EXCL_LINE */
2448 : }
2449 :
2450 : /* return x / d. x: rational; d,result: integral. */
2451 : static GEN
2452 146477484 : Q_divi_to_int(GEN x, GEN d)
2453 : {
2454 146477484 : switch(typ(x))
2455 : {
2456 123196649 : case t_INT:
2457 123196649 : return diviiexact(x,d);
2458 :
2459 17605775 : case t_VEC:
2460 : case t_COL:
2461 136342561 : case t_MAT: pari_APPLY_same(Q_divi_to_int(gel(x,i), d));
2462 21182533 : case t_POL: pari_APPLY_pol_normalized(Q_divi_to_int(gel(x,i), d));
2463 :
2464 0 : case t_RFRAC:
2465 0 : return gdiv(x,d);
2466 :
2467 5082 : case t_POLMOD:
2468 5082 : retmkpolmod(Q_divi_to_int(gel(x,2), d), RgX_copy(gel(x,1)));
2469 : }
2470 0 : pari_err_TYPE("Q_divi_to_int",x);
2471 : return NULL; /* LCOV_EXCL_LINE */
2472 : }
2473 : /* c t_FRAC */
2474 : static GEN
2475 8665349 : Q_divq_to_int(GEN x, GEN c)
2476 : {
2477 8665349 : GEN n = gel(c,1), d = gel(c,2);
2478 8665349 : if (is_pm1(n)) {
2479 6547021 : GEN y = Q_muli_to_int(x,d);
2480 6547021 : if (signe(n) < 0) y = gneg(y);
2481 6547021 : return y;
2482 : }
2483 2118328 : return Q_divmuli_to_int(x, n,d);
2484 : }
2485 :
2486 : /* return y = x / c, assuming x,c rational, and y integral */
2487 : GEN
2488 433539 : Q_div_to_int(GEN x, GEN c)
2489 : {
2490 433539 : switch(typ(c))
2491 : {
2492 432829 : case t_INT: return Q_divi_to_int(x, c);
2493 710 : case t_FRAC: return Q_divq_to_int(x, c);
2494 : }
2495 0 : pari_err_TYPE("Q_div_to_int",c);
2496 : return NULL; /* LCOV_EXCL_LINE */
2497 : }
2498 : /* return y = x * c, assuming x,c rational, and y integral */
2499 : GEN
2500 0 : Q_mul_to_int(GEN x, GEN c)
2501 : {
2502 : GEN d, n;
2503 0 : switch(typ(c))
2504 : {
2505 0 : case t_INT: return Q_muli_to_int(x, c);
2506 0 : case t_FRAC:
2507 0 : n = gel(c,1);
2508 0 : d = gel(c,2);
2509 0 : return Q_divmuli_to_int(x, d,n);
2510 : }
2511 0 : pari_err_TYPE("Q_mul_to_int",c);
2512 : return NULL; /* LCOV_EXCL_LINE */
2513 : }
2514 :
2515 : GEN
2516 74171730 : Q_primitive_part(GEN x, GEN *ptc)
2517 : {
2518 74171730 : pari_sp av = avma;
2519 74171730 : GEN c = Q_content_safe(x);
2520 74171730 : if (c)
2521 : {
2522 74171484 : if (typ(c) == t_INT)
2523 : {
2524 65506845 : if (equali1(c)) { set_avma(av); c = NULL; }
2525 12110898 : else if (signe(c)) x = Q_divi_to_int(x, c);
2526 : }
2527 8664639 : else x = Q_divq_to_int(x, c);
2528 : }
2529 74171730 : if (ptc) *ptc = c;
2530 74171730 : return x;
2531 : }
2532 : GEN
2533 8334535 : Q_primpart(GEN x) { return Q_primitive_part(x, NULL); }
2534 : GEN
2535 99520 : vec_Q_primpart(GEN x)
2536 568089 : { pari_APPLY_same(Q_primpart(gel(x,i))) }
2537 : GEN
2538 53166 : row_Q_primpart(GEN M)
2539 53166 : { return shallowtrans(vec_Q_primpart(shallowtrans(M))); }
2540 :
2541 : /*******************************************************************/
2542 : /* */
2543 : /* SUBRESULTANT */
2544 : /* */
2545 : /*******************************************************************/
2546 : /* for internal use */
2547 : GEN
2548 19954043 : gdivexact(GEN x, GEN y)
2549 : {
2550 : long i,lx;
2551 : GEN z;
2552 19954043 : if (gequal1(y)) return x;
2553 19952254 : if (typ(y) == t_POLMOD) return gmul(x, ginv(y));
2554 19952184 : switch(typ(x))
2555 : {
2556 16632235 : case t_INT:
2557 16632235 : if (typ(y)==t_INT) return diviiexact(x,y);
2558 30 : if (!signe(x)) return gen_0;
2559 0 : break;
2560 6476 : case t_INTMOD:
2561 : case t_FFELT:
2562 6476 : case t_POLMOD: return gmul(x,ginv(y));
2563 3309219 : case t_POL:
2564 3309219 : switch(typ(y))
2565 : {
2566 612 : case t_INTMOD:
2567 : case t_FFELT:
2568 612 : case t_POLMOD: return gmul(x,ginv(y));
2569 137272 : case t_POL: { /* not stack-clean */
2570 : long v;
2571 137272 : if (varn(x)!=varn(y)) break;
2572 136582 : v = RgX_valrem(y,&y);
2573 136582 : if (v) x = RgX_shift_shallow(x,-v);
2574 136582 : if (!degpol(y)) { y = gel(y,2); break; }
2575 134659 : return RgX_div(x,y);
2576 : }
2577 0 : case t_RFRAC:
2578 0 : if (varn(gel(y,2)) != varn(x)) break;
2579 0 : return gdiv(x, y);
2580 : }
2581 3173948 : return RgX_Rg_divexact(x, y);
2582 4254 : case t_VEC: case t_COL: case t_MAT:
2583 4254 : lx = lg(x); z = new_chunk(lx);
2584 46602 : for (i=1; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
2585 4254 : z[0] = x[0]; return z;
2586 : }
2587 0 : if (DEBUGLEVEL) pari_warn(warner,"missing case in gdivexact");
2588 0 : return gdiv(x,y);
2589 : }
2590 :
2591 : static GEN
2592 1131055 : init_resultant(GEN x, GEN y)
2593 : {
2594 1131055 : long tx = typ(x), ty = typ(y), vx, vy;
2595 1131055 : if (is_scalar_t(tx) || is_scalar_t(ty))
2596 : {
2597 10 : if (gequal0(x) || gequal0(y)) return gmul(x,y); /* keep type info */
2598 10 : if (tx==t_POL) return gpowgs(y, degpol(x));
2599 0 : if (ty==t_POL) return gpowgs(x, degpol(y));
2600 0 : return gen_1;
2601 : }
2602 1131045 : if (tx!=t_POL) pari_err_TYPE("resultant",x);
2603 1131045 : if (ty!=t_POL) pari_err_TYPE("resultant",y);
2604 1131045 : if (!signe(x) || !signe(y)) return gmul(Rg_get_0(x),Rg_get_0(y)); /*type*/
2605 1131040 : vx = varn(x);
2606 1131040 : vy = varn(y); if (vx == vy) return NULL;
2607 5 : return (varncmp(vx,vy) < 0)? gpowgs(y,degpol(x)): gpowgs(x,degpol(y));
2608 : }
2609 :
2610 : /* x an RgX, y a scalar */
2611 : static GEN
2612 5 : scalar_res(GEN x, GEN y, GEN *U, GEN *V)
2613 : {
2614 5 : *V = gpowgs(y,degpol(x)-1);
2615 5 : *U = gen_0; return gmul(y, *V);
2616 : }
2617 :
2618 : /* return 0 if the subresultant chain can be interrupted.
2619 : * Set u = NULL if the resultant is 0. */
2620 : static int
2621 10059 : subres_step(GEN *u, GEN *v, GEN *g, GEN *h, GEN *uze, GEN *um1, long *signh)
2622 : {
2623 10059 : GEN u0, c, r, q = RgX_pseudodivrem(*u,*v, &r);
2624 : long du, dv, dr, degq;
2625 :
2626 10059 : if (gequal0(leading_coeff(r))) r = RgX_renormalize(r);
2627 10059 : dr = lg(r); if (!signe(r)) { *u = NULL; return 0; }
2628 9849 : du = degpol(*u);
2629 9849 : dv = degpol(*v);
2630 9849 : degq = du - dv;
2631 9849 : if (*um1 == gen_1)
2632 5475 : u0 = gpowgs(gel(*v,dv+2),degq+1);
2633 4374 : else if (*um1 == gen_0)
2634 1976 : u0 = gen_0;
2635 : else /* except in those 2 cases, um1 is an RgX */
2636 2398 : u0 = RgX_Rg_mul(*um1, gpowgs(gel(*v,dv+2),degq+1));
2637 :
2638 9849 : if (*uze == gen_0) /* except in that case, uze is an RgX */
2639 5475 : u0 = scalarpol(u0, varn(*u)); /* now an RgX */
2640 : else
2641 4374 : u0 = gsub(u0, gmul(q,*uze));
2642 :
2643 9849 : *um1 = *uze;
2644 9849 : *uze = u0; /* uze <- lead(v)^(degq + 1) * um1 - q * uze */
2645 :
2646 9849 : *u = *v; c = *g; *g = leading_coeff(*u);
2647 9849 : switch(degq)
2648 : {
2649 1421 : case 0: break;
2650 6877 : case 1:
2651 6877 : c = gmul(*h,c); *h = *g; break;
2652 1551 : default:
2653 1551 : c = gmul(gpowgs(*h,degq), c);
2654 1551 : *h = gdivexact(gpowgs(*g,degq), gpowgs(*h,degq-1));
2655 : }
2656 9849 : if (typ(c) == t_POLMOD)
2657 : {
2658 761 : c = ginv(c);
2659 761 : *v = RgX_Rg_mul(r,c);
2660 761 : *uze = RgX_Rg_mul(*uze,c);
2661 : }
2662 : else
2663 : {
2664 9088 : *v = RgX_Rg_divexact(r,c);
2665 9088 : *uze= RgX_Rg_divexact(*uze,c);
2666 : }
2667 9849 : if (both_odd(du, dv)) *signh = -*signh;
2668 9849 : return (dr > 3);
2669 : }
2670 :
2671 : /* compute U, V s.t Ux + Vy = resultant(x,y) */
2672 : static GEN
2673 1994 : subresext_i(GEN x, GEN y, GEN *U, GEN *V)
2674 : {
2675 : pari_sp av, av2;
2676 1994 : long dx, dy, du, signh, tx = typ(x), ty = typ(y);
2677 : GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze;
2678 :
2679 1994 : if (!is_extscalar_t(tx)) pari_err_TYPE("subresext",x);
2680 1994 : if (!is_extscalar_t(ty)) pari_err_TYPE("subresext",y);
2681 1994 : if (gequal0(x) || gequal0(y)) { *U = *V = gen_0; return gen_0; }
2682 1994 : if (tx != t_POL) {
2683 5 : if (ty != t_POL) { *U = ginv(x); *V = gen_0; return gen_1; }
2684 5 : return scalar_res(y,x,V,U);
2685 : }
2686 1989 : if (ty != t_POL) return scalar_res(x,y,U,V);
2687 1989 : if (varn(x) != varn(y))
2688 0 : return varncmp(varn(x), varn(y)) < 0? scalar_res(x,y,U,V)
2689 0 : : scalar_res(y,x,V,U);
2690 1989 : if (gequal0(leading_coeff(x))) x = RgX_renormalize(x);
2691 1989 : if (gequal0(leading_coeff(y))) y = RgX_renormalize(y);
2692 1989 : dx = degpol(x);
2693 1989 : dy = degpol(y);
2694 1989 : signh = 1;
2695 1989 : if (dx < dy)
2696 : {
2697 713 : pswap(U,V); lswap(dx,dy); swap(x,y);
2698 713 : if (both_odd(dx, dy)) signh = -signh;
2699 : }
2700 1989 : if (dy == 0)
2701 : {
2702 0 : *V = gpowgs(gel(y,2),dx-1);
2703 0 : *U = gen_0; return gmul(*V,gel(y,2));
2704 : }
2705 1989 : av = avma;
2706 1989 : u = x = primitive_part(x, &cu);
2707 1989 : v = y = primitive_part(y, &cv);
2708 1989 : g = h = gen_1; av2 = avma;
2709 1989 : um1 = gen_1; uze = gen_0;
2710 : for(;;)
2711 : {
2712 5973 : if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
2713 3984 : if (gc_needed(av2,1))
2714 : {
2715 0 : if(DEBUGMEM>1) pari_warn(warnmem,"subresext, dr = %ld", degpol(v));
2716 0 : (void)gc_all(av2,6, &u,&v,&g,&h,&uze,&um1);
2717 : }
2718 : }
2719 : /* uze an RgX */
2720 1989 : if (!u) { *U = *V = gen_0; return gc_const(av, gen_0); }
2721 1983 : z = gel(v,2); du = degpol(u);
2722 1983 : if (du > 1)
2723 : { /* z = gdivexact(gpowgs(z,du), gpowgs(h,du-1)); */
2724 207 : p1 = gpowgs(gdiv(z,h),du-1);
2725 207 : z = gmul(z,p1);
2726 207 : uze = RgX_Rg_mul(uze, p1);
2727 : }
2728 1983 : if (signh < 0) { z = gneg_i(z); uze = RgX_neg(uze); }
2729 :
2730 1983 : vze = RgX_divrem(Rg_RgX_sub(z, RgX_mul(uze,x)), y, &r);
2731 1983 : if (signe(r)) pari_warn(warner,"inexact computation in subresext");
2732 : /* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */
2733 1983 : p1 = gen_1;
2734 1983 : if (cu) p1 = gmul(p1, gpowgs(cu,dy));
2735 1983 : if (cv) p1 = gmul(p1, gpowgs(cv,dx));
2736 1983 : cu = cu? gdiv(p1,cu): p1;
2737 1983 : cv = cv? gdiv(p1,cv): p1;
2738 1983 : z = gmul(z,p1);
2739 1983 : *U = RgX_Rg_mul(uze,cu);
2740 1983 : *V = RgX_Rg_mul(vze,cv);
2741 1983 : return z;
2742 : }
2743 : GEN
2744 0 : subresext(GEN x, GEN y, GEN *U, GEN *V)
2745 : {
2746 0 : pari_sp av = avma;
2747 0 : GEN z = subresext_i(x, y, U, V);
2748 0 : return gc_all(av, 3, &z, U, V);
2749 : }
2750 :
2751 : static GEN
2752 372 : zero_extgcd(GEN y, GEN *U, GEN *V, long vx)
2753 : {
2754 372 : GEN x=content(y);
2755 372 : *U=pol_0(vx); *V = scalarpol(ginv(x), vx); return gmul(y,*V);
2756 : }
2757 :
2758 : static int
2759 3736 : must_negate(GEN x)
2760 : {
2761 3736 : GEN t = leading_coeff(x);
2762 3736 : switch(typ(t))
2763 : {
2764 3652 : case t_INT: case t_REAL:
2765 3652 : return (signe(t) < 0);
2766 0 : case t_FRAC:
2767 0 : return (signe(gel(t,1)) < 0);
2768 : }
2769 84 : return 0;
2770 : }
2771 :
2772 : static GEN
2773 184 : gc_gcdext(pari_sp av, GEN r, GEN *u, GEN *v)
2774 : {
2775 184 : if (!u && !v) return gc_upto(av, r);
2776 184 : if (u && v) return gc_all(av, 3, &r, u, v);
2777 0 : return gc_all(av, 2, &r, u ? u: v);
2778 : }
2779 :
2780 : static GEN
2781 113 : RgX_extgcd_FpX(GEN x, GEN y, GEN p, GEN *u, GEN *v)
2782 : {
2783 113 : pari_sp av = avma;
2784 113 : GEN r = FpX_extgcd(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, u, v);
2785 113 : if (u) *u = FpX_to_mod(*u, p);
2786 113 : if (v) *v = FpX_to_mod(*v, p);
2787 113 : return gc_gcdext(av, FpX_to_mod(r, p), u, v);
2788 : }
2789 :
2790 : static GEN
2791 5 : RgX_extgcd_FpXQX(GEN x, GEN y, GEN pol, GEN p, GEN *U, GEN *V)
2792 : {
2793 5 : pari_sp av = avma;
2794 5 : GEN r, T = RgX_to_FpX(pol, p);
2795 5 : r = FpXQX_extgcd(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p, U, V);
2796 5 : return gc_gcdext(av, FpXQX_to_mod(r, T, p), U, V);
2797 : }
2798 :
2799 : static GEN
2800 3872 : RgX_extgcd_fast(GEN x, GEN y, GEN *U, GEN *V)
2801 : {
2802 : GEN p, pol;
2803 : long pa;
2804 3872 : long t = RgX_type(x, &p,&pol,&pa);
2805 3872 : switch(t)
2806 : {
2807 16 : case t_FFELT: return FFX_extgcd(x, y, pol, U, V);
2808 113 : case t_INTMOD: return RgX_extgcd_FpX(x, y, p, U, V);
2809 5 : case RgX_type_code(t_POLMOD, t_INTMOD):
2810 5 : return RgX_extgcd_FpXQX(x, y, pol, p, U, V);
2811 3738 : default: return NULL;
2812 : }
2813 : }
2814 :
2815 : /* compute U, V s.t Ux + Vy = GCD(x,y) using subresultant */
2816 : GEN
2817 4250 : RgX_extgcd(GEN x, GEN y, GEN *U, GEN *V)
2818 : {
2819 : pari_sp av, av2, tetpil;
2820 : long signh; /* junk */
2821 4250 : long dx, dy, vx, tx = typ(x), ty = typ(y);
2822 : GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze;
2823 :
2824 4250 : if (tx!=t_POL) pari_err_TYPE("RgX_extgcd",x);
2825 4250 : if (ty!=t_POL) pari_err_TYPE("RgX_extgcd",y);
2826 4250 : if ( varncmp(varn(x),varn(y))) pari_err_VAR("RgX_extgcd",x,y);
2827 4250 : vx=varn(x);
2828 4250 : if (!signe(x))
2829 : {
2830 12 : if (signe(y)) return zero_extgcd(y,U,V,vx);
2831 6 : *U = pol_0(vx); *V = pol_0(vx);
2832 6 : return pol_0(vx);
2833 : }
2834 4238 : if (!signe(y)) return zero_extgcd(x,V,U,vx);
2835 3872 : r = RgX_extgcd_fast(x, y, U, V);
2836 3872 : if (r) return r;
2837 3738 : dx = degpol(x); dy = degpol(y);
2838 3738 : if (dx < dy) { pswap(U,V); lswap(dx,dy); swap(x,y); }
2839 3738 : if (dy==0) { *U=pol_0(vx); *V=ginv(y); return pol_1(vx); }
2840 :
2841 3576 : av = avma;
2842 3576 : u = x = primitive_part(x, &cu);
2843 3576 : v = y = primitive_part(y, &cv);
2844 3576 : g = h = gen_1; av2 = avma;
2845 3576 : um1 = gen_1; uze = gen_0;
2846 : for(;;)
2847 : {
2848 3738 : if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
2849 162 : if (gc_needed(av2,1))
2850 : {
2851 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_extgcd, dr = %ld",degpol(v));
2852 0 : (void)gc_all(av2,6,&u,&v,&g,&h,&uze,&um1);
2853 : }
2854 : }
2855 3576 : if (uze != gen_0) {
2856 : GEN r;
2857 3396 : vze = RgX_divrem(RgX_sub(v, RgX_mul(uze,x)), y, &r);
2858 3396 : if (signe(r)) pari_warn(warner,"inexact computation in RgX_extgcd");
2859 3396 : if (cu) uze = RgX_Rg_div(uze,cu);
2860 3396 : if (cv) vze = RgX_Rg_div(vze,cv);
2861 3396 : p1 = ginv(content(v));
2862 : }
2863 : else /* y | x */
2864 : {
2865 180 : vze = cv ? RgX_Rg_div(pol_1(vx),cv): pol_1(vx);
2866 180 : uze = pol_0(vx);
2867 180 : p1 = gen_1;
2868 : }
2869 3576 : if (must_negate(v)) p1 = gneg(p1);
2870 3576 : tetpil = avma;
2871 3576 : z = RgX_Rg_mul(v,p1);
2872 3576 : *U = RgX_Rg_mul(uze,p1);
2873 3576 : *V = RgX_Rg_mul(vze,p1);
2874 3576 : return gc_all_unsafe(av,tetpil, 3, &z, U, V);
2875 : }
2876 :
2877 : static GEN
2878 12 : RgX_halfgcd_all_i(GEN a, GEN b, GEN *pa, GEN *pb)
2879 : {
2880 12 : pari_sp av=avma;
2881 12 : long m = degpol(a), va = varn(a);
2882 : GEN R, u,u1,v,v1;
2883 12 : u1 = v = pol_0(va);
2884 12 : u = v1 = pol_1(va);
2885 12 : if (degpol(a)<degpol(b))
2886 : {
2887 0 : swap(a,b);
2888 0 : swap(u,v); swap(u1,v1);
2889 : }
2890 36 : while (2*degpol(b) >= m)
2891 : {
2892 24 : GEN r, q = RgX_pseudodivrem(a,b,&r);
2893 24 : GEN l = gpowgs(leading_coeff(b), degpol(a)-degpol(b)+1);
2894 24 : GEN g = ggcd(l, content(r));
2895 24 : q = RgX_Rg_div(q, g);
2896 24 : r = RgX_Rg_div(r, g);
2897 24 : l = gdiv(l, g);
2898 24 : a = b; b = r; swap(u,v); swap(u1,v1);
2899 24 : v = RgX_sub(gmul(l,v), RgX_mul(u, q));
2900 24 : v1 = RgX_sub(gmul(l,v1), RgX_mul(u1, q));
2901 24 : if (gc_needed(av,2))
2902 : {
2903 0 : if (DEBUGMEM>1) pari_warn(warnmem,"halfgcd (d = %ld)",degpol(b));
2904 0 : (void)gc_all(av,6, &a,&b,&u1,&v1,&u,&v);
2905 : }
2906 : }
2907 12 : if (pa) *pa = a;
2908 12 : if (pb) *pb = b;
2909 12 : R = mkmat22(u,u1,v,v1);
2910 12 : return !pa && pb ? gc_all(av, 2, &R, pb): gc_all(av, 1+!!pa+!!pb, &R, pa, pb);
2911 : }
2912 :
2913 : static GEN
2914 24 : RgX_halfgcd_all_FpX(GEN x, GEN y, GEN p, GEN *a, GEN *b)
2915 : {
2916 24 : pari_sp av = avma;
2917 : GEN M;
2918 24 : if (lgefint(p) == 3)
2919 : {
2920 12 : ulong pp = uel(p, 2);
2921 12 : GEN xp = RgX_to_Flx(x, pp), yp = RgX_to_Flx(y, pp);
2922 12 : M = Flx_halfgcd_all(xp, yp, pp, a, b);
2923 12 : M = FlxM_to_ZXM(M); *a = Flx_to_ZX(*a); *b = Flx_to_ZX(*b);
2924 : }
2925 : else
2926 : {
2927 12 : x = RgX_to_FpX(x, p); y = RgX_to_FpX(y, p);
2928 12 : M = FpX_halfgcd_all(x, y, p, a, b);
2929 : }
2930 24 : return !a && b ? gc_all(av, 2, &M, b): gc_all(av, 1+!!a+!!b, &M, a, b);
2931 : }
2932 :
2933 : static GEN
2934 0 : RgX_halfgcd_all_FpXQX(GEN x, GEN y, GEN pol, GEN p, GEN *a, GEN *b)
2935 : {
2936 0 : pari_sp av = avma;
2937 0 : GEN M, T = RgX_to_FpX(pol, p);
2938 0 : if (signe(T)==0) pari_err_OP("halfgcd", x, y);
2939 0 : x = RgX_to_FpXQX(x, T, p); y = RgX_to_FpXQX(y, T, p);
2940 0 : M = FpXQX_halfgcd_all(x, y, T, p, a, b);
2941 0 : if (a) *a = FqX_to_mod(*a, T, p);
2942 0 : if (b) *b = FqX_to_mod(*b, T, p);
2943 0 : M = FqXM_to_mod(M, T, p);
2944 0 : return !a && b ? gc_all(av, 2, &M, b): gc_all(av, 1+!!a+!!b, &M, a, b);
2945 : }
2946 :
2947 : static GEN
2948 54 : RgX_halfgcd_all_fast(GEN x, GEN y, GEN *a, GEN *b)
2949 : {
2950 : GEN p, pol;
2951 : long pa;
2952 54 : long t = RgX_type2(x,y, &p,&pol,&pa);
2953 54 : switch(t)
2954 : {
2955 18 : case t_FFELT: return FFX_halfgcd_all(x, y, pol, a, b);
2956 24 : case t_INTMOD: return RgX_halfgcd_all_FpX(x, y, p, a, b);
2957 0 : case RgX_type_code(t_POLMOD, t_INTMOD):
2958 0 : return RgX_halfgcd_all_FpXQX(x, y, pol, p, a, b);
2959 12 : default: return NULL;
2960 : }
2961 : }
2962 :
2963 : GEN
2964 54 : RgX_halfgcd_all(GEN x, GEN y, GEN *a, GEN *b)
2965 : {
2966 54 : GEN z = RgX_halfgcd_all_fast(x, y, a, b);
2967 54 : if (z) return z;
2968 12 : return RgX_halfgcd_all_i(x, y, a, b);
2969 : }
2970 :
2971 : GEN
2972 0 : RgX_halfgcd(GEN x, GEN y)
2973 0 : { return RgX_halfgcd_all(x, y, NULL, NULL); }
2974 :
2975 : int
2976 96 : RgXQ_ratlift(GEN x, GEN T, long amax, long bmax, GEN *P, GEN *Q)
2977 : {
2978 96 : pari_sp av = avma, av2, tetpil;
2979 : long signh; /* junk */
2980 : long vx;
2981 : GEN g, h, p1, cu, cv, u, v, um1, uze;
2982 :
2983 96 : if (typ(x)!=t_POL) pari_err_TYPE("RgXQ_ratlift",x);
2984 96 : if (typ(T)!=t_POL) pari_err_TYPE("RgXQ_ratlift",T);
2985 96 : if ( varncmp(varn(x),varn(T)) ) pari_err_VAR("RgXQ_ratlift",x,T);
2986 96 : if (bmax < 0) pari_err_DOMAIN("ratlift", "bmax", "<", gen_0, stoi(bmax));
2987 96 : if (!signe(T)) {
2988 0 : if (degpol(x) <= amax) {
2989 0 : *P = RgX_copy(x);
2990 0 : *Q = pol_1(varn(x));
2991 0 : return 1;
2992 : }
2993 0 : return 0;
2994 : }
2995 96 : if (amax+bmax >= degpol(T))
2996 0 : pari_err_DOMAIN("ratlift", "amax+bmax", ">=", stoi(degpol(T)),
2997 : mkvec3(stoi(amax), stoi(bmax), T));
2998 96 : vx = varn(T);
2999 96 : u = x = primitive_part(x, &cu);
3000 96 : v = T = primitive_part(T, &cv);
3001 96 : g = h = gen_1; av2 = avma;
3002 96 : um1 = gen_1; uze = gen_0;
3003 : for(;;)
3004 : {
3005 348 : (void) subres_step(&u, &v, &g, &h, &uze, &um1, &signh);
3006 348 : if (!u || (typ(uze)==t_POL && degpol(uze)>bmax)) return gc_bool(av,0);
3007 348 : if (typ(v)!=t_POL || degpol(v)<=amax) break;
3008 252 : if (gc_needed(av2,1))
3009 : {
3010 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXQ_ratlift, dr = %ld", degpol(v));
3011 0 : (void)gc_all(av2,6,&u,&v,&g,&h,&uze,&um1);
3012 : }
3013 : }
3014 96 : if (uze == gen_0)
3015 : {
3016 0 : set_avma(av); *P = pol_0(vx); *Q = pol_1(vx);
3017 0 : return 1;
3018 : }
3019 96 : if (cu) uze = RgX_Rg_div(uze,cu);
3020 96 : p1 = ginv(content(v));
3021 96 : if (must_negate(v)) p1 = gneg(p1);
3022 96 : tetpil = avma;
3023 96 : *P = RgX_Rg_mul(v,p1);
3024 96 : *Q = RgX_Rg_mul(uze,p1);
3025 96 : (void)gc_all_unsafe(av,tetpil,2,P,Q); return 1;
3026 : }
3027 :
3028 : GEN
3029 0 : RgX_chinese_coprime(GEN x, GEN y, GEN Tx, GEN Ty, GEN Tz)
3030 : {
3031 0 : pari_sp av = avma;
3032 0 : GEN ax = RgX_mul(RgXQ_inv(Tx,Ty), Tx);
3033 0 : GEN p1 = RgX_mul(ax, RgX_sub(y,x));
3034 0 : p1 = RgX_add(x,p1);
3035 0 : if (!Tz) Tz = RgX_mul(Tx,Ty);
3036 0 : p1 = RgX_rem(p1, Tz);
3037 0 : return gc_upto(av,p1);
3038 : }
3039 :
3040 : /*******************************************************************/
3041 : /* */
3042 : /* RESULTANT USING DUCOS VARIANT */
3043 : /* */
3044 : /*******************************************************************/
3045 : /* x^n / y^(n-1), assume n > 0 */
3046 : static GEN
3047 109979 : Lazard(GEN x, GEN y, long n)
3048 : {
3049 : long a;
3050 : GEN c;
3051 :
3052 109979 : if (n == 1) return x;
3053 630 : a = 1 << expu(n); /* a = 2^k <= n < 2^(k+1) */
3054 630 : c=x; n-=a;
3055 1365 : while (a>1)
3056 : {
3057 735 : a>>=1; c=gdivexact(gsqr(c),y);
3058 735 : if (n>=a) { c=gdivexact(gmul(c,x),y); n -= a; }
3059 : }
3060 630 : return c;
3061 : }
3062 :
3063 : /* F (x/y)^(n-1), assume n >= 1 */
3064 : static GEN
3065 244592 : Lazard2(GEN F, GEN x, GEN y, long n)
3066 : {
3067 244592 : if (n == 1) return F;
3068 1594 : return RgX_Rg_divexact(RgX_Rg_mul(F, Lazard(x,y,n-1)), y);
3069 : }
3070 :
3071 : static GEN
3072 244592 : RgX_neg_i(GEN x, long lx)
3073 : {
3074 : long i;
3075 244592 : GEN y = cgetg(lx, t_POL); y[1] = x[1];
3076 839064 : for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
3077 244592 : return y;
3078 : }
3079 : static GEN
3080 733083 : RgX_Rg_mul_i(GEN y, GEN x, long ly)
3081 : {
3082 : long i;
3083 : GEN z;
3084 733083 : if (isrationalzero(x)) return pol_0(varn(y));
3085 733065 : z = cgetg(ly,t_POL); z[1] = y[1];
3086 2512741 : for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
3087 733065 : return z;
3088 : }
3089 : static long
3090 731577 : reductum_lg(GEN x, long lx)
3091 : {
3092 731577 : long i = lx-2;
3093 736666 : while (i > 1 && gequal0(gel(x,i))) i--;
3094 731577 : return i+1;
3095 : }
3096 :
3097 : #define addshift(x,y) RgX_addmulXn_shallow((x),(y),1)
3098 : /* delta = deg(P) - deg(Q) > 0, deg(Q) > 0, P,Q,Z t_POL in the same variable,
3099 : * s "scalar". Return prem(P, -Q) / s^delta lc(P) */
3100 : static GEN
3101 244592 : nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
3102 : {
3103 244592 : GEN p0, q0, h0, TMP, H, A, z0 = leading_coeff(Z);
3104 : long p, q, j, lP, lQ;
3105 : pari_sp av;
3106 :
3107 244592 : p = degpol(P); p0 = gel(P,p+2); lP = reductum_lg(P,lg(P));
3108 244592 : q = degpol(Q); q0 = gel(Q,q+2); lQ = reductum_lg(Q,lg(Q));
3109 : /* p > q. Very often p - 1 = q */
3110 244592 : av = avma;
3111 : /* H = RgX_neg(reductum(Z)) optimized, using Q ~ Z */
3112 244592 : H = RgX_neg_i(Z, lQ); /* deg H < q */
3113 :
3114 244592 : A = (q+2 < lP)? RgX_Rg_mul_i(H, gel(P,q+2), lQ): NULL;
3115 247265 : for (j = q+1; j < p; j++)
3116 : {
3117 2673 : if (degpol(H) == q-1)
3118 : { /* h0 = coeff of degree q-1 = leading coeff */
3119 2046 : h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1);
3120 2046 : H = addshift(H, RgX_Rg_divexact(RgX_Rg_mul_i(Q, gneg(h0), lQ), q0));
3121 : }
3122 : else
3123 627 : H = RgX_shift_shallow(H, 1);
3124 2673 : if (j+2 < lP)
3125 : {
3126 1784 : TMP = RgX_Rg_mul(H, gel(P,j+2));
3127 1784 : A = A? RgX_add(A, TMP): TMP;
3128 : }
3129 2673 : if (gc_needed(av,1))
3130 : {
3131 102 : if(DEBUGMEM>1) pari_warn(warnmem,"nextSousResultant j = %ld/%ld",j,p);
3132 102 : (void)gc_all(av,A?2:1,&H,&A);
3133 : }
3134 : }
3135 244592 : if (q+2 < lP) lP = reductum_lg(P, q+3);
3136 244592 : TMP = RgX_Rg_mul_i(P, z0, lP);
3137 244592 : A = A? RgX_add(A, TMP): TMP;
3138 244592 : A = RgX_Rg_divexact(A, p0);
3139 244592 : if (degpol(H) == q-1)
3140 : {
3141 244052 : h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1); /* destroy old H */
3142 244052 : A = RgX_add(RgX_Rg_mul(addshift(H,A),q0), RgX_Rg_mul_i(Q, gneg(h0), lQ));
3143 : }
3144 : else
3145 540 : A = RgX_Rg_mul(addshift(H,A), q0);
3146 244592 : return RgX_Rg_divexact(A, s);
3147 : }
3148 : #undef addshift
3149 :
3150 : static GEN
3151 216770 : RgX_pseudodenom(GEN x)
3152 : {
3153 216770 : GEN m = NULL;
3154 216770 : long l = lg(x), i;
3155 1254692 : for (i = 2; i < l; i++)
3156 : {
3157 1037922 : GEN xi = gel(x, i);
3158 1037922 : if (typ(xi) == t_RFRAC)
3159 : {
3160 32 : GEN d = denom_i(xi);
3161 32 : if (!m || signe(RgX_pseudorem(m, d)))
3162 32 : m = m ? gmul(m, d): d;
3163 : }
3164 : }
3165 216770 : return m;
3166 : }
3167 :
3168 : /* Ducos's subresultant */
3169 : GEN
3170 249528 : RgX_resultant_all(GEN P, GEN Q, GEN *sol)
3171 : {
3172 : pari_sp av, av2;
3173 249528 : long dP, dQ, delta, sig = 1;
3174 : GEN DP, DQ, cP, cQ, Z, s;
3175 :
3176 249528 : dP = degpol(P);
3177 249528 : dQ = degpol(Q); delta = dP - dQ;
3178 249528 : if (delta < 0)
3179 : {
3180 1571 : if (both_odd(dP, dQ)) sig = -1;
3181 1571 : swap(P,Q); lswap(dP, dQ); delta = -delta;
3182 : }
3183 249528 : if (sol) *sol = gen_0;
3184 249528 : av = avma;
3185 249528 : if (dQ <= 0)
3186 : {
3187 924 : if (dQ < 0) return Rg_get_0(P);
3188 924 : s = gpowgs(gel(Q,2), dP);
3189 924 : if (sig == -1) s = gc_upto(av, gneg(s));
3190 924 : return s;
3191 : }
3192 248604 : if (dQ == 1)
3193 : {
3194 140219 : if (sol) *sol = Q;
3195 140219 : s = RgX_homogenous_evalpow(P, gel(Q,2), gpowers(gneg(gel(Q,3)), dP));
3196 140219 : if (sig==-1) s = gneg(s);
3197 140219 : return gc_all(av, sol ? 2: 1, &s, sol);
3198 : }
3199 : /* primitive_part is also possible here, but possibly very costly,
3200 : * and hardly ever worth it */
3201 :
3202 108385 : DP = RgX_pseudodenom(P); if (DP) P = gmul(P,DP);
3203 108385 : DQ = RgX_pseudodenom(Q); if (DQ) Q = gmul(Q,DQ);
3204 108385 : P = Q_primitive_part(P, &cP);
3205 108385 : Q = Q_primitive_part(Q, &cQ);
3206 108385 : av2 = avma;
3207 108385 : s = gpowgs(leading_coeff(Q),delta);
3208 108385 : if (both_odd(dP, dQ)) sig = -sig;
3209 108385 : Z = Q;
3210 108385 : Q = RgX_pseudorem(P, Q);
3211 108385 : P = Z;
3212 352977 : while(degpol(Q) > 0)
3213 : {
3214 244592 : delta = degpol(P) - degpol(Q); /* > 0 */
3215 244592 : Z = Lazard2(Q, leading_coeff(Q), s, delta);
3216 244592 : if (both_odd(degpol(P), degpol(Q))) sig = -sig;
3217 244592 : Q = nextSousResultant(P, Q, Z, s);
3218 244592 : P = Z;
3219 244592 : if (gc_needed(av,1))
3220 : {
3221 9 : if(DEBUGMEM>1) pari_warn(warnmem,"resultant_all, degpol Q = %ld",degpol(Q));
3222 9 : (void)gc_all(av2,2,&P,&Q);
3223 : }
3224 244592 : s = leading_coeff(P);
3225 : }
3226 108385 : if (!signe(Q)) { set_avma(av); return Rg_get_0(Q); }
3227 108385 : s = Lazard(leading_coeff(Q), s, degpol(P));
3228 108385 : if (sig == -1) s = gneg(s);
3229 108385 : if (DP) s = gdiv(s, gpowgs(DP,dQ));
3230 108385 : if (DQ) s = gdiv(s, gpowgs(DQ,dP));
3231 108385 : if (cP) s = gmul(s, gpowgs(cP,dQ));
3232 108385 : if (cQ) s = gmul(s, gpowgs(cQ,dP));
3233 108385 : if (!sol) return gc_GEN(av, s);
3234 2269 : *sol = P; return gc_all(av, 2, &s, sol);
3235 : }
3236 :
3237 : static GEN
3238 20 : RgX_resultant_FpX(GEN x, GEN y, GEN p)
3239 : {
3240 20 : pari_sp av = avma;
3241 : GEN r;
3242 20 : if (lgefint(p) == 3)
3243 : {
3244 10 : ulong pp = uel(p, 2);
3245 10 : r = utoi(Flx_resultant(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), pp));
3246 : }
3247 : else
3248 10 : r = FpX_resultant(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3249 20 : return gc_upto(av, Fp_to_mod(r, p));
3250 : }
3251 :
3252 : static GEN
3253 15 : RgX_resultant_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3254 : {
3255 15 : pari_sp av = avma;
3256 15 : GEN r, T = RgX_to_FpX(pol, p);
3257 15 : r = FpXQX_resultant(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
3258 15 : return gc_upto(av, FpX_to_mod(r, p));
3259 : }
3260 :
3261 : static GEN
3262 1131032 : resultant_fast(GEN x, GEN y)
3263 : {
3264 : GEN p, pol;
3265 : long pa, t;
3266 1131032 : p = init_resultant(x,y);
3267 1131032 : if (p) return p;
3268 1131012 : t = RgX_type2(x,y, &p,&pol,&pa);
3269 1131012 : switch(t)
3270 : {
3271 226 : case t_INT: return ZX_resultant(x,y);
3272 40 : case t_FRAC: return QX_resultant(x,y);
3273 15 : case t_FFELT: return FFX_resultant(x,y,pol);
3274 20 : case t_INTMOD: return RgX_resultant_FpX(x, y, p);
3275 15 : case RgX_type_code(t_POLMOD, t_INTMOD):
3276 15 : return RgX_resultant_FpXQX(x, y, pol, p);
3277 807644 : case RgX_type_code(t_POL, t_INT):
3278 : {
3279 807644 : long v = -1;
3280 807644 : if (varn(x)==varn(y) && RgX_is_ZXX(x, &v) && RgX_is_ZXX(y, &v) && v>=0)
3281 807481 : return ZXX_resultant(x,y,v);
3282 : } /* FALL THROUGH */
3283 323215 : default: return NULL;
3284 : }
3285 : }
3286 :
3287 : static GEN
3288 145647 : RgX_resultant_sylvester(GEN x, GEN y)
3289 : {
3290 145647 : pari_sp av = avma;
3291 145647 : return gc_upto(av, det(RgX_sylvestermatrix(x,y)));
3292 : }
3293 :
3294 : /* Return resultant(P,Q).
3295 : * Uses Sylvester's matrix if P or Q inexact, a modular algorithm if they
3296 : * are in Q[X], and Ducos/Lazard optimization of the subresultant algorithm
3297 : * in the "generic" case. */
3298 : GEN
3299 1131032 : resultant(GEN P, GEN Q)
3300 : {
3301 1131032 : GEN z = resultant_fast(P,Q);
3302 1131032 : if (z) return z;
3303 323215 : if (isinexact(P) || isinexact(Q)) return RgX_resultant_sylvester(P,Q);
3304 177591 : return RgX_resultant_all(P, Q, NULL);
3305 : }
3306 :
3307 : /*******************************************************************/
3308 : /* */
3309 : /* RESULTANT USING SYLVESTER MATRIX */
3310 : /* */
3311 : /*******************************************************************/
3312 : static GEN
3313 318594 : syl_RgC(GEN x, long j, long d, long D, long cp)
3314 : {
3315 318594 : GEN c = cgetg(d+1,t_COL);
3316 : long i;
3317 848615 : for (i=1; i< j; i++) gel(c,i) = gen_0;
3318 1836220 : for ( ; i<=D; i++) { GEN t = gel(x,D-i+2); gel(c,i) = cp? gcopy(t): t; }
3319 848615 : for ( ; i<=d; i++) gel(c,i) = gen_0;
3320 318594 : return c;
3321 : }
3322 : static GEN
3323 145652 : syl_RgM(GEN x, GEN y, long cp)
3324 : {
3325 145652 : long j, d, dx = degpol(x), dy = degpol(y);
3326 : GEN M;
3327 145652 : if (dx < 0) return dy < 0? cgetg(1,t_MAT): zeromat(dy,dy);
3328 145652 : if (dy < 0) return zeromat(dx,dx);
3329 145652 : d = dx+dy; M = cgetg(d+1,t_MAT);
3330 378938 : for (j=1; j<=dy; j++) gel(M,j) = syl_RgC(x,j,d,j+dx, cp);
3331 230960 : for (j=1; j<=dx; j++) gel(M,j+dy) = syl_RgC(y,j,d,j+dy, cp);
3332 145652 : return M;
3333 : }
3334 : GEN
3335 145647 : RgX_sylvestermatrix(GEN x, GEN y) { return syl_RgM(x,y,0); }
3336 : GEN
3337 5 : sylvestermatrix(GEN x, GEN y)
3338 : {
3339 5 : if (typ(x)!=t_POL) pari_err_TYPE("sylvestermatrix",x);
3340 5 : if (typ(y)!=t_POL) pari_err_TYPE("sylvestermatrix",y);
3341 5 : if (varn(x) != varn(y)) pari_err_VAR("sylvestermatrix",x,y);
3342 5 : return syl_RgM(x,y,1);
3343 : }
3344 :
3345 : GEN
3346 23 : resultant2(GEN x, GEN y)
3347 : {
3348 23 : GEN r = init_resultant(x,y);
3349 23 : return r? r: RgX_resultant_sylvester(x,y);
3350 : }
3351 :
3352 : /* let vx = main variable of x, v0 a variable of highest priority;
3353 : * return a t_POL in variable v0:
3354 : * if vx <= v, return subst(x, v, pol_x(v0))
3355 : * if vx > v, return scalarpol(x, v0) */
3356 : static GEN
3357 272 : fix_pol(GEN x, long v, long v0)
3358 : {
3359 272 : long vx, tx = typ(x);
3360 272 : if (tx != t_POL)
3361 31 : vx = gvar(x);
3362 : else
3363 : { /* shortcut: almost nothing to do */
3364 241 : vx = varn(x);
3365 241 : if (v == vx)
3366 : {
3367 97 : if (v0 != v) { x = leafcopy(x); setvarn(x, v0); }
3368 97 : return x;
3369 : }
3370 : }
3371 175 : if (varncmp(v, vx) > 0)
3372 : {
3373 170 : x = gsubst(x, v, pol_x(v0));
3374 170 : if (typ(x) != t_POL) vx = gvar(x);
3375 : else
3376 : {
3377 164 : vx = varn(x);
3378 164 : if (vx == v0) return x;
3379 : }
3380 : }
3381 37 : if (varncmp(vx, v0) <= 0) pari_err_TYPE("polresultant", x);
3382 31 : return scalarpol_shallow(x, v0);
3383 : }
3384 :
3385 : /* resultant of x and y with respect to variable v, or with respect to their
3386 : * main variable if v < 0. */
3387 : GEN
3388 373 : polresultant0(GEN x, GEN y, long v, long flag)
3389 : {
3390 373 : pari_sp av = avma;
3391 :
3392 373 : if (v >= 0)
3393 : {
3394 117 : long v0 = fetch_var_higher();
3395 117 : x = fix_pol(x,v, v0);
3396 117 : y = fix_pol(y,v, v0);
3397 : }
3398 373 : switch(flag)
3399 : {
3400 368 : case 0: x=resultant(x,y); break;
3401 5 : case 1: x=resultant2(x,y); break;
3402 0 : case 2: x=RgX_resultant_all(x,y,NULL); break;
3403 0 : default: pari_err_FLAG("polresultant");
3404 : }
3405 373 : if (v >= 0) (void)delete_var();
3406 373 : return gc_upto(av,x);
3407 : }
3408 :
3409 : static GEN
3410 66 : RgX_extresultant_FpX(GEN x, GEN y, GEN p, GEN *u, GEN *v)
3411 : {
3412 66 : pari_sp av = avma;
3413 66 : GEN r = FpX_extresultant(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, u, v);
3414 66 : if (signe(r) == 0) { *u = gen_0; *v = gen_0; return gc_const(av, gen_0); }
3415 66 : if (u) *u = FpX_to_mod(*u, p);
3416 66 : if (v) *v = FpX_to_mod(*v, p);
3417 66 : return gc_gcdext(av, Fp_to_mod(r, p), u, v);
3418 : }
3419 :
3420 : static GEN
3421 1340 : RgX_extresultant_fast(GEN x, GEN y, GEN *U, GEN *V)
3422 : {
3423 : GEN p, pol;
3424 : long pa;
3425 1340 : long t = RgX_type2(x, y, &p,&pol,&pa);
3426 1340 : switch(t)
3427 : {
3428 66 : case t_INTMOD: return RgX_extresultant_FpX(x, y, p, U, V);
3429 1274 : default: return NULL;
3430 : }
3431 : }
3432 :
3433 : GEN
3434 1345 : polresultantext0(GEN x, GEN y, long v)
3435 : {
3436 1345 : GEN R = NULL, U, V;
3437 1345 : pari_sp av = avma;
3438 :
3439 1345 : if (v >= 0)
3440 : {
3441 10 : long v0 = fetch_var_higher();
3442 10 : x = fix_pol(x,v, v0);
3443 10 : y = fix_pol(y,v, v0);
3444 : }
3445 1345 : if (typ(x)==t_POL && typ(y)==t_POL)
3446 1340 : R = RgX_extresultant_fast(x, y, &U, &V);
3447 1345 : if (!R)
3448 1279 : R = subresext_i(x,y, &U,&V);
3449 1345 : if (v >= 0)
3450 : {
3451 10 : (void)delete_var();
3452 10 : if (typ(U) == t_POL && varn(U) != v) U = poleval(U, pol_x(v));
3453 10 : if (typ(V) == t_POL && varn(V) != v) V = poleval(V, pol_x(v));
3454 : }
3455 1345 : return gc_GEN(av, mkvec3(U,V,R));
3456 : }
3457 : GEN
3458 1254 : polresultantext(GEN x, GEN y) { return polresultantext0(x,y,-1); }
3459 :
3460 : /*******************************************************************/
3461 : /* */
3462 : /* CHARACTERISTIC POLYNOMIAL USING RESULTANT */
3463 : /* */
3464 : /*******************************************************************/
3465 :
3466 : static GEN
3467 12 : RgXQ_charpoly_FpXQ(GEN x, GEN T, GEN p, long v)
3468 : {
3469 12 : pari_sp av = avma;
3470 : GEN r;
3471 12 : if (lgefint(p)==3)
3472 : {
3473 0 : ulong pp = p[2];
3474 0 : r = Flx_to_ZX(Flxq_charpoly(RgX_to_Flx(x, pp), RgX_to_Flx(T, pp), pp));
3475 : }
3476 : else
3477 12 : r = FpXQ_charpoly(RgX_to_FpX(x, p), RgX_to_FpX(T, p), p);
3478 12 : r = FpX_to_mod(r, p); setvarn(r, v);
3479 12 : return gc_upto(av, r);
3480 : }
3481 :
3482 : static GEN
3483 10955 : RgXQ_charpoly_fast(GEN x, GEN T, long v)
3484 : {
3485 : GEN p, pol;
3486 10955 : long pa, t = RgX_type2(x,T, &p,&pol,&pa);
3487 10955 : switch(t)
3488 : {
3489 7962 : case t_INT: return ZXQ_charpoly(x, T, v);
3490 1871 : case t_FRAC:
3491 : {
3492 1871 : pari_sp av = avma;
3493 : GEN cT;
3494 1871 : T = Q_primitive_part(T, &cT);
3495 1871 : T = QXQ_charpoly(x, T, v);
3496 1871 : if (cT) T = gc_upto(av, T); /* silly rare case */
3497 1871 : return T;
3498 : }
3499 12 : case t_INTMOD: return RgXQ_charpoly_FpXQ(x, T, p, v);
3500 1110 : default: return NULL;
3501 : }
3502 : }
3503 :
3504 : /* (v - x)^d */
3505 : static GEN
3506 96 : caract_const(pari_sp av, GEN x, long v, long d)
3507 96 : { return gc_upto(av, gpowgs(gsub(pol_x(v), x), d)); }
3508 :
3509 : GEN
3510 983925 : RgXQ_charpoly_i(GEN x, GEN T, long v)
3511 : {
3512 983925 : pari_sp av = avma;
3513 983925 : long d = degpol(T), dx = degpol(x), v0;
3514 : GEN ch, L;
3515 983925 : if (dx >= degpol(T)) { x = RgX_rem(x, T); dx = degpol(x); }
3516 983925 : if (dx <= 0) return dx? pol_xn(d, v): caract_const(av, gel(x,2), v, d);
3517 :
3518 983867 : v0 = fetch_var_higher();
3519 983867 : x = RgX_neg(x);
3520 983867 : gel(x,2) = gadd(gel(x,2), pol_x(v));
3521 983867 : setvarn(x, v0);
3522 983867 : T = leafcopy(T); setvarn(T, v0);
3523 983867 : ch = resultant(T, x);
3524 983867 : (void)delete_var();
3525 : /* test for silly input: x mod (deg 0 polynomial) */
3526 983867 : if (typ(ch) != t_POL)
3527 6 : pari_err_PRIORITY("RgXQ_charpoly", pol_x(v), "<", gvar(ch));
3528 983861 : L = leading_coeff(ch);
3529 983861 : if (!gequal1(L)) ch = RgX_Rg_div(ch, L);
3530 983861 : return gc_upto(av, ch);
3531 : }
3532 :
3533 : /* return caract(Mod(x,T)) in variable v */
3534 : GEN
3535 10955 : RgXQ_charpoly(GEN x, GEN T, long v)
3536 : {
3537 10955 : GEN ch = RgXQ_charpoly_fast(x, T, v);
3538 10955 : if (ch) return ch;
3539 1110 : return RgXQ_charpoly_i(x, T, v);
3540 : }
3541 :
3542 : /* characteristic polynomial (in v) of x over nf, where x is an element of the
3543 : * algebra nf[t]/(Q(t)) */
3544 : GEN
3545 160 : rnfcharpoly(GEN nf, GEN Q, GEN x, long v)
3546 : {
3547 160 : const char *f = "rnfcharpoly";
3548 160 : long dQ = degpol(Q);
3549 160 : pari_sp av = avma;
3550 : GEN T;
3551 :
3552 160 : if (v < 0) v = 0;
3553 160 : nf = checknf(nf); T = nf_get_pol(nf);
3554 160 : Q = RgX_nffix(f, T,Q,0);
3555 160 : switch(typ(x))
3556 : {
3557 20 : case t_INT:
3558 20 : case t_FRAC: return caract_const(av, x, v, dQ);
3559 65 : case t_POLMOD:
3560 65 : x = polmod_nffix2(f,T,Q, x,0);
3561 40 : break;
3562 40 : case t_POL:
3563 40 : x = varn(x) == varn(T)? Rg_nffix(f,T,x,0): RgX_nffix(f, T,x,0);
3564 30 : break;
3565 35 : default: pari_err_TYPE(f,x);
3566 : }
3567 70 : if (typ(x) != t_POL) return caract_const(av, x, v, dQ);
3568 : /* x a t_POL in variable vQ */
3569 40 : if (degpol(x) >= dQ) x = RgX_rem(x, Q);
3570 40 : if (dQ <= 1) return caract_const(av, constant_coeff(x), v, 1);
3571 40 : return gc_GEN(av, lift_if_rational( RgXQ_charpoly(x, Q, v) ));
3572 : }
3573 :
3574 : /*******************************************************************/
3575 : /* */
3576 : /* GCD USING SUBRESULTANT */
3577 : /* */
3578 : /*******************************************************************/
3579 : static int inexact(GEN x, int *simple);
3580 : static int
3581 1832 : isinexactall(GEN x, int *simple)
3582 : {
3583 1832 : long i, lx = lg(x);
3584 11204 : for (i=2; i<lx; i++)
3585 9384 : if (inexact(gel(x,i), simple)) return 1;
3586 1820 : return 0;
3587 : }
3588 : /* return 1 if coeff explosion is not possible */
3589 : static int
3590 9600 : inexact(GEN x, int *simple)
3591 : {
3592 9600 : int junk = 0;
3593 9600 : switch(typ(x))
3594 : {
3595 6474 : case t_INT: case t_FRAC: return 0;
3596 :
3597 6 : case t_REAL: case t_PADIC: case t_SER: return 1;
3598 :
3599 1758 : case t_INTMOD:
3600 : case t_FFELT:
3601 1758 : if (!*simple) *simple = 1;
3602 1758 : return 0;
3603 :
3604 66 : case t_COMPLEX:
3605 66 : return inexact(gel(x,1), simple)
3606 66 : || inexact(gel(x,2), simple);
3607 0 : case t_QUAD:
3608 0 : *simple = 0;
3609 0 : return inexact(gel(x,2), &junk)
3610 0 : || inexact(gel(x,3), &junk);
3611 :
3612 666 : case t_POLMOD:
3613 666 : return isinexactall(gel(x,1), simple);
3614 588 : case t_POL:
3615 588 : *simple = -1;
3616 588 : return isinexactall(x, &junk);
3617 42 : case t_RFRAC:
3618 42 : *simple = -1;
3619 42 : return inexact(gel(x,1), &junk)
3620 42 : || inexact(gel(x,2), &junk);
3621 : }
3622 0 : *simple = -1; return 0;
3623 : }
3624 :
3625 : /* x monomial, y t_POL in the same variable */
3626 : static GEN
3627 2657 : gcdmonome(GEN x, GEN y)
3628 : {
3629 2657 : pari_sp av = avma;
3630 2657 : long dx = degpol(x), e = RgX_valrem(y, &y);
3631 2657 : long i, l = lg(y);
3632 2657 : GEN t, v = cgetg(l, t_VEC);
3633 2657 : gel(v,1) = gel(x,dx+2);
3634 5509 : for (i = 2; i < l; i++) gel(v,i) = gel(y,i);
3635 2657 : t = content(v); /* gcd(lc(x), cont(y)) */
3636 2657 : t = simplify_shallow(t);
3637 2657 : if (dx < e) e = dx;
3638 2657 : return gc_upto(av, monomialcopy(t, e, varn(x)));
3639 : }
3640 :
3641 : static GEN
3642 93260 : RgX_gcd_FpX(GEN x, GEN y, GEN p)
3643 : {
3644 93260 : pari_sp av = avma;
3645 : GEN r;
3646 93260 : if (lgefint(p) == 3)
3647 : {
3648 93248 : ulong pp = uel(p, 2);
3649 93248 : r = Flx_to_ZX_inplace(Flx_gcd(RgX_to_Flx(x, pp),
3650 : RgX_to_Flx(y, pp), pp));
3651 : }
3652 : else
3653 12 : r = FpX_gcd(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3654 93260 : return gc_upto(av, FpX_to_mod(r, p));
3655 : }
3656 :
3657 : static GEN
3658 5 : RgX_gcd_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3659 : {
3660 5 : pari_sp av = avma;
3661 5 : GEN r, T = RgX_to_FpX(pol, p);
3662 5 : if (signe(T)==0) pari_err_OP("gcd", x, y);
3663 5 : r = FpXQX_gcd(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
3664 5 : return gc_upto(av, FpXQX_to_mod(r, T, p));
3665 : }
3666 :
3667 : static GEN
3668 516 : RgX_gcd_FpXk(GEN x, GEN y, GEN p)
3669 : {
3670 516 : pari_sp av = avma;
3671 516 : GEN r = FpXk_gcd(Rg_to_FpXk(x, p), Rg_to_FpXk(y, p), p);
3672 516 : return gc_upto(av, gmul(r, gmodulsg(1,p)));
3673 : }
3674 :
3675 : static GEN
3676 9014 : RgX_liftred(GEN x, GEN T)
3677 9014 : { return RgXQX_red(liftpol_shallow(x), T); }
3678 :
3679 : static GEN
3680 1937 : RgX_gcd_ZXQX(GEN x, GEN y, GEN T)
3681 : {
3682 1937 : pari_sp av = avma;
3683 1937 : GEN r = ZXQX_gcd(RgX_liftred(x, T), RgX_liftred(y, T), T);
3684 1937 : return gc_GEN(av, QXQX_to_mod_shallow(r, T));
3685 : }
3686 :
3687 : static GEN
3688 2570 : RgX_gcd_QXQX(GEN x, GEN y, GEN T)
3689 : {
3690 2570 : pari_sp av = avma;
3691 2570 : GEN r = QXQX_gcd(RgX_liftred(x, T), RgX_liftred(y, T), T);
3692 2570 : return gc_GEN(av, QXQX_to_mod_shallow(r, T));
3693 : }
3694 :
3695 : static GEN
3696 8223574 : RgX_gcd_fast(GEN x, GEN y)
3697 : {
3698 : GEN p, pol;
3699 : long pa;
3700 8223574 : long t = RgX_type2(x,y, &p,&pol,&pa);
3701 8223574 : switch(t)
3702 : {
3703 6645366 : case t_INT: return ZX_gcd(x, y);
3704 6777 : case t_FRAC: return QX_gcd(x, y);
3705 2208 : case t_FFELT: return FFX_gcd(x, y, pol);
3706 93260 : case t_INTMOD: return RgX_gcd_FpX(x, y, p);
3707 5 : case RgX_type_code(t_POLMOD, t_INTMOD):
3708 5 : return RgX_gcd_FpXQX(x, y, pol, p);
3709 1943 : case RgX_type_code(t_POLMOD, t_INT):
3710 1943 : return ZX_is_monic(pol)? RgX_gcd_ZXQX(x,y,pol): NULL;
3711 2580 : case RgX_type_code(t_POLMOD, t_FRAC):
3712 5160 : return RgX_is_ZX(pol) && ZX_is_monic(pol) ?
3713 5160 : RgX_gcd_QXQX(x,y,pol): NULL;
3714 1467826 : case RgX_type_code(t_POL, t_INT):
3715 1467826 : return ZXk_gcd(x,y);
3716 160 : case RgX_type_code(t_POL, t_FRAC):
3717 160 : return QXk_gcd(x,y);
3718 516 : case RgX_type_code(t_POL, t_INTMOD):
3719 516 : return RgX_gcd_FpXk(x,y,p);
3720 2933 : default: return NULL;
3721 : }
3722 : }
3723 :
3724 : /* x, y are t_POL in the same variable */
3725 : GEN
3726 8223574 : RgX_gcd(GEN x, GEN y)
3727 : {
3728 : long dx, dy;
3729 : pari_sp av, av1;
3730 : GEN d, g, h, p1, p2, u, v;
3731 8223574 : int simple = 0;
3732 8223574 : GEN z = RgX_gcd_fast(x, y);
3733 8223574 : if (z) return z;
3734 2949 : if (isexactzero(y)) return RgX_copy(x);
3735 2949 : if (isexactzero(x)) return RgX_copy(y);
3736 2949 : if (RgX_is_monomial(x)) return gcdmonome(x,y);
3737 357 : if (RgX_is_monomial(y)) return gcdmonome(y,x);
3738 292 : if (isinexactall(x,&simple) || isinexactall(y,&simple))
3739 : {
3740 6 : av = avma; u = ggcd(content(x), content(y));
3741 6 : return gc_upto(av, scalarpol(u, varn(x)));
3742 : }
3743 :
3744 286 : av = avma;
3745 286 : if (simple > 0) x = RgX_gcd_simple(x,y);
3746 : else
3747 : {
3748 286 : dx = lg(x); dy = lg(y);
3749 286 : if (dx < dy) { swap(x,y); lswap(dx,dy); }
3750 286 : if (dy==3)
3751 : {
3752 0 : d = ggcd(gel(y,2), content(x));
3753 0 : return gc_upto(av, scalarpol(d, varn(x)));
3754 : }
3755 286 : u = primitive_part(x, &p1); if (!p1) p1 = gen_1;
3756 286 : v = primitive_part(y, &p2); if (!p2) p2 = gen_1;
3757 286 : d = ggcd(p1,p2);
3758 286 : av1 = avma;
3759 286 : g = h = gen_1;
3760 : for(;;)
3761 256 : {
3762 542 : GEN r = RgX_pseudorem(u,v);
3763 542 : long degq, du, dv, dr = lg(r);
3764 :
3765 542 : if (!signe(r)) break;
3766 478 : if (dr <= 3)
3767 : {
3768 222 : set_avma(av1);
3769 222 : return gc_upto(av, scalarpol(d, varn(x)));
3770 : }
3771 256 : du = lg(u); dv = lg(v); degq = du-dv;
3772 256 : u = v; p1 = g; g = leading_coeff(u);
3773 256 : switch(degq)
3774 : {
3775 10 : case 0: break;
3776 240 : case 1:
3777 240 : p1 = gmul(h,p1); h = g; break;
3778 6 : default:
3779 6 : p1 = gmul(gpowgs(h,degq), p1);
3780 6 : h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
3781 : }
3782 256 : v = RgX_Rg_div(r,p1);
3783 256 : if (gc_needed(av1,1))
3784 : {
3785 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd, dr = %ld", degpol(r));
3786 0 : (void)gc_all(av1,4, &u,&v,&g,&h);
3787 : }
3788 : }
3789 64 : x = RgX_Rg_mul(primpart(v), d);
3790 : }
3791 64 : if (must_negate(x)) x = RgX_neg(x);
3792 64 : return gc_upto(av,x);
3793 : }
3794 :
3795 : /* disc P = (-1)^(n(n-1)/2) lc(P)^(n - deg P' - 2) Res(P,P'), n = deg P */
3796 : static GEN
3797 337 : RgX_disc_i(GEN P)
3798 : {
3799 337 : long n = degpol(P), dd;
3800 : GEN N, D, L, y;
3801 337 : if (!signe(P) || !n) return Rg_get_0(P);
3802 331 : if (n == 1) return Rg_get_1(P);
3803 331 : if (n == 2) {
3804 120 : GEN a = gel(P,4), b = gel(P,3), c = gel(P,2);
3805 120 : return gsub(gsqr(b), gmul2n(gmul(a,c),2));
3806 : }
3807 211 : y = RgX_deriv(P);
3808 211 : N = characteristic(P);
3809 211 : if (signe(N)) y = gmul(y, mkintmod(gen_1,N));
3810 211 : if (!signe(y)) return Rg_get_0(y);
3811 211 : dd = n - 2 - degpol(y);
3812 211 : if (isinexact(P))
3813 18 : D = resultant2(P,y);
3814 : else
3815 : {
3816 193 : D = RgX_resultant_all(P, y, NULL);
3817 193 : if (D == gen_0) return Rg_get_0(y);
3818 : }
3819 211 : L = leading_coeff(P);
3820 211 : if (dd && !gequal1(L)) D = (dd == -1)? gdiv(D, L): gmul(D, gpowgs(L, dd));
3821 211 : if (n & 2) D = gneg(D);
3822 211 : return D;
3823 : }
3824 :
3825 : static GEN
3826 33 : RgX_disc_FpX(GEN x, GEN p)
3827 : {
3828 33 : pari_sp av = avma;
3829 33 : GEN r = FpX_disc(RgX_to_FpX(x, p), p);
3830 33 : return gc_upto(av, Fp_to_mod(r, p));
3831 : }
3832 :
3833 : static GEN
3834 23 : RgX_disc_FpXQX(GEN x, GEN pol, GEN p)
3835 : {
3836 23 : pari_sp av = avma;
3837 23 : GEN r, T = RgX_to_FpX(pol, p);
3838 23 : r = FpXQX_disc(RgX_to_FpXQX(x, T, p), T, p);
3839 23 : return gc_upto(av, FpX_to_mod(r, p));
3840 : }
3841 :
3842 : static GEN
3843 105376 : RgX_disc_fast(GEN x)
3844 : {
3845 : GEN p, pol;
3846 : long pa;
3847 105376 : long t = RgX_type(x, &p,&pol,&pa);
3848 105376 : switch(t)
3849 : {
3850 104950 : case t_INT: return ZX_disc(x);
3851 5 : case t_FRAC: return QX_disc(x);
3852 28 : case t_FFELT: return FFX_disc(x, pol);
3853 33 : case t_INTMOD: return RgX_disc_FpX(x, p);
3854 23 : case RgX_type_code(t_POLMOD, t_INTMOD):
3855 23 : return RgX_disc_FpXQX(x, pol, p);
3856 337 : default: return NULL;
3857 : }
3858 : }
3859 :
3860 : GEN
3861 105376 : RgX_disc(GEN x)
3862 : {
3863 : pari_sp av;
3864 105376 : GEN z = RgX_disc_fast(x);
3865 105376 : if (z) return z;
3866 337 : av = avma;
3867 337 : return gc_upto(av, RgX_disc_i(x));
3868 : }
3869 :
3870 : GEN
3871 2968 : poldisc0(GEN x, long v)
3872 : {
3873 2968 : long v0, tx = typ(x);
3874 : pari_sp av;
3875 : GEN D;
3876 2968 : if (tx == t_POL && (v < 0 || v == varn(x))) return RgX_disc(x);
3877 23 : switch(tx)
3878 : {
3879 0 : case t_QUAD:
3880 0 : return quad_disc(x);
3881 0 : case t_POLMOD:
3882 0 : if (v >= 0 && varn(gel(x,1)) != v) break;
3883 0 : return RgX_disc(gel(x,1));
3884 5 : case t_QFB:
3885 5 : return icopy(qfb_disc(x));
3886 0 : case t_VEC: case t_COL: case t_MAT:
3887 0 : pari_APPLY_same(poldisc0(gel(x,i), v));
3888 : }
3889 18 : if (v < 0) pari_err_TYPE("poldisc",x);
3890 18 : av = avma; v0 = fetch_var_higher();
3891 18 : x = fix_pol(x,v, v0);
3892 12 : D = RgX_disc(x); (void)delete_var();
3893 12 : return gc_upto(av, D);
3894 : }
3895 :
3896 : GEN
3897 5 : reduceddiscsmith(GEN x)
3898 : {
3899 5 : long j, n = degpol(x);
3900 5 : pari_sp av = avma;
3901 : GEN xp, M;
3902 :
3903 5 : if (typ(x) != t_POL) pari_err_TYPE("poldiscreduced",x);
3904 5 : if (n<=0) pari_err_CONSTPOL("poldiscreduced");
3905 5 : RgX_check_ZX(x,"poldiscreduced");
3906 5 : if (!gequal1(gel(x,n+2)))
3907 0 : pari_err_IMPL("nonmonic polynomial in poldiscreduced");
3908 5 : M = cgetg(n+1,t_MAT);
3909 5 : xp = ZX_deriv(x);
3910 20 : for (j=1; j<=n; j++)
3911 : {
3912 15 : gel(M,j) = RgX_to_RgC(xp, n);
3913 15 : if (j<n) xp = RgX_rem(RgX_shift_shallow(xp, 1), x);
3914 : }
3915 5 : return gc_upto(av, ZM_snf(M));
3916 : }
3917 :
3918 : /***********************************************************************/
3919 : /** **/
3920 : /** STURM ALGORITHM **/
3921 : /** (number of real roots of x in [a,b]) **/
3922 : /** **/
3923 : /***********************************************************************/
3924 : static GEN
3925 380 : R_to_Q_up(GEN x)
3926 : {
3927 : long e;
3928 380 : switch(typ(x))
3929 : {
3930 380 : case t_INT: case t_FRAC: case t_INFINITY: return x;
3931 0 : case t_REAL:
3932 0 : x = mantissa_real(x,&e);
3933 0 : return gmul2n(addiu(x,1), -e);
3934 0 : default: pari_err_TYPE("R_to_Q_up", x);
3935 : return NULL; /* LCOV_EXCL_LINE */
3936 : }
3937 : }
3938 : static GEN
3939 380 : R_to_Q_down(GEN x)
3940 : {
3941 : long e;
3942 380 : switch(typ(x))
3943 : {
3944 380 : case t_INT: case t_FRAC: case t_INFINITY: return x;
3945 0 : case t_REAL:
3946 0 : x = mantissa_real(x,&e);
3947 0 : return gmul2n(subiu(x,1), -e);
3948 0 : default: pari_err_TYPE("R_to_Q_down", x);
3949 : return NULL; /* LCOV_EXCL_LINE */
3950 : }
3951 : }
3952 :
3953 : static long
3954 850 : sturmpart_i(GEN x, GEN ab)
3955 : {
3956 850 : long tx = typ(x);
3957 850 : if (gequal0(x)) pari_err_ROOTS0("sturm");
3958 850 : if (tx != t_POL)
3959 : {
3960 0 : if (is_real_t(tx)) return 0;
3961 0 : pari_err_TYPE("sturm",x);
3962 : }
3963 850 : if (lg(x) == 3) return 0;
3964 850 : if (!RgX_is_ZX(x)) x = RgX_rescale_to_int(x);
3965 850 : (void)ZX_gcd_all(x, ZX_deriv(x), &x);
3966 850 : if (ab)
3967 : {
3968 : GEN A, B;
3969 380 : if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("RgX_sturmpart", ab);
3970 380 : A = R_to_Q_down(gel(ab,1));
3971 380 : B = R_to_Q_up(gel(ab,2));
3972 380 : ab = mkvec2(A, B);
3973 : }
3974 850 : return ZX_sturmpart(x, ab);
3975 : }
3976 : /* Deprecated: RgX_sturmpart() should be preferred */
3977 : long
3978 280 : sturmpart(GEN x, GEN a, GEN b)
3979 : {
3980 280 : pari_sp av = avma;
3981 280 : if (!b && a && typ(a) == t_VEC) return RgX_sturmpart(x, a);
3982 280 : if (!a) a = mkmoo();
3983 280 : if (!b) b = mkoo();
3984 280 : return gc_long(av, sturmpart_i(x, mkvec2(a,b)));
3985 : }
3986 : long
3987 570 : RgX_sturmpart(GEN x, GEN ab)
3988 570 : { pari_sp av = avma; return gc_long(av, sturmpart_i(x, ab)); }
3989 :
3990 : /***********************************************************************/
3991 : /** **/
3992 : /** GENERIC EXTENDED GCD **/
3993 : /** **/
3994 : /***********************************************************************/
3995 : /* assume typ(x) = typ(y) = t_POL */
3996 : static GEN
3997 715 : RgXQ_inv_i(GEN x, GEN y)
3998 : {
3999 715 : long vx=varn(x), vy=varn(y);
4000 : pari_sp av;
4001 : GEN u, v, d;
4002 :
4003 715 : while (vx != vy)
4004 : {
4005 0 : if (varncmp(vx,vy) > 0)
4006 : {
4007 0 : d = (vx == NO_VARIABLE)? ginv(x): gred_rfrac_simple(gen_1, x);
4008 0 : return scalarpol(d, vy);
4009 : }
4010 0 : if (lg(x)!=3) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
4011 0 : x = gel(x,2); vx = gvar(x);
4012 : }
4013 715 : av = avma; d = subresext_i(x,y,&u,&v/*junk*/);
4014 715 : if (gequal0(d)) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
4015 715 : d = gdiv(u,d);
4016 715 : if (typ(d) != t_POL || varn(d) != vy) d = scalarpol(d, vy);
4017 715 : return gc_upto(av, d);
4018 : }
4019 :
4020 : /*Assume x is a polynomial and y is not */
4021 : static GEN
4022 96 : scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
4023 : {
4024 96 : long vx = varn(x);
4025 96 : int xis0 = signe(x)==0, yis0 = gequal0(y);
4026 96 : if (xis0 && yis0) { *U = *V = pol_0(vx); return pol_0(vx); }
4027 72 : if (yis0) { *U=pol_1(vx); *V = pol_0(vx); return RgX_copy(x);}
4028 48 : *U=pol_0(vx); *V= ginv(y); return pol_1(vx);
4029 : }
4030 : /* Assume x==0, y!=0 */
4031 : static GEN
4032 54 : zero_bezout(GEN y, GEN *U, GEN *V)
4033 : {
4034 54 : *U=gen_0; *V = ginv(y); return gen_1;
4035 : }
4036 :
4037 : GEN
4038 355 : gbezout(GEN x, GEN y, GEN *u, GEN *v)
4039 : {
4040 355 : long tx=typ(x), ty=typ(y), vx;
4041 355 : if (tx == t_INT && ty == t_INT) return bezout(x,y,u,v);
4042 326 : if (tx != t_POL)
4043 : {
4044 120 : if (ty == t_POL)
4045 48 : return scalar_bezout(y,x,v,u);
4046 : else
4047 : {
4048 72 : int xis0 = gequal0(x), yis0 = gequal0(y);
4049 72 : if (xis0 && yis0) { *u = *v = gen_0; return gen_0; }
4050 54 : if (xis0) return zero_bezout(y,u,v);
4051 36 : else return zero_bezout(x,v,u);
4052 : }
4053 : }
4054 206 : else if (ty != t_POL) return scalar_bezout(x,y,u,v);
4055 158 : vx = varn(x);
4056 158 : if (vx != varn(y))
4057 0 : return varncmp(vx, varn(y)) < 0? scalar_bezout(x,y,u,v)
4058 0 : : scalar_bezout(y,x,v,u);
4059 158 : return RgX_extgcd(x,y,u,v);
4060 : }
4061 :
4062 : GEN
4063 355 : gcdext0(GEN x, GEN y)
4064 : {
4065 355 : GEN z=cgetg(4,t_VEC);
4066 355 : gel(z,3) = gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
4067 355 : return z;
4068 : }
4069 :
4070 : /*******************************************************************/
4071 : /* */
4072 : /* GENERIC (modular) INVERSE */
4073 : /* */
4074 : /*******************************************************************/
4075 :
4076 : GEN
4077 29845 : ginvmod(GEN x, GEN y)
4078 : {
4079 29845 : long tx=typ(x);
4080 :
4081 29845 : switch(typ(y))
4082 : {
4083 29845 : case t_POL:
4084 29845 : if (tx==t_POL) return RgXQ_inv(x,y);
4085 11317 : if (is_scalar_t(tx)) return ginv(x);
4086 0 : break;
4087 0 : case t_INT:
4088 0 : if (tx==t_INT) return Fp_inv(x,y);
4089 0 : if (tx==t_POL) return gen_0;
4090 : }
4091 0 : pari_err_TYPE2("ginvmod",x,y);
4092 : return NULL; /* LCOV_EXCL_LINE */
4093 : }
4094 :
4095 : /***********************************************************************/
4096 : /** **/
4097 : /** NEWTON POLYGON **/
4098 : /** **/
4099 : /***********************************************************************/
4100 :
4101 : /* assume leading coeff of x is nonzero */
4102 : GEN
4103 20 : newtonpoly(GEN x, GEN p)
4104 : {
4105 20 : pari_sp av = avma;
4106 : long n, ind, a, b;
4107 : GEN y, vval;
4108 :
4109 20 : if (typ(x) != t_POL) pari_err_TYPE("newtonpoly",x);
4110 20 : n = degpol(x); if (n<=0) return cgetg(1,t_VEC);
4111 20 : vval = new_chunk(n+1);
4112 20 : y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
4113 120 : for (a = 0; a <= n; a++) vval[a] = gvaluation(gel(x,a),p);
4114 30 : for (a = 0, ind = 1; a < n; a++)
4115 : {
4116 30 : if (vval[a] != LONG_MAX) break;
4117 10 : gel(y,ind++) = mkoo();
4118 : }
4119 60 : for (b = a+1; b <= n; a = b, b = a+1)
4120 : {
4121 : long u1, u2, c;
4122 50 : while (vval[b] == LONG_MAX) b++;
4123 40 : u1 = vval[a] - vval[b];
4124 40 : u2 = b - a;
4125 110 : for (c = b+1; c <= n; c++)
4126 : {
4127 : long r1, r2;
4128 70 : if (vval[c] == LONG_MAX) continue;
4129 50 : r1 = vval[a] - vval[c];
4130 50 : r2 = c - a;
4131 50 : if (u1*r2 <= u2*r1) { u1 = r1; u2 = r2; b = c; }
4132 : }
4133 110 : while (ind <= b) gel(y,ind++) = sstoQ(u1,u2);
4134 : }
4135 20 : stackdummy((pari_sp)vval, av); return y;
4136 : }
4137 :
4138 : static GEN
4139 235117 : RgXQ_mul_FpXQ(GEN x, GEN y, GEN T, GEN p)
4140 : {
4141 235117 : pari_sp av = avma;
4142 : GEN r;
4143 235117 : if (lgefint(p) == 3)
4144 : {
4145 130538 : ulong pp = uel(p, 2);
4146 130538 : r = Flx_to_ZX_inplace(Flxq_mul(RgX_to_Flx(x, pp),
4147 : RgX_to_Flx(y, pp), RgX_to_Flx(T, pp), pp));
4148 : }
4149 : else
4150 104579 : r = FpXQ_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), RgX_to_FpX(T, p), p);
4151 235117 : return gc_upto(av, FpX_to_mod(r, p));
4152 : }
4153 :
4154 : static GEN
4155 10 : RgXQ_sqr_FpXQ(GEN x, GEN y, GEN p)
4156 : {
4157 10 : pari_sp av = avma;
4158 : GEN r;
4159 10 : if (lgefint(p) == 3)
4160 : {
4161 5 : ulong pp = uel(p, 2);
4162 5 : r = Flx_to_ZX_inplace(Flxq_sqr(RgX_to_Flx(x, pp),
4163 : RgX_to_Flx(y, pp), pp));
4164 : }
4165 : else
4166 5 : r = FpXQ_sqr(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
4167 10 : return gc_upto(av, FpX_to_mod(r, p));
4168 : }
4169 :
4170 : static GEN
4171 10327 : RgXQ_inv_FpXQ(GEN x, GEN y, GEN p)
4172 : {
4173 10327 : pari_sp av = avma;
4174 : GEN r;
4175 10327 : if (lgefint(p) == 3)
4176 : {
4177 5212 : ulong pp = uel(p, 2);
4178 5212 : r = Flx_to_ZX_inplace(Flxq_inv(RgX_to_Flx(x, pp),
4179 : RgX_to_Flx(y, pp), pp));
4180 : }
4181 : else
4182 5115 : r = FpXQ_inv(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
4183 10327 : return gc_upto(av, FpX_to_mod(r, p));
4184 : }
4185 :
4186 : static GEN
4187 329 : RgXQ_mul_FpXQXQ(GEN x, GEN y, GEN S, GEN pol, GEN p)
4188 : {
4189 329 : pari_sp av = avma;
4190 : GEN r;
4191 329 : GEN T = RgX_to_FpX(pol, p);
4192 329 : if (signe(T)==0) pari_err_OP("*",x,y);
4193 329 : if (lgefint(p) == 3)
4194 : {
4195 203 : ulong pp = uel(p, 2);
4196 203 : GEN Tp = ZX_to_Flx(T, pp);
4197 203 : r = FlxX_to_ZXX(FlxqXQ_mul(RgX_to_FlxqX(x, Tp, pp),
4198 : RgX_to_FlxqX(y, Tp, pp),
4199 : RgX_to_FlxqX(S, Tp, pp), Tp, pp));
4200 : }
4201 : else
4202 126 : r = FpXQXQ_mul(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p),
4203 : RgX_to_FpXQX(S, T, p), T, p);
4204 329 : return gc_upto(av, FpXQX_to_mod(r, T, p));
4205 : }
4206 :
4207 : static GEN
4208 0 : RgXQ_sqr_FpXQXQ(GEN x, GEN y, GEN pol, GEN p)
4209 : {
4210 0 : pari_sp av = avma;
4211 : GEN r;
4212 0 : GEN T = RgX_to_FpX(pol, p);
4213 0 : if (signe(T)==0) pari_err_OP("*",x,x);
4214 0 : if (lgefint(p) == 3)
4215 : {
4216 0 : ulong pp = uel(p, 2);
4217 0 : GEN Tp = ZX_to_Flx(T, pp);
4218 0 : r = FlxX_to_ZXX(FlxqXQ_sqr(RgX_to_FlxqX(x, Tp, pp),
4219 : RgX_to_FlxqX(y, Tp, pp), Tp, pp));
4220 : }
4221 : else
4222 0 : r = FpXQXQ_sqr(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
4223 0 : return gc_upto(av, FpXQX_to_mod(r, T, p));
4224 : }
4225 :
4226 : static GEN
4227 5 : RgXQ_inv_FpXQXQ(GEN x, GEN y, GEN pol, GEN p)
4228 : {
4229 5 : pari_sp av = avma;
4230 : GEN r;
4231 5 : GEN T = RgX_to_FpX(pol, p);
4232 5 : if (signe(T)==0) pari_err_OP("^",x,gen_m1);
4233 5 : if (lgefint(p) == 3)
4234 : {
4235 5 : ulong pp = uel(p, 2);
4236 5 : GEN Tp = ZX_to_Flx(T, pp);
4237 5 : r = FlxX_to_ZXX(FlxqXQ_inv(RgX_to_FlxqX(x, Tp, pp),
4238 : RgX_to_FlxqX(y, Tp, pp), Tp, pp));
4239 : }
4240 : else
4241 0 : r = FpXQXQ_inv(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
4242 5 : return gc_upto(av, FpXQX_to_mod(r, T, p));
4243 : }
4244 :
4245 : static GEN
4246 1298917 : RgXQ_mul_fast(GEN x, GEN y, GEN T)
4247 : {
4248 : GEN p, pol;
4249 : long pa;
4250 1298917 : long t = RgX_type3(x,y,T, &p,&pol,&pa);
4251 1298917 : switch(t)
4252 : {
4253 484469 : case t_INT: return ZX_is_monic(T) ? ZXQ_mul(x,y,T): NULL;
4254 545699 : case t_FRAC: return RgX_is_ZX(T) && ZX_is_monic(T) ? QXQ_mul(x,y,T): NULL;
4255 88 : case t_FFELT: return FFXQ_mul(x, y, T, pol);
4256 235117 : case t_INTMOD: return RgXQ_mul_FpXQ(x, y, T, p);
4257 329 : case RgX_type_code(t_POLMOD, t_INTMOD):
4258 329 : return RgXQ_mul_FpXQXQ(x, y, T, pol, p);
4259 33215 : default: return NULL;
4260 : }
4261 : }
4262 :
4263 : GEN
4264 1298917 : RgXQ_mul(GEN x, GEN y, GEN T)
4265 : {
4266 1298917 : GEN z = RgXQ_mul_fast(x, y, T);
4267 1298917 : if (!z) z = RgX_rem(RgX_mul(x, y), T);
4268 1298917 : return z;
4269 : }
4270 :
4271 : static GEN
4272 396985 : RgXQ_sqr_fast(GEN x, GEN T)
4273 : {
4274 : GEN p, pol;
4275 : long pa;
4276 396985 : long t = RgX_type2(x, T, &p,&pol,&pa);
4277 396985 : switch(t)
4278 : {
4279 93520 : case t_INT: return ZX_is_monic(T) ? ZXQ_sqr(x,T): NULL;
4280 297605 : case t_FRAC: return RgX_is_ZX(T) && ZX_is_monic(T) ? QXQ_sqr(x,T): NULL;
4281 6 : case t_FFELT: return FFXQ_sqr(x, T, pol);
4282 10 : case t_INTMOD: return RgXQ_sqr_FpXQ(x, T, p);
4283 0 : case RgX_type_code(t_POLMOD, t_INTMOD):
4284 0 : return RgXQ_sqr_FpXQXQ(x, T, pol, p);
4285 5844 : default: return NULL;
4286 : }
4287 : }
4288 :
4289 : GEN
4290 396985 : RgXQ_sqr(GEN x, GEN T)
4291 : {
4292 396985 : GEN z = RgXQ_sqr_fast(x, T);
4293 396985 : if (!z) z = RgX_rem(RgX_sqr(x), T);
4294 396985 : return z;
4295 : }
4296 :
4297 : static GEN
4298 115506 : RgXQ_inv_fast(GEN x, GEN y)
4299 : {
4300 : GEN p, pol;
4301 : long pa;
4302 115506 : long t = RgX_type2(x,y, &p,&pol,&pa);
4303 115506 : switch(t)
4304 : {
4305 77115 : case t_INT: return QXQ_inv(x,y);
4306 27339 : case t_FRAC: return RgX_is_ZX(y)? QXQ_inv(x,y): NULL;
4307 10 : case t_FFELT: return FFXQ_inv(x, y, pol);
4308 10327 : case t_INTMOD: return RgXQ_inv_FpXQ(x, y, p);
4309 5 : case RgX_type_code(t_POLMOD, t_INTMOD):
4310 5 : return RgXQ_inv_FpXQXQ(x, y, pol, p);
4311 710 : default: return NULL;
4312 : }
4313 : }
4314 :
4315 : GEN
4316 115506 : RgXQ_inv(GEN x, GEN y)
4317 : {
4318 115506 : GEN z = RgXQ_inv_fast(x, y);
4319 115495 : if (!z) z = RgXQ_inv_i(x, y);
4320 115495 : return z;
4321 : }
|
