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