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