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 127982 : polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N)
33 : {
34 127982 : long dP=degpol(P), i, k, m;
35 : pari_sp av1, av2;
36 : GEN s,y,P_lead;
37 :
38 127982 : if (n<0) pari_err_IMPL("polsym of a negative n");
39 127982 : if (typ(P) != t_POL) pari_err_TYPE("polsym",P);
40 127982 : if (!signe(P)) pari_err_ROOTS0("polsym");
41 127982 : y = cgetg(n+2,t_COL);
42 127982 : 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 114745 : m = 1;
51 114745 : gel(y,1) = stoi(dP);
52 : }
53 127983 : P += 2; /* strip codewords */
54 :
55 127983 : P_lead = gel(P,dP); if (gequal1(P_lead)) P_lead = NULL;
56 127983 : 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 380104 : for (k=m; k<=n; k++)
62 : {
63 252120 : av1 = avma; s = (dP>=k)? gmulsg(k,gel(P,dP-k)): gen_0;
64 725537 : for (i=1; i<k && i<=dP; i++)
65 473433 : s = gadd(s, gmul(gel(y,k-i+1),gel(P,dP-i)));
66 252104 : 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 233344 : 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 233344 : if (P_lead) s = gdiv(s, P_lead);
78 252104 : av2 = avma; gel(y,k+1) = gerepile(av1,av2, gneg(s));
79 : }
80 127984 : return y;
81 : }
82 :
83 : GEN
84 110567 : polsym(GEN x, long n)
85 : {
86 110567 : 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 85008097 : centermodii(GEN x, GEN p, GEN po2)
92 : {
93 85008097 : GEN y = remii(x, p);
94 85181398 : switch(signe(y))
95 : {
96 10524046 : case 0: break;
97 51485772 : case 1: if (po2 && abscmpii(y,po2) > 0) y = subii(y, p);
98 51346177 : break;
99 23504771 : case -1: if (!po2 || abscmpii(y,po2) > 0) y = addii(y, p);
100 23304761 : break;
101 : }
102 84841793 : 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 12648404 : centermod_i(GEN x, GEN p, GEN ps2)
116 : {
117 : long i, lx;
118 : pari_sp av;
119 : GEN y;
120 :
121 12648404 : if (!ps2) ps2 = shifti(p,-1);
122 12646122 : switch(typ(x))
123 : {
124 1123830 : case t_INT: return centermodii(x,p,ps2);
125 :
126 6143289 : case t_POL: lx = lg(x);
127 6143289 : y = cgetg(lx,t_POL); y[1] = x[1];
128 42161020 : for (i=2; i<lx; i++)
129 : {
130 36015603 : av = avma;
131 36015603 : gel(y,i) = gerepileuptoint(av, centermodii(gel(x,i),p,ps2));
132 : }
133 6145417 : return normalizepol_lg(y, lx);
134 :
135 5377932 : case t_COL: lx = lg(x);
136 5377932 : y = cgetg(lx,t_COL);
137 25451437 : for (i=1; i<lx; i++) gel(y,i) = centermodii(gel(x,i),p,ps2);
138 5376867 : 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 9069201 : 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 9121 : RgX_cmbf(GEN p, long i, GEN BLOC, GEN Lmod, GEN Lfac, GEN *F)
228 : {
229 : GEN q;
230 9121 : if (i == lg(Lmod)) return 0;
231 4704 : if (RgX_cmbf(p, i+1, BLOC, Lmod, Lfac, F) && p) return 1;
232 4501 : if (!gel(Lmod,i)) return 0;
233 4487 : p = p? RgX_mul(p, gel(Lmod,i)): gel(Lmod,i);
234 4487 : q = RgV_to_RgX(RgX_digits(p, BLOC), varn(*F));
235 4487 : if (degpol(q))
236 : {
237 4130 : GEN R, Q = RgX_divrem(*F, q, &R);
238 4130 : if (signe(R)==0) { vectrunc_append(Lfac, q); *F = Q; return 1; }
239 : }
240 4158 : if (RgX_cmbf(p, i+1, BLOC, Lmod, Lfac, F)) { gel(Lmod,i) = NULL; return 1; }
241 3892 : return 0;
242 : }
243 :
244 : static GEN factor_domain(GEN x, GEN flag);
245 :
246 : static GEN
247 434 : ok_bloc(GEN f, GEN BLOC, ulong c)
248 : {
249 434 : GEN F = poleval(f, BLOC);
250 434 : return issquarefree(c ? gmul(F,mkintmodu(1,c)): F)? F: NULL;
251 : }
252 : static GEN
253 112 : random_FpX_monic(long n, long v, GEN p)
254 : {
255 112 : long i, d = n + 2;
256 112 : GEN y = cgetg(d + 1, t_POL); y[1] = evalsigne(1) | evalvarn(v);
257 364 : for (i = 2; i < d; i++) gel(y,i) = randomi(p);
258 112 : gel(y,i) = gen_1; return y;
259 : }
260 : static GEN
261 273 : RgXY_factor_squarefree(GEN f, GEN dom)
262 : {
263 273 : pari_sp av = avma;
264 273 : ulong i, c = itou_or_0(residual_characteristic(f));
265 273 : long vy = gvar2(f), val = RgX_valrem(f, &f), n = RgXY_degreex(f);
266 273 : GEN y, Lmod, F = NULL, BLOC = NULL, Lfac = coltrunc_init(degpol(f)+2);
267 273 : GEN gc = c? utoipos(c): NULL;
268 273 : if (val)
269 : {
270 35 : GEN x = pol_x(varn(f));
271 35 : if (dom)
272 : {
273 14 : GEN c = Rg_get_1(dom);
274 14 : if (typ(c) != t_INT) x = RgX_Rg_mul(x, c);
275 : }
276 35 : vectrunc_append(Lfac, x); if (!degpol(f)) return Lfac;
277 : }
278 259 : y = pol_x(vy);
279 : for(;;)
280 : {
281 322 : for (i = 0; !c || i < c; i++)
282 : {
283 322 : BLOC = gpowgs(gaddgs(y, i), n+1);
284 322 : if ((F = ok_bloc(f, BLOC, c))) break;
285 147 : if (c)
286 : {
287 112 : BLOC = random_FpX_monic(n, vy, gc);
288 112 : if ((F = ok_bloc(f, BLOC, c))) break;
289 : }
290 : }
291 259 : if (!c || i < c) break;
292 0 : n++;
293 : }
294 259 : if (DEBUGLEVEL >= 2)
295 0 : err_printf("bifactor: bloc:(x+%ld)^%ld, deg f=%ld\n",i,n,RgXY_degreex(f));
296 259 : Lmod = gel(factor_domain(F,dom),1);
297 259 : if (DEBUGLEVEL >= 2)
298 0 : err_printf("bifactor: %ld local factors\n",lg(Lmod)-1);
299 259 : (void)RgX_cmbf(NULL, 1, BLOC, Lmod, Lfac, &f);
300 259 : if (degpol(f)) vectrunc_append(Lfac, f);
301 259 : return gerepilecopy(av, Lfac);
302 : }
303 :
304 : static GEN
305 245 : FE_matconcat(GEN F, GEN E, long l)
306 : {
307 245 : setlg(E,l); E = shallowconcat1(E);
308 245 : setlg(F,l); F = shallowconcat1(F); return mkmat2(F,E);
309 : }
310 :
311 : static int
312 385 : gen_cmp_RgXY(void *data, GEN x, GEN y)
313 : {
314 385 : long vx = varn(x), vy = varn(y);
315 385 : return (vx == vy)? gen_cmp_RgX(data, x, y): -varncmp(vx, vy);
316 : }
317 : static GEN
318 245 : RgXY_factor(GEN f, GEN dom)
319 : {
320 245 : pari_sp av = avma;
321 : GEN C, F, E, cf, V;
322 : long i, j, l;
323 245 : if (dom) { GEN c = Rg_get_1(dom); if (typ(c) != t_INT) f = RgX_Rg_mul(f,c); }
324 245 : cf = content(f);
325 245 : V = RgXY_squff(gdiv(f, cf)); l = lg(V);
326 245 : C = factor_domain(cf, dom);
327 245 : F = cgetg(l+1, t_VEC); gel(F,1) = gel(C,1);
328 245 : E = cgetg(l+1, t_VEC); gel(E,1) = gel(C,2);
329 756 : for (i=1, j=2; i < l; i++)
330 : {
331 511 : GEN v = gel(V,i);
332 511 : if (degpol(v))
333 : {
334 273 : gel(F,j) = v = RgXY_factor_squarefree(v, dom);
335 273 : gel(E,j) = const_col(lg(v)-1, utoipos(i));
336 273 : j++;
337 : }
338 : }
339 245 : f = FE_matconcat(F,E,j);
340 245 : (void)sort_factor(f,(void*)cmp_universal, &gen_cmp_RgXY);
341 245 : return gerepilecopy(av, f);
342 : }
343 :
344 : /***********************************************************************/
345 : /** **/
346 : /** FACTORIZATION **/
347 : /** **/
348 : /***********************************************************************/
349 : static long RgX_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var);
350 : #define assign_or_fail(x,y) { GEN __x = x;\
351 : if (!*y) *y=__x; else if (!gequal(__x,*y)) return 0;\
352 : }
353 : #define update_prec(x,y) { long __x = x; if (__x < *y) *y=__x; }
354 :
355 : static const long tsh = 6;
356 : #define code(t1,t2) ((t1 << 6) | t2)
357 : void
358 89456 : RgX_type_decode(long x, long *t1, long *t2)
359 : {
360 89456 : *t1 = x >> tsh;
361 89456 : *t2 = (x & ((1L<<tsh)-1));
362 89456 : }
363 : int
364 9112494 : RgX_type_is_composite(long t) { return t >= tsh; }
365 :
366 : static int
367 2139029264 : settype(GEN c, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
368 : {
369 : long j;
370 2139029264 : switch(typ(c))
371 : {
372 1603450618 : case t_INT:
373 1603450618 : break;
374 25922044 : case t_FRAC:
375 25922044 : t[1]=1; break;
376 : break;
377 234776474 : case t_REAL:
378 234776474 : update_prec(precision(c), pa);
379 234775668 : t[2]=1; break;
380 23503244 : case t_INTMOD:
381 23503244 : assign_or_fail(gel(c,1),p);
382 23503244 : t[3]=1; break;
383 1870168 : case t_FFELT:
384 1870168 : if (!*ff) *ff=c; else if (!FF_samefield(c,*ff)) return 0;
385 1870168 : assign_or_fail(FF_p_i(c),p);
386 1870168 : t[5]=1; break;
387 215557010 : case t_COMPLEX:
388 646668530 : for (j=1; j<=2; j++)
389 : {
390 431111641 : GEN d = gel(c,j);
391 431111641 : switch(typ(d))
392 : {
393 1643416 : case t_INT: case t_FRAC:
394 1643416 : if (!*t2) *t2 = t_COMPLEX;
395 1643416 : t[1]=1; break;
396 429468190 : case t_REAL:
397 429468190 : update_prec(precision(d), pa);
398 429468076 : if (!*t2) *t2 = t_COMPLEX;
399 429468076 : t[2]=1; break;
400 14 : case t_INTMOD:
401 14 : assign_or_fail(gel(d,1),p);
402 14 : if (!signe(*p) || mod4(*p) != 3) return 0;
403 7 : if (!*t2) *t2 = t_COMPLEX;
404 7 : t[3]=1; break;
405 21 : case t_PADIC:
406 21 : update_prec(precp(d)+valp(d), pa);
407 21 : assign_or_fail(gel(d,2),p);
408 21 : if (!*t2) *t2 = t_COMPLEX;
409 21 : t[7]=1; break;
410 0 : default: return 0;
411 : }
412 : }
413 215556889 : if (!t[2]) assign_or_fail(mkpoln(3, gen_1,gen_0,gen_1), pol); /*x^2+1*/
414 215556889 : break;
415 3296524 : case t_PADIC:
416 3296524 : update_prec(precp(c)+valp(c), pa);
417 3296524 : assign_or_fail(gel(c,2),p);
418 3296524 : t[7]=1; break;
419 1491 : case t_QUAD:
420 1491 : assign_or_fail(gel(c,1),pol);
421 4473 : for (j=2; j<=3; j++)
422 : {
423 2982 : GEN d = gel(c,j);
424 2982 : switch(typ(d))
425 : {
426 2947 : case t_INT: case t_FRAC:
427 2947 : t[8]=1; break;
428 28 : case t_INTMOD:
429 28 : assign_or_fail(gel(d,1),p);
430 28 : if (*t2 != t_POLMOD) *t2 = t_QUAD;
431 28 : t[3]=1; break;
432 7 : case t_PADIC:
433 7 : update_prec(precp(d)+valp(d), pa);
434 7 : assign_or_fail(gel(d,2),p);
435 7 : if (*t2 != t_POLMOD) *t2 = t_QUAD;
436 7 : t[7]=1; break;
437 0 : default: return 0;
438 : }
439 : }
440 1491 : break;
441 3694434 : case t_POLMOD:
442 3694434 : assign_or_fail(gel(c,1),pol);
443 3694219 : if (typ(gel(c,2))==t_POL && varn(gel(c,2))!=varn(gel(c,1))) return 0;
444 11055280 : for (j=1; j<=2; j++)
445 : {
446 : GEN pbis, polbis;
447 : long pabis;
448 7376858 : *t2 = t_POLMOD;
449 7376858 : switch(Rg_type(gel(c,j),&pbis,&polbis,&pabis))
450 : {
451 3944226 : case t_INT: break;
452 844426 : case t_FRAC: t[1]=1; break;
453 2576875 : case t_INTMOD: t[3]=1; break;
454 7 : case t_PADIC: t[7]=1; update_prec(pabis,pa); break;
455 11337 : default: return 0;
456 : }
457 7365534 : if (pbis) assign_or_fail(pbis,p);
458 7365534 : if (polbis) assign_or_fail(polbis,pol);
459 : }
460 3678422 : break;
461 5325148 : case t_RFRAC: t[10] = 1;
462 5325148 : if (!settype(gel(c,1),t,p,pol,pa,ff,t2,var)) return 0;
463 5325148 : c = gel(c,2); /* fall through */
464 26955556 : case t_POL: t[10] = 1;
465 26955556 : if (!RgX_settype(c,t,p,pol,pa,ff,t2,var)) return 0;
466 26968484 : if (*var == NO_VARIABLE) { *var = varn(c); break; }
467 : /* if more than one free var, ensure varn() == *var fails. FIXME: should
468 : * keep the list of all variables, later t_POLMOD may cancel them */
469 19264296 : if (*var != varn(c)) *var = MAXVARN+1;
470 19264296 : break;
471 1701 : default: return 0;
472 : }
473 2139023552 : return 1;
474 : }
475 : /* t[0] unused. Other values, if set, indicate a coefficient of type
476 : * t[1] : t_FRAC
477 : * t[2] : t_REAL
478 : * t[3] : t_INTMOD
479 : * t[4] : Unused
480 : * t[5] : t_FFELT
481 : * t[6] : Unused
482 : * t[7] : t_PADIC
483 : * t[8] : t_QUAD of rationals (t_INT/t_FRAC)
484 : * t[9]: Unused
485 : * t[10]: t_POL (recursive factorisation) */
486 : /* if t2 != 0: t_POLMOD/t_QUAD/t_COMPLEX of modular (t_INTMOD/t_PADIC,
487 : * given by t) */
488 : static long
489 155400011 : choosetype(long *t, long t2, GEN ff, GEN *pol, long var)
490 : {
491 155400011 : if (t[10] && (!*pol || var!=varn(*pol))) return t_POL;
492 147707717 : if (t2) /* polmod/quad/complex of intmod/padic */
493 : {
494 20875416 : if (t[2] && (t[3]||t[7])) return 0;
495 20875416 : if (t[3]) return code(t2,t_INTMOD);
496 20846163 : if (t[7]) return code(t2,t_PADIC);
497 20846114 : if (t[2]) return t_COMPLEX;
498 487747 : if (t[1]) return code(t2,t_FRAC);
499 223391 : return code(t2,t_INT);
500 : }
501 126832301 : if (t[5]) /* ffelt */
502 : {
503 219745 : if (t[2]||t[8]||t[9]) return 0;
504 219745 : *pol=ff; return t_FFELT;
505 : }
506 126612556 : if (t[2]) /* inexact, real */
507 : {
508 41580476 : if (t[3]||t[7]||t[9]) return 0;
509 41580508 : return t_REAL;
510 : }
511 85032080 : if (t[10]) return t_POL;
512 85032080 : if (t[8]) return code(t_QUAD,t_INT);
513 85031443 : if (t[3]) return t_INTMOD;
514 81339728 : if (t[7]) return t_PADIC;
515 80743370 : if (t[1]) return t_FRAC;
516 74884609 : return t_INT;
517 : }
518 :
519 : static long
520 229484839 : RgX_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
521 : {
522 229484839 : long i, lx = lg(x);
523 883317205 : for (i=2; i<lx; i++)
524 653854533 : if (!settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
525 229462672 : return 1;
526 : }
527 :
528 : static long
529 243040370 : RgC_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
530 : {
531 243040370 : long i, l = lg(x);
532 1722290471 : for (i = 1; i<l; i++)
533 1479252631 : if (!settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
534 243037840 : return 1;
535 : }
536 :
537 : static long
538 49477089 : RgM_settype(GEN x, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2, long *var)
539 : {
540 49477089 : long i, l = lg(x);
541 250928309 : for (i = 1; i < l; i++)
542 201453673 : if (!RgC_settype(gel(x,i),t,p,pol,pa,ff,t2,var)) return 0;
543 49474636 : return 1;
544 : }
545 :
546 : long
547 16489358 : Rg_type(GEN x, GEN *p, GEN *pol, long *pa)
548 : {
549 16489358 : long t[] = {0,0,0,0,0,0,0,0,0,0,0};
550 16489358 : long t2 = 0, var = NO_VARIABLE;
551 16489358 : GEN ff = NULL;
552 16489358 : *p = *pol = NULL; *pa = LONG_MAX;
553 16489358 : switch(typ(x))
554 : {
555 602080 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_FFELT:
556 : case t_COMPLEX: case t_PADIC: case t_QUAD:
557 602080 : if (!settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
558 602082 : break;
559 15304640 : case t_POL: case t_SER:
560 15304640 : if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
561 15304650 : break;
562 14 : case t_VEC: case t_COL:
563 14 : if(!RgC_settype(x, t, p, pol, pa, &ff, &t2, &var)) return 0;
564 14 : break;
565 126 : case t_MAT:
566 126 : if(!RgM_settype(x, t, p, pol, pa, &ff, &t2, &var)) return 0;
567 126 : break;
568 582498 : default: return 0;
569 : }
570 15906872 : return choosetype(t,t2,ff,pol,var);
571 : }
572 :
573 : long
574 1256868 : RgX_type(GEN x, GEN *p, GEN *pol, long *pa)
575 : {
576 1256868 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
577 1256868 : long t2 = 0, var = NO_VARIABLE;
578 1256868 : GEN ff = NULL;
579 1256868 : *p = *pol = NULL; *pa = LONG_MAX;
580 1256868 : if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
581 1256820 : return choosetype(t,t2,ff,pol,var);
582 : }
583 :
584 : long
585 294 : RgX_Rg_type(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
586 : {
587 294 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
588 294 : long t2 = 0, var = NO_VARIABLE;
589 294 : GEN ff = NULL;
590 294 : *p = *pol = NULL; *pa = LONG_MAX;
591 294 : if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
592 294 : if (!settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
593 294 : return choosetype(t,t2,ff,pol,var);
594 : }
595 :
596 : long
597 91722912 : RgX_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
598 : {
599 91722912 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
600 91722912 : long t2 = 0, var = NO_VARIABLE;
601 91722912 : GEN ff = NULL;
602 91722912 : *p = *pol = NULL; *pa = LONG_MAX;
603 183434474 : if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
604 91723922 : !RgX_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
605 91711444 : return choosetype(t,t2,ff,pol,var);
606 : }
607 :
608 : long
609 839214 : RgX_type3(GEN x, GEN y, GEN z, GEN *p, GEN *pol, long *pa)
610 : {
611 839214 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
612 839214 : long t2 = 0, var = NO_VARIABLE;
613 839214 : GEN ff = NULL;
614 839214 : *p = *pol = NULL; *pa = LONG_MAX;
615 1675937 : if (!RgX_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
616 1673445 : !RgX_settype(y,t,p,pol,pa,&ff,&t2,&var) ||
617 839215 : !RgX_settype(z,t,p,pol,pa,&ff,&t2,&var)) return 0;
618 836720 : return choosetype(t,t2,ff,pol,var);
619 : }
620 :
621 : long
622 316430 : RgM_type(GEN x, GEN *p, GEN *pol, long *pa)
623 : {
624 316430 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
625 316430 : long t2 = 0, var = NO_VARIABLE;
626 316430 : GEN ff = NULL;
627 316430 : *p = *pol = NULL; *pa = LONG_MAX;
628 316430 : if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
629 315459 : return choosetype(t,t2,ff,pol,var);
630 : }
631 :
632 : long
633 70 : RgV_type(GEN x, GEN *p, GEN *pol, long *pa)
634 : {
635 70 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
636 70 : long t2 = 0, var = NO_VARIABLE;
637 70 : GEN ff = NULL;
638 70 : *p = *pol = NULL; *pa = LONG_MAX;
639 70 : if (!RgC_settype(x,t,p,pol,pa,&ff,&t2,&var)) return 0;
640 70 : return choosetype(t,t2,ff,pol,var);
641 : }
642 :
643 : long
644 203 : RgV_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
645 : {
646 203 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
647 203 : long t2 = 0, var = NO_VARIABLE;
648 203 : GEN ff = NULL;
649 203 : *p = *pol = NULL; *pa = LONG_MAX;
650 406 : if (!RgC_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
651 203 : !RgC_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
652 203 : return choosetype(t,t2,ff,pol,var);
653 : }
654 :
655 : long
656 41587557 : RgM_RgC_type(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
657 : {
658 41587557 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
659 41587557 : long t2 = 0, var = NO_VARIABLE;
660 41587557 : GEN ff = NULL;
661 41587557 : *p = *pol = NULL; *pa = LONG_MAX;
662 83174773 : if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
663 41588439 : !RgC_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
664 41586417 : return choosetype(t,t2,ff,pol,var);
665 : }
666 :
667 : long
668 3786675 : RgM_type2(GEN x, GEN y, GEN *p, GEN *pol, long *pa)
669 : {
670 3786675 : long t[] = {0,0,0,0,0,0,0,0,0,0,0,0};
671 3786675 : long t2 = 0, var = NO_VARIABLE;
672 3786675 : GEN ff = NULL;
673 3786675 : *p = *pol = NULL; *pa = LONG_MAX;
674 7572893 : if (!RgM_settype(x,t,p,pol,pa,&ff,&t2,&var) ||
675 3786801 : !RgM_settype(y,t,p,pol,pa,&ff,&t2,&var)) return 0;
676 3786113 : return choosetype(t,t2,ff,pol,var);
677 : }
678 :
679 : GEN
680 59326 : factor0(GEN x, GEN flag)
681 : {
682 : ulong B;
683 59326 : long tx = typ(x);
684 59326 : if (!flag) return factor(x);
685 224 : if ((tx != t_INT && tx!=t_FRAC) || typ(flag) != t_INT)
686 175 : return factor_domain(x, flag);
687 49 : if (signe(flag) < 0) pari_err_FLAG("factor");
688 49 : switch(lgefint(flag))
689 : {
690 7 : case 2: B = 0; break;
691 42 : case 3: B = flag[2]; break;
692 0 : default: pari_err_OVERFLOW("factor [large prime bound]");
693 : return NULL; /*LCOV_EXCL_LINE*/
694 : }
695 49 : return boundfact(x, B);
696 : }
697 :
698 : GEN
699 134185 : deg1_from_roots(GEN L, long v)
700 : {
701 134185 : long i, l = lg(L);
702 134185 : GEN z = cgetg(l,t_COL);
703 411278 : for (i=1; i<l; i++)
704 277094 : gel(z,i) = deg1pol_shallow(gen_1, gneg(gel(L,i)), v);
705 134184 : return z;
706 : }
707 : GEN
708 63852 : roots_from_deg1(GEN x)
709 : {
710 63852 : long i,l = lg(x);
711 63852 : GEN r = cgetg(l,t_VEC);
712 392640 : for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = gneg(gel(P,2)); }
713 63853 : return r;
714 : }
715 :
716 : static GEN
717 42 : gauss_factor_p(GEN p)
718 : {
719 42 : GEN a, b; (void)cornacchia(gen_1, p, &a,&b);
720 42 : return mkcomplex(a, b);
721 : }
722 :
723 : static GEN
724 49 : gauss_primpart(GEN x, GEN *c)
725 : {
726 49 : GEN a = real_i(x), b = imag_i(x), n = gcdii(a, b);
727 49 : *c = n; if (n == gen_1) return x;
728 49 : retmkcomplex(diviiexact(a,n), diviiexact(b,n));
729 : }
730 :
731 : static GEN
732 70 : gauss_primpart_try(GEN x, GEN c)
733 : {
734 : GEN r, y;
735 70 : if (typ(x) == t_INT)
736 : {
737 42 : y = dvmdii(x, c, &r); if (r != gen_0) return NULL;
738 : }
739 : else
740 : {
741 28 : GEN a = gel(x,1), b = gel(x,2); y = cgetg(3, t_COMPLEX);
742 28 : gel(y,1) = dvmdii(a, c, &r); if (r != gen_0) return NULL;
743 14 : gel(y,2) = dvmdii(b, c, &r); if (r != gen_0) return NULL;
744 : }
745 56 : return y;
746 : }
747 :
748 : static int
749 91 : gauss_cmp(GEN x, GEN y)
750 : {
751 : int v;
752 91 : if (typ(x) != t_COMPLEX)
753 0 : return (typ(y) == t_COMPLEX)? -1: gcmp(x, y);
754 91 : if (typ(y) != t_COMPLEX) return 1;
755 63 : v = cmpii(gel(x,2), gel(y,2));
756 63 : if (v) return v;
757 28 : return gcmp(gel(x,1), gel(y,1));
758 : }
759 :
760 : /* 0 or canonical representative in Z[i]^* / <i> (impose imag(x) >= 0) */
761 : static GEN
762 469 : gauss_normal(GEN x)
763 : {
764 469 : if (typ(x) != t_COMPLEX) return (signe(x) < 0)? absi(x): x;
765 392 : if (signe(gel(x,1)) < 0) x = gneg(x);
766 392 : if (signe(gel(x,2)) < 0) x = mulcxI(x);
767 392 : return x;
768 : }
769 :
770 : static GEN
771 49 : gauss_factor(GEN x)
772 : {
773 49 : pari_sp av = avma;
774 49 : GEN a = real_i(x), b = imag_i(x), d = gen_1, n, y, fa, P, E, P2, E2;
775 49 : long t1 = typ(a);
776 49 : long t2 = typ(b), i, j, l, exp = 0;
777 49 : if (t1 == t_FRAC) d = gel(a,2);
778 49 : if (t2 == t_FRAC) d = lcmii(d, gel(b,2));
779 49 : if (d == gen_1) y = x;
780 : else
781 : {
782 21 : y = gmul(x, d);
783 21 : a = real_i(y); t1 = typ(a);
784 21 : b = imag_i(y); t2 = typ(b);
785 : }
786 49 : if (t1 != t_INT || t2 != t_INT) return NULL;
787 49 : y = gauss_primpart(y, &n);
788 49 : fa = factor(cxnorm(y));
789 49 : P = gel(fa,1);
790 49 : E = gel(fa,2); l = lg(P);
791 49 : P2 = cgetg(l, t_COL);
792 49 : E2 = cgetg(l, t_COL);
793 105 : for (j = 1, i = l-1; i > 0; i--) /* remove largest factors first */
794 : { /* either p = 2 (ramified) or those factors split in Q(i) */
795 56 : GEN p = gel(P,i), w, w2, t, we, pe;
796 56 : long v, e = itos(gel(E,i));
797 56 : int is2 = absequaliu(p, 2);
798 56 : w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
799 56 : w2 = gauss_normal( conj_i(w) );
800 : /* w * w2 * I^3 = p, w2 = conj(w) * I */
801 56 : pe = powiu(p, e);
802 56 : we = gpowgs(w, e);
803 56 : t = gauss_primpart_try( gmul(y, conj_i(we)), pe );
804 56 : if (t) y = t; /* y /= w^e */
805 : else {
806 : /* y /= conj(w)^e, should be y /= w2^e */
807 14 : y = gauss_primpart_try( gmul(y, we), pe );
808 14 : swap(w, w2); exp -= e; /* += 3*e mod 4 */
809 : }
810 56 : gel(P,i) = w;
811 56 : v = Z_pvalrem(n, p, &n);
812 56 : if (v) {
813 7 : exp -= v; /* += 3*v mod 4 */
814 7 : if (is2) v <<= 1; /* 2 = w^2 I^3 */
815 : else {
816 0 : gel(P2,j) = w2;
817 0 : gel(E2,j) = utoipos(v); j++;
818 : }
819 7 : gel(E,i) = stoi(e + v);
820 : }
821 56 : v = Z_pvalrem(d, p, &d);
822 56 : if (v) {
823 7 : exp += v; /* -= 3*v mod 4 */
824 7 : if (is2) v <<= 1; /* 2 is ramified */
825 : else {
826 7 : gel(P2,j) = w2;
827 7 : gel(E2,j) = utoineg(v); j++;
828 : }
829 7 : gel(E,i) = stoi(e - v);
830 : }
831 56 : exp &= 3;
832 : }
833 49 : if (j > 1) {
834 7 : long k = 1;
835 7 : GEN P1 = cgetg(l, t_COL);
836 7 : GEN E1 = cgetg(l, t_COL);
837 : /* remove factors with exponent 0 */
838 14 : for (i = 1; i < l; i++)
839 7 : if (signe(gel(E,i)))
840 : {
841 0 : gel(P1,k) = gel(P,i);
842 0 : gel(E1,k) = gel(E,i);
843 0 : k++;
844 : }
845 7 : setlg(P1, k); setlg(E1, k);
846 7 : setlg(P2, j); setlg(E2, j);
847 7 : fa = famat_mul_shallow(mkmat2(P1,E1), mkmat2(P2,E2));
848 : }
849 49 : if (!equali1(n) || !equali1(d))
850 : {
851 28 : GEN Fa = factor(Qdivii(n, d));
852 28 : P = gel(Fa,1); l = lg(P);
853 28 : E = gel(Fa,2);
854 70 : for (i = 1; i < l; i++)
855 : {
856 42 : GEN w, p = gel(P,i);
857 : long e;
858 : int is2;
859 42 : switch(mod4(p))
860 : {
861 14 : case 3: continue;
862 14 : case 2: is2 = 1; break;
863 14 : default:is2 = 0; break;
864 : }
865 28 : e = itos(gel(E,i));
866 28 : w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
867 28 : gel(P,i) = w;
868 28 : if (is2)
869 14 : gel(E,i) = stoi(2*e);
870 : else
871 : {
872 14 : P = vec_append(P, gauss_normal( conj_i(w) ));
873 14 : E = vec_append(E, gel(E,i));
874 : }
875 28 : exp -= e; /* += 3*e mod 4 */
876 28 : exp &= 3;
877 : }
878 28 : gel(Fa,1) = P;
879 28 : gel(Fa,2) = E;
880 28 : fa = famat_mul_shallow(fa, Fa);
881 : }
882 49 : fa = sort_factor(fa, (void*)&gauss_cmp, &cmp_nodata);
883 :
884 49 : y = gmul(y, powIs(exp));
885 49 : if (!gequal1(y)) {
886 35 : gel(fa,1) = vec_prepend(gel(fa,1), y);
887 35 : gel(fa,2) = vec_prepend(gel(fa,2), gen_1);
888 : }
889 49 : return gerepilecopy(av, fa);
890 : }
891 :
892 : GEN
893 9754 : Q_factor_limit(GEN x, ulong lim)
894 : {
895 9754 : pari_sp av = avma;
896 : GEN a, b;
897 9754 : if (typ(x) == t_INT) return Z_factor_limit(x, lim);
898 5292 : a = Z_factor_limit(gel(x,1), lim);
899 5292 : b = Z_factor_limit(gel(x,2), lim); gel(b,2) = ZC_neg(gel(b,2));
900 5292 : return gerepilecopy(av, merge_factor(a,b,(void*)&cmpii,cmp_nodata));
901 : }
902 : GEN
903 20812 : Q_factor(GEN x)
904 : {
905 20812 : pari_sp av = avma;
906 : GEN a, b;
907 20812 : if (typ(x) == t_INT) return Z_factor(x);
908 33 : a = Z_factor(gel(x,1));
909 35 : b = Z_factor(gel(x,2)); gel(b,2) = ZC_neg(gel(b,2));
910 35 : return gerepilecopy(av, merge_factor(a,b,(void*)&cmpii,cmp_nodata));
911 : }
912 : /* replace t_QUAD/t_COMPLEX coeffs by t_POLMOD in T */
913 : static GEN
914 126 : RgX_fix_quad(GEN x, GEN T)
915 : {
916 126 : long i, l, v = varn(T);
917 126 : GEN y = cgetg_copy(x,&l);
918 630 : for (i = 2; i < l; i++)
919 : {
920 504 : GEN c = gel(x,i);
921 504 : switch(typ(c))
922 : {
923 56 : case t_QUAD: c++;/* fall through */
924 98 : case t_COMPLEX: c = deg1pol_shallow(gel(c,2),gel(c,1),v);
925 : }
926 504 : gel(y,i) = c;
927 : }
928 126 : y[1] = x[1]; return y;
929 : }
930 :
931 : static GEN
932 13629 : RgX_factor(GEN x, GEN dom)
933 : {
934 : pari_sp av;
935 : long pa, v, lx, r1, i;
936 : GEN p, pol, y, p1, p2;
937 13629 : long tx = dom ? RgX_Rg_type(x,dom,&p,&pol,&pa): RgX_type(x,&p,&pol,&pa);
938 13629 : switch(tx)
939 : {
940 7 : case 0: pari_err_IMPL("factor for general polynomials");
941 245 : case t_POL: return RgXY_factor(x, dom);
942 12600 : case t_INT: return ZX_factor(x);
943 7 : case t_FRAC: return QX_factor(x);
944 329 : case t_INTMOD: return factmod(x,p);
945 42 : case t_PADIC: return factorpadic(x,p,pa);
946 98 : case t_FFELT: return FFX_factor(x,pol);
947 :
948 21 : case t_COMPLEX: y = cgetg(3,t_MAT);
949 21 : av = avma; p1 = deg1_from_roots(roots(x,pa), varn(x));
950 21 : gel(y,1) = p1 = gerepileupto(av, p1);
951 21 : gel(y,2) = const_col(lg(p1)-1, gen_1); return y;
952 :
953 28 : case t_REAL: y=cgetg(3,t_MAT); v=varn(x);
954 28 : av=avma; p1=cleanroots(x,pa);
955 28 : lx = lg(p1);
956 70 : for (r1 = 1; r1 < lx; r1++)
957 49 : if (typ(gel(p1,r1)) == t_COMPLEX) break;
958 28 : lx=(r1+lx)>>1; p2=cgetg(lx,t_COL);
959 70 : for (i = 1; i < r1; i++)
960 42 : gel(p2,i) = deg1pol_shallow(gen_1, negr(gel(p1,i)), v);
961 35 : for ( ; i < lx; i++)
962 : {
963 7 : GEN a = gel(p1,2*i-r1);
964 7 : p = cgetg(5, t_POL); gel(p2,i) = p;
965 7 : p[1] = x[1];
966 7 : gel(p,2) = gnorm(a);
967 7 : gel(p,3) = gmul2n(gel(a,1),1); togglesign(gel(p,3));
968 7 : gel(p,4) = gen_1;
969 : }
970 28 : gel(y,1) = gerepileupto(av,p2);
971 28 : gel(y,2) = const_col(lx-1, gen_1); return y;
972 :
973 252 : default:
974 : {
975 252 : GEN w = NULL, T = pol;
976 : long t1, t2;
977 252 : av = avma;
978 252 : RgX_type_decode(tx, &t1, &t2);
979 252 : if (t1 == t_COMPLEX) w = gen_I();
980 203 : else if (t1 == t_QUAD) w = mkquad(pol,gen_0,gen_1);
981 252 : if (w)
982 : { /* substitute I or w by t_POLMOD */
983 126 : T = leafcopy(pol); setvarn(T, fetch_var());
984 126 : x = RgX_fix_quad(x, T);
985 : }
986 252 : switch (t2)
987 : {
988 161 : case t_INT: case t_FRAC: p1 = nffactor(T,x); break;
989 56 : case t_INTMOD:
990 56 : T = RgX_to_FpX(T,p);
991 56 : if (FpX_is_irred(T,p)) { p1 = factmod(x,mkvec2(p,T)); break; }
992 : /*fall through*/
993 : default:
994 56 : if (w) (void)delete_var();
995 56 : pari_err_IMPL("factor for general polynomial");
996 : return NULL; /* LCOV_EXCL_LINE */
997 : }
998 196 : if (t1 == t_POLMOD) return gerepileupto(av, p1);
999 : /* substitute back I or w */
1000 98 : gel(p1,1) = gsubst(liftpol_shallow(gel(p1,1)), varn(T), w);
1001 98 : (void)delete_var(); return gerepilecopy(av, p1);
1002 : }
1003 : }
1004 : }
1005 :
1006 : static GEN
1007 62791 : factor_domain(GEN x, GEN dom)
1008 : {
1009 62791 : long tx = typ(x);
1010 62791 : long tdom = dom ? typ(dom): 0;
1011 : pari_sp av;
1012 :
1013 62791 : if (gequal0(x))
1014 63 : switch(tx)
1015 : {
1016 63 : case t_INT:
1017 : case t_COMPLEX:
1018 : case t_POL:
1019 63 : case t_RFRAC: return prime_fact(x);
1020 0 : default: pari_err_TYPE("factor",x);
1021 : }
1022 62728 : av = avma;
1023 62728 : switch(tx)
1024 : {
1025 2429 : case t_POL: return RgX_factor(x, dom);
1026 35 : case t_RFRAC: {
1027 35 : GEN a = gel(x,1), b = gel(x,2);
1028 35 : GEN y = famat_inv_shallow(RgX_factor(b, dom));
1029 35 : if (typ(a)==t_POL) y = famat_mul_shallow(RgX_factor(a, dom), y);
1030 35 : return gerepilecopy(av, sort_factor_pol(y, cmp_universal));
1031 : }
1032 60194 : case t_INT: if (tdom==0 || tdom==t_INT) return Z_factor(x);
1033 28 : case t_FRAC: if (tdom==0 || tdom==t_INT) return Q_factor(x);
1034 : case t_COMPLEX: /* fall through */
1035 49 : if (tdom==0 || tdom==t_COMPLEX)
1036 : {
1037 49 : GEN y = gauss_factor(x); if (y) return y;
1038 : }
1039 : /* fall through */
1040 : }
1041 0 : pari_err_TYPE("factor",x);
1042 : return NULL; /* LCOV_EXCL_LINE */
1043 : }
1044 :
1045 : GEN
1046 62112 : factor(GEN x) { return factor_domain(x, NULL); }
1047 :
1048 : /*******************************************************************/
1049 : /* */
1050 : /* ROOTS --> MONIC POLYNOMIAL */
1051 : /* */
1052 : /*******************************************************************/
1053 : static GEN
1054 1684321 : normalized_mul(void *E, GEN x, GEN y)
1055 : {
1056 1684321 : long a = gel(x,1)[1], b = gel(y,1)[1];
1057 : (void) E;
1058 1684308 : return mkvec2(mkvecsmall(a + b),
1059 1684321 : RgX_mul_normalized(gel(x,2),a, gel(y,2),b));
1060 : }
1061 : /* L = [Vecsmall([a]), A], with a > 0, A an RgX, deg(A) < a; return X^a + A */
1062 : static GEN
1063 931554 : normalized_to_RgX(GEN L)
1064 : {
1065 931554 : long i, a = gel(L,1)[1];
1066 931554 : GEN A = gel(L,2);
1067 931554 : GEN z = cgetg(a + 3, t_POL);
1068 931554 : z[1] = evalsigne(1) | evalvarn(varn(A));
1069 5439245 : for (i = 2; i < lg(A); i++) gel(z,i) = gcopy(gel(A,i));
1070 934012 : for ( ; i < a+2; i++) gel(z,i) = gen_0;
1071 931555 : gel(z,i) = gen_1; return z;
1072 : }
1073 :
1074 : /* compute prod (x - a[i]) */
1075 : GEN
1076 671978 : roots_to_pol(GEN a, long v)
1077 : {
1078 671978 : pari_sp av = avma;
1079 671978 : long i, k, lx = lg(a);
1080 : GEN L;
1081 671978 : if (lx == 1) return pol_1(v);
1082 671930 : L = cgetg(lx, t_VEC);
1083 1426080 : for (k=1,i=1; i<lx-1; i+=2)
1084 : {
1085 754150 : GEN s = gel(a,i), t = gel(a,i+1);
1086 754150 : GEN x0 = gmul(s,t);
1087 754149 : GEN x1 = gneg(gadd(s,t));
1088 754150 : gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
1089 : }
1090 1321654 : if (i < lx) gel(L,k++) = mkvec2(mkvecsmall(1),
1091 649724 : scalarpol_shallow(gneg(gel(a,i)), v));
1092 671930 : setlg(L, k); L = gen_product(L, NULL, normalized_mul);
1093 671930 : return gerepileupto(av, normalized_to_RgX(L));
1094 : }
1095 :
1096 : /* prod_{i=1..r1} (x - a[i]) prod_{i=1..r2} (x - a[i])(x - conj(a[i]))*/
1097 : GEN
1098 259625 : roots_to_pol_r1(GEN a, long v, long r1)
1099 : {
1100 259625 : pari_sp av = avma;
1101 259625 : long i, k, lx = lg(a);
1102 : GEN L;
1103 259625 : if (lx == 1) return pol_1(v);
1104 259625 : L = cgetg(lx, t_VEC);
1105 775223 : for (k=1,i=1; i<r1; i+=2)
1106 : {
1107 515596 : GEN s = gel(a,i), t = gel(a,i+1);
1108 515596 : GEN x0 = gmul(s,t);
1109 515597 : GEN x1 = gneg(gadd(s,t));
1110 515598 : gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
1111 : }
1112 331499 : if (i < r1+1) gel(L,k++) = mkvec2(mkvecsmall(1),
1113 71872 : scalarpol_shallow(gneg(gel(a,i)), v));
1114 884135 : for (i=r1+1; i<lx; i++)
1115 : {
1116 624507 : GEN s = gel(a,i);
1117 624507 : GEN x0 = gnorm(s);
1118 624496 : GEN x1 = gneg(gtrace(s));
1119 624500 : gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
1120 : }
1121 259628 : setlg(L, k); L = gen_product(L, NULL, normalized_mul);
1122 259625 : return gerepileupto(av, normalized_to_RgX(L));
1123 : }
1124 :
1125 : /*******************************************************************/
1126 : /* */
1127 : /* FACTORBACK */
1128 : /* */
1129 : /*******************************************************************/
1130 : static GEN
1131 53921832 : mul(void *a, GEN x, GEN y) { (void)a; return gmul(x,y);}
1132 : static GEN
1133 78371063 : powi(void *a, GEN x, GEN y) { (void)a; return powgi(x,y);}
1134 : static GEN
1135 25240512 : Fpmul(void *a, GEN x, GEN y) { return Fp_mul(x,y,(GEN)a); }
1136 : static GEN
1137 205261 : Fppow(void *a, GEN x, GEN n) { return Fp_pow(x,n,(GEN)a); }
1138 :
1139 : /* [L,e] = [fa, NULL] or [elts, NULL] or [elts, exponents] */
1140 : GEN
1141 31982723 : gen_factorback(GEN L, GEN e, void *data, GEN (*_mul)(void*,GEN,GEN),
1142 : GEN (*_pow)(void*,GEN,GEN), GEN (*_one)(void*))
1143 : {
1144 31982723 : pari_sp av = avma;
1145 : long k, l, lx;
1146 : GEN p,x;
1147 :
1148 31982723 : if (e) /* supplied vector of exponents */
1149 1309823 : p = L;
1150 : else
1151 : {
1152 30672900 : switch(typ(L)) {
1153 6951031 : case t_VEC:
1154 : case t_COL: /* product of the L[i] */
1155 6951031 : if (lg(L)==1) return _one? _one(data): gen_1;
1156 6877482 : return gerepileupto(av, gen_product(L, data, _mul));
1157 23721881 : case t_MAT: /* genuine factorization */
1158 23721881 : l = lg(L);
1159 23721881 : if (l == 3) break;
1160 : /*fall through*/
1161 : default:
1162 6 : pari_err_TYPE("factorback [not a factorization]", L);
1163 : }
1164 23721874 : p = gel(L,1);
1165 23721874 : e = gel(L,2);
1166 : }
1167 : /* p = elts, e = expo */
1168 25031697 : lx = lg(p);
1169 : /* check whether e is an integral vector of correct length */
1170 25031697 : switch(typ(e))
1171 : {
1172 94584 : case t_VECSMALL:
1173 94584 : if (lx != lg(e))
1174 0 : pari_err_TYPE("factorback [not an exponent vector]", e);
1175 94584 : if (lx == 1) return _one? _one(data): gen_1;
1176 87241 : x = cgetg(lx,t_VEC);
1177 1149868 : for (l=1,k=1; k<lx; k++)
1178 1062627 : if (e[k]) gel(x,l++) = _pow(data, gel(p,k), stoi(e[k]));
1179 87241 : break;
1180 24937111 : case t_VEC: case t_COL:
1181 24937111 : if (lx != lg(e) || !RgV_is_ZV(e))
1182 7 : pari_err_TYPE("factorback [not an exponent vector]", e);
1183 24937104 : if (lx == 1) return _one? _one(data): gen_1;
1184 24836120 : x = cgetg(lx,t_VEC);
1185 104368292 : for (l=1,k=1; k<lx; k++)
1186 79532175 : if (signe(gel(e,k))) gel(x,l++) = _pow(data, gel(p,k), gel(e,k));
1187 24836117 : break;
1188 2 : default:
1189 2 : pari_err_TYPE("factorback [not an exponent vector]", e);
1190 : return NULL;/*LCOV_EXCL_LINE*/
1191 : }
1192 24923358 : if (l==1) return gerepileupto(av, _one? _one(data): gen_1);
1193 24831065 : x[0] = evaltyp(t_VEC) | _evallg(l);
1194 24831065 : return gerepileupto(av, gen_product(x, data, _mul));
1195 : }
1196 :
1197 : GEN
1198 7104765 : FpV_factorback(GEN L, GEN e, GEN p)
1199 7104765 : { return gen_factorback(L, e, (void*)p, &Fpmul, &Fppow, NULL); }
1200 :
1201 : ulong
1202 89887 : Flv_factorback(GEN L, GEN e, ulong p)
1203 : {
1204 89887 : long i, l = lg(e);
1205 89887 : ulong r = 1UL, ri = 1UL;
1206 437549 : for (i = 1; i < l; i++)
1207 : {
1208 347662 : long c = e[i];
1209 347662 : if (!c) continue;
1210 124383 : if (c < 0)
1211 0 : ri = Fl_mul(ri, Fl_powu(L[i],-c,p), p);
1212 : else
1213 124383 : r = Fl_mul(r, Fl_powu(L[i],c,p), p);
1214 : }
1215 89887 : if (ri != 1UL) r = Fl_div(r, ri, p);
1216 89887 : return r;
1217 : }
1218 : GEN
1219 2467 : FlxqV_factorback(GEN L, GEN e, GEN Tp, ulong p)
1220 : {
1221 2467 : pari_sp av = avma;
1222 2467 : GEN Hi = NULL, H = NULL;
1223 2467 : long i, l = lg(L), v = get_Flx_var(Tp);
1224 170338 : for (i = 1; i < l; i++)
1225 : {
1226 167805 : GEN x, ei = gel(e,i);
1227 167805 : long s = signe(ei);
1228 167805 : if (!s) continue;
1229 159572 : x = Flxq_pow(gel(L,i), s > 0? ei: negi(ei), Tp, p);
1230 159564 : if (s > 0)
1231 82484 : H = H? Flxq_mul(H, x, Tp, p): x;
1232 : else
1233 77080 : Hi = Hi? Flxq_mul(Hi, x, Tp, p): x;
1234 : }
1235 2533 : if (!Hi)
1236 : {
1237 0 : if (!H) { set_avma(av); return mkvecsmall2(v,1); }
1238 0 : return gerepileuptoleaf(av, H);
1239 : }
1240 2533 : Hi = Flxq_inv(Hi, Tp, p);
1241 2467 : return gerepileuptoleaf(av, H? Flxq_mul(H,Hi,Tp,p): Hi);
1242 : }
1243 : GEN
1244 14 : FqV_factorback(GEN L, GEN e, GEN Tp, GEN p)
1245 : {
1246 14 : pari_sp av = avma;
1247 14 : GEN Hi = NULL, H = NULL;
1248 14 : long i, l = lg(L), small = typ(e) == t_VECSMALL;
1249 1744 : for (i = 1; i < l; i++)
1250 : {
1251 : GEN x;
1252 : long s;
1253 1730 : if (small)
1254 : {
1255 0 : s = e[i]; if (!s) continue;
1256 0 : x = Fq_powu(gel(L,i), labs(s), Tp, p);
1257 : }
1258 : else
1259 : {
1260 1730 : GEN ei = gel(e,i);
1261 1730 : s = signe(ei); if (!s) continue;
1262 1730 : x = Fq_pow(gel(L,i), s > 0? ei: negi(ei), Tp, p);
1263 : }
1264 1730 : if (s > 0)
1265 846 : H = H? Fq_mul(H, x, Tp, p): x;
1266 : else
1267 884 : Hi = Hi? Fq_mul(Hi, x, Tp, p): x;
1268 : }
1269 14 : if (Hi)
1270 : {
1271 7 : Hi = Fq_inv(Hi, Tp, p);
1272 7 : H = H? Fq_mul(H,Hi,Tp,p): Hi;
1273 : }
1274 7 : else if (!H) return gc_const(av, gen_1);
1275 14 : return gerepileupto(av, H);
1276 : }
1277 :
1278 : GEN
1279 24490718 : factorback2(GEN L, GEN e) { return gen_factorback(L, e, NULL, &mul, &powi, NULL); }
1280 : GEN
1281 1163779 : factorback(GEN fa) { return factorback2(fa, NULL); }
1282 :
1283 : GEN
1284 70 : vecprod(GEN v)
1285 : {
1286 70 : pari_sp av = avma;
1287 70 : if (!is_vec_t(typ(v)))
1288 0 : pari_err_TYPE("vecprod", v);
1289 70 : if (lg(v) == 1) return gen_1;
1290 63 : return gerepilecopy(av, gen_product(v, NULL, mul));
1291 : }
1292 :
1293 : static int
1294 11165 : RgX_is_irred_i(GEN x)
1295 : {
1296 : GEN y, p, pol;
1297 11165 : long l = lg(x), pa;
1298 :
1299 11165 : if (!signe(x) || l <= 3) return 0;
1300 11165 : switch(RgX_type(x,&p,&pol,&pa))
1301 : {
1302 21 : case t_INTMOD: return FpX_is_irred(RgX_to_FpX(x,p), p);
1303 0 : case t_COMPLEX: return l == 4;
1304 0 : case t_REAL:
1305 0 : if (l == 4) return 1;
1306 0 : if (l > 5) return 0;
1307 0 : return gsigne(RgX_disc(x)) > 0;
1308 : }
1309 11144 : y = RgX_factor(x, NULL);
1310 11144 : return (lg(gcoeff(y,1,1))==l);
1311 : }
1312 : static int
1313 11165 : RgX_is_irred(GEN x)
1314 11165 : { pari_sp av = avma; return gc_bool(av, RgX_is_irred_i(x)); }
1315 : long
1316 11165 : polisirreducible(GEN x)
1317 : {
1318 11165 : long tx = typ(x);
1319 11165 : if (tx == t_POL) return RgX_is_irred(x);
1320 0 : if (!is_scalar_t(tx)) pari_err_TYPE("polisirreducible",x);
1321 0 : return 0;
1322 : }
1323 :
1324 : /*******************************************************************/
1325 : /* */
1326 : /* GENERIC GCD */
1327 : /* */
1328 : /*******************************************************************/
1329 : /* x is a COMPLEX or a QUAD */
1330 : static GEN
1331 2555 : triv_cont_gcd(GEN x, GEN y)
1332 : {
1333 2555 : pari_sp av = avma;
1334 : GEN c;
1335 2555 : if (typ(x)==t_COMPLEX)
1336 : {
1337 2233 : GEN a = gel(x,1), b = gel(x,2);
1338 2233 : if (typ(a) == t_REAL || typ(b) == t_REAL) return gen_1;
1339 21 : c = ggcd(a,b);
1340 : }
1341 : else
1342 322 : c = ggcd(gel(x,2),gel(x,3));
1343 343 : return gerepileupto(av, ggcd(c,y));
1344 : }
1345 :
1346 : /* y is a PADIC, x a rational number or an INTMOD */
1347 : static GEN
1348 1407 : padic_gcd(GEN x, GEN y)
1349 : {
1350 1407 : GEN p = gel(y,2);
1351 1407 : long v = gvaluation(x,p), w = valp(y);
1352 1407 : if (w < v) v = w;
1353 1407 : return powis(p, v);
1354 : }
1355 :
1356 : /* x,y in Z[i], at least one of which is t_COMPLEX */
1357 : static GEN
1358 399 : gauss_gcd(GEN x, GEN y)
1359 : {
1360 399 : pari_sp av = avma;
1361 : GEN dx, dy;
1362 399 : dx = denom_i(x); x = gmul(x, dx);
1363 399 : dy = denom_i(y); y = gmul(y, dy);
1364 847 : while (!gequal0(y))
1365 : {
1366 448 : GEN z = gsub(x, gmul(ground(gdiv(x,y)), y));
1367 448 : x = y; y = z;
1368 : }
1369 399 : x = gauss_normal(x);
1370 399 : if (typ(x) == t_COMPLEX)
1371 : {
1372 196 : if (gequal0(gel(x,2))) x = gel(x,1);
1373 196 : else if (gequal0(gel(x,1))) x = gel(x,2);
1374 : }
1375 399 : return gerepileupto(av, gdiv(x, lcmii(dx, dy)));
1376 : }
1377 :
1378 : static int
1379 2590 : c_is_rational(GEN x)
1380 2590 : { return is_rational_t(typ(gel(x,1))) && is_rational_t(typ(gel(x,2))); }
1381 : static GEN
1382 1267 : c_zero_gcd(GEN c)
1383 : {
1384 1267 : GEN x = gel(c,1), y = gel(c,2);
1385 1267 : long tx = typ(x), ty = typ(y);
1386 1267 : if (tx == t_REAL || ty == t_REAL) return gen_1;
1387 42 : if (tx == t_PADIC || tx == t_INTMOD
1388 42 : || ty == t_PADIC || ty == t_INTMOD) return ggcd(x, y);
1389 35 : return gauss_gcd(c, gen_0);
1390 : }
1391 :
1392 : /* gcd(x, 0) */
1393 : static GEN
1394 8244858 : zero_gcd(GEN x)
1395 : {
1396 : pari_sp av;
1397 8244858 : switch(typ(x))
1398 : {
1399 42322 : case t_INT: return absi(x);
1400 13850 : case t_FRAC: return absfrac(x);
1401 1267 : case t_COMPLEX: return c_zero_gcd(x);
1402 616 : case t_REAL: return gen_1;
1403 728 : case t_PADIC: return powis(gel(x,2), valp(x));
1404 245 : case t_SER: return pol_xnall(valp(x), varn(x));
1405 3172 : case t_POLMOD: {
1406 3172 : GEN d = gel(x,2);
1407 3172 : if (typ(d) == t_POL && varn(d) == varn(gel(x,1))) return content(d);
1408 406 : return isinexact(d)? zero_gcd(d): gcopy(d);
1409 : }
1410 7920140 : case t_POL:
1411 7920140 : if (!isinexact(x)) break;
1412 14 : av = avma;
1413 14 : return gerepileupto(av, monomialcopy(content(x), RgX_val(x), varn(x)));
1414 :
1415 217067 : case t_RFRAC:
1416 217067 : if (!isinexact(x)) break;
1417 0 : av = avma;
1418 0 : return gerepileupto(av, gdiv(zero_gcd(gel(x,1)), gel(x,2)));
1419 : }
1420 8182644 : return gcopy(x);
1421 : }
1422 : /* z is an exact zero, t_INT, t_INTMOD or t_FFELT */
1423 : static GEN
1424 9357573 : zero_gcd2(GEN y, GEN z)
1425 : {
1426 : pari_sp av;
1427 9357573 : switch(typ(z))
1428 : {
1429 8229817 : case t_INT: return zero_gcd(y);
1430 1126706 : case t_INTMOD:
1431 1126706 : av = avma;
1432 1126706 : return gerepileupto(av, gmul(y, mkintmod(gen_1,gel(z,1))));
1433 1050 : case t_FFELT:
1434 1050 : av = avma;
1435 1050 : return gerepileupto(av, gmul(y, FF_1(z)));
1436 0 : default:
1437 0 : pari_err_TYPE("zero_gcd", z);
1438 : return NULL;/*LCOV_EXCL_LINE*/
1439 : }
1440 : }
1441 : static GEN
1442 2786273 : cont_gcd_pol_i(GEN x, GEN y) { return scalarpol(ggcd(content(x),y), varn(x));}
1443 : /* tx = t_POL, y considered as constant */
1444 : static GEN
1445 2786266 : cont_gcd_pol(GEN x, GEN y)
1446 2786266 : { pari_sp av = avma; return gerepileupto(av, cont_gcd_pol_i(x,y)); }
1447 : /* tx = t_RFRAC, y considered as constant */
1448 : static GEN
1449 864932 : cont_gcd_rfrac(GEN x, GEN y)
1450 : {
1451 864932 : pari_sp av = avma;
1452 864932 : GEN cx; x = primitive_part(x, &cx);
1453 : /* e.g. Mod(1,2) / (2*y+1) => primitive_part = Mod(1,2)*y^0 */
1454 864932 : if (typ(x) != t_RFRAC) x = cont_gcd_pol_i(x, y);
1455 864925 : else x = gred_rfrac_simple(ggcd(cx? cx: gen_1, y), gel(x,2));
1456 864932 : return gerepileupto(av, x);
1457 : }
1458 : /* !is_const_t(tx), tx != t_POL,t_RFRAC, y considered as constant */
1459 : static GEN
1460 2492 : cont_gcd_gen(GEN x, GEN y)
1461 : {
1462 2492 : pari_sp av = avma;
1463 2492 : return gerepileupto(av, ggcd(content(x),y));
1464 : }
1465 : /* !is_const(tx), y considered as constant */
1466 : static GEN
1467 2799678 : cont_gcd(GEN x, long tx, GEN y)
1468 : {
1469 2799678 : switch(tx)
1470 : {
1471 10934 : case t_RFRAC: return cont_gcd_rfrac(x,y);
1472 2786252 : case t_POL: return cont_gcd_pol(x,y);
1473 2492 : default: return cont_gcd_gen(x,y);
1474 : }
1475 : }
1476 : static GEN
1477 11202039 : gcdiq(GEN x, GEN y)
1478 : {
1479 : GEN z;
1480 11202039 : if (!signe(x)) return Q_abs(y);
1481 4113473 : z = cgetg(3,t_FRAC);
1482 4113495 : gel(z,1) = gcdii(x,gel(y,1));
1483 4113453 : gel(z,2) = icopy(gel(y,2));
1484 4113478 : return z;
1485 : }
1486 : static GEN
1487 25534678 : gcdqq(GEN x, GEN y)
1488 : {
1489 25534678 : GEN z = cgetg(3,t_FRAC);
1490 25534679 : gel(z,1) = gcdii(gel(x,1), gel(y,1));
1491 25534264 : gel(z,2) = lcmii(gel(x,2), gel(y,2));
1492 25534604 : return z;
1493 : }
1494 : /* assume x,y t_INT or t_FRAC */
1495 : GEN
1496 239025833 : Q_gcd(GEN x, GEN y)
1497 : {
1498 239025833 : long tx = typ(x), ty = typ(y);
1499 239025833 : if (tx == t_INT)
1500 204890426 : { return (ty == t_INT)? gcdii(x,y): gcdiq(x,y); }
1501 : else
1502 34135407 : { return (ty == t_INT)? gcdiq(y,x): gcdqq(x,y); }
1503 : }
1504 :
1505 : GEN
1506 32492512 : ggcd(GEN x, GEN y)
1507 : {
1508 32492512 : long vx, vy, tx = typ(x), ty = typ(y);
1509 : pari_sp av, tetpil;
1510 : GEN p1,z;
1511 :
1512 64985024 : if (is_noncalc_t(tx) || is_matvec_t(tx) ||
1513 64985024 : is_noncalc_t(ty) || is_matvec_t(ty)) pari_err_TYPE2("gcd",x,y);
1514 32492513 : if (tx>ty) { swap(x,y); lswap(tx,ty); }
1515 : /* tx <= ty */
1516 32492513 : z = gisexactzero(x); if (z) return zero_gcd2(y,z);
1517 29413611 : z = gisexactzero(y); if (z) return zero_gcd2(x,z);
1518 23134940 : if (is_const_t(tx))
1519 : {
1520 15171318 : if (ty == tx) switch(tx)
1521 : {
1522 8132995 : case t_INT:
1523 8132995 : return gcdii(x,y);
1524 :
1525 4084710 : case t_INTMOD: z=cgetg(3,t_INTMOD);
1526 4084710 : if (equalii(gel(x,1),gel(y,1)))
1527 4084703 : gel(z,1) = icopy(gel(x,1));
1528 : else
1529 7 : gel(z,1) = gcdii(gel(x,1),gel(y,1));
1530 4084710 : if (gequal1(gel(z,1))) gel(z,2) = gen_0;
1531 : else
1532 : {
1533 4084710 : av = avma; p1 = gcdii(gel(z,1),gel(x,2));
1534 4084710 : if (!equali1(p1))
1535 : {
1536 7 : p1 = gcdii(p1,gel(y,2));
1537 7 : if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
1538 7 : else p1 = gerepileuptoint(av, p1);
1539 : }
1540 4084710 : gel(z,2) = p1;
1541 : }
1542 4084710 : return z;
1543 :
1544 158716 : case t_FRAC:
1545 158716 : return gcdqq(x,y);
1546 :
1547 5397 : case t_FFELT:
1548 5397 : if (!FF_samefield(x,y)) pari_err_OP("gcd",x,y);
1549 5397 : return FF_equal0(x) && FF_equal0(y)? FF_zero(y): FF_1(y);
1550 :
1551 21 : case t_COMPLEX:
1552 21 : if (c_is_rational(x) && c_is_rational(y)) return gauss_gcd(x,y);
1553 7 : return triv_cont_gcd(y,x);
1554 :
1555 14 : case t_PADIC:
1556 14 : if (!equalii(gel(x,2),gel(y,2))) return gen_1;
1557 7 : return powis(gel(y,2), minss(valp(x), valp(y)));
1558 :
1559 133 : case t_QUAD:
1560 133 : av=avma; p1=gdiv(x,y);
1561 133 : if (gequal0(gel(p1,3)))
1562 : {
1563 14 : p1=gel(p1,2);
1564 14 : if (typ(p1)==t_INT) { set_avma(av); return gcopy(y); }
1565 7 : tetpil=avma; return gerepile(av,tetpil, gdiv(y,gel(p1,2)));
1566 : }
1567 119 : if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) {set_avma(av); return gcopy(y);}
1568 112 : p1 = ginv(p1); set_avma(av);
1569 112 : if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) return gcopy(x);
1570 105 : return triv_cont_gcd(y,x);
1571 :
1572 0 : default: return gen_1; /* t_REAL */
1573 : }
1574 2789332 : if (is_const_t(ty)) switch(tx)
1575 : {
1576 16786 : case t_INT:
1577 16786 : switch(ty)
1578 : {
1579 70 : case t_INTMOD: z = cgetg(3,t_INTMOD);
1580 70 : gel(z,1) = icopy(gel(y,1)); av = avma;
1581 70 : p1 = gcdii(gel(y,1),gel(y,2));
1582 70 : if (!equali1(p1)) {
1583 14 : p1 = gcdii(x,p1);
1584 14 : if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
1585 : else
1586 14 : p1 = gerepileuptoint(av, p1);
1587 : }
1588 70 : gel(z,2) = p1; return z;
1589 :
1590 8162 : case t_REAL: return gen_1;
1591 :
1592 4452 : case t_FRAC:
1593 4452 : return gcdiq(x,y);
1594 :
1595 2471 : case t_COMPLEX:
1596 2471 : if (c_is_rational(y)) return gauss_gcd(x,y);
1597 2128 : return triv_cont_gcd(y,x);
1598 :
1599 112 : case t_FFELT:
1600 112 : if (!FF_equal0(y)) return FF_1(y);
1601 0 : return dvdii(x, gel(y,4))? FF_zero(y): FF_1(y);
1602 :
1603 1393 : case t_PADIC:
1604 1393 : return padic_gcd(x,y);
1605 :
1606 126 : case t_QUAD:
1607 126 : return triv_cont_gcd(y,x);
1608 0 : default:
1609 0 : pari_err_TYPE2("gcd",x,y);
1610 : }
1611 :
1612 14 : case t_REAL:
1613 14 : switch(ty)
1614 : {
1615 14 : case t_INTMOD:
1616 : case t_FFELT:
1617 14 : case t_PADIC: pari_err_TYPE2("gcd",x,y);
1618 0 : default: return gen_1;
1619 : }
1620 :
1621 49 : case t_INTMOD:
1622 49 : switch(ty)
1623 : {
1624 14 : case t_FRAC:
1625 14 : av = avma; p1=gcdii(gel(x,1),gel(y,2)); set_avma(av);
1626 14 : if (!equali1(p1)) pari_err_OP("gcd",x,y);
1627 7 : return ggcd(gel(y,1), x);
1628 :
1629 14 : case t_FFELT:
1630 : {
1631 14 : GEN p = gel(y,4);
1632 14 : if (!dvdii(gel(x,1), p)) pari_err_OP("gcd",x,y);
1633 7 : if (!FF_equal0(y)) return FF_1(y);
1634 0 : return dvdii(gel(x,2),p)? FF_zero(y): FF_1(y);
1635 : }
1636 :
1637 14 : case t_COMPLEX: case t_QUAD:
1638 14 : return triv_cont_gcd(y,x);
1639 :
1640 7 : case t_PADIC:
1641 7 : return padic_gcd(x,y);
1642 :
1643 0 : default: pari_err_TYPE2("gcd",x,y);
1644 : }
1645 :
1646 210 : case t_FRAC:
1647 210 : switch(ty)
1648 : {
1649 84 : case t_COMPLEX:
1650 84 : if (c_is_rational(y)) return gauss_gcd(x,y);
1651 : case t_QUAD:
1652 154 : return triv_cont_gcd(y,x);
1653 42 : case t_FFELT:
1654 : {
1655 42 : GEN p = gel(y,4);
1656 42 : if (dvdii(gel(x,2), p)) pari_err_OP("gcd",x,y);
1657 21 : if (!FF_equal0(y)) return FF_1(y);
1658 0 : return dvdii(gel(x,1),p)? FF_zero(y): FF_1(y);
1659 : }
1660 :
1661 7 : case t_PADIC:
1662 7 : return padic_gcd(x,y);
1663 :
1664 0 : default: pari_err_TYPE2("gcd",x,y);
1665 : }
1666 70 : case t_FFELT:
1667 70 : switch(ty)
1668 : {
1669 42 : case t_PADIC:
1670 : {
1671 42 : GEN p = gel(y,2);
1672 42 : long v = valp(y);
1673 42 : if (!equalii(p, gel(x,4)) || v < 0) pari_err_OP("gcd",x,y);
1674 14 : return (v && FF_equal0(x))? FF_zero(x): FF_1(x);
1675 : }
1676 28 : default: pari_err_TYPE2("gcd",x,y);
1677 : }
1678 :
1679 14 : case t_COMPLEX:
1680 14 : switch(ty)
1681 : {
1682 14 : case t_PADIC:
1683 14 : case t_QUAD: return triv_cont_gcd(x,y);
1684 0 : default: pari_err_TYPE2("gcd",x,y);
1685 : }
1686 :
1687 7 : case t_PADIC:
1688 7 : switch(ty)
1689 : {
1690 7 : case t_QUAD: return triv_cont_gcd(y,x);
1691 0 : default: pari_err_TYPE2("gcd",x,y);
1692 : }
1693 :
1694 0 : default: return gen_1; /* tx = t_REAL */
1695 : }
1696 2772182 : return cont_gcd(y,ty, x);
1697 : }
1698 :
1699 7963622 : if (tx == t_POLMOD)
1700 : {
1701 6119 : if (ty == t_POLMOD)
1702 : {
1703 6063 : GEN T = gel(x,1);
1704 6063 : z = cgetg(3,t_POLMOD);
1705 6063 : T = RgX_equal_var(T,gel(y,1))? RgX_copy(T): RgX_gcd(T, gel(y,1));
1706 6063 : gel(z,1) = T;
1707 6063 : if (degpol(T) <= 0) gel(z,2) = gen_0;
1708 : else
1709 : {
1710 : GEN X, Y, d;
1711 6063 : av = avma; X = gel(x,2); Y = gel(y,2);
1712 6063 : d = ggcd(content(X), content(Y));
1713 6063 : if (!gequal1(d)) { X = gdiv(X,d); Y = gdiv(Y,d); }
1714 6063 : p1 = ggcd(T, X);
1715 6063 : gel(z,2) = gerepileupto(av, gmul(d, ggcd(p1, Y)));
1716 : }
1717 6063 : return z;
1718 : }
1719 56 : vx = varn(gel(x,1));
1720 56 : switch(ty)
1721 : {
1722 28 : case t_POL:
1723 28 : vy = varn(y);
1724 28 : if (varncmp(vy,vx) < 0) return cont_gcd_pol(y, x);
1725 14 : z = cgetg(3,t_POLMOD);
1726 14 : gel(z,1) = RgX_copy(gel(x,1));
1727 14 : av = avma; p1 = ggcd(gel(x,1),gel(x,2));
1728 14 : gel(z,2) = gerepileupto(av, ggcd(p1,y));
1729 14 : return z;
1730 :
1731 28 : case t_RFRAC:
1732 28 : vy = varn(gel(y,2));
1733 28 : if (varncmp(vy,vx) < 0) return cont_gcd_rfrac(y, x);
1734 28 : av = avma;
1735 28 : p1 = ggcd(gel(x,1),gel(y,2));
1736 28 : if (degpol(p1)) pari_err_OP("gcd",x,y);
1737 21 : set_avma(av); return gdiv(ggcd(gel(y,1),x), content(gel(y,2)));
1738 : }
1739 : }
1740 :
1741 7957503 : vx = gvar(x);
1742 7957503 : vy = gvar(y);
1743 7957503 : if (varncmp(vy, vx) < 0) return cont_gcd(y,ty, x);
1744 7945484 : if (varncmp(vy, vx) > 0) return cont_gcd(x,tx, y);
1745 :
1746 : /* vx = vy: same main variable */
1747 7930007 : switch(tx)
1748 : {
1749 7769073 : case t_POL:
1750 : switch(ty)
1751 : {
1752 6915054 : case t_POL: return RgX_gcd(x,y);
1753 21 : case t_SER:
1754 21 : z = ggcd(content(x), content(y));
1755 21 : return monomialcopy(z, minss(valp(y),gval(x,vx)), vx);
1756 853998 : case t_RFRAC: return cont_gcd_rfrac(y, x);
1757 : }
1758 0 : break;
1759 :
1760 14 : case t_SER:
1761 14 : z = ggcd(content(x), content(y));
1762 : switch(ty)
1763 : {
1764 7 : case t_SER: return monomialcopy(z, minss(valp(x),valp(y)), vx);
1765 7 : case t_RFRAC: return monomialcopy(z, minss(valp(x),gval(y,vx)), vx);
1766 : }
1767 0 : break;
1768 :
1769 160920 : case t_RFRAC:
1770 : {
1771 160920 : GEN xd = gel(x,2), yd = gel(y,2);
1772 160920 : if (ty != t_RFRAC) pari_err_TYPE2("gcd",x,y);
1773 160920 : z = cgetg(3,t_RFRAC); av = avma;
1774 160920 : gel(z,2) = gerepileupto(av, RgX_mul(xd, RgX_div(yd, RgX_gcd(xd, yd))));
1775 160920 : gel(z,1) = ggcd(gel(x,1), gel(y,1)); return z;
1776 : }
1777 : }
1778 0 : pari_err_TYPE2("gcd",x,y);
1779 : return NULL; /* LCOV_EXCL_LINE */
1780 : }
1781 : GEN
1782 5203504 : ggcd0(GEN x, GEN y) { return y? ggcd(x,y): content(x); }
1783 :
1784 : static GEN
1785 3568 : fix_lcm(GEN x)
1786 : {
1787 : GEN t;
1788 3568 : switch(typ(x))
1789 : {
1790 3470 : case t_INT: if (signe(x)<0) x = negi(x);
1791 3470 : break;
1792 98 : case t_POL:
1793 98 : if (lg(x) <= 2) break;
1794 98 : t = leading_coeff(x);
1795 98 : if (typ(t) == t_INT && signe(t) < 0) x = gneg(x);
1796 : }
1797 3568 : return x;
1798 : }
1799 : GEN
1800 2877 : glcm0(GEN x, GEN y)
1801 : {
1802 2877 : if (!y) return fix_lcm(gassoc_proto(glcm,x,y));
1803 2828 : return glcm(x,y);
1804 : }
1805 : GEN
1806 3470 : ZV_lcm(GEN x) { return fix_lcm(gassoc_proto(lcmii,x,NULL)); }
1807 : GEN
1808 3171 : glcm(GEN x, GEN y)
1809 : {
1810 : pari_sp av;
1811 : GEN z;
1812 3171 : if (typ(x)==t_INT && typ(y)==t_INT) return lcmii(x,y);
1813 63 : av = avma; z = ggcd(x,y);
1814 63 : if (!gequal1(z))
1815 : {
1816 63 : if (gequal0(z)) { set_avma(av); return gmul(x,y); }
1817 49 : y = gdiv(y,z);
1818 : }
1819 49 : return gerepileupto(av, fix_lcm(gmul(x,y)));
1820 : }
1821 :
1822 : /* x + r ~ x ? Assume x,r are t_POL, deg(r) <= deg(x) */
1823 : static int
1824 0 : pol_approx0(GEN r, GEN x, int exact)
1825 : {
1826 : long i, l;
1827 0 : if (exact) return !signe(r);
1828 0 : l = minss(lg(x), lg(r));
1829 0 : for (i = 2; i < l; i++)
1830 0 : if (!cx_approx0(gel(r,i), gel(x,i))) return 0;
1831 0 : return 1;
1832 : }
1833 :
1834 : GEN
1835 0 : RgX_gcd_simple(GEN x, GEN y)
1836 : {
1837 0 : pari_sp av1, av = avma;
1838 0 : GEN r, yorig = y;
1839 0 : int exact = !(isinexactreal(x) || isinexactreal(y));
1840 :
1841 : for(;;)
1842 : {
1843 0 : av1 = avma; r = RgX_rem(x,y);
1844 0 : if (pol_approx0(r, x, exact))
1845 : {
1846 0 : set_avma(av1);
1847 0 : if (y == yorig) return RgX_copy(y);
1848 0 : y = normalizepol_approx(y, lg(y));
1849 0 : if (lg(y) == 3) { set_avma(av); return pol_1(varn(x)); }
1850 0 : return gerepileupto(av,y);
1851 : }
1852 0 : x = y; y = r;
1853 0 : if (gc_needed(av,1)) {
1854 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd_simple");
1855 0 : gerepileall(av,2, &x,&y);
1856 : }
1857 : }
1858 : }
1859 : GEN
1860 0 : RgX_extgcd_simple(GEN a, GEN b, GEN *pu, GEN *pv)
1861 : {
1862 0 : pari_sp av = avma;
1863 : GEN q, r, d, d1, u, v, v1;
1864 0 : int exact = !(isinexactreal(a) || isinexactreal(b));
1865 :
1866 0 : d = a; d1 = b; v = gen_0; v1 = gen_1;
1867 : for(;;)
1868 : {
1869 0 : if (pol_approx0(d1, a, exact)) break;
1870 0 : q = poldivrem(d,d1, &r);
1871 0 : v = gsub(v, gmul(q,v1));
1872 0 : u=v; v=v1; v1=u;
1873 0 : u=r; d=d1; d1=u;
1874 : }
1875 0 : u = gsub(d, gmul(b,v));
1876 0 : u = RgX_div(u,a);
1877 :
1878 0 : gerepileall(av, 3, &u,&v,&d);
1879 0 : *pu = u;
1880 0 : *pv = v; return d;
1881 : }
1882 :
1883 : GEN
1884 42 : ghalfgcd(GEN x, GEN y)
1885 : {
1886 42 : long tx = typ(x), ty = typ(y);
1887 42 : if (tx==t_INT && ty==t_INT) return halfgcdii(x, y);
1888 35 : if (tx==t_POL && ty==t_POL && varn(x)==varn(y)) return RgX_halfgcd(x, y);
1889 0 : pari_err_OP("halfgcd", x, y);
1890 : return NULL; /* LCOV_EXCL_LINE */
1891 : }
1892 :
1893 : /*******************************************************************/
1894 : /* */
1895 : /* CONTENT / PRIMITIVE PART */
1896 : /* */
1897 : /*******************************************************************/
1898 :
1899 : GEN
1900 74074595 : content(GEN x)
1901 : {
1902 74074595 : long lx, i, t, tx = typ(x);
1903 74074595 : pari_sp av = avma;
1904 : GEN c;
1905 :
1906 74074595 : if (is_scalar_t(tx)) return zero_gcd(x);
1907 74062355 : switch(tx)
1908 : {
1909 864967 : case t_RFRAC:
1910 : {
1911 864967 : GEN n = gel(x,1), d = gel(x,2);
1912 : /* -- varncmp(vn, vd) < 0 can't happen
1913 : * -- if n is POLMOD, its main variable (in the sense of gvar2)
1914 : * has lower priority than denominator */
1915 864967 : if (typ(n) == t_POLMOD || varncmp(gvar(n), varn(d)) > 0)
1916 824622 : n = isinexact(n)? zero_gcd(n): gcopy(n);
1917 : else
1918 40345 : n = content(n);
1919 864967 : return gerepileupto(av, gdiv(n, content(d)));
1920 : }
1921 :
1922 2568326 : case t_VEC: case t_COL:
1923 2568326 : lx = lg(x); if (lx==1) return gen_0;
1924 2568319 : break;
1925 :
1926 1652 : case t_MAT:
1927 : {
1928 : long hx, j;
1929 1652 : lx = lg(x);
1930 1652 : if (lx == 1) return gen_0;
1931 1645 : hx = lgcols(x);
1932 1645 : if (hx == 1) return gen_0;
1933 1638 : if (lx == 2) { x = gel(x,1); lx = lg(x); break; }
1934 1638 : if (hx == 2) { x = row_i(x, 1, 1, lx-1); break; }
1935 1638 : c = content(gel(x,1));
1936 3276 : for (j=2; j<lx; j++)
1937 4914 : for (i=1; i<hx; i++) c = ggcd(c,gcoeff(x,i,j));
1938 1638 : if (typ(c) == t_INTMOD || isinexact(c)) return gc_const(av, gen_1);
1939 1638 : return gerepileupto(av,c);
1940 : }
1941 :
1942 70627207 : case t_POL: case t_SER:
1943 70627207 : lx = lg(x); if (lx == 2) return gen_0;
1944 70604023 : break;
1945 21 : case t_VECSMALL: return utoi(zv_content(x));
1946 182 : case t_QFB:
1947 182 : lx = 4; break;
1948 :
1949 0 : default: pari_err_TYPE("content",x);
1950 : return NULL; /* LCOV_EXCL_LINE */
1951 : }
1952 215793778 : for (i=lontyp[tx]; i<lx; i++)
1953 155812306 : if (typ(gel(x,i)) != t_INT) break;
1954 73172525 : lx--; c = gel(x,lx);
1955 73172525 : t = typ(c); if (is_matvec_t(t)) c = content(c);
1956 73172528 : if (i > lx)
1957 : { /* integer coeffs */
1958 62875072 : while (lx-- > lontyp[tx])
1959 : {
1960 61381162 : c = gcdii(c, gel(x,lx));
1961 61381129 : if (equali1(c)) return gc_const(av, gen_1);
1962 : }
1963 : }
1964 : else
1965 : {
1966 13191053 : if (isinexact(c)) c = zero_gcd(c);
1967 36639912 : while (lx-- > lontyp[tx])
1968 : {
1969 23448859 : GEN d = gel(x,lx);
1970 23448859 : t = typ(d); if (is_matvec_t(t)) d = content(d);
1971 23448859 : c = ggcd(c, d);
1972 : }
1973 13191053 : if (isinexact(c)) return gc_const(av, gen_1);
1974 : }
1975 14684963 : switch(typ(c))
1976 : {
1977 1498825 : case t_INT:
1978 1498825 : if (signe(c) < 0) c = negi(c);
1979 1498825 : break;
1980 0 : case t_VEC: case t_COL: case t_MAT:
1981 0 : pari_err_TYPE("content",x);
1982 : }
1983 :
1984 14684977 : return av==avma? gcopy(c): gerepileupto(av,c);
1985 : }
1986 :
1987 : GEN
1988 1700922 : primitive_part(GEN x, GEN *ptc)
1989 : {
1990 1700922 : pari_sp av = avma;
1991 1700922 : GEN c = content(x);
1992 1700908 : if (gequal1(c)) { set_avma(av); c = NULL; }
1993 99807 : else if (!gequal0(c)) x = gdiv(x,c);
1994 1700913 : if (ptc) *ptc = c;
1995 1700913 : return x;
1996 : }
1997 : GEN
1998 2233 : primpart(GEN x) { return primitive_part(x, NULL); }
1999 :
2000 : static GEN
2001 71608248 : Q_content_v(GEN x, long imin, long l)
2002 : {
2003 71608248 : pari_sp av = avma;
2004 71608248 : long i = l-1;
2005 71608248 : GEN d = Q_content_safe(gel(x,i));
2006 71606833 : if (!d) return NULL;
2007 310577513 : for (i--; i >= imin; i--)
2008 : {
2009 239001496 : GEN c = Q_content_safe(gel(x,i));
2010 239012440 : if (!c) return NULL;
2011 239012426 : d = Q_gcd(d, c);
2012 238972050 : if (gc_needed(av,1)) d = gerepileupto(av, d);
2013 : }
2014 71576017 : return gerepileupto(av, d);
2015 : }
2016 : /* As content(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
2017 : * of Q(x2,...,xn)[x1] */
2018 : GEN
2019 359467160 : Q_content_safe(GEN x)
2020 : {
2021 : long l;
2022 359467160 : switch(typ(x))
2023 : {
2024 253822740 : case t_INT: return absi(x);
2025 33686428 : case t_FRAC: return absfrac(x);
2026 35839147 : case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
2027 35839147 : l = lg(x); return l==1? gen_1: Q_content_v(x, 1, l);
2028 36100585 : case t_POL:
2029 36100585 : l = lg(x); return l==2? gen_0: Q_content_v(x, 2, l);
2030 32683 : case t_POLMOD: return Q_content_safe(gel(x,2));
2031 70 : case t_RFRAC:
2032 : {
2033 : GEN a, b;
2034 70 : a = Q_content_safe(gel(x,1)); if (!a) return NULL;
2035 63 : b = Q_content_safe(gel(x,2)); if (!b) return NULL;
2036 63 : return gdiv(a, b);
2037 : }
2038 : }
2039 366 : return NULL;
2040 : }
2041 : GEN
2042 308318 : Q_content(GEN x)
2043 : {
2044 308318 : GEN c = Q_content_safe(x);
2045 308319 : if (!c) pari_err_TYPE("Q_content",x);
2046 308319 : return c;
2047 : }
2048 :
2049 : GEN
2050 13146 : ZX_content(GEN x)
2051 : {
2052 13146 : long i, l = lg(x);
2053 : GEN d;
2054 : pari_sp av;
2055 :
2056 13146 : if (l == 2) return gen_0;
2057 13146 : d = gel(x,2);
2058 13146 : if (l == 3) return absi(d);
2059 9170 : av = avma;
2060 18977 : for (i=3; !is_pm1(d) && i<l; i++) d = gcdii(d, gel(x,i));
2061 9170 : if (signe(d) < 0) d = negi(d);
2062 9170 : return gerepileuptoint(av, d);
2063 : }
2064 :
2065 : static GEN
2066 43169 : Z_content_v(GEN x, long i, long l)
2067 : {
2068 43169 : pari_sp av = avma;
2069 43169 : GEN d = Z_content(gel(x,i));
2070 43169 : if (!d) return NULL;
2071 273098 : for (i++; i<l; i++)
2072 : {
2073 235060 : GEN c = Z_content(gel(x,i));
2074 235060 : if (!c) return NULL;
2075 234878 : d = gcdii(d, c); if (equali1(d)) return NULL;
2076 234752 : if ((i & 255) == 0) d = gerepileuptoint(av, d);
2077 : }
2078 38038 : return gerepileuptoint(av, d);
2079 : }
2080 : /* return NULL for 1 */
2081 : GEN
2082 280049 : Z_content(GEN x)
2083 : {
2084 : long l;
2085 280049 : switch(typ(x))
2086 : {
2087 217875 : case t_INT:
2088 217875 : if (is_pm1(x)) return NULL;
2089 216692 : return absi(x);
2090 9002 : case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
2091 9002 : l = lg(x); return l==1? NULL: Z_content_v(x, 1, l);
2092 53172 : case t_POL:
2093 53172 : l = lg(x); return l==2? gen_0: Z_content_v(x, 2, l);
2094 0 : case t_POLMOD: return Z_content(gel(x,2));
2095 : }
2096 0 : pari_err_TYPE("Z_content", x);
2097 : return NULL; /* LCOV_EXCL_LINE */
2098 : }
2099 :
2100 : static GEN
2101 31737440 : Q_denom_v(GEN x, long i, long l)
2102 : {
2103 31737440 : pari_sp av = avma;
2104 31737440 : GEN d = Q_denom_safe(gel(x,i));
2105 31737283 : if (!d) return NULL;
2106 123542005 : for (i++; i<l; i++)
2107 : {
2108 91804885 : GEN D = Q_denom_safe(gel(x,i));
2109 91804651 : if (!D) return NULL;
2110 91804651 : if (D != gen_1) d = lcmii(d, D);
2111 91804563 : if ((i & 255) == 0) d = gerepileuptoint(av, d);
2112 : }
2113 31737120 : return gerepileuptoint(av, d);
2114 : }
2115 : /* NOT MEMORY CLEAN (because of t_FRAC).
2116 : * As denom(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
2117 : * of Q(x2,...,xn)[x1] */
2118 : GEN
2119 161606191 : Q_denom_safe(GEN x)
2120 : {
2121 : long l;
2122 161606191 : switch(typ(x))
2123 : {
2124 109233985 : case t_INT: return gen_1;
2125 28 : case t_PADIC: l = valp(x); return l < 0? powiu(gel(x,2), -l): gen_1;
2126 20447020 : case t_FRAC: return gel(x,2);
2127 504 : case t_QUAD: return Q_denom_v(x, 2, 4);
2128 21426506 : case t_COMPLEX: case t_VEC: case t_COL: case t_MAT:
2129 21426506 : l = lg(x); return l==1? gen_1: Q_denom_v(x, 1, l);
2130 10401308 : case t_POL: case t_SER:
2131 10401308 : l = lg(x); return l==2? gen_1: Q_denom_v(x, 2, l);
2132 93282 : case t_POLMOD: return Q_denom(gel(x,2));
2133 8134 : case t_RFRAC:
2134 : {
2135 : GEN a, b;
2136 8134 : a = Q_content(gel(x,1)); if (!a) return NULL;
2137 8134 : b = Q_content(gel(x,2)); if (!b) return NULL;
2138 8134 : return Q_denom(gdiv(a, b));
2139 : }
2140 : }
2141 66 : return NULL;
2142 : }
2143 : GEN
2144 2160891 : Q_denom(GEN x)
2145 : {
2146 2160891 : GEN d = Q_denom_safe(x);
2147 2160882 : if (!d) pari_err_TYPE("Q_denom",x);
2148 2160882 : return d;
2149 : }
2150 :
2151 : GEN
2152 35906154 : Q_remove_denom(GEN x, GEN *ptd)
2153 : {
2154 35906154 : GEN d = Q_denom_safe(x);
2155 35905642 : if (d) { if (d == gen_1) d = NULL; else x = Q_muli_to_int(x,d); }
2156 35905592 : if (ptd) *ptd = d;
2157 35905592 : return x;
2158 : }
2159 :
2160 : /* return y = x * d, assuming x rational, and d,y integral */
2161 : GEN
2162 106140852 : Q_muli_to_int(GEN x, GEN d)
2163 : {
2164 : long i, l;
2165 : GEN y, xn, xd;
2166 : pari_sp av;
2167 :
2168 106140852 : if (typ(d) != t_INT) pari_err_TYPE("Q_muli_to_int",d);
2169 106145830 : switch (typ(x))
2170 : {
2171 34458058 : case t_INT:
2172 34458058 : return mulii(x,d);
2173 :
2174 50555480 : case t_FRAC:
2175 50555480 : xn = gel(x,1);
2176 50555480 : xd = gel(x,2); av = avma;
2177 50555480 : y = mulii(xn, diviiexact(d, xd));
2178 50546708 : return gerepileuptoint(av, y);
2179 7 : case t_COMPLEX:
2180 7 : y = cgetg(3,t_COMPLEX);
2181 7 : gel(y,1) = Q_muli_to_int(gel(x,1),d);
2182 7 : gel(y,2) = Q_muli_to_int(gel(x,2),d);
2183 7 : return y;
2184 14 : case t_PADIC:
2185 14 : y = gcopy(x); if (!isint1(d)) setvalp(y, 0);
2186 14 : return y;
2187 175 : case t_QUAD:
2188 175 : y = cgetg(4,t_QUAD);
2189 175 : gel(y,1) = ZX_copy(gel(x,1));
2190 175 : gel(y,2) = Q_muli_to_int(gel(x,2),d);
2191 175 : gel(y,3) = Q_muli_to_int(gel(x,3),d); return y;
2192 :
2193 12222528 : case t_VEC: case t_COL: case t_MAT:
2194 12222528 : y = cgetg_copy(x, &l);
2195 70559820 : for (i=1; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
2196 12221941 : return y;
2197 :
2198 8859315 : case t_POL: case t_SER:
2199 8859315 : y = cgetg_copy(x, &l); y[1] = x[1];
2200 41950708 : for (i=2; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
2201 8858117 : return y;
2202 :
2203 50232 : case t_POLMOD:
2204 50232 : retmkpolmod(Q_muli_to_int(gel(x,2), d), RgX_copy(gel(x,1)));
2205 21 : case t_RFRAC:
2206 21 : return gmul(x, d);
2207 : }
2208 0 : pari_err_TYPE("Q_muli_to_int",x);
2209 : return NULL; /* LCOV_EXCL_LINE */
2210 : }
2211 :
2212 : static void
2213 5866092 : rescale_init(GEN c, int *exact, long *emin, GEN *D)
2214 : {
2215 : long e, i;
2216 5866092 : switch(typ(c))
2217 : {
2218 4214257 : case t_REAL:
2219 4214257 : *exact = 0;
2220 4214257 : if (!signe(c)) return;
2221 4187362 : e = expo(c) + 1 - bit_prec(c);
2222 4215344 : for (i = lg(c)-1; i > 2; i--, e += BITS_IN_LONG)
2223 3976083 : if (c[i]) break;
2224 4187363 : e += vals(c[i]); break; /* e[2] != 0 */
2225 1647315 : case t_INT:
2226 1647315 : if (!signe(c)) return;
2227 726583 : e = expi(c);
2228 726583 : break;
2229 4494 : case t_FRAC:
2230 4494 : e = expi(gel(c,1)) - expi(gel(c,2));
2231 4494 : if (*exact) *D = lcmii(*D, gel(c,2));
2232 4494 : break;
2233 26 : default:
2234 26 : pari_err_TYPE("rescale_to_int",c);
2235 : return; /* LCOV_EXCL_LINE */
2236 : }
2237 4918438 : if (e < *emin) *emin = e;
2238 : }
2239 : GEN
2240 561797 : RgM_rescale_to_int(GEN x)
2241 : {
2242 561797 : long lx = lg(x), i,j, hx, emin;
2243 : GEN D;
2244 : int exact;
2245 :
2246 561797 : if (lx == 1) return cgetg(1,t_MAT);
2247 561797 : hx = lgcols(x);
2248 561798 : exact = 1;
2249 561798 : emin = HIGHEXPOBIT;
2250 561798 : D = gen_1;
2251 2168066 : for (j = 1; j < lx; j++)
2252 7280648 : for (i = 1; i < hx; i++) rescale_init(gcoeff(x,i,j), &exact, &emin, &D);
2253 561769 : if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
2254 561699 : return grndtoi(gmul2n(x, -emin), NULL);
2255 : }
2256 : GEN
2257 37387 : RgX_rescale_to_int(GEN x)
2258 : {
2259 37387 : long lx = lg(x), i, emin;
2260 : GEN D;
2261 : int exact;
2262 37387 : if (lx == 2) return gcopy(x); /* rare */
2263 37387 : exact = 1;
2264 37387 : emin = HIGHEXPOBIT;
2265 37387 : D = gen_1;
2266 229100 : for (i = 2; i < lx; i++) rescale_init(gel(x,i), &exact, &emin, &D);
2267 37387 : if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
2268 36281 : return grndtoi(gmul2n(x, -emin), NULL);
2269 : }
2270 :
2271 : /* return x * n/d. x: rational; d,n,result: integral; d,n coprime */
2272 : static GEN
2273 9532643 : Q_divmuli_to_int(GEN x, GEN d, GEN n)
2274 : {
2275 : long i, l;
2276 : GEN y, xn, xd;
2277 : pari_sp av;
2278 :
2279 9532643 : switch(typ(x))
2280 : {
2281 2781146 : case t_INT:
2282 2781146 : av = avma; y = diviiexact(x,d);
2283 2781146 : return gerepileuptoint(av, mulii(y,n));
2284 :
2285 4586763 : case t_FRAC:
2286 4586763 : xn = gel(x,1);
2287 4586763 : xd = gel(x,2); av = avma;
2288 4586763 : y = mulii(diviiexact(xn, d), diviiexact(n, xd));
2289 4586763 : return gerepileuptoint(av, y);
2290 :
2291 453057 : case t_VEC: case t_COL: case t_MAT:
2292 453057 : y = cgetg_copy(x, &l);
2293 4095019 : for (i=1; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
2294 453057 : return y;
2295 :
2296 1711677 : case t_POL:
2297 1711677 : y = cgetg_copy(x, &l); y[1] = x[1];
2298 5808911 : for (i=2; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
2299 1711677 : return y;
2300 :
2301 0 : case t_RFRAC:
2302 0 : av = avma;
2303 0 : return gerepileupto(av, gmul(x,mkfrac(n,d)));
2304 :
2305 0 : case t_POLMOD:
2306 0 : retmkpolmod(Q_divmuli_to_int(gel(x,2), d,n), RgX_copy(gel(x,1)));
2307 : }
2308 0 : pari_err_TYPE("Q_divmuli_to_int",x);
2309 : return NULL; /* LCOV_EXCL_LINE */
2310 : }
2311 :
2312 : /* return x / d. x: rational; d,result: integral. */
2313 : static GEN
2314 38174797 : Q_divi_to_int(GEN x, GEN d)
2315 : {
2316 : long i, l;
2317 : GEN y;
2318 :
2319 38174797 : switch(typ(x))
2320 : {
2321 30440551 : case t_INT:
2322 30440551 : return diviiexact(x,d);
2323 :
2324 4044844 : case t_VEC: case t_COL: case t_MAT:
2325 4044844 : y = cgetg_copy(x, &l);
2326 23393518 : for (i=1; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
2327 4043675 : return y;
2328 :
2329 3684820 : case t_POL:
2330 3684820 : y = cgetg_copy(x, &l); y[1] = x[1];
2331 16868963 : for (i=2; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
2332 3684771 : return y;
2333 :
2334 7 : case t_RFRAC:
2335 7 : return gdiv(x,d);
2336 :
2337 5915 : case t_POLMOD:
2338 5915 : retmkpolmod(Q_divi_to_int(gel(x,2), d), RgX_copy(gel(x,1)));
2339 : }
2340 0 : pari_err_TYPE("Q_divi_to_int",x);
2341 : return NULL; /* LCOV_EXCL_LINE */
2342 : }
2343 : /* c t_FRAC */
2344 : static GEN
2345 8310854 : Q_divq_to_int(GEN x, GEN c)
2346 : {
2347 8310854 : GEN n = gel(c,1), d = gel(c,2);
2348 8310854 : if (is_pm1(n)) {
2349 6517404 : GEN y = Q_muli_to_int(x,d);
2350 6517382 : if (signe(n) < 0) y = gneg(y);
2351 6517382 : return y;
2352 : }
2353 1793447 : return Q_divmuli_to_int(x, n,d);
2354 : }
2355 :
2356 : /* return y = x / c, assuming x,c rational, and y integral */
2357 : GEN
2358 4746 : Q_div_to_int(GEN x, GEN c)
2359 : {
2360 4746 : switch(typ(c))
2361 : {
2362 3283 : case t_INT: return Q_divi_to_int(x, c);
2363 1463 : case t_FRAC: return Q_divq_to_int(x, c);
2364 : }
2365 0 : pari_err_TYPE("Q_div_to_int",c);
2366 : return NULL; /* LCOV_EXCL_LINE */
2367 : }
2368 : /* return y = x * c, assuming x,c rational, and y integral */
2369 : GEN
2370 0 : Q_mul_to_int(GEN x, GEN c)
2371 : {
2372 : GEN d, n;
2373 0 : switch(typ(c))
2374 : {
2375 0 : case t_INT: return Q_muli_to_int(x, c);
2376 0 : case t_FRAC:
2377 0 : n = gel(c,1);
2378 0 : d = gel(c,2);
2379 0 : return Q_divmuli_to_int(x, d,n);
2380 : }
2381 0 : pari_err_TYPE("Q_mul_to_int",c);
2382 : return NULL; /* LCOV_EXCL_LINE */
2383 : }
2384 :
2385 : GEN
2386 48512047 : Q_primitive_part(GEN x, GEN *ptc)
2387 : {
2388 48512047 : pari_sp av = avma;
2389 48512047 : GEN c = Q_content_safe(x);
2390 48510386 : if (c)
2391 : {
2392 48510234 : if (typ(c) == t_INT)
2393 : {
2394 40200842 : if (equali1(c)) { set_avma(av); c = NULL; }
2395 6061636 : else if (signe(c)) x = Q_divi_to_int(x, c);
2396 : }
2397 8309392 : else x = Q_divq_to_int(x, c);
2398 : }
2399 48509237 : if (ptc) *ptc = c;
2400 48509237 : return x;
2401 : }
2402 : GEN
2403 5006522 : Q_primpart(GEN x) { return Q_primitive_part(x, NULL); }
2404 : GEN
2405 110801 : vec_Q_primpart(GEN x)
2406 626427 : { pari_APPLY_same(Q_primpart(gel(x,i))) }
2407 : GEN
2408 16128 : row_Q_primpart(GEN M)
2409 16128 : { return shallowtrans(vec_Q_primpart(shallowtrans(M))); }
2410 :
2411 : /*******************************************************************/
2412 : /* */
2413 : /* SUBRESULTANT */
2414 : /* */
2415 : /*******************************************************************/
2416 : /* for internal use */
2417 : GEN
2418 20027971 : gdivexact(GEN x, GEN y)
2419 : {
2420 : long i,lx;
2421 : GEN z;
2422 20027971 : if (gequal1(y)) return x;
2423 20060893 : if (typ(y) == t_POLMOD) return gmul(x, ginv(y));
2424 20060816 : switch(typ(x))
2425 : {
2426 17187252 : case t_INT:
2427 17187252 : if (typ(y)==t_INT) return diviiexact(x,y);
2428 75 : if (!signe(x)) return gen_0;
2429 40 : break;
2430 8316 : case t_INTMOD:
2431 : case t_FFELT:
2432 8316 : case t_POLMOD: return gmul(x,ginv(y));
2433 2867736 : case t_POL:
2434 2867736 : switch(typ(y))
2435 : {
2436 714 : case t_INTMOD:
2437 : case t_FFELT:
2438 714 : case t_POLMOD: return gmul(x,ginv(y));
2439 36386 : case t_POL: { /* not stack-clean */
2440 : long v;
2441 36386 : if (varn(x)!=varn(y)) break;
2442 35420 : v = RgX_valrem(y,&y);
2443 35420 : if (v) x = RgX_shift_shallow(x,-v);
2444 35420 : if (!degpol(y)) { y = gel(y,2); break; }
2445 31564 : return RgX_div(x,y);
2446 : }
2447 42 : case t_RFRAC:
2448 42 : if (varn(gel(y,2)) != varn(x)) break;
2449 42 : return gdiv(x, y);
2450 : }
2451 2835416 : return RgX_Rg_divexact(x, y);
2452 4844 : case t_VEC: case t_COL: case t_MAT:
2453 4844 : lx = lg(x); z = new_chunk(lx);
2454 53060 : for (i=1; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
2455 4844 : z[0] = x[0]; return z;
2456 : }
2457 36 : if (DEBUGLEVEL) pari_warn(warner,"missing case in gdivexact");
2458 36 : return gdiv(x,y);
2459 : }
2460 :
2461 : static GEN
2462 399603 : init_resultant(GEN x, GEN y)
2463 : {
2464 399603 : long tx = typ(x), ty = typ(y), vx, vy;
2465 399603 : if (is_scalar_t(tx) || is_scalar_t(ty))
2466 : {
2467 14 : if (gequal0(x) || gequal0(y)) return gmul(x,y); /* keep type info */
2468 14 : if (tx==t_POL) return gpowgs(y, degpol(x));
2469 0 : if (ty==t_POL) return gpowgs(x, degpol(y));
2470 0 : return gen_1;
2471 : }
2472 399588 : if (tx!=t_POL) pari_err_TYPE("resultant",x);
2473 399588 : if (ty!=t_POL) pari_err_TYPE("resultant",y);
2474 399588 : if (!signe(x) || !signe(y)) return gmul(Rg_get_0(x),Rg_get_0(y)); /*type*/
2475 399511 : vx = varn(x);
2476 399511 : vy = varn(y); if (vx == vy) return NULL;
2477 7 : return (varncmp(vx,vy) < 0)? gpowgs(y,degpol(x)): gpowgs(x,degpol(y));
2478 : }
2479 :
2480 : /* x an RgX, y a scalar */
2481 : static GEN
2482 0 : scalar_res(GEN x, GEN y, GEN *U, GEN *V)
2483 : {
2484 0 : *V = gpowgs(y,degpol(x)-1);
2485 0 : *U = gen_0; return gmul(y, *V);
2486 : }
2487 :
2488 : /* return 0 if the subresultant chain can be interrupted.
2489 : * Set u = NULL if the resultant is 0. */
2490 : static int
2491 11439 : subres_step(GEN *u, GEN *v, GEN *g, GEN *h, GEN *uze, GEN *um1, long *signh)
2492 : {
2493 11439 : GEN u0, c, r, q = RgX_pseudodivrem(*u,*v, &r);
2494 : long du, dv, dr, degq;
2495 :
2496 11439 : if (gequal0(leading_coeff(r))) r = RgX_renormalize(r);
2497 11439 : dr = lg(r); if (!signe(r)) { *u = NULL; return 0; }
2498 11292 : du = degpol(*u);
2499 11292 : dv = degpol(*v);
2500 11292 : degq = du - dv;
2501 11292 : if (*um1 == gen_1)
2502 6266 : u0 = gpowgs(gel(*v,dv+2),degq+1);
2503 5026 : else if (*um1 == gen_0)
2504 2275 : u0 = gen_0;
2505 : else /* except in those 2 cases, um1 is an RgX */
2506 2751 : u0 = RgX_Rg_mul(*um1, gpowgs(gel(*v,dv+2),degq+1));
2507 :
2508 11292 : if (*uze == gen_0) /* except in that case, uze is an RgX */
2509 6266 : u0 = scalarpol(u0, varn(*u)); /* now an RgX */
2510 : else
2511 5026 : u0 = gsub(u0, gmul(q,*uze));
2512 :
2513 11292 : *um1 = *uze;
2514 11292 : *uze = u0; /* uze <- lead(v)^(degq + 1) * um1 - q * uze */
2515 :
2516 11292 : *u = *v; c = *g; *g = leading_coeff(*u);
2517 11292 : switch(degq)
2518 : {
2519 1666 : case 0: break;
2520 7833 : case 1:
2521 7833 : c = gmul(*h,c); *h = *g; break;
2522 1793 : default:
2523 1793 : c = gmul(gpowgs(*h,degq), c);
2524 1793 : *h = gdivexact(gpowgs(*g,degq), gpowgs(*h,degq-1));
2525 : }
2526 11292 : if (typ(c) == t_POLMOD)
2527 : {
2528 918 : c = ginv(c);
2529 918 : *v = RgX_Rg_mul(r,c);
2530 918 : *uze = RgX_Rg_mul(*uze,c);
2531 : }
2532 : else
2533 : {
2534 10374 : *v = RgX_Rg_divexact(r,c);
2535 10374 : *uze= RgX_Rg_divexact(*uze,c);
2536 : }
2537 11292 : if (both_odd(du, dv)) *signh = -*signh;
2538 11292 : return (dr > 3);
2539 : }
2540 :
2541 : /* compute U, V s.t Ux + Vy = resultant(x,y) */
2542 : static GEN
2543 2395 : subresext_i(GEN x, GEN y, GEN *U, GEN *V)
2544 : {
2545 : pari_sp av, av2;
2546 2395 : long dx, dy, du, signh, tx = typ(x), ty = typ(y);
2547 : GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze;
2548 :
2549 2395 : if (!is_extscalar_t(tx)) pari_err_TYPE("subresext",x);
2550 2395 : if (!is_extscalar_t(ty)) pari_err_TYPE("subresext",y);
2551 2395 : if (gequal0(x) || gequal0(y)) { *U = *V = gen_0; return gen_0; }
2552 2395 : if (tx != t_POL) {
2553 0 : if (ty != t_POL) { *U = ginv(x); *V = gen_0; return gen_1; }
2554 0 : return scalar_res(y,x,V,U);
2555 : }
2556 2395 : if (ty != t_POL) return scalar_res(x,y,U,V);
2557 2395 : if (varn(x) != varn(y))
2558 0 : return varncmp(varn(x), varn(y)) < 0? scalar_res(x,y,U,V)
2559 0 : : scalar_res(y,x,V,U);
2560 2395 : if (gequal0(leading_coeff(x))) x = RgX_renormalize(x);
2561 2395 : if (gequal0(leading_coeff(y))) y = RgX_renormalize(y);
2562 2395 : dx = degpol(x);
2563 2395 : dy = degpol(y);
2564 2395 : signh = 1;
2565 2395 : if (dx < dy)
2566 : {
2567 904 : pswap(U,V); lswap(dx,dy); swap(x,y);
2568 904 : if (both_odd(dx, dy)) signh = -signh;
2569 : }
2570 2395 : if (dy == 0)
2571 : {
2572 0 : *V = gpowgs(gel(y,2),dx-1);
2573 0 : *U = gen_0; return gmul(*V,gel(y,2));
2574 : }
2575 2395 : av = avma;
2576 2395 : u = x = primitive_part(x, &cu);
2577 2395 : v = y = primitive_part(y, &cv);
2578 2395 : g = h = gen_1; av2 = avma;
2579 2395 : um1 = gen_1; uze = gen_0;
2580 : for(;;)
2581 : {
2582 7064 : if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
2583 4669 : if (gc_needed(av2,1))
2584 : {
2585 0 : if(DEBUGMEM>1) pari_warn(warnmem,"subresext, dr = %ld", degpol(v));
2586 0 : gerepileall(av2,6, &u,&v,&g,&h,&uze,&um1);
2587 : }
2588 : }
2589 : /* uze an RgX */
2590 2395 : if (!u) { *U = *V = gen_0; return gc_const(av, gen_0); }
2591 2388 : z = gel(v,2); du = degpol(u);
2592 2388 : if (du > 1)
2593 : { /* z = gdivexact(gpowgs(z,du), gpowgs(h,du-1)); */
2594 245 : p1 = gpowgs(gdiv(z,h),du-1);
2595 245 : z = gmul(z,p1);
2596 245 : uze = RgX_Rg_mul(uze, p1);
2597 : }
2598 2388 : if (signh < 0) { z = gneg_i(z); uze = RgX_neg(uze); }
2599 :
2600 2388 : vze = RgX_divrem(Rg_RgX_sub(z, RgX_mul(uze,x)), y, &r);
2601 2388 : if (signe(r)) pari_warn(warner,"inexact computation in subresext");
2602 : /* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */
2603 2388 : p1 = gen_1;
2604 2388 : if (cu) p1 = gmul(p1, gpowgs(cu,dy));
2605 2388 : if (cv) p1 = gmul(p1, gpowgs(cv,dx));
2606 2388 : cu = cu? gdiv(p1,cu): p1;
2607 2388 : cv = cv? gdiv(p1,cv): p1;
2608 2388 : z = gmul(z,p1);
2609 2388 : *U = RgX_Rg_mul(uze,cu);
2610 2388 : *V = RgX_Rg_mul(vze,cv);
2611 2388 : return z;
2612 : }
2613 : GEN
2614 0 : subresext(GEN x, GEN y, GEN *U, GEN *V)
2615 : {
2616 0 : pari_sp av = avma;
2617 0 : GEN z = subresext_i(x, y, U, V);
2618 0 : return gc_all(av, 3, &z, U, V);
2619 : }
2620 :
2621 : static GEN
2622 28 : zero_extgcd(GEN y, GEN *U, GEN *V, long vx)
2623 : {
2624 28 : GEN x=content(y);
2625 28 : *U=pol_0(vx); *V = scalarpol(ginv(x), vx); return gmul(y,*V);
2626 : }
2627 :
2628 : static int
2629 6167 : must_negate(GEN x)
2630 : {
2631 6167 : GEN t = leading_coeff(x);
2632 6167 : switch(typ(t))
2633 : {
2634 4263 : case t_INT: case t_REAL:
2635 4263 : return (signe(t) < 0);
2636 0 : case t_FRAC:
2637 0 : return (signe(gel(t,1)) < 0);
2638 : }
2639 1904 : return 0;
2640 : }
2641 :
2642 : static GEN
2643 63 : gc_gcdext(pari_sp av, GEN r, GEN *u, GEN *v)
2644 : {
2645 63 : if (!u && !v) return gerepileupto(av, r);
2646 63 : if (u && v) return gc_all(av, 3, &r, u, v);
2647 0 : return gc_all(av, 2, &r, u ? u: v);
2648 : }
2649 :
2650 : static GEN
2651 56 : RgX_extgcd_FpX(GEN x, GEN y, GEN p, GEN *u, GEN *v)
2652 : {
2653 56 : pari_sp av = avma;
2654 56 : GEN r = FpX_extgcd(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p, u, v);
2655 56 : if (u) *u = FpX_to_mod(*u, p);
2656 56 : if (v) *v = FpX_to_mod(*v, p);
2657 56 : return gc_gcdext(av, FpX_to_mod(r, p), u, v);
2658 : }
2659 :
2660 : static GEN
2661 7 : RgX_extgcd_FpXQX(GEN x, GEN y, GEN pol, GEN p, GEN *U, GEN *V)
2662 : {
2663 7 : pari_sp av = avma;
2664 7 : GEN r, T = RgX_to_FpX(pol, p);
2665 7 : r = FpXQX_extgcd(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p, U, V);
2666 7 : return gc_gcdext(av, FpXQX_to_mod(r, T, p), U, V);
2667 : }
2668 :
2669 : static GEN
2670 4221 : RgX_extgcd_fast(GEN x, GEN y, GEN *U, GEN *V)
2671 : {
2672 : GEN p, pol;
2673 : long pa;
2674 4221 : long t = RgX_type(x, &p,&pol,&pa);
2675 4221 : switch(t)
2676 : {
2677 21 : case t_FFELT: return FFX_extgcd(x, y, pol, U, V);
2678 56 : case t_INTMOD: return RgX_extgcd_FpX(x, y, p, U, V);
2679 7 : case code(t_POLMOD, t_INTMOD):
2680 7 : return RgX_extgcd_FpXQX(x, y, pol, p, U, V);
2681 4137 : default: return NULL;
2682 : }
2683 : }
2684 :
2685 : /* compute U, V s.t Ux + Vy = GCD(x,y) using subresultant */
2686 : GEN
2687 4256 : RgX_extgcd(GEN x, GEN y, GEN *U, GEN *V)
2688 : {
2689 : pari_sp av, av2, tetpil;
2690 : long signh; /* junk */
2691 4256 : long dx, dy, vx, tx = typ(x), ty = typ(y);
2692 : GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze, *gptr[3];
2693 :
2694 4256 : if (tx!=t_POL) pari_err_TYPE("RgX_extgcd",x);
2695 4256 : if (ty!=t_POL) pari_err_TYPE("RgX_extgcd",y);
2696 4256 : if ( varncmp(varn(x),varn(y))) pari_err_VAR("RgX_extgcd",x,y);
2697 4256 : vx=varn(x);
2698 4256 : if (!signe(x))
2699 : {
2700 14 : if (signe(y)) return zero_extgcd(y,U,V,vx);
2701 7 : *U = pol_0(vx); *V = pol_0(vx);
2702 7 : return pol_0(vx);
2703 : }
2704 4242 : if (!signe(y)) return zero_extgcd(x,V,U,vx);
2705 4221 : r = RgX_extgcd_fast(x, y, U, V);
2706 4221 : if (r) return r;
2707 4137 : dx = degpol(x); dy = degpol(y);
2708 4137 : if (dx < dy) { pswap(U,V); lswap(dx,dy); swap(x,y); }
2709 4137 : if (dy==0) { *U=pol_0(vx); *V=ginv(y); return pol_1(vx); }
2710 :
2711 3948 : av = avma;
2712 3948 : u = x = primitive_part(x, &cu);
2713 3948 : v = y = primitive_part(y, &cv);
2714 3948 : g = h = gen_1; av2 = avma;
2715 3948 : um1 = gen_1; uze = gen_0;
2716 : for(;;)
2717 : {
2718 4088 : if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
2719 140 : if (gc_needed(av2,1))
2720 : {
2721 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_extgcd, dr = %ld",degpol(v));
2722 0 : gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
2723 : }
2724 : }
2725 3948 : if (uze != gen_0) {
2726 : GEN r;
2727 3808 : vze = RgX_divrem(RgX_sub(v, RgX_mul(uze,x)), y, &r);
2728 3808 : if (signe(r)) pari_warn(warner,"inexact computation in RgX_extgcd");
2729 3808 : if (cu) uze = RgX_Rg_div(uze,cu);
2730 3808 : if (cv) vze = RgX_Rg_div(vze,cv);
2731 3808 : p1 = ginv(content(v));
2732 : }
2733 : else /* y | x */
2734 : {
2735 140 : vze = cv ? RgX_Rg_div(pol_1(vx),cv): pol_1(vx);
2736 140 : uze = pol_0(vx);
2737 140 : p1 = gen_1;
2738 : }
2739 3948 : if (must_negate(v)) p1 = gneg(p1);
2740 3948 : tetpil = avma;
2741 3948 : z = RgX_Rg_mul(v,p1);
2742 3948 : *U = RgX_Rg_mul(uze,p1);
2743 3948 : *V = RgX_Rg_mul(vze,p1);
2744 3948 : gptr[0] = &z;
2745 3948 : gptr[1] = U;
2746 3948 : gptr[2] = V;
2747 3948 : gerepilemanysp(av,tetpil,gptr,3); return z;
2748 : }
2749 :
2750 : static GEN
2751 14 : RgX_halfgcd_i(GEN a, GEN b)
2752 : {
2753 14 : pari_sp av=avma;
2754 14 : long m = degpol(a), va = varn(a);
2755 : GEN u,u1,v,v1;
2756 14 : u1 = v = pol_0(va);
2757 14 : u = v1 = pol_1(va);
2758 14 : if (degpol(a)<degpol(b))
2759 : {
2760 0 : swap(a,b);
2761 0 : swap(u,v); swap(u1,v1);
2762 : }
2763 42 : while (2*degpol(b) >= m)
2764 : {
2765 28 : GEN r, q = RgX_pseudodivrem(a,b,&r);
2766 28 : GEN l = gpowgs(leading_coeff(b), degpol(a)-degpol(b)+1);
2767 28 : GEN g = ggcd(l, content(r));
2768 28 : q = RgX_Rg_div(q, g);
2769 28 : r = RgX_Rg_div(r, g);
2770 28 : l = gdiv(l, g);
2771 28 : a = b; b = r; swap(u,v); swap(u1,v1);
2772 28 : v = RgX_sub(gmul(l,v), RgX_mul(u, q));
2773 28 : v1 = RgX_sub(gmul(l,v1), RgX_mul(u1, q));
2774 28 : if (gc_needed(av,2))
2775 : {
2776 0 : if (DEBUGMEM>1) pari_warn(warnmem,"halfgcd (d = %ld)",degpol(b));
2777 0 : gerepileall(av,6, &a,&b,&u1,&v1,&u,&v);
2778 : }
2779 : }
2780 14 : return gerepilecopy(av, mkvec2(mkmat22(u,u1,v,v1), mkcol2(a, b)));
2781 : }
2782 :
2783 : static GEN
2784 7 : RgX_halfgcd_FFX(GEN x, GEN y, GEN fa)
2785 : {
2786 7 : pari_sp av = avma;
2787 7 : GEN M = FFX_halfgcd(x, y, fa);
2788 7 : GEN a = FFX_add(FFX_mul(gcoeff(M,1,1), x, fa),
2789 7 : FFX_mul(gcoeff(M,1,2), y, fa), fa);
2790 7 : GEN b = FFX_add(FFX_mul(gcoeff(M,2,1), x, fa),
2791 7 : FFX_mul(gcoeff(M,2,2), y, fa), fa);
2792 7 : return gerepilecopy(av, mkvec2(M, mkcol2(a, b)));
2793 : }
2794 :
2795 : static GEN
2796 14 : RgX_halfgcd_FpX(GEN x, GEN y, GEN p)
2797 : {
2798 14 : pari_sp av = avma;
2799 : GEN M, V, a, b;
2800 14 : if (lgefint(p) == 3)
2801 : {
2802 7 : ulong pp = uel(p, 2);
2803 7 : GEN xp = RgX_to_Flx(x, pp), yp = RgX_to_Flx(y, pp);
2804 7 : M = Flx_halfgcd(xp, yp, pp);
2805 7 : a = Flx_add(Flx_mul(gcoeff(M,1,1), xp, pp),
2806 7 : Flx_mul(gcoeff(M,1,2), yp, pp), pp);
2807 7 : b = Flx_add(Flx_mul(gcoeff(M,2,1), xp, pp),
2808 7 : Flx_mul(gcoeff(M,2,2), yp, pp), pp);
2809 7 : M = FlxM_to_ZXM(M); a = Flx_to_ZX(a); b = Flx_to_ZX(b);
2810 : }
2811 : else
2812 : {
2813 7 : x = RgX_to_FpX(x, p); y = RgX_to_FpX(y, p);
2814 7 : M = FpX_halfgcd(x, y, p);
2815 7 : a = FpX_add(FpX_mul(gcoeff(M,1,1), x, p),
2816 7 : FpX_mul(gcoeff(M,1,2), y, p), p);
2817 7 : b = FpX_add(FpX_mul(gcoeff(M,2,1), x, p),
2818 7 : FpX_mul(gcoeff(M,2,2), y, p), p);
2819 : }
2820 14 : V = mkcol2(a, b);
2821 14 : return gerepilecopy(av, mkvec2(FpXM_to_mod(M, p), FpXC_to_mod(V, p)));
2822 : }
2823 :
2824 : static GEN
2825 0 : RgX_halfgcd_FpXQX(GEN x, GEN y, GEN pol, GEN p)
2826 : {
2827 0 : pari_sp av = avma;
2828 0 : GEN M, V, a, b, T = RgX_to_FpX(pol, p);
2829 0 : if (signe(T)==0) pari_err_OP("halfgcd", x, y);
2830 0 : x = RgX_to_FpXQX(x, T, p); y = RgX_to_FpXQX(y, T, p);
2831 0 : M = FpXQX_halfgcd(x, y, T, p);
2832 0 : a = FpXX_add(FpXQX_mul(gcoeff(M,1,1), x, T, p), FpXQX_mul(gcoeff(M,1,2), y, T, p), p);
2833 0 : b = FpXX_add(FpXQX_mul(gcoeff(M,2,1), x, T, p), FpXQX_mul(gcoeff(M,2,2), y, T, p), p);
2834 0 : V = mkcol2(a, b);
2835 0 : return gerepilecopy(av, mkvec2(FqXM_to_mod(M, T, p), FqXC_to_mod(V, T, p)));
2836 : }
2837 :
2838 : static GEN
2839 35 : RgX_halfgcd_fast(GEN x, GEN y)
2840 : {
2841 : GEN p, pol;
2842 : long pa;
2843 35 : long t = RgX_type2(x,y, &p,&pol,&pa);
2844 35 : switch(t)
2845 : {
2846 7 : case t_FFELT: return RgX_halfgcd_FFX(x, y, pol);
2847 14 : case t_INTMOD: return RgX_halfgcd_FpX(x, y, p);
2848 0 : case code(t_POLMOD, t_INTMOD):
2849 0 : return RgX_halfgcd_FpXQX(x, y, pol, p);
2850 14 : default: return NULL;
2851 : }
2852 : }
2853 :
2854 : GEN
2855 35 : RgX_halfgcd(GEN a, GEN b)
2856 : {
2857 35 : GEN z = RgX_halfgcd_fast(a,b);
2858 35 : if (z) return z;
2859 14 : return RgX_halfgcd_i(a, b);
2860 : }
2861 :
2862 : int
2863 70 : RgXQ_ratlift(GEN x, GEN T, long amax, long bmax, GEN *P, GEN *Q)
2864 : {
2865 70 : pari_sp av = avma, av2, tetpil;
2866 : long signh; /* junk */
2867 : long vx;
2868 : GEN g, h, p1, cu, cv, u, v, um1, uze, *gptr[2];
2869 :
2870 70 : if (typ(x)!=t_POL) pari_err_TYPE("RgXQ_ratlift",x);
2871 70 : if (typ(T)!=t_POL) pari_err_TYPE("RgXQ_ratlift",T);
2872 70 : if ( varncmp(varn(x),varn(T)) ) pari_err_VAR("RgXQ_ratlift",x,T);
2873 70 : if (bmax < 0) pari_err_DOMAIN("ratlift", "bmax", "<", gen_0, stoi(bmax));
2874 70 : if (!signe(T)) {
2875 0 : if (degpol(x) <= amax) {
2876 0 : *P = RgX_copy(x);
2877 0 : *Q = pol_1(varn(x));
2878 0 : return 1;
2879 : }
2880 0 : return 0;
2881 : }
2882 70 : if (amax+bmax >= degpol(T))
2883 0 : pari_err_DOMAIN("ratlift", "amax+bmax", ">=", stoi(degpol(T)),
2884 : mkvec3(stoi(amax), stoi(bmax), T));
2885 70 : vx = varn(T);
2886 70 : u = x = primitive_part(x, &cu);
2887 70 : v = T = primitive_part(T, &cv);
2888 70 : g = h = gen_1; av2 = avma;
2889 70 : um1 = gen_1; uze = gen_0;
2890 : for(;;)
2891 : {
2892 287 : (void) subres_step(&u, &v, &g, &h, &uze, &um1, &signh);
2893 287 : if (!u || (typ(uze)==t_POL && degpol(uze)>bmax)) return gc_bool(av,0);
2894 287 : if (typ(v)!=t_POL || degpol(v)<=amax) break;
2895 217 : if (gc_needed(av2,1))
2896 : {
2897 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXQ_ratlift, dr = %ld", degpol(v));
2898 0 : gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
2899 : }
2900 : }
2901 70 : if (uze == gen_0)
2902 : {
2903 0 : set_avma(av); *P = pol_0(vx); *Q = pol_1(vx);
2904 0 : return 1;
2905 : }
2906 70 : if (cu) uze = RgX_Rg_div(uze,cu);
2907 70 : p1 = ginv(content(v));
2908 70 : if (must_negate(v)) p1 = gneg(p1);
2909 70 : tetpil = avma;
2910 70 : *P = RgX_Rg_mul(v,p1);
2911 70 : *Q = RgX_Rg_mul(uze,p1);
2912 70 : gptr[0] = P;
2913 70 : gptr[1] = Q;
2914 70 : gerepilemanysp(av,tetpil,gptr,2); return 1;
2915 : }
2916 :
2917 : GEN
2918 0 : RgX_chinese_coprime(GEN x, GEN y, GEN Tx, GEN Ty, GEN Tz)
2919 : {
2920 0 : pari_sp av = avma;
2921 0 : GEN ax = RgX_mul(RgXQ_inv(Tx,Ty), Tx);
2922 0 : GEN p1 = RgX_mul(ax, RgX_sub(y,x));
2923 0 : p1 = RgX_add(x,p1);
2924 0 : if (!Tz) Tz = RgX_mul(Tx,Ty);
2925 0 : p1 = RgX_rem(p1, Tz);
2926 0 : return gerepileupto(av,p1);
2927 : }
2928 :
2929 : /*******************************************************************/
2930 : /* */
2931 : /* RESULTANT USING DUCOS VARIANT */
2932 : /* */
2933 : /*******************************************************************/
2934 : /* x^n / y^(n-1), assume n > 0 */
2935 : static GEN
2936 129121 : Lazard(GEN x, GEN y, long n)
2937 : {
2938 : long a;
2939 : GEN c;
2940 :
2941 129121 : if (n == 1) return x;
2942 8147 : a = 1 << expu(n); /* a = 2^k <= n < 2^(k+1) */
2943 8147 : c=x; n-=a;
2944 16918 : while (a>1)
2945 : {
2946 8770 : a>>=1; c=gdivexact(gsqr(c),y);
2947 8771 : if (n>=a) { c=gdivexact(gmul(c,x),y); n -= a; }
2948 : }
2949 8148 : return c;
2950 : }
2951 :
2952 : /* F (x/y)^(n-1), assume n >= 1 */
2953 : static GEN
2954 99226 : Lazard2(GEN F, GEN x, GEN y, long n)
2955 : {
2956 99226 : if (n == 1) return F;
2957 1228 : return RgX_Rg_divexact(RgX_Rg_mul(F, Lazard(x,y,n-1)), y);
2958 : }
2959 :
2960 : static GEN
2961 99226 : RgX_neg_i(GEN x, long lx)
2962 : {
2963 : long i;
2964 99226 : GEN y = cgetg(lx, t_POL); y[1] = x[1];
2965 220511 : for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
2966 99226 : return y;
2967 : }
2968 : static GEN
2969 259913 : RgX_Rg_mul_i(GEN y, GEN x, long ly)
2970 : {
2971 : long i;
2972 : GEN z;
2973 259913 : if (isrationalzero(x)) return pol_0(varn(y));
2974 259853 : z = cgetg(ly,t_POL); z[1] = y[1];
2975 578408 : for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
2976 259846 : return z;
2977 : }
2978 : static long
2979 256593 : reductum_lg(GEN x, long lx)
2980 : {
2981 256593 : long i = lx-2;
2982 305535 : while (i > 1 && gequal0(gel(x,i))) i--;
2983 256593 : return i+1;
2984 : }
2985 :
2986 : #define addshift(x,y) RgX_addmulXn_shallow((x),(y),1)
2987 : /* delta = deg(P) - deg(Q) > 0, deg(Q) > 0, P,Q,Z t_POL in the same variable,
2988 : * s "scalar". Return prem(P, -Q) / s^delta lc(P) */
2989 : static GEN
2990 99226 : nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
2991 : {
2992 99226 : GEN p0, q0, h0, TMP, H, A, z0 = leading_coeff(Z);
2993 : long p, q, j, lP, lQ;
2994 : pari_sp av;
2995 :
2996 99226 : p = degpol(P); p0 = gel(P,p+2); lP = reductum_lg(P,lg(P));
2997 99225 : q = degpol(Q); q0 = gel(Q,q+2); lQ = reductum_lg(Q,lg(Q));
2998 : /* p > q. Very often p - 1 = q */
2999 99226 : av = avma;
3000 : /* H = RgX_neg(reductum(Z)) optimized, using Q ~ Z */
3001 99226 : H = RgX_neg_i(Z, lQ); /* deg H < q */
3002 :
3003 99226 : A = (q+2 < lP)? RgX_Rg_mul_i(H, gel(P,q+2), lQ): NULL;
3004 103561 : for (j = q+1; j < p; j++)
3005 : {
3006 4336 : if (degpol(H) == q-1)
3007 : { /* h0 = coeff of degree q-1 = leading coeff */
3008 3713 : h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1);
3009 3713 : H = addshift(H, RgX_Rg_divexact(RgX_Rg_mul_i(Q, gneg(h0), lQ), q0));
3010 : }
3011 : else
3012 623 : H = RgX_shift_shallow(H, 1);
3013 4336 : if (j+2 < lP)
3014 : {
3015 3249 : TMP = RgX_Rg_mul(H, gel(P,j+2));
3016 3249 : A = A? RgX_add(A, TMP): TMP;
3017 : }
3018 4336 : if (gc_needed(av,1))
3019 : {
3020 147 : if(DEBUGMEM>1) pari_warn(warnmem,"nextSousResultant j = %ld/%ld",j,p);
3021 147 : gerepileall(av,A?2:1,&H,&A);
3022 : }
3023 : }
3024 99225 : if (q+2 < lP) lP = reductum_lg(P, q+3);
3025 99225 : TMP = RgX_Rg_mul_i(P, z0, lP);
3026 99223 : A = A? RgX_add(A, TMP): TMP;
3027 99222 : A = RgX_Rg_divexact(A, p0);
3028 99225 : if (degpol(H) == q-1)
3029 : {
3030 98833 : h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1); /* destroy old H */
3031 98832 : A = RgX_add(RgX_Rg_mul(addshift(H,A),q0), RgX_Rg_mul_i(Q, gneg(h0), lQ));
3032 : }
3033 : else
3034 392 : A = RgX_Rg_mul(addshift(H,A), q0);
3035 99226 : return RgX_Rg_divexact(A, s);
3036 : }
3037 : #undef addshift
3038 :
3039 : /* Ducos's subresultant */
3040 : GEN
3041 128541 : RgX_resultant_all(GEN P, GEN Q, GEN *sol)
3042 : {
3043 : pari_sp av, av2;
3044 128541 : long dP, dQ, delta, sig = 1;
3045 : GEN cP, cQ, Z, s;
3046 :
3047 128541 : dP = degpol(P);
3048 128541 : dQ = degpol(Q); delta = dP - dQ;
3049 128541 : if (delta < 0)
3050 : {
3051 595 : if (both_odd(dP, dQ)) sig = -1;
3052 595 : swap(P,Q); lswap(dP, dQ); delta = -delta;
3053 : }
3054 128541 : if (sol) *sol = gen_0;
3055 128541 : av = avma;
3056 128541 : if (dQ <= 0)
3057 : {
3058 651 : if (dQ < 0) return Rg_get_0(P);
3059 651 : s = gpowgs(gel(Q,2), dP);
3060 651 : if (sig == -1) s = gerepileupto(av, gneg(s));
3061 651 : return s;
3062 : }
3063 : /* primitive_part is also possible here, but possibly very costly,
3064 : * and hardly ever worth it */
3065 127890 : P = Q_primitive_part(P, &cP);
3066 127893 : Q = Q_primitive_part(Q, &cQ);
3067 127890 : av2 = avma;
3068 127890 : s = gpowgs(leading_coeff(Q),delta);
3069 127892 : if (both_odd(dP, dQ)) sig = -sig;
3070 127892 : Z = Q;
3071 127892 : Q = RgX_pseudorem(P, Q);
3072 127893 : P = Z;
3073 227119 : while(degpol(Q) > 0)
3074 : {
3075 99225 : delta = degpol(P) - degpol(Q); /* > 0 */
3076 99225 : Z = Lazard2(Q, leading_coeff(Q), s, delta);
3077 99226 : if (both_odd(degpol(P), degpol(Q))) sig = -sig;
3078 99226 : Q = nextSousResultant(P, Q, Z, s);
3079 99226 : P = Z;
3080 99226 : if (gc_needed(av,1))
3081 : {
3082 13 : if(DEBUGMEM>1) pari_warn(warnmem,"resultant_all, degpol Q = %ld",degpol(Q));
3083 13 : gerepileall(av2,2,&P,&Q);
3084 : }
3085 99226 : s = leading_coeff(P);
3086 : }
3087 127894 : if (!signe(Q)) { set_avma(av); return Rg_get_0(Q); }
3088 127894 : s = Lazard(leading_coeff(Q), s, degpol(P));
3089 127894 : if (sig == -1) s = gneg(s);
3090 127894 : if (cP) s = gmul(s, gpowgs(cP,dQ));
3091 127894 : if (cQ) s = gmul(s, gpowgs(cQ,dP));
3092 127894 : if (!sol) return gerepilecopy(av, s);
3093 9499 : *sol = P; return gc_all(av, 2, &s, sol);
3094 : }
3095 :
3096 : static GEN
3097 14 : RgX_resultant_FpX(GEN x, GEN y, GEN p)
3098 : {
3099 14 : pari_sp av = avma;
3100 : GEN r;
3101 14 : if (lgefint(p) == 3)
3102 : {
3103 7 : ulong pp = uel(p, 2);
3104 7 : r = utoi(Flx_resultant(RgX_to_Flx(x, pp), RgX_to_Flx(y, pp), pp));
3105 : }
3106 : else
3107 7 : r = FpX_resultant(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3108 14 : return gerepileupto(av, Fp_to_mod(r, p));
3109 : }
3110 :
3111 : static GEN
3112 21 : RgX_resultant_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3113 : {
3114 21 : pari_sp av = avma;
3115 21 : GEN r, T = RgX_to_FpX(pol, p);
3116 21 : r = FpXQX_resultant(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
3117 21 : return gerepileupto(av, FpX_to_mod(r, p));
3118 : }
3119 :
3120 : static GEN
3121 399581 : resultant_fast(GEN x, GEN y)
3122 : {
3123 : GEN p, pol;
3124 : long pa, t;
3125 399581 : p = init_resultant(x,y);
3126 399581 : if (p) return p;
3127 399483 : t = RgX_type2(x,y, &p,&pol,&pa);
3128 399484 : switch(t)
3129 : {
3130 2660 : case t_INT: return ZX_resultant(x,y);
3131 56 : case t_FRAC: return QX_resultant(x,y);
3132 21 : case t_FFELT: return FFX_resultant(x,y,pol);
3133 14 : case t_INTMOD: return RgX_resultant_FpX(x, y, p);
3134 21 : case code(t_POLMOD, t_INTMOD):
3135 21 : return RgX_resultant_FpXQX(x, y, pol, p);
3136 396712 : default: return NULL;
3137 : }
3138 : }
3139 :
3140 : static GEN
3141 278334 : RgX_resultant_sylvester(GEN x, GEN y)
3142 : {
3143 278334 : pari_sp av = avma;
3144 278334 : return gerepileupto(av, det(RgX_sylvestermatrix(x,y)));
3145 : }
3146 :
3147 : /* Return resultant(P,Q).
3148 : * Uses Sylvester's matrix if P or Q inexact, a modular algorithm if they
3149 : * are in Q[X], and Ducos/Lazard optimization of the subresultant algorithm
3150 : * in the "generic" case. */
3151 : GEN
3152 399580 : resultant(GEN P, GEN Q)
3153 : {
3154 399580 : GEN z = resultant_fast(P,Q);
3155 399582 : if (z) return z;
3156 396712 : if (isinexact(P) || isinexact(Q)) return RgX_resultant_sylvester(P,Q);
3157 118399 : return RgX_resultant_all(P, Q, NULL);
3158 : }
3159 :
3160 : /*******************************************************************/
3161 : /* */
3162 : /* RESULTANT USING SYLVESTER MATRIX */
3163 : /* */
3164 : /*******************************************************************/
3165 : static GEN
3166 531622 : syl_RgC(GEN x, long j, long d, long D, long cp)
3167 : {
3168 531622 : GEN c = cgetg(d+1,t_COL);
3169 : long i;
3170 1380827 : for (i=1; i< j; i++) gel(c,i) = gen_0;
3171 2978318 : for ( ; i<=D; i++) { GEN t = gel(x,D-i+2); gel(c,i) = cp? gcopy(t): t; }
3172 1380827 : for ( ; i<=d; i++) gel(c,i) = gen_0;
3173 531622 : return c;
3174 : }
3175 : static GEN
3176 278341 : syl_RgM(GEN x, GEN y, long cp)
3177 : {
3178 278341 : long j, d, dx = degpol(x), dy = degpol(y);
3179 : GEN M;
3180 278341 : if (dx < 0) return dy < 0? cgetg(1,t_MAT): zeromat(dy,dy);
3181 278341 : if (dy < 0) return zeromat(dx,dx);
3182 278341 : d = dx+dy; M = cgetg(d+1,t_MAT);
3183 684537 : for (j=1; j<=dy; j++) gel(M,j) = syl_RgC(x,j,d,j+dx, cp);
3184 403767 : for (j=1; j<=dx; j++) gel(M,j+dy) = syl_RgC(y,j,d,j+dy, cp);
3185 278341 : return M;
3186 : }
3187 : GEN
3188 278334 : RgX_sylvestermatrix(GEN x, GEN y) { return syl_RgM(x,y,0); }
3189 : GEN
3190 7 : sylvestermatrix(GEN x, GEN y)
3191 : {
3192 7 : if (typ(x)!=t_POL) pari_err_TYPE("sylvestermatrix",x);
3193 7 : if (typ(y)!=t_POL) pari_err_TYPE("sylvestermatrix",y);
3194 7 : if (varn(x) != varn(y)) pari_err_VAR("sylvestermatrix",x,y);
3195 7 : return syl_RgM(x,y,1);
3196 : }
3197 :
3198 : GEN
3199 21 : resultant2(GEN x, GEN y)
3200 : {
3201 21 : GEN r = init_resultant(x,y);
3202 21 : return r? r: RgX_resultant_sylvester(x,y);
3203 : }
3204 :
3205 : /* let vx = main variable of x, v0 a variable of highest priority;
3206 : * return a t_POL in variable v0:
3207 : * if vx <= v, return subst(x, v, pol_x(v0))
3208 : * if vx > v, return scalarpol(x, v0) */
3209 : static GEN
3210 329 : fix_pol(GEN x, long v, long v0)
3211 : {
3212 329 : long vx, tx = typ(x);
3213 329 : if (tx != t_POL)
3214 35 : vx = gvar(x);
3215 : else
3216 : { /* shortcut: almost nothing to do */
3217 294 : vx = varn(x);
3218 294 : if (v == vx)
3219 : {
3220 119 : if (v0 != v) { x = leafcopy(x); setvarn(x, v0); }
3221 119 : return x;
3222 : }
3223 : }
3224 210 : if (varncmp(v, vx) > 0)
3225 : {
3226 203 : x = gsubst(x, v, pol_x(v0));
3227 203 : if (typ(x) != t_POL) vx = gvar(x);
3228 : else
3229 : {
3230 196 : vx = varn(x);
3231 196 : if (vx == v0) return x;
3232 : }
3233 : }
3234 49 : if (varncmp(vx, v0) <= 0) pari_err_TYPE("polresultant", x);
3235 42 : return scalarpol_shallow(x, v0);
3236 : }
3237 :
3238 : /* resultant of x and y with respect to variable v, or with respect to their
3239 : * main variable if v < 0. */
3240 : GEN
3241 476 : polresultant0(GEN x, GEN y, long v, long flag)
3242 : {
3243 476 : pari_sp av = avma;
3244 :
3245 476 : if (v >= 0)
3246 : {
3247 140 : long v0 = fetch_var_higher();
3248 140 : x = fix_pol(x,v, v0);
3249 140 : y = fix_pol(y,v, v0);
3250 : }
3251 476 : switch(flag)
3252 : {
3253 469 : case 2:
3254 469 : case 0: x=resultant(x,y); break;
3255 7 : case 1: x=resultant2(x,y); break;
3256 0 : default: pari_err_FLAG("polresultant");
3257 : }
3258 476 : if (v >= 0) (void)delete_var();
3259 476 : return gerepileupto(av,x);
3260 : }
3261 : GEN
3262 1491 : polresultantext0(GEN x, GEN y, long v)
3263 : {
3264 : GEN R, U, V;
3265 1491 : pari_sp av = avma;
3266 :
3267 1491 : if (v >= 0)
3268 : {
3269 14 : long v0 = fetch_var_higher();
3270 14 : x = fix_pol(x,v, v0);
3271 14 : y = fix_pol(y,v, v0);
3272 : }
3273 1491 : R = subresext_i(x,y, &U,&V);
3274 1491 : if (v >= 0)
3275 : {
3276 14 : (void)delete_var();
3277 14 : if (typ(U) == t_POL && varn(U) != v) U = poleval(U, pol_x(v));
3278 14 : if (typ(V) == t_POL && varn(V) != v) V = poleval(V, pol_x(v));
3279 : }
3280 1491 : return gerepilecopy(av, mkvec3(U,V,R));
3281 : }
3282 : GEN
3283 1463 : polresultantext(GEN x, GEN y) { return polresultantext0(x,y,-1); }
3284 :
3285 : /*******************************************************************/
3286 : /* */
3287 : /* CHARACTERISTIC POLYNOMIAL USING RESULTANT */
3288 : /* */
3289 : /*******************************************************************/
3290 :
3291 : /* (v - x)^d */
3292 : static GEN
3293 140 : caract_const(pari_sp av, GEN x, long v, long d)
3294 140 : { return gerepileupto(av, gpowgs(deg1pol_shallow(gen_1, gneg_i(x), v), d)); }
3295 :
3296 : /* return caract(Mod(x,T)) in variable v */
3297 : GEN
3298 116851 : RgXQ_charpoly(GEN x, GEN T, long v)
3299 : {
3300 116851 : pari_sp av = avma;
3301 116851 : long d = degpol(T), dx, vx, vp, v0;
3302 : GEN ch, L;
3303 :
3304 116852 : if (typ(x) != t_POL) return caract_const(av, x, v, d);
3305 116838 : vx = varn(x);
3306 116838 : vp = varn(T);
3307 116838 : if (varncmp(vx, vp) > 0) return caract_const(av, x, v, d);
3308 116838 : if (varncmp(vx, vp) < 0) pari_err_PRIORITY("RgXQ_charpoly", x, "<", vp);
3309 116838 : dx = degpol(x);
3310 116838 : if (dx >= degpol(T)) { x = RgX_rem(x, T); dx = degpol(x); }
3311 116838 : if (dx <= 0) return dx? pol_xn(d, v): caract_const(av, gel(x,2), v, d);
3312 :
3313 116768 : v0 = fetch_var_higher();
3314 116768 : x = RgX_neg(x);
3315 116767 : gel(x,2) = gadd(gel(x,2), pol_x(v));
3316 116766 : setvarn(x, v0);
3317 116766 : T = leafcopy(T); setvarn(T, v0);
3318 116766 : ch = resultant(T, x);
3319 116771 : (void)delete_var();
3320 : /* test for silly input: x mod (deg 0 polynomial) */
3321 116771 : if (typ(ch) != t_POL)
3322 7 : pari_err_PRIORITY("RgXQ_charpoly", pol_x(v), "<", gvar(ch));
3323 116764 : L = leading_coeff(ch);
3324 116764 : if (!gequal1(L)) ch = RgX_Rg_div(ch, L);
3325 116764 : return gerepileupto(av, ch);
3326 : }
3327 :
3328 : /* characteristic polynomial (in v) of x over nf, where x is an element of the
3329 : * algebra nf[t]/(Q(t)) */
3330 : GEN
3331 224 : rnfcharpoly(GEN nf, GEN Q, GEN x, long v)
3332 : {
3333 224 : const char *f = "rnfcharpoly";
3334 224 : long dQ = degpol(Q);
3335 224 : pari_sp av = avma;
3336 : GEN T;
3337 :
3338 224 : if (v < 0) v = 0;
3339 224 : nf = checknf(nf); T = nf_get_pol(nf);
3340 224 : Q = RgX_nffix(f, T,Q,0);
3341 224 : switch(typ(x))
3342 : {
3343 28 : case t_INT:
3344 28 : case t_FRAC: return caract_const(av, x, v, dQ);
3345 91 : case t_POLMOD:
3346 91 : x = polmod_nffix2(f,T,Q, x,0);
3347 56 : break;
3348 56 : case t_POL:
3349 56 : x = varn(x) == varn(T)? Rg_nffix(f,T,x,0): RgX_nffix(f, T,x,0);
3350 42 : break;
3351 49 : default: pari_err_TYPE(f,x);
3352 : }
3353 98 : if (typ(x) != t_POL) return caract_const(av, x, v, dQ);
3354 : /* x a t_POL in variable vQ */
3355 56 : if (degpol(x) >= dQ) x = RgX_rem(x, Q);
3356 56 : if (dQ <= 1) return caract_const(av, constant_coeff(x), v, 1);
3357 56 : return gerepilecopy(av, lift_if_rational( RgXQ_charpoly(x, Q, v) ));
3358 : }
3359 :
3360 : /*******************************************************************/
3361 : /* */
3362 : /* GCD USING SUBRESULTANT */
3363 : /* */
3364 : /*******************************************************************/
3365 : static int inexact(GEN x, int *simple);
3366 : static int
3367 38367 : isinexactall(GEN x, int *simple)
3368 : {
3369 38367 : long i, lx = lg(x);
3370 139020 : for (i=2; i<lx; i++)
3371 100667 : if (inexact(gel(x,i), simple)) return 1;
3372 38353 : return 0;
3373 : }
3374 : /* return 1 if coeff explosion is not possible */
3375 : static int
3376 100919 : inexact(GEN x, int *simple)
3377 : {
3378 100919 : int junk = 0;
3379 100919 : switch(typ(x))
3380 : {
3381 64729 : case t_INT: case t_FRAC: return 0;
3382 :
3383 7 : case t_REAL: case t_PADIC: case t_SER: return 1;
3384 :
3385 4011 : case t_INTMOD:
3386 : case t_FFELT:
3387 4011 : if (!*simple) *simple = 1;
3388 4011 : return 0;
3389 :
3390 77 : case t_COMPLEX:
3391 77 : return inexact(gel(x,1), simple)
3392 77 : || inexact(gel(x,2), simple);
3393 0 : case t_QUAD:
3394 0 : *simple = 0;
3395 0 : return inexact(gel(x,2), &junk)
3396 0 : || inexact(gel(x,3), &junk);
3397 :
3398 819 : case t_POLMOD:
3399 819 : return isinexactall(gel(x,1), simple);
3400 31227 : case t_POL:
3401 31227 : *simple = -1;
3402 31227 : return isinexactall(x, &junk);
3403 49 : case t_RFRAC:
3404 49 : *simple = -1;
3405 49 : return inexact(gel(x,1), &junk)
3406 49 : || inexact(gel(x,2), &junk);
3407 : }
3408 0 : *simple = -1; return 0;
3409 : }
3410 :
3411 : /* x monomial, y t_POL in the same variable */
3412 : static GEN
3413 1706468 : gcdmonome(GEN x, GEN y)
3414 : {
3415 1706468 : pari_sp av = avma;
3416 1706468 : long dx = degpol(x), e = RgX_valrem(y, &y);
3417 1706468 : long i, l = lg(y);
3418 1706468 : GEN t, v = cgetg(l, t_VEC);
3419 1706468 : gel(v,1) = gel(x,dx+2);
3420 3426099 : for (i = 2; i < l; i++) gel(v,i) = gel(y,i);
3421 1706468 : t = content(v); /* gcd(lc(x), cont(y)) */
3422 1706468 : t = simplify_shallow(t);
3423 1706468 : if (dx < e) e = dx;
3424 1706468 : return gerepileupto(av, monomialcopy(t, e, varn(x)));
3425 : }
3426 :
3427 : static GEN
3428 198912 : RgX_gcd_FpX(GEN x, GEN y, GEN p)
3429 : {
3430 198912 : pari_sp av = avma;
3431 : GEN r;
3432 198912 : if (lgefint(p) == 3)
3433 : {
3434 198898 : ulong pp = uel(p, 2);
3435 198898 : r = Flx_to_ZX_inplace(Flx_gcd(RgX_to_Flx(x, pp),
3436 : RgX_to_Flx(y, pp), pp));
3437 : }
3438 : else
3439 14 : r = FpX_gcd(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3440 198912 : return gerepileupto(av, FpX_to_mod(r, p));
3441 : }
3442 :
3443 : static GEN
3444 7 : RgX_gcd_FpXQX(GEN x, GEN y, GEN pol, GEN p)
3445 : {
3446 7 : pari_sp av = avma;
3447 7 : GEN r, T = RgX_to_FpX(pol, p);
3448 7 : if (signe(T)==0) pari_err_OP("gcd", x, y);
3449 7 : r = FpXQX_gcd(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
3450 7 : return gerepileupto(av, FpXQX_to_mod(r, T, p));
3451 : }
3452 :
3453 : static GEN
3454 7392 : RgX_liftred(GEN x, GEN T)
3455 7392 : { return RgXQX_red(liftpol_shallow(x), T); }
3456 :
3457 : static GEN
3458 2268 : RgX_gcd_ZXQX(GEN x, GEN y, GEN T)
3459 : {
3460 2268 : pari_sp av = avma;
3461 2268 : GEN r = ZXQX_gcd(RgX_liftred(x, T), RgX_liftred(y, T), T);
3462 2268 : return gerepilecopy(av, QXQX_to_mod_shallow(r, T));
3463 : }
3464 :
3465 : static GEN
3466 1428 : RgX_gcd_QXQX(GEN x, GEN y, GEN T)
3467 : {
3468 1428 : pari_sp av = avma;
3469 1428 : GEN r = QXQX_gcd(RgX_liftred(x, T), RgX_liftred(y, T), T);
3470 1428 : return gerepilecopy(av, QXQX_to_mod_shallow(r, T));
3471 : }
3472 :
3473 : static GEN
3474 11324206 : RgX_gcd_fast(GEN x, GEN y)
3475 : {
3476 : GEN p, pol;
3477 : long pa;
3478 11324206 : long t = RgX_type2(x,y, &p,&pol,&pa);
3479 11324206 : switch(t)
3480 : {
3481 9406842 : case t_INT: return ZX_gcd(x, y);
3482 2499 : case t_FRAC: return QX_gcd(x, y);
3483 2520 : case t_FFELT: return FFX_gcd(x, y, pol);
3484 198912 : case t_INTMOD: return RgX_gcd_FpX(x, y, p);
3485 7 : case code(t_POLMOD, t_INTMOD):
3486 7 : return RgX_gcd_FpXQX(x, y, pol, p);
3487 2275 : case code(t_POLMOD, t_INT):
3488 2275 : return ZX_is_monic(pol)? RgX_gcd_ZXQX(x,y,pol): NULL;
3489 1442 : case code(t_POLMOD, t_FRAC):
3490 2884 : return RgX_is_ZX(pol) && ZX_is_monic(pol) ?
3491 2884 : RgX_gcd_QXQX(x,y,pol): NULL;
3492 1709709 : default: return NULL;
3493 : }
3494 : }
3495 :
3496 : /* x, y are t_POL in the same variable */
3497 : GEN
3498 11324206 : RgX_gcd(GEN x, GEN y)
3499 : {
3500 : long dx, dy;
3501 : pari_sp av, av1;
3502 : GEN d, g, h, p1, p2, u, v;
3503 11324206 : int simple = 0;
3504 11324206 : GEN z = RgX_gcd_fast(x, y);
3505 11324206 : if (z) return z;
3506 1709730 : if (isexactzero(y)) return RgX_copy(x);
3507 1709632 : if (isexactzero(x)) return RgX_copy(y);
3508 1709632 : if (RgX_is_monomial(x)) return gcdmonome(x,y);
3509 4270 : if (RgX_is_monomial(y)) return gcdmonome(y,x);
3510 3164 : if (isinexactall(x,&simple) || isinexactall(y,&simple))
3511 : {
3512 7 : av = avma; u = ggcd(content(x), content(y));
3513 7 : return gerepileupto(av, scalarpol(u, varn(x)));
3514 : }
3515 :
3516 3157 : av = avma;
3517 3157 : if (simple > 0) x = RgX_gcd_simple(x,y);
3518 : else
3519 : {
3520 3157 : dx = lg(x); dy = lg(y);
3521 3157 : if (dx < dy) { swap(x,y); lswap(dx,dy); }
3522 3157 : if (dy==3)
3523 : {
3524 0 : d = ggcd(gel(y,2), content(x));
3525 0 : return gerepileupto(av, scalarpol(d, varn(x)));
3526 : }
3527 3157 : u = primitive_part(x, &p1); if (!p1) p1 = gen_1;
3528 3157 : v = primitive_part(y, &p2); if (!p2) p2 = gen_1;
3529 3157 : d = ggcd(p1,p2);
3530 3157 : av1 = avma;
3531 3157 : g = h = gen_1;
3532 : for(;;)
3533 1197 : {
3534 4354 : GEN r = RgX_pseudorem(u,v);
3535 4354 : long degq, du, dv, dr = lg(r);
3536 :
3537 4354 : if (!signe(r)) break;
3538 2205 : if (dr <= 3)
3539 : {
3540 1008 : set_avma(av1);
3541 1008 : return gerepileupto(av, scalarpol(d, varn(x)));
3542 : }
3543 1197 : du = lg(u); dv = lg(v); degq = du-dv;
3544 1197 : u = v; p1 = g; g = leading_coeff(u);
3545 1197 : switch(degq)
3546 : {
3547 189 : case 0: break;
3548 917 : case 1:
3549 917 : p1 = gmul(h,p1); h = g; break;
3550 91 : default:
3551 91 : p1 = gmul(gpowgs(h,degq), p1);
3552 91 : h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
3553 : }
3554 1197 : v = RgX_Rg_div(r,p1);
3555 1197 : if (gc_needed(av1,1))
3556 : {
3557 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd, dr = %ld", degpol(r));
3558 0 : gerepileall(av1,4, &u,&v,&g,&h);
3559 : }
3560 : }
3561 2149 : x = RgX_Rg_mul(primpart(v), d);
3562 : }
3563 2149 : if (must_negate(x)) x = RgX_neg(x);
3564 2149 : return gerepileupto(av,x);
3565 : }
3566 :
3567 : /* disc P = (-1)^(n(n-1)/2) lc(P)^(n - deg P' - 2) Res(P,P'), n = deg P */
3568 : static GEN
3569 308 : RgX_disc_i(GEN P)
3570 : {
3571 308 : long n = degpol(P), dd;
3572 : GEN N, D, L, y;
3573 308 : if (!signe(P) || !n) return Rg_get_0(P);
3574 301 : if (n == 1) return Rg_get_1(P);
3575 301 : if (n == 2) {
3576 105 : GEN a = gel(P,4), b = gel(P,3), c = gel(P,2);
3577 105 : return gsub(gsqr(b), gmul2n(gmul(a,c),2));
3578 : }
3579 196 : y = RgX_deriv(P);
3580 196 : N = characteristic(P);
3581 196 : if (signe(N)) y = gmul(y, mkintmod(gen_1,N));
3582 196 : if (!signe(y)) return Rg_get_0(y);
3583 196 : dd = n - 2 - degpol(y);
3584 196 : if (isinexact(P))
3585 14 : D = resultant2(P,y);
3586 : else
3587 : {
3588 182 : D = RgX_resultant_all(P, y, NULL);
3589 182 : if (D == gen_0) return Rg_get_0(y);
3590 : }
3591 196 : L = leading_coeff(P);
3592 196 : if (dd && !gequal1(L)) D = (dd == -1)? gdiv(D, L): gmul(D, gpowgs(L, dd));
3593 196 : if (n & 2) D = gneg(D);
3594 196 : return D;
3595 : }
3596 :
3597 : static GEN
3598 42 : RgX_disc_FpX(GEN x, GEN p)
3599 : {
3600 42 : pari_sp av = avma;
3601 42 : GEN r = FpX_disc(RgX_to_FpX(x, p), p);
3602 42 : return gerepileupto(av, Fp_to_mod(r, p));
3603 : }
3604 :
3605 : static GEN
3606 28 : RgX_disc_FpXQX(GEN x, GEN pol, GEN p)
3607 : {
3608 28 : pari_sp av = avma;
3609 28 : GEN r, T = RgX_to_FpX(pol, p);
3610 28 : r = FpXQX_disc(RgX_to_FpXQX(x, T, p), T, p);
3611 28 : return gerepileupto(av, FpX_to_mod(r, p));
3612 : }
3613 :
3614 : static GEN
3615 122902 : RgX_disc_fast(GEN x)
3616 : {
3617 : GEN p, pol;
3618 : long pa;
3619 122902 : long t = RgX_type(x, &p,&pol,&pa);
3620 122902 : switch(t)
3621 : {
3622 122482 : case t_INT: return ZX_disc(x);
3623 7 : case t_FRAC: return QX_disc(x);
3624 35 : case t_FFELT: return FFX_disc(x, pol);
3625 42 : case t_INTMOD: return RgX_disc_FpX(x, p);
3626 28 : case code(t_POLMOD, t_INTMOD):
3627 28 : return RgX_disc_FpXQX(x, pol, p);
3628 308 : default: return NULL;
3629 : }
3630 : }
3631 :
3632 : GEN
3633 122902 : RgX_disc(GEN x)
3634 : {
3635 : pari_sp av;
3636 122902 : GEN z = RgX_disc_fast(x);
3637 122902 : if (z) return z;
3638 308 : av = avma;
3639 308 : return gerepileupto(av, RgX_disc_i(x));
3640 : }
3641 :
3642 : GEN
3643 4672 : poldisc0(GEN x, long v)
3644 : {
3645 4672 : long v0, tx = typ(x);
3646 : pari_sp av;
3647 : GEN D;
3648 4672 : if (tx == t_POL && (v < 0 || v == varn(x))) return RgX_disc(x);
3649 35 : switch(tx)
3650 : {
3651 0 : case t_QUAD:
3652 0 : return quad_disc(x);
3653 0 : case t_POLMOD:
3654 0 : if (v >= 0 && varn(gel(x,1)) != v) break;
3655 0 : return RgX_disc(gel(x,1));
3656 14 : case t_QFB:
3657 14 : return icopy(qfb_disc(x));
3658 0 : case t_VEC: case t_COL: case t_MAT:
3659 : {
3660 : long i;
3661 0 : GEN z = cgetg_copy(x, &i);
3662 0 : for (i--; i; i--) gel(z,i) = poldisc0(gel(x,i), v);
3663 0 : return z;
3664 : }
3665 : }
3666 21 : if (v < 0) pari_err_TYPE("poldisc",x);
3667 21 : av = avma; v0 = fetch_var_higher();
3668 21 : x = fix_pol(x,v, v0);
3669 14 : D = RgX_disc(x); (void)delete_var();
3670 14 : return gerepileupto(av, D);
3671 : }
3672 :
3673 : GEN
3674 7 : reduceddiscsmith(GEN x)
3675 : {
3676 7 : long j, n = degpol(x);
3677 7 : pari_sp av = avma;
3678 : GEN xp, M;
3679 :
3680 7 : if (typ(x) != t_POL) pari_err_TYPE("poldiscreduced",x);
3681 7 : if (n<=0) pari_err_CONSTPOL("poldiscreduced");
3682 7 : RgX_check_ZX(x,"poldiscreduced");
3683 7 : if (!gequal1(gel(x,n+2)))
3684 0 : pari_err_IMPL("nonmonic polynomial in poldiscreduced");
3685 7 : M = cgetg(n+1,t_MAT);
3686 7 : xp = ZX_deriv(x);
3687 28 : for (j=1; j<=n; j++)
3688 : {
3689 21 : gel(M,j) = RgX_to_RgC(xp, n);
3690 21 : if (j<n) xp = RgX_rem(RgX_shift_shallow(xp, 1), x);
3691 : }
3692 7 : return gerepileupto(av, ZM_snf(M));
3693 : }
3694 :
3695 : /***********************************************************************/
3696 : /** **/
3697 : /** STURM ALGORITHM **/
3698 : /** (number of real roots of x in [a,b]) **/
3699 : /** **/
3700 : /***********************************************************************/
3701 : static GEN
3702 1169 : R_to_Q_up(GEN x)
3703 : {
3704 : long e;
3705 1169 : switch(typ(x))
3706 : {
3707 1169 : case t_INT: case t_FRAC: case t_INFINITY: return x;
3708 0 : case t_REAL:
3709 0 : x = mantissa_real(x,&e);
3710 0 : return gmul2n(addiu(x,1), -e);
3711 0 : default: pari_err_TYPE("R_to_Q_up", x);
3712 : return NULL; /* LCOV_EXCL_LINE */
3713 : }
3714 : }
3715 : static GEN
3716 1169 : R_to_Q_down(GEN x)
3717 : {
3718 : long e;
3719 1169 : switch(typ(x))
3720 : {
3721 1155 : case t_INT: case t_FRAC: case t_INFINITY: return x;
3722 14 : case t_REAL:
3723 14 : x = mantissa_real(x,&e);
3724 14 : return gmul2n(subiu(x,1), -e);
3725 0 : default: pari_err_TYPE("R_to_Q_down", x);
3726 : return NULL; /* LCOV_EXCL_LINE */
3727 : }
3728 : }
3729 :
3730 : static long
3731 1169 : sturmpart_i(GEN x, GEN ab)
3732 : {
3733 1169 : long tx = typ(x);
3734 1169 : if (gequal0(x)) pari_err_ROOTS0("sturm");
3735 1169 : if (tx != t_POL)
3736 : {
3737 0 : if (is_real_t(tx)) return 0;
3738 0 : pari_err_TYPE("sturm",x);
3739 : }
3740 1169 : if (lg(x) == 3) return 0;
3741 1169 : if (!RgX_is_ZX(x)) x = RgX_rescale_to_int(x);
3742 1169 : (void)ZX_gcd_all(x, ZX_deriv(x), &x);
3743 1169 : if (ab)
3744 : {
3745 : GEN A, B;
3746 1169 : if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("RgX_sturmpart", ab);
3747 1169 : A = R_to_Q_down(gel(ab,1));
3748 1169 : B = R_to_Q_up(gel(ab,2));
3749 1169 : ab = mkvec2(A, B);
3750 : }
3751 1169 : return ZX_sturmpart(x, ab);
3752 : }
3753 : /* Deprecated: RgX_sturmpart() should be preferred */
3754 : long
3755 1169 : sturmpart(GEN x, GEN a, GEN b)
3756 : {
3757 1169 : pari_sp av = avma;
3758 1169 : if (!b && a && typ(a) == t_VEC) return RgX_sturmpart(x, a);
3759 1029 : if (!a) a = mkmoo();
3760 1029 : if (!b) b = mkoo();
3761 1029 : return gc_long(av, sturmpart_i(x, mkvec2(a,b)));
3762 : }
3763 : long
3764 140 : RgX_sturmpart(GEN x, GEN ab)
3765 140 : { pari_sp av = avma; return gc_long(av, sturmpart_i(x, ab)); }
3766 :
3767 : /***********************************************************************/
3768 : /** **/
3769 : /** GENERIC EXTENDED GCD **/
3770 : /** **/
3771 : /***********************************************************************/
3772 : /* assume typ(x) = typ(y) = t_POL */
3773 : static GEN
3774 904 : RgXQ_inv_i(GEN x, GEN y)
3775 : {
3776 904 : long vx=varn(x), vy=varn(y);
3777 : pari_sp av;
3778 : GEN u, v, d;
3779 :
3780 904 : while (vx != vy)
3781 : {
3782 0 : if (varncmp(vx,vy) > 0)
3783 : {
3784 0 : d = (vx == NO_VARIABLE)? ginv(x): gred_rfrac_simple(gen_1, x);
3785 0 : return scalarpol(d, vy);
3786 : }
3787 0 : if (lg(x)!=3) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
3788 0 : x = gel(x,2); vx = gvar(x);
3789 : }
3790 904 : av = avma; d = subresext_i(x,y,&u,&v/*junk*/);
3791 904 : if (gequal0(d)) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
3792 904 : d = gdiv(u,d);
3793 904 : if (typ(d) != t_POL || varn(d) != vy) d = scalarpol(d, vy);
3794 904 : return gerepileupto(av, d);
3795 : }
3796 :
3797 : /*Assume x is a polynomial and y is not */
3798 : static GEN
3799 112 : scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
3800 : {
3801 112 : long vx = varn(x);
3802 112 : int xis0 = signe(x)==0, yis0 = gequal0(y);
3803 112 : if (xis0 && yis0) { *U = *V = pol_0(vx); return pol_0(vx); }
3804 84 : if (yis0) { *U=pol_1(vx); *V = pol_0(vx); return RgX_copy(x);}
3805 56 : *U=pol_0(vx); *V= ginv(y); return pol_1(vx);
3806 : }
3807 : /* Assume x==0, y!=0 */
3808 : static GEN
3809 63 : zero_bezout(GEN y, GEN *U, GEN *V)
3810 : {
3811 63 : *U=gen_0; *V = ginv(y); return gen_1;
3812 : }
3813 :
3814 : GEN
3815 350 : gbezout(GEN x, GEN y, GEN *u, GEN *v)
3816 : {
3817 350 : long tx=typ(x), ty=typ(y), vx;
3818 350 : if (tx == t_INT && ty == t_INT) return bezout(x,y,u,v);
3819 315 : if (tx != t_POL)
3820 : {
3821 140 : if (ty == t_POL)
3822 56 : return scalar_bezout(y,x,v,u);
3823 : else
3824 : {
3825 84 : int xis0 = gequal0(x), yis0 = gequal0(y);
3826 84 : if (xis0 && yis0) { *u = *v = gen_0; return gen_0; }
3827 63 : if (xis0) return zero_bezout(y,u,v);
3828 42 : else return zero_bezout(x,v,u);
3829 : }
3830 : }
3831 175 : else if (ty != t_POL) return scalar_bezout(x,y,u,v);
3832 119 : vx = varn(x);
3833 119 : if (vx != varn(y))
3834 0 : return varncmp(vx, varn(y)) < 0? scalar_bezout(x,y,u,v)
3835 0 : : scalar_bezout(y,x,v,u);
3836 119 : return RgX_extgcd(x,y,u,v);
3837 : }
3838 :
3839 : GEN
3840 350 : gcdext0(GEN x, GEN y)
3841 : {
3842 350 : GEN z=cgetg(4,t_VEC);
3843 350 : gel(z,3) = gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
3844 350 : return z;
3845 : }
3846 :
3847 : /*******************************************************************/
3848 : /* */
3849 : /* GENERIC (modular) INVERSE */
3850 : /* */
3851 : /*******************************************************************/
3852 :
3853 : GEN
3854 27251 : ginvmod(GEN x, GEN y)
3855 : {
3856 27251 : long tx=typ(x);
3857 :
3858 27251 : switch(typ(y))
3859 : {
3860 27251 : case t_POL:
3861 27251 : if (tx==t_POL) return RgXQ_inv(x,y);
3862 13727 : if (is_scalar_t(tx)) return ginv(x);
3863 0 : break;
3864 0 : case t_INT:
3865 0 : if (tx==t_INT) return Fp_inv(x,y);
3866 0 : if (tx==t_POL) return gen_0;
3867 : }
3868 0 : pari_err_TYPE2("ginvmod",x,y);
3869 : return NULL; /* LCOV_EXCL_LINE */
3870 : }
3871 :
3872 : /***********************************************************************/
3873 : /** **/
3874 : /** NEWTON POLYGON **/
3875 : /** **/
3876 : /***********************************************************************/
3877 :
3878 : /* assume leading coeff of x is nonzero */
3879 : GEN
3880 28 : newtonpoly(GEN x, GEN p)
3881 : {
3882 28 : pari_sp av = avma;
3883 : long n, ind, a, b;
3884 : GEN y, vval;
3885 :
3886 28 : if (typ(x) != t_POL) pari_err_TYPE("newtonpoly",x);
3887 28 : n = degpol(x); if (n<=0) return cgetg(1,t_VEC);
3888 28 : vval = new_chunk(n+1);
3889 28 : y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
3890 168 : for (a = 0; a <= n; a++) vval[a] = gvaluation(gel(x,a),p);
3891 42 : for (a = 0, ind = 1; a < n; a++)
3892 : {
3893 42 : if (vval[a] != LONG_MAX) break;
3894 14 : gel(y,ind++) = mkoo();
3895 : }
3896 84 : for (b = a+1; b <= n; a = b, b = a+1)
3897 : {
3898 : long u1, u2, c;
3899 70 : while (vval[b] == LONG_MAX) b++;
3900 56 : u1 = vval[a] - vval[b];
3901 56 : u2 = b - a;
3902 154 : for (c = b+1; c <= n; c++)
3903 : {
3904 : long r1, r2;
3905 98 : if (vval[c] == LONG_MAX) continue;
3906 70 : r1 = vval[a] - vval[c];
3907 70 : r2 = c - a;
3908 70 : if (u1*r2 <= u2*r1) { u1 = r1; u2 = r2; b = c; }
3909 : }
3910 154 : while (ind <= b) gel(y,ind++) = sstoQ(u1,u2);
3911 : }
3912 28 : stackdummy((pari_sp)vval, av); return y;
3913 : }
3914 :
3915 : static GEN
3916 274288 : RgXQ_mul_FpXQ(GEN x, GEN y, GEN T, GEN p)
3917 : {
3918 274288 : pari_sp av = avma;
3919 : GEN r;
3920 274288 : if (lgefint(p) == 3)
3921 : {
3922 152381 : ulong pp = uel(p, 2);
3923 152381 : r = Flx_to_ZX_inplace(Flxq_mul(RgX_to_Flx(x, pp),
3924 : RgX_to_Flx(y, pp), RgX_to_Flx(T, pp), pp));
3925 : }
3926 : else
3927 121907 : r = FpXQ_mul(RgX_to_FpX(x, p), RgX_to_FpX(y, p), RgX_to_FpX(T, p), p);
3928 274288 : return gerepileupto(av, FpX_to_mod(r, p));
3929 : }
3930 :
3931 : static GEN
3932 14 : RgXQ_sqr_FpXQ(GEN x, GEN y, GEN p)
3933 : {
3934 14 : pari_sp av = avma;
3935 : GEN r;
3936 14 : if (lgefint(p) == 3)
3937 : {
3938 7 : ulong pp = uel(p, 2);
3939 7 : r = Flx_to_ZX_inplace(Flxq_sqr(RgX_to_Flx(x, pp),
3940 : RgX_to_Flx(y, pp), pp));
3941 : }
3942 : else
3943 7 : r = FpXQ_sqr(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3944 14 : return gerepileupto(av, FpX_to_mod(r, p));
3945 : }
3946 :
3947 : static GEN
3948 12054 : RgXQ_inv_FpXQ(GEN x, GEN y, GEN p)
3949 : {
3950 12054 : pari_sp av = avma;
3951 : GEN r;
3952 12054 : if (lgefint(p) == 3)
3953 : {
3954 6088 : ulong pp = uel(p, 2);
3955 6088 : r = Flx_to_ZX_inplace(Flxq_inv(RgX_to_Flx(x, pp),
3956 : RgX_to_Flx(y, pp), pp));
3957 : }
3958 : else
3959 5966 : r = FpXQ_inv(RgX_to_FpX(x, p), RgX_to_FpX(y, p), p);
3960 12054 : return gerepileupto(av, FpX_to_mod(r, p));
3961 : }
3962 :
3963 : static GEN
3964 385 : RgXQ_mul_FpXQXQ(GEN x, GEN y, GEN S, GEN pol, GEN p)
3965 : {
3966 385 : pari_sp av = avma;
3967 : GEN r;
3968 385 : GEN T = RgX_to_FpX(pol, p);
3969 385 : if (signe(T)==0) pari_err_OP("*",x,y);
3970 385 : if (lgefint(p) == 3)
3971 : {
3972 241 : ulong pp = uel(p, 2);
3973 241 : GEN Tp = ZX_to_Flx(T, pp);
3974 241 : r = FlxX_to_ZXX(FlxqXQ_mul(RgX_to_FlxqX(x, Tp, pp),
3975 : RgX_to_FlxqX(y, Tp, pp),
3976 : RgX_to_FlxqX(S, Tp, pp), Tp, pp));
3977 : }
3978 : else
3979 144 : r = FpXQXQ_mul(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p),
3980 : RgX_to_FpXQX(S, T, p), T, p);
3981 385 : return gerepileupto(av, FpXQX_to_mod(r, T, p));
3982 : }
3983 :
3984 : static GEN
3985 0 : RgXQ_sqr_FpXQXQ(GEN x, GEN y, GEN pol, GEN p)
3986 : {
3987 0 : pari_sp av = avma;
3988 : GEN r;
3989 0 : GEN T = RgX_to_FpX(pol, p);
3990 0 : if (signe(T)==0) pari_err_OP("*",x,x);
3991 0 : if (lgefint(p) == 3)
3992 : {
3993 0 : ulong pp = uel(p, 2);
3994 0 : GEN Tp = ZX_to_Flx(T, pp);
3995 0 : r = FlxX_to_ZXX(FlxqXQ_sqr(RgX_to_FlxqX(x, Tp, pp),
3996 : RgX_to_FlxqX(y, Tp, pp), Tp, pp));
3997 : }
3998 : else
3999 0 : r = FpXQXQ_sqr(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
4000 0 : return gerepileupto(av, FpXQX_to_mod(r, T, p));
4001 : }
4002 :
4003 : static GEN
4004 7 : RgXQ_inv_FpXQXQ(GEN x, GEN y, GEN pol, GEN p)
4005 : {
4006 7 : pari_sp av = avma;
4007 : GEN r;
4008 7 : GEN T = RgX_to_FpX(pol, p);
4009 7 : if (signe(T)==0) pari_err_OP("^",x,gen_m1);
4010 7 : if (lgefint(p) == 3)
4011 : {
4012 7 : ulong pp = uel(p, 2);
4013 7 : GEN Tp = ZX_to_Flx(T, pp);
4014 7 : r = FlxX_to_ZXX(FlxqXQ_inv(RgX_to_FlxqX(x, Tp, pp),
4015 : RgX_to_FlxqX(y, Tp, pp), Tp, pp));
4016 : }
4017 : else
4018 0 : r = FpXQXQ_inv(RgX_to_FpXQX(x, T, p), RgX_to_FpXQX(y, T, p), T, p);
4019 7 : return gerepileupto(av, FpXQX_to_mod(r, T, p));
4020 : }
4021 :
4022 : static GEN
4023 839214 : RgXQ_mul_fast(GEN x, GEN y, GEN T)
4024 : {
4025 : GEN p, pol;
4026 : long pa;
4027 839214 : long t = RgX_type3(x,y,T, &p,&pol,&pa);
4028 839212 : switch(t)
4029 : {
4030 252004 : case t_INT: return ZX_is_monic(T) ? ZXQ_mul(x,y,T): NULL;
4031 277755 : case t_FRAC: return RgX_is_ZX(T) && ZX_is_monic(T) ? QXQ_mul(x,y,T): NULL;
4032 105 : case t_FFELT: return FFXQ_mul(x, y, T, pol);
4033 274288 : case t_INTMOD: return RgXQ_mul_FpXQ(x, y, T, p);
4034 385 : case code(t_POLMOD, t_INTMOD):
4035 385 : return RgXQ_mul_FpXQXQ(x, y, T, pol, p);
4036 34675 : default: return NULL;
4037 : }
4038 : }
4039 :
4040 : GEN
4041 839215 : RgXQ_mul(GEN x, GEN y, GEN T)
4042 : {
4043 839215 : GEN z = RgXQ_mul_fast(x, y, T);
4044 839216 : if (!z) z = RgX_rem(RgX_mul(x, y), T);
4045 839216 : return z;
4046 : }
4047 :
4048 : static GEN
4049 104391 : RgXQ_sqr_fast(GEN x, GEN T)
4050 : {
4051 : GEN p, pol;
4052 : long pa;
4053 104391 : long t = RgX_type2(x, T, &p,&pol,&pa);
4054 104391 : switch(t)
4055 : {
4056 78176 : case t_INT: return ZX_is_monic(T) ? ZXQ_sqr(x,T): NULL;
4057 19308 : case t_FRAC: return RgX_is_ZX(T) && ZX_is_monic(T) ? QXQ_sqr(x,T): NULL;
4058 7 : case t_FFELT: return FFXQ_sqr(x, T, pol);
4059 14 : case t_INTMOD: return RgXQ_sqr_FpXQ(x, T, p);
4060 0 : case code(t_POLMOD, t_INTMOD):
4061 0 : return RgXQ_sqr_FpXQXQ(x, T, pol, p);
4062 6886 : default: return NULL;
4063 : }
4064 : }
4065 :
4066 : GEN
4067 104391 : RgXQ_sqr(GEN x, GEN T)
4068 : {
4069 104391 : GEN z = RgXQ_sqr_fast(x, T);
4070 104391 : if (!z) z = RgX_rem(RgX_sqr(x), T);
4071 104391 : return z;
4072 : }
4073 :
4074 : static GEN
4075 19709 : RgXQ_inv_fast(GEN x, GEN y)
4076 : {
4077 : GEN p, pol;
4078 : long pa;
4079 19709 : long t = RgX_type2(x,y, &p,&pol,&pa);
4080 19709 : switch(t)
4081 : {
4082 5869 : case t_INT: return QXQ_inv(x,y);
4083 868 : case t_FRAC: return RgX_is_ZX(y)? QXQ_inv(x,y): NULL;
4084 14 : case t_FFELT: return FFXQ_inv(x, y, pol);
4085 12054 : case t_INTMOD: return RgXQ_inv_FpXQ(x, y, p);
4086 7 : case code(t_POLMOD, t_INTMOD):
4087 7 : return RgXQ_inv_FpXQXQ(x, y, pol, p);
4088 897 : default: return NULL;
4089 : }
4090 : }
4091 :
4092 : GEN
4093 19709 : RgXQ_inv(GEN x, GEN y)
4094 : {
4095 19709 : GEN z = RgXQ_inv_fast(x, y);
4096 19695 : if (!z) z = RgXQ_inv_i(x, y);
4097 19695 : return z;
4098 : }
|
