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