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