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