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